diff --git a/RNAlien.cabal b/RNAlien.cabal --- a/RNAlien.cabal +++ b/RNAlien.cabal @@ -1,5 +1,5 @@ name: RNAlien -version: 1.3.1 +version: 1.3.7 synopsis: Unsupervized construction of RNA family models description: RNAlien is a tool for automatic construction of RNAfamily models from a single sequence. . @@ -50,42 +50,36 @@ source-repository this type: git - location: https://github.com/eggzilla/RNAlien/tree/1.3.1 - tag: 1.3.1 + location: https://github.com/eggzilla/RNAlien/tree/1.3.7 + tag: 1.3.7 executable RNAlien Hs-Source-Dirs: ./src/Bio/ main-is: RNAlien.hs - ghc-options: -Wall -O2 + ghc-options: -Wall build-depends: base >=4.5 && <5, cmdargs, directory, biofasta, random, biocore, containers, RNAlien, time, either-unwrap, filepath executable RNAlienStatistics Hs-Source-Dirs: ./src/Bio/ main-is: RNAlienStatistics.hs - ghc-options: -Wall -O2 + ghc-options: -Wall build-depends: base >=4.5 && <5, cmdargs, cassava, vector, process, bytestring, either-unwrap, RNAlien, directory, biofasta, biocore, split, filepath, ViennaRNAParser>=1.3.2 -executable SelectSequences - Hs-Source-Dirs: ./src/Bio/ - main-is: SelectSequences.hs - ghc-options: -Wall -O2 - build-depends: base >=4.5 && <5, cmdargs, either-unwrap, RNAlien - executable cmsearchToBed Hs-Source-Dirs: ./src/Bio/ main-is: cmsearchToBED.hs - ghc-options: -Wall -O2 + ghc-options: -Wall build-depends: base >=4.5 && <5, cmdargs, either-unwrap, RNAlien, bytestring, text executable RNAcentralHTTPRequest Hs-Source-Dirs: ./src/Bio/ main-is: RNAcentralHTTPRequest.hs - ghc-options: -Wall -O2 + ghc-options: -Wall build-depends: base >=4.5 && <5, cmdargs, either-unwrap, RNAlien Library Hs-Source-Dirs: ./src/ - ghc-options: -Wall -O2 -fno-warn-unused-do-bind + ghc-options: -Wall -fno-warn-unused-do-bind build-depends: base >=4.5 && <5, cmdargs, ViennaRNAParser>=1.3.2, process, directory, blastxml>=0.3.2, biofasta, parsec, random, BlastHTTP>=1.2.1, biocore, bytestring, Taxonomy >= 1.0.2, either-unwrap, containers, ClustalParser>=1.2.1, EntrezHTTP>=1.0.3, vector, edit-distance, cassava, matrix, hierarchical-clustering, filepath, HTTP, http-conduit, hxt, network, aeson, text, transformers, pureMD5, http-types, text-metrics Exposed-Modules: Bio.RNAlienData Bio.RNAlienLibrary diff --git a/src/Bio/InfernalParser.hs b/src/Bio/InfernalParser.hs --- a/src/Bio/InfernalParser.hs +++ b/src/Bio/InfernalParser.hs @@ -1,7 +1,7 @@ -- | This module contains parsing functions for Infernal programs module Bio.InfernalParser ( - module Bio.RNAlienData, + module Bio.RNAlienData, readCMSearch, readCMSearches, parseCMSearch, @@ -11,30 +11,30 @@ ) where -import Text.ParserCombinators.Parsec +import Text.ParserCombinators.Parsec import Bio.RNAlienData import qualified Data.ByteString.Lazy.Char8 as L import qualified Control.Exception.Base as CE -- | parse from input filePath parseCMSearch :: String -> Either ParseError CMsearch -parseCMSearch = parse genParserCMSearch "parseCMsearch" +parseCMSearch = parse genParserCMSearch "parseCMsearch" -- | parse from input filePath parseCMSearches :: String -> Either ParseError CMsearch parseCMSearches = parse genParserCMSearches "parseCMsearch" -- | parse from input filePath -readCMSearch :: String -> IO (Either ParseError CMsearch) -readCMSearch filePath = do +readCMSearch :: String -> IO (Either ParseError CMsearch) +readCMSearch filePath = do parsedFile <- parseFromFile genParserCMSearch filePath - CE.evaluate parsedFile + CE.evaluate parsedFile -- | parse from input filePath -readCMSearches :: String -> IO (Either ParseError CMsearch) -readCMSearches filePath = do +readCMSearches :: String -> IO (Either ParseError CMsearch) +readCMSearches filePath = do parsedFile <- parseFromFile genParserCMSearches filePath - CE.evaluate parsedFile + CE.evaluate parsedFile genParserCMSearches :: GenParser Char st CMsearch genParserCMSearches = do @@ -42,12 +42,12 @@ newline string "# INFERNAL " many1 (noneOf "\n") - newline + newline string "# Copyright (C) 201" many1 (noneOf "\n") - newline + newline string "# Freely distributed under the GNU General Public License (GPLv3)." - newline + newline string "# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -" newline string "# query CM file:" @@ -55,7 +55,7 @@ queryCMfile' <- many1 (noneOf "\n") newline string "# target sequence database:" - many1 space + many1 space targetSequenceDatabase' <- many1 (noneOf "\n") newline optional (try (genParserCMsearchHeaderField "# CM configuration")) @@ -72,19 +72,19 @@ optional (string "[ok]\n") eof return $ CMsearch queryCMfile' targetSequenceDatabase' numberOfWorkerThreads' (concat cmSearchesHits) - + genParserCMSearch :: GenParser Char st CMsearch genParserCMSearch = do string "# cmsearch :: search CM(s) against a sequence database" newline string "# INFERNAL " many1 (noneOf "\n") - newline + newline string "# Copyright (C) 201" many1 (noneOf "\n") - newline + newline string "# Freely distributed under the GNU General Public License (GPLv3)." - newline + newline string "# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -" newline string "# query CM file:" @@ -92,7 +92,7 @@ queryCMfile' <- many1 (noneOf "\n") newline string "# target sequence database:" - many1 space + many1 space targetSequenceDatabase' <- many1 (noneOf "\n") newline optional (try (genParserCMsearchHeaderField "# CM configuration")) @@ -106,32 +106,32 @@ newline optional newline string "Query:" - many1 (noneOf "\n") + many1 (noneOf "\n") newline optional (try (genParserCMsearchHeaderField "Accession")) optional (try (genParserCMsearchHeaderField "Description")) string "Hit scores:" newline choice [try (string " rank"), try (string " rank") , try (string " rank"), try (string " rank"),try (string " rank"),try (string " rank")] - many1 space + many1 space string "E-value" - many1 space + many1 space string "score" - many1 space + many1 space string "bias" - many1 space + many1 space string "sequence" - many1 space + many1 space string "start" - many1 space + many1 space string "end" - many1 space + many1 space string "mdl" - many1 space + many1 space string "trunc" - many1 space + many1 space string "gc" - many1 space + many1 space string "description" newline string " -" @@ -152,32 +152,32 @@ --optional newline --optional string "//" string "Query:" - many1 (noneOf "\n") + many1 (noneOf "\n") newline optional (try (genParserCMsearchHeaderField "Accession")) optional (try (genParserCMsearchHeaderField "Description")) string "Hit scores:" newline choice [try (string " rank"), try (string " rank") , try (string " rank"), try (string " rank"),try (string " rank"),try (string " rank")] - many1 space + many1 space string "E-value" - many1 space + many1 space string "score" - many1 space + many1 space string "bias" - many1 space + many1 space string "sequence" - many1 space + many1 space string "start" - many1 space + many1 space string "end" - many1 space + many1 space string "mdl" - many1 space + many1 space string "trunc" - many1 space + many1 space string "gc" - many1 space + many1 space string "description" newline string " -" @@ -210,33 +210,33 @@ genParserCMsearchHit :: GenParser Char st CMsearchHit genParserCMsearchHit = do many1 space - string "(" + string "(" hitRank' <- many1 digit string ")" many1 space hitSignificant' <- choice [char '!', char '?'] - many1 space + many1 space hitEValue' <- many1 (oneOf "0123456789.e-") - many1 space + many1 space hitScore' <- many1 (oneOf "0123456789.e-") - many1 space + many1 space hitBias' <- many1 (oneOf "0123456789.e-") many1 space hitSequenceHeader' <- many1 (noneOf " ") - many1 space + many1 space hitStart' <- many1 digit many1 space hitEnd' <- many1 digit - many1 space + many1 space hitStrand' <- choice [char '+', char '-', char '.'] - many1 space + many1 space hitModel' <- many1 letter - many1 space + many1 space hitTruncation' <- many1 (choice [alphaNum, char '\'']) - many1 space + many1 space hitGCcontent' <- many1 (oneOf "0123456789.e-") - many1 space - hitDescription' <- many1 (noneOf "\n") + many1 space + hitDescription' <- many1 (noneOf "\n") newline optional (try (string " ------ inclusion threshold ------")) optional (try newline) @@ -247,23 +247,23 @@ parseCMstat = parse genParserCMstat "parseCMstat" -- | parse from input filePath -readCMstat :: String -> IO (Either ParseError CMstat) -readCMstat filePath = do +readCMstat :: String -> IO (Either ParseError CMstat) +readCMstat filePath = do parsedFile <- parseFromFile genParserCMstat filePath - CE.evaluate parsedFile - + CE.evaluate parsedFile + genParserCMstat :: GenParser Char st CMstat genParserCMstat = do string "# cmstat :: display summary statistics for CMs" newline string "# INFERNAL " many1 (noneOf "\n") - newline + newline string "# Copyright (C) 201" many1 (noneOf "\n") - newline + newline string "# Freely distributed under the GNU General Public License (GPLv3)." - newline + newline string "# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -" newline char '#' @@ -275,27 +275,27 @@ many1 (char '-') newline char '#' - many1 space + many1 space string "idx" - many1 space + many1 space string "name" - many1 space + many1 space string "accession" - many1 space + many1 space string "nseq" - many1 space + many1 space string "eff_nseq" - many1 space + many1 space string "clen" - many1 space + many1 space string "W" - many1 space + many1 space string "bps" - many1 space + many1 space string "bifs" - many1 space + many1 space string "model" - many1 space + many1 space string "cm" many1 space string "hmm" @@ -303,34 +303,34 @@ string "#" many1 (try (oneOf " -")) newline - many1 space + many1 space _statIndex <- many1 digit many1 space _statName <- many1 letter - many1 space + many1 space _statAccession <- many1 (noneOf " ") - many1 space + many1 space _statSequenceNumber <- many1 digit - many1 space + many1 space _statEffectiveSequences <- many1 (oneOf "0123456789.e-") many1 space _statConsensusLength <- many digit - many1 space + many1 space _statW <- many1 digit many1 space _statBasepaires <- many1 digit - many1 space + many1 space _statBifurcations <- many1 digit - many1 space + many1 space _statModel <- many1 letter - many1 space + many1 space _relativeEntropyCM <- many1 (oneOf "0123456789.e-") - many1 space + many1 space _relativeEntropyHMM <- many1 (oneOf "0123456789.e-") newline char '#' newline - eof + eof return $ CMstat (readInt _statIndex) _statName _statAccession (readInt _statSequenceNumber) (readDouble _statEffectiveSequences) (readInt _statConsensusLength) (readInt _statW) (readInt _statBasepaires) (readInt _statBifurcations) _statModel (readDouble _relativeEntropyCM) (readDouble _relativeEntropyHMM) -- readInt :: String -> Int diff --git a/src/Bio/RNAcentralHTTP.hs b/src/Bio/RNAcentralHTTP.hs --- a/src/Bio/RNAcentralHTTP.hs +++ b/src/Bio/RNAcentralHTTP.hs @@ -1,27 +1,27 @@ {-# LANGUAGE OverloadedStrings #-} -{-# LANGUAGE Arrows #-} + {-# LANGUAGE DeriveGeneric #-} -- | Interface for the RNAcentral REST webservice. -- module Bio.RNAcentralHTTP (rnaCentralHTTP, buildSequenceViaMD5Query, - buildStringViaMD5Query, + buildStringViaMD5Query, getRNACentralEntries, showRNAcentralAlienEvaluation, RNAcentralEntryResponse(..), RNAcentralEntry(..) ) where -import Network.HTTP.Conduit -import qualified Data.ByteString.Lazy.Char8 as L8 +import Network.HTTP.Conduit +import qualified Data.ByteString.Lazy.Char8 as L8 import Network import Control.Concurrent import Data.Text import Data.Aeson import GHC.Generics import qualified Data.Digest.Pure.MD5 as M -import Bio.Core.Sequence +import Bio.Core.Sequence import Bio.Sequence.Fasta import Data.Either @@ -40,7 +40,7 @@ toJSON = genericToJSON defaultOptions --toEncoding = genericToEncoding defaultOptions -instance FromJSON RNAcentralEntryResponse +instance FromJSON RNAcentralEntryResponse data RNAcentralEntry = RNAcentralEntry { @@ -58,7 +58,7 @@ toJSON = genericToJSON defaultOptions --toEncoding = genericToEncoding defaultOptions -instance FromJSON RNAcentralEntry +instance FromJSON RNAcentralEntry -- | Send query and parse return XML startSession :: String -> IO (Either String RNAcentralEntryResponse) @@ -68,7 +68,7 @@ --putStr (L8.unpack requestXml) let eitherErrorResponse = eitherDecode requestXml :: Either String RNAcentralEntryResponse return eitherErrorResponse - + -- | Send query and return response XML sendQuery :: String -> IO L8.ByteString sendQuery query' = do @@ -76,10 +76,10 @@ let request = address ++ query' --putStrLn request simpleHttp request - + -- | Function for querying the RNAcentral REST interface. rnaCentralHTTP :: String -> IO (Either String RNAcentralEntryResponse) -rnaCentralHTTP query' = do +rnaCentralHTTP query' = startSession query' -- | Function for delayed queries to the RNAcentral REST interface. Enforces the maximum 20 requests per second policy. @@ -88,17 +88,16 @@ threadDelay 55000 startSession query' -getRNACentralEntries :: [String] -> IO [(Either String RNAcentralEntryResponse)] +getRNACentralEntries :: [String] -> IO [Either String RNAcentralEntryResponse] getRNACentralEntries queries = do - responses <- mapM delayedRNACentralHTTP queries - return responses + mapM delayedRNACentralHTTP queries --Build a query from a input sequence buildSequenceViaMD5Query :: Sequence -> String buildSequenceViaMD5Query s = qString where querySequence = unSD (seqdata s) querySequenceUreplacedwithT = L8.map bsreplaceUT querySequence - querySequenceU2Twolb = L8.filter (\a -> a /= '\n') querySequenceUreplacedwithT + querySequenceU2Twolb = L8.filter ((/= '\n')) querySequenceUreplacedwithT md5Sequence = M.md5 querySequenceU2Twolb qString = "?md5=" ++ show md5Sequence @@ -106,17 +105,17 @@ buildStringViaMD5Query :: String -> String buildStringViaMD5Query s = qString where querySequenceUreplacedwithT = L8.map bsreplaceUT (L8.pack s) - querySequenceU2Twolb = L8.filter (\a -> a /= '\n') querySequenceUreplacedwithT + querySequenceU2Twolb = L8.filter ((/= '\n')) querySequenceUreplacedwithT md5Sequence = M.md5 querySequenceU2Twolb qString = "?md5=" ++ show md5Sequence - -showRNAcentralAlienEvaluation :: [(Either String RNAcentralEntryResponse)] -> String + +showRNAcentralAlienEvaluation :: [Either String RNAcentralEntryResponse] -> String showRNAcentralAlienEvaluation responses = output where resultEntries = Prelude.concatMap results (rights responses) resulthead = "rnacentral_id\tmd5\tlength\n" resultentries = Prelude.concatMap showRNAcentralAlienEvaluationLine resultEntries - output = if resultentries == [] then "No matching sequences found in RNAcentral\n" else resulthead ++ resultentries - + output = if Prelude.null resultentries then "No matching sequences found in RNAcentral\n" else resulthead ++ resultentries + showRNAcentralAlienEvaluationLine :: RNAcentralEntry -> String showRNAcentralAlienEvaluationLine entry = unpack (rnacentral_id entry) ++ "\t" ++ unpack (md5 entry) ++ "\t" ++ show (Bio.RNAcentralHTTP.length entry) ++"\n" diff --git a/src/Bio/RNAcentralHTTPRequest.hs b/src/Bio/RNAcentralHTTPRequest.hs --- a/src/Bio/RNAcentralHTTPRequest.hs +++ b/src/Bio/RNAcentralHTTPRequest.hs @@ -4,23 +4,23 @@ -- | RNAcentralHTTPRequest -- Testcommand: dist/build/RNAcentralHTTPRequest/RNAcentralHTTPRequest -i ATACTTACCTGGCACAGGGGATACCACGATCACCAAGGTGGTTCCCCCAAGACGAGGCTCACCATTGCACTCCGGTGGCGCTGACCCTTGCAATGACCCCAAATGTGGGTTACTCGGGTGTGTAATTTCTGTTAGCTGGGGACTGCGTTCGCGCTTTCCCCTT module Main where - -import System.Console.CmdArgs + +import System.Console.CmdArgs import Bio.RNAcentralHTTP -data Options = Options +data Options = Options { inputSequence :: String } deriving (Show,Data,Typeable) options :: Options options = Options { inputSequence = def &= name "i" &= help "input sequence" - } &= summary "RNAcentralHTTPRequest" &= help "Florian Eggenhofer 2016" &= verbosity - + } &= summary "RNAcentralHTTPRequest" &= help "Florian Eggenhofer 2016" &= verbosity + main :: IO () main = do Options{..} <- cmdArgs options - let query = buildStringViaMD5Query inputSequence + let query = buildStringViaMD5Query inputSequence rnacentralentries <- getRNACentralEntries [query] print rnacentralentries - + diff --git a/src/Bio/RNAlien.hs b/src/Bio/RNAlien.hs --- a/src/Bio/RNAlien.hs +++ b/src/Bio/RNAlien.hs @@ -5,10 +5,10 @@ -- For more information on RNA family models consult -- Testcommand: dist/build/RNAlien/RNAlien -i ~egg/initialfasta/RybB.fa -c 3 -o /scr/kronos/egg/temp/ > ~egg/Desktop/alieninitialtest module Main where - -import System.Console.CmdArgs -import System.Directory -import Bio.Sequence.Fasta + +import System.Console.CmdArgs +import System.Directory +import Bio.Sequence.Fasta import Bio.RNAlienData import Bio.RNAlienLibrary import Data.Maybe @@ -18,8 +18,8 @@ import Paths_RNAlien (version) import Data.Version (showVersion) -data Options = Options - { inputFastaFilePath :: String, +data Options = Options + { inputFastaFilePath :: String, outputPath :: String, inputTaxId :: Maybe Int, inputnSCICutoff :: Maybe Double, @@ -34,17 +34,18 @@ threads :: Int, taxonomyRestriction :: Maybe String, sessionIdentificator :: Maybe String, - performEvaluation :: Bool + performEvaluation :: Bool, + checkSetup :: Bool } deriving (Show,Data,Typeable) options :: Options options = Options - { inputFastaFilePath = def &= name "i" &= help "Path to input fasta file", - outputPath = def &= name "o" &= help "Path to output directory", + { inputFastaFilePath = def &= name "i" &= help "Path to input fasta file", + outputPath = def &= name "o" &= help "Path to output directory. Default: current working directory", inputTaxId = Nothing &= name "t" &= help "NCBI taxonomy ID number of input RNA organism", - inputnSCICutoff = (Just (1 :: Double)) &= name "z" &= help "Only candidate sequences with a normalized structure conservation index (nSCI) higher than this value are accepted. Default: 1", - inputEvalueCutoff = (Just (0.001 :: Double)) &= name "e" &= help "Evalue cutoff for cmsearch filtering. Default: 0.001", - inputBlastDatabase = Just "nt" &= name "b" &= help "Specify name of blast database to use. Default: nt", + inputnSCICutoff = Just (1 :: Double) &= name "z" &= help "Only candidate sequences with a normalized structure conservation index (nSCI) higher than this value are accepted. Default: 1", + inputEvalueCutoff = Just (0.001 :: Double) &= name "e" &= help "Evalue cutoff for cmsearch filtering. Default: 0.001", + inputBlastDatabase = Just "nt" &= name "b" &= help "Specify name of blast database to use. Default: nt", lengthFilter = True &= name "l" &= help "Filter blast hits per genomic length. Default: True", coverageFilter = True &= name "a" &= help "Filter blast hits by coverage of at least 80%. Default: True", singleHitperTax = False &= name "s" &= help "Only the best blast hit per taxonomic entry is considered. Default: False", @@ -54,64 +55,74 @@ threads = 1 &= name "c" &= help "Number of available cpu slots/cores. Default: 1", taxonomyRestriction = Nothing &= name "r" &= help "Restrict search space to taxonomic kingdom (bacteria,archea,eukaryia). Default: not set", sessionIdentificator = Nothing &= name "d" &= help "Optional session id that is used instead of automatically generated one.", - performEvaluation = True &= name "x" &= help "Perform evaluation step. Default: True" - } &= summary ("RNAlien " ++ alienVersion) &= help "Florian Eggenhofer, Ivo L. Hofacker, Christian Höner zu Siederdissen - 2013 - 2017" &= verbosity - + performEvaluation = True &= name "x" &= help "Perform evaluation step. Default: True", + checkSetup = False &= name "g" &= help "Just prints installed tool versions and performs connection check. Default: False" + } &= summary ("RNAlien " ++ alienVersion) &= help "Florian Eggenhofer, Ivo L. Hofacker, Christian Hoener zu Siederdissen - 2013 - 2017" &= verbosity + main :: IO () main = do Options{..} <- cmdArgs options verboseLevel <- getVerbosity + let tools = if inputQuerySelectionMethod == "clustering" then ["clustalo","mlocarna","RNAfold","RNAalifold","cmcalibrate","cmstat","cmbuild","RNAz","RNAcode"] else ["mlocarna","RNAfold","RNAalifold","cmcalibrate","cmstat","cmbuild","RNAz","RNAcode"] -- Generate SessionID sessionId <- createSessionID sessionIdentificator timestamp <- getCurrentTime - let temporaryDirectoryPath = FP.addTrailingPathSeparator outputPath ++ sessionId ++ "/" + currentWorkDirectory <- getCurrentDirectory + let selectedOutputPath = if null outputPath then currentWorkDirectory else outputPath + let temporaryDirectoryPath = FP.addTrailingPathSeparator selectedOutputPath ++ sessionId ++ "/" createDirectoryIfMissing False temporaryDirectoryPath - createDirectoryIfMissing False (temporaryDirectoryPath ++ "log") - -- Create Log files - writeFile (temporaryDirectoryPath ++ "Log") ("RNAlien " ++ alienVersion ++ "\n") - writeFile (temporaryDirectoryPath ++ "log/warnings") ("") - logMessage ("Timestamp: " ++ (show timestamp) ++ "\n") temporaryDirectoryPath - logMessage ("Temporary Directory: " ++ temporaryDirectoryPath ++ "\n") temporaryDirectoryPath - inputFasta <- readFasta inputFastaFilePath networkCheck <- checkNCBIConnection - if isLeft networkCheck + if checkSetup then do - putStrLn ("Error - Could not contact NCBI server: " ++ fromLeft networkCheck ++ "\n") - logMessage ("Error - Could not contact NCBI server: " ++ fromLeft networkCheck ++ "\n") temporaryDirectoryPath - else do - if null inputFasta + toolsCheck <- checkTools tools inputQuerySelectionMethod temporaryDirectoryPath + let setupCheckPath = temporaryDirectoryPath ++ "setupCheck" + let toolCheckResult = either id id toolsCheck + let networkCheckResult = either id id networkCheck + writeFile setupCheckPath (toolCheckResult ++ "\n" ++ networkCheckResult ++ "\n") + else + if isLeft networkCheck then do - putStrLn "Error: Input fasta file is empty." - logMessage "Error: Input fasta file is empty.\n" temporaryDirectoryPath - else do - let iterationNumber = 0 - let tools = if inputQuerySelectionMethod == "clustering" then ["clustalo","mlocarna","RNAfold","RNAalifold","cmcalibrate","cmstat","cmbuild","RNAz","RNAcode"] else ["mlocarna","RNAfold","RNAalifold","cmcalibrate","cmstat","cmbuild","RNAz","RNAcode"] - toolsCheck <- checkTools tools inputQuerySelectionMethod temporaryDirectoryPath - -- Check required commandline tools - if isLeft toolsCheck - then do - putStrLn ("Error - Not all required tools could be found in $PATH: " ++ fromLeft toolsCheck ++ "\n") - logMessage ("Error - Not all required tools could be found in $PATH: " ++ fromLeft toolsCheck ++ "\n") temporaryDirectoryPath - else do - logToolVersions inputQuerySelectionMethod temporaryDirectoryPath - let inputSequence = reformatFasta (head inputFasta) - initialTaxId <- setInitialTaxId inputBlastDatabase temporaryDirectoryPath inputTaxId inputSequence - let checkedTaxonomyRestriction = checkTaxonomyRestriction taxonomyRestriction - let staticOptions = StaticOptions temporaryDirectoryPath sessionId (fromJust inputnSCICutoff) inputTaxId singleHitperTax inputQuerySelectionMethod inputQueryNumber lengthFilter coverageFilter blastSoftmasking threads inputBlastDatabase checkedTaxonomyRestriction (setVerbose verboseLevel) - let initialization = ModelConstruction iterationNumber inputSequence [] initialTaxId Nothing (fromJust inputEvalueCutoff) False [] [] - logMessage (show initialization) temporaryDirectoryPath - modelConstructionResults <- modelConstructer staticOptions initialization - let resultTaxonomyRecordsCSVTable = constructTaxonomyRecordsCSVTable modelConstructionResults - writeFile (temporaryDirectoryPath ++ "result.csv") resultTaxonomyRecordsCSVTable - if performEvaluation - then do - resultEvaluation <- evaluateConstructionResult staticOptions modelConstructionResults - appendFile (temporaryDirectoryPath ++ "Log") resultEvaluation - resultSummary modelConstructionResults staticOptions - writeFile (temporaryDirectoryPath ++ "done") "" - else do - resultSummary modelConstructionResults staticOptions - writeFile (temporaryDirectoryPath ++ "done") "" + putStrLn ("Error - Could not contact NCBI server: " ++ fromLeft networkCheck ++ "\n") + logMessage ("Error - Could not contact NCBI server: " ++ fromLeft networkCheck ++ "\n") temporaryDirectoryPath + else do + createDirectoryIfMissing False (temporaryDirectoryPath ++ "log") + -- Create Log files + writeFile (temporaryDirectoryPath ++ "Log") ("RNAlien " ++ alienVersion ++ "\n") + writeFile (temporaryDirectoryPath ++ "log/warnings") ("") + logMessage ("Timestamp: " ++ (show timestamp) ++ "\n") temporaryDirectoryPath + logMessage ("Temporary Directory: " ++ temporaryDirectoryPath ++ "\n") temporaryDirectoryPath + inputFasta <- readFasta inputFastaFilePath + if null inputFasta + then do + putStrLn "Error: Input fasta file is empty." + logMessage "Error: Input fasta file is empty.\n" temporaryDirectoryPath + else do + let iterationNumber = 0 + toolsCheck <- checkTools tools inputQuerySelectionMethod temporaryDirectoryPath + if isLeft toolsCheck + then do + putStrLn ("Error - Not all required tools could be found in $PATH: " ++ fromLeft toolsCheck ++ "\n") + logMessage ("Error - Not all required tools could be found in $PATH: " ++ fromLeft toolsCheck ++ "\n") temporaryDirectoryPath + else do + logToolVersions inputQuerySelectionMethod temporaryDirectoryPath + let inputSequence = reformatFasta (head inputFasta) + initialTaxId <- setInitialTaxId inputBlastDatabase temporaryDirectoryPath inputTaxId inputSequence + let checkedTaxonomyRestriction = checkTaxonomyRestriction taxonomyRestriction + let staticOptions = StaticOptions temporaryDirectoryPath sessionId (fromJust inputnSCICutoff) inputTaxId singleHitperTax inputQuerySelectionMethod inputQueryNumber lengthFilter coverageFilter blastSoftmasking threads inputBlastDatabase checkedTaxonomyRestriction (setVerbose verboseLevel) + let initialization = ModelConstruction iterationNumber inputSequence [] initialTaxId Nothing (fromJust inputEvalueCutoff) False [] [] + logMessage (show initialization) temporaryDirectoryPath + modelConstructionResults <- modelConstructer staticOptions initialization + let resultTaxonomyRecordsCSVTable = constructTaxonomyRecordsCSVTable modelConstructionResults + writeFile (temporaryDirectoryPath ++ "result.csv") resultTaxonomyRecordsCSVTable + if performEvaluation + then do + resultEvaluation <- evaluateConstructionResult staticOptions modelConstructionResults + appendFile (temporaryDirectoryPath ++ "Log") resultEvaluation + resultSummary modelConstructionResults staticOptions + writeFile (temporaryDirectoryPath ++ "done") "" + else do + resultSummary modelConstructionResults staticOptions + writeFile (temporaryDirectoryPath ++ "done") "" alienVersion :: String alienVersion = showVersion version diff --git a/src/Bio/RNAlienData.hs b/src/Bio/RNAlienData.hs --- a/src/Bio/RNAlienData.hs +++ b/src/Bio/RNAlienData.hs @@ -4,7 +4,7 @@ import qualified Data.ByteString.Lazy.Char8 as L import Bio.Core.Sequence -import Bio.Sequence.Fasta +import Bio.Sequence.Fasta import Bio.Taxonomy -- | Static construction options @@ -23,26 +23,26 @@ blastDatabase :: Maybe String, taxRestriction :: Maybe String, verbositySwitch :: Bool - } deriving (Show) + } deriving (Show) -- | Keeps track of model construction data ModelConstruction = ModelConstruction { iterationNumber :: Int, - inputFasta :: Sequence, + inputFasta :: Sequence, taxRecords :: [TaxonomyRecord], --Taxonomy ID of the highest node in taxonomic subtree used in search upperTaxonomyLimit :: Maybe Int, taxonomicContext :: Maybe Taxon, - evalueThreshold :: Double, + evalueThreshold :: Double, alignmentModeInfernal :: Bool, selectedQueries :: [Sequence], potentialMembers :: [SearchResult] - } + } instance Show ModelConstruction where show (ModelConstruction _iterationNumber _inputFasta _taxRecords _upperTaxonomyLimit _taxonomicContext _evalueThreshold _alignmentModeInfernal _selectedQueries _potentialMembers) = a ++ b ++ c ++ d ++ e ++ g ++ h ++ i - where a = "Modelconstruction iteration: " ++ show _iterationNumber ++ "\n" - b = "Input fasta:\n" ++ L.unpack (unSL (seqheader _inputFasta)) ++ "\n" ++ L.unpack (unSD (seqdata _inputFasta)) ++ "\n" + where a = "Modelconstruction iteration: " ++ show _iterationNumber ++ "\n" + b = "Input fasta:\n" ++ L.unpack (unSL (seqheader _inputFasta)) ++ "\n" ++ L.unpack (unSD (seqdata _inputFasta)) ++ "\n" c = show _taxRecords d = "Upper taxonomy limit: " ++ maybe "not set" show _upperTaxonomyLimit ++ "\n" e = "Taxonomic Context: " ++ maybe "not set" show _taxonomicContext ++ "\n" @@ -57,7 +57,7 @@ instance Show TaxonomyRecord where show (TaxonomyRecord _recordTaxonomyId _sequenceRecords) = a ++ b - where a = "TaxonomyRecord TaxonomyId: " ++ show _recordTaxonomyId ++ "\n" + where a = "TaxonomyRecord TaxonomyId: " ++ show _recordTaxonomyId ++ "\n" b = show _sequenceRecords data SequenceRecord = SequenceRecord @@ -65,14 +65,14 @@ nucleotideSequence :: Sequence, -- 0 is unaligned, number is the iteration the sequence has been included into the alignment aligned :: Int, - recordDescription :: L.ByteString - } + recordDescription :: L.ByteString + } instance Show SequenceRecord where - show (SequenceRecord _nucleotideSequence _aligned _recordDescription) = a ++ b ++ c - where a = "Record Description: " ++ (L.unpack _recordDescription) ++ "\n" - b = "Aligned in iteration: " ++ show _aligned ++ "\n" - c = "Sequence:" ++ show _nucleotideSequence ++ "\n" + show (SequenceRecord _nucleotideSequence _aligned _recordDescription) = a ++ b ++ c + where a = "Record Description: " ++ L.unpack _recordDescription ++ "\n" + b = "Aligned in iteration: " ++ show _aligned ++ "\n" + c = "Sequence:" ++ show _nucleotideSequence ++ "\n" -- | data CMsearch = CMsearch { queryCMfile :: String, @@ -81,7 +81,7 @@ cmsearchHits :: [CMsearchHit] -- hitAlignments :: [CMsearchHitAlignment] -- internalCMPipelineStatisticsSummary - } deriving (Show, Eq, Read) + } deriving (Show, Eq, Read) -- | data CMsearchHit = CMsearchHit @@ -98,7 +98,7 @@ hitTruncation :: L.ByteString, hitGCContent :: Double, hitDescription :: L.ByteString - } deriving (Show, Eq, Read) + } deriving (Show, Eq, Read) data SearchResult = SearchResult { candidates :: [(Sequence,Int,L.ByteString)], @@ -106,9 +106,9 @@ } instance Show SearchResult where - show (SearchResult _candidates _blastDatabaseSize) = a ++ b - where a = "SearchResults :\n " ++ concatMap show _candidates ++ "\n" - b = "BlastDb Size: " ++ show _blastDatabaseSize ++ "\n" + show (SearchResult _candidates _blastDatabaseSize) = a ++ b + where a = "SearchResults :\n " ++ concatMap show _candidates ++ "\n" + b = "BlastDb Size: " ++ show _blastDatabaseSize ++ "\n" -- | data CMstat = CMstat @@ -125,12 +125,12 @@ statModel :: String, relativeEntropyCM :: Double, relativeEntropyHMM :: Double - } deriving (Eq, Read) + } deriving (Eq, Read) instance Show CMstat where show (CMstat _statIndex _statName _statAccession _statSequenceNumber _statEffectiveSequences _statConsensusLength _statW _statBasepairs _statBifurcations _statModel _relativeEntropyCM _relativeEntropyHMM) = a ++ b ++ c ++ d ++ e ++ f ++ g ++ h ++ i ++ j ++ k ++ l - where a = "CMstat - covariance model statistics:\nIndex: " ++ show _statIndex ++ "\n" - b = "Name: " ++ show _statName ++ "\n" + where a = "CMstat - covariance model statistics:\nIndex: " ++ show _statIndex ++ "\n" + b = "Name: " ++ show _statName ++ "\n" c = "Accession: " ++ show _statAccession ++ "\n" d = "Sequence Number: " ++ show _statSequenceNumber ++ "\n" e = "Effective Sequences: " ++ show _statEffectiveSequences ++ "\n" @@ -141,4 +141,4 @@ j = "Modeltype: " ++ show _statModel ++ "\n" k = "Relative Entropy CM: " ++ show _relativeEntropyCM ++ "\n" l = "Relative Entropy HMM: " ++ show _relativeEntropyHMM ++ "\n" - + diff --git a/src/Bio/RNAlienLibrary.hs b/src/Bio/RNAlienLibrary.hs --- a/src/Bio/RNAlienLibrary.hs +++ b/src/Bio/RNAlienLibrary.hs @@ -31,33 +31,33 @@ evaluePartitionTrimCMsearchHits ) where - -import System.Process + +import System.Process import qualified System.FilePath as FP -import Text.ParserCombinators.Parsec +import Text.ParserCombinators.Parsec import Data.List import Data.Char -import Bio.Core.Sequence -import Bio.Sequence.Fasta +import Bio.Core.Sequence +import Bio.Sequence.Fasta import Bio.BlastXML import Bio.ClustalParser import Data.Int (Int16) import Bio.RNAlienData import qualified Data.ByteString.Lazy.Char8 as L import qualified Data.ByteString.Char8 as B -import Bio.Taxonomy +import Bio.Taxonomy import Data.Either.Unwrap import Data.Maybe -import Bio.EntrezHTTP +import Bio.EntrezHTTP import System.Exit import Data.Either (lefts,rights,Either) -import qualified Text.EditDistance as ED +import qualified Text.EditDistance as ED import qualified Data.Vector as V -import Control.Concurrent +import Control.Concurrent import System.Random import Data.Csv import Data.Matrix -import Bio.BlastHTTP +import Bio.BlastHTTP import Data.Clustering.Hierarchical import System.Directory import System.Console.CmdArgs @@ -78,55 +78,57 @@ import qualified Data.Text.Lazy.IO as TIO import Text.Printf import qualified Data.Text.Metrics as TM +import Control.Monad +import Control.Arrow -- | Initial RNA family model construction - generates iteration number, seed alignment and model modelConstructer :: StaticOptions -> ModelConstruction -> IO ModelConstruction modelConstructer staticOptions modelConstruction = do logMessage ("Iteration: " ++ show (iterationNumber modelConstruction) ++ "\n") (tempDirPath staticOptions) iterationSummary modelConstruction staticOptions - let currentIterationNumber = (iterationNumber modelConstruction) + let currentIterationNumber = iterationNumber modelConstruction let foundSequenceNumber = length (concatMap sequenceRecords (taxRecords modelConstruction)) --extract queries let queries = extractQueries foundSequenceNumber modelConstruction logVerboseMessage (verbositySwitch staticOptions) ("Queries:" ++ show queries ++ "\n") (tempDirPath staticOptions) - let iterationDirectory = (tempDirPath staticOptions) ++ (show currentIterationNumber) ++ "/" + let iterationDirectory = tempDirPath staticOptions ++ show currentIterationNumber ++ "/" let maybeLastTaxId = extractLastTaxId (taxonomicContext modelConstruction) - if (isNothing maybeLastTaxId) then logMessage ("Lineage: Could not extract last tax id \n") (tempDirPath staticOptions) else (return ()) + Control.Monad.when (isNothing maybeLastTaxId) $ logMessage ("Lineage: Could not extract last tax id \n") (tempDirPath staticOptions) --If highest node in linage was used as upper taxonomy limit, taxonomic tree is exhausted - if (maybe True (\uppertaxlimit -> maybe True (\lastTaxId -> uppertaxlimit /= lastTaxId) maybeLastTaxId) (upperTaxonomyLimit modelConstruction)) + if maybe True (\uppertaxlimit -> maybe True (\lastTaxId -> uppertaxlimit /= lastTaxId) maybeLastTaxId) (upperTaxonomyLimit modelConstruction) then do - createDirectory (iterationDirectory) + createDirectory iterationDirectory let (upperTaxLimit,lowerTaxLimit) = setTaxonomicContextEntrez currentIterationNumber (taxonomicContext modelConstruction) (upperTaxonomyLimit modelConstruction) - logVerboseMessage (verbositySwitch staticOptions) ("Upper taxonomy limit: " ++ (show upperTaxLimit) ++ "\n " ++ "Lower taxonomy limit: "++ show lowerTaxLimit ++ "\n") (tempDirPath staticOptions) + logVerboseMessage (verbositySwitch staticOptions) ("Upper taxonomy limit: " ++ show upperTaxLimit ++ "\n " ++ "Lower taxonomy limit: "++ show lowerTaxLimit ++ "\n") (tempDirPath staticOptions) --search queries let expectThreshold = setBlastExpectThreshold modelConstruction - searchResults <- catchAll (searchCandidates staticOptions Nothing currentIterationNumber upperTaxLimit lowerTaxLimit expectThreshold queries) + searchResults <- catchAll (searchCandidates staticOptions Nothing currentIterationNumber upperTaxLimit lowerTaxLimit expectThreshold queries) (\e -> do logWarning ("Warning: Search results iteration" ++ show (iterationNumber modelConstruction) ++ " - exception: " ++ show e) (tempDirPath staticOptions) return (SearchResult [] Nothing)) currentTaxonomicContext <- getTaxonomicContextEntrez upperTaxLimit (taxonomicContext modelConstruction) if null (candidates searchResults) - then do + then alignmentConstructionWithoutCandidates currentTaxonomicContext upperTaxLimit staticOptions modelConstruction - else do + else alignmentConstructionWithCandidates currentTaxonomicContext upperTaxLimit searchResults staticOptions modelConstruction else do - logMessage ("Message: Modelconstruction complete: Out of queries or taxonomic tree exhausted\n") (tempDirPath staticOptions) + logMessage "Message: Modelconstruction complete: Out of queries or taxonomic tree exhausted\n" (tempDirPath staticOptions) modelConstructionResult staticOptions modelConstruction catchAll :: IO a -> (CE.SomeException -> IO a) -> IO a catchAll = CE.catch setInitialTaxId :: Maybe String -> String -> Maybe Int -> Sequence -> IO (Maybe Int) -setInitialTaxId inputBlastDatabase tempdir inputTaxId inputSequence = do +setInitialTaxId inputBlastDatabase tempdir inputTaxId inputSequence = if (isNothing inputTaxId) then do initialTaxId <- findTaxonomyStart inputBlastDatabase tempdir inputSequence return (Just initialTaxId) - else do + else do return inputTaxId extractLastTaxId :: Maybe Taxon -> Maybe Int -extractLastTaxId taxon +extractLastTaxId taxon | isNothing taxon = Nothing | V.null lineageExVector = Nothing | otherwise = Just (lineageTaxId (V.head lineageExVector)) @@ -144,8 +146,8 @@ let queries = extractQueries foundSequenceNumber modelConstruction --- --let alignedSequences' = map nucleotideSequence (concatMap sequenceRecords (taxRecords modelConstruction)) --- logVerboseMessage (verbositySwitch staticOptions) ("Queries:" ++ show queries ++ "\n") outputDirectory - let iterationDirectory = outputDirectory ++ (show currentIterationNumber) ++ "/" - createDirectory (iterationDirectory) + let iterationDirectory = outputDirectory ++ show currentIterationNumber ++ "/" + createDirectory iterationDirectory let logFileDirectoryPath = iterationDirectory ++ "log" createDirectoryIfMissing False logFileDirectoryPath let expectThreshold = setBlastExpectThreshold modelConstruction @@ -168,7 +170,7 @@ candidates1 <- catchAll (searchCandidates staticOptions (Just "archea") currentIterationNumber upperTaxLimit1 lowerTaxLimit1 expectThreshold queries) (\e -> do logWarning ("Warning: Search results iteration" ++ show currentIterationNumber ++ " - exception: " ++ show e) outputDirectory return (SearchResult [] Nothing)) - let uniqueCandidates1 = filterDuplicates modelConstruction candidates1 + let uniqueCandidates1 = filterDuplicates modelConstruction candidates1 (alignmentResults1,potentialMembers1) <- catchAll (alignCandidates staticOptions modelConstruction "archea" uniqueCandidates1) (\e -> do logWarning ("Warning: Alignment results iteration" ++ show currentIterationNumber ++ " - exception: " ++ show e) outputDirectory return ([],[])) @@ -198,10 +200,10 @@ let preliminaryAlignmentPath = iterationDirectory ++ "model.stockholm" let preliminaryCMLogPath = iterationDirectory ++ "model.cm.log" let nextModelConstructionInput = constructNext currentIterationNumber modelConstruction alignmentResults Nothing Nothing [] currentPotentialMembers (alignmentModeInfernal modelConstruction) - if (length alignmentResults == 0) && (not (alignmentModeInfernal modelConstruction)) + if (null alignmentResults) && not (alignmentModeInfernal modelConstruction) then do - logVerboseMessage (verbositySwitch staticOptions) ("Alignment result initial mode\n") outputDirectory - logMessage ("Message: No sequences found that statisfy filters. Try to reconstruct model with less strict cutoff parameters.") outputDirectory + logVerboseMessage (verbositySwitch staticOptions) "Alignment result initial mode\n" outputDirectory + logMessage "Message: No sequences found that statisfy filters. Try to reconstruct model with less strict cutoff parameters." outputDirectory let alignedSequences = extractAlignedSequences (iterationNumber modelConstruction) modelConstruction let alignmentSequences = map snd (V.toList (V.concat [alignedSequences])) writeFasta preliminaryFastaPath alignmentSequences @@ -216,9 +218,8 @@ writeFile preliminaryAlignmentPath stockholAlignment _ <- systemCMbuild cmBuildOptions preliminaryAlignmentPath preliminaryCMPath cmBuildFilepath _ <- systemCMcalibrate "fast" (cpuThreads staticOptions) preliminaryCMPath preliminaryCMLogPath - resultModelConstruction <- reevaluatePotentialMembers staticOptions nextModelConstructionInput - return resultModelConstruction - else do + reevaluatePotentialMembers staticOptions nextModelConstructionInput + else if (alignmentModeInfernal modelConstruction) then do logVerboseMessage (verbositySwitch staticOptions) ("Alignment construction with candidates - infernal mode\n") outputDirectory @@ -245,8 +246,8 @@ let outputDirectory = tempDirPath staticOptions iterationSummary modelConstruction staticOptions logMessage ("Reevaluation of potential members iteration: " ++ show currentIterationNumber ++ "\n") outputDirectory - let iterationDirectory = outputDirectory ++ (show currentIterationNumber) ++ "/" - createDirectory (iterationDirectory) + let iterationDirectory = outputDirectory ++ show currentIterationNumber ++ "/" + createDirectory iterationDirectory let indexedPotentialMembers = V.indexed (V.fromList (potentialMembers modelConstruction)) potentialMembersAlignmentResultVector <- V.mapM (\(i,searchresult) -> (alignCandidates staticOptions modelConstruction (show i ++ "_") searchresult)) indexedPotentialMembers let potentialMembersAlignmentResults = V.toList potentialMembersAlignmentResultVector @@ -256,6 +257,7 @@ let resultFastaPath = outputDirectory ++ "result.fa" let resultCMPath = outputDirectory ++ "result.cm" let resultAlignmentPath = outputDirectory ++ "result.stockholm" + let resultClustalFilepath = outputDirectory ++ "result.clustal" let resultCMLogPath = outputDirectory ++ "log/result.cm.log" if null alignmentResults then do @@ -267,13 +269,14 @@ copyFile lastIterationFastaPath resultFastaPath copyFile lastIterationAlignmentPath resultAlignmentPath _ <- systemCMcalibrate "standard" (cpuThreads staticOptions) resultCMPath resultCMLogPath + systemCMalign ("--outformat=Clustal --cpu " ++ show (cpuThreads staticOptions)) resultCMPath resultFastaPath resultClustalFilepath writeFile (iterationDirectory ++ "done") "" return modelConstruction else do let lastIterationFastaPath = outputDirectory ++ show currentIterationNumber ++ "/model.fa" let lastIterationAlignmentPath = outputDirectory ++ show currentIterationNumber ++ "/model.stockholm" let lastIterationCMPath = outputDirectory ++ show currentIterationNumber ++ "/model.cm" - logVerboseMessage (verbositySwitch staticOptions) ("Alignment construction with candidates - infernal mode\n") outputDirectory + logVerboseMessage (verbositySwitch staticOptions) "Alignment construction with candidates - infernal mode\n" outputDirectory let nextModelConstructionInput = constructNext currentIterationNumber modelConstruction alignmentResults Nothing Nothing [] [] (alignmentModeInfernal modelConstruction) constructModel nextModelConstructionInput staticOptions copyFile lastIterationCMPath resultCMPath @@ -284,42 +287,42 @@ logMessage (iterationSummaryLog nextModelConstructionInput) outputDirectory logVerboseMessage (verbositySwitch staticOptions) (show nextModelConstructionInput) outputDirectory _ <- systemCMcalibrate "standard" (cpuThreads staticOptions) resultCMPath resultCMLogPath + systemCMalign ("--outformat=Clustal --cpu " ++ show (cpuThreads staticOptions)) resultCMPath resultFastaPath resultClustalFilepath writeFile (iterationDirectory ++ "done") "" return nextModelConstructionInput - + --------------------------------------------------------- - + alignmentConstructionWithCandidates :: Maybe Taxon -> Maybe Int -> SearchResult -> StaticOptions -> ModelConstruction -> IO ModelConstruction alignmentConstructionWithCandidates currentTaxonomicContext currentUpperTaxonomyLimit searchResults staticOptions modelConstruction = do --candidates usedUpperTaxonomyLimit blastDatabaseSize - let currentIterationNumber = (iterationNumber modelConstruction) - let iterationDirectory = (tempDirPath staticOptions) ++ (show currentIterationNumber) ++ "/" + let currentIterationNumber = iterationNumber modelConstruction + let iterationDirectory = tempDirPath staticOptions ++ show currentIterationNumber ++ "/" --let usedUpperTaxonomyLimit = (snd (head candidates)) --align search result (alignmentResults,potentialMemberEntries) <- catchAll (alignCandidates staticOptions modelConstruction "" searchResults) (\e -> do logWarning ("Warning: Alignment results iteration" ++ show (iterationNumber modelConstruction) ++ " - exception: " ++ show e) (tempDirPath staticOptions) return ([],[])) let currentPotentialMembers = [SearchResult potentialMemberEntries (blastDatabaseSize searchResults)] - if (length alignmentResults == 0) && (not (alignmentModeInfernal modelConstruction)) + if (null alignmentResults) && not (alignmentModeInfernal modelConstruction) then do logVerboseMessage (verbositySwitch staticOptions) ("Alignment construction with candidates - length 1 - inital mode" ++ "\n") (tempDirPath staticOptions) --too few sequences for alignment. because of lack in sequences no cm was constructed before --reusing previous modelconstruction with increased upperTaxonomyLimit but include found sequence --prepare next iteration - let newTaxEntries = (taxRecords modelConstruction) ++ (buildTaxRecords alignmentResults currentIterationNumber) - let nextModelConstructionInputWithThreshold = modelConstruction {iterationNumber = (currentIterationNumber + 1),upperTaxonomyLimit = currentUpperTaxonomyLimit, taxRecords = newTaxEntries,taxonomicContext = currentTaxonomicContext} + let newTaxEntries = taxRecords modelConstruction ++ buildTaxRecords alignmentResults currentIterationNumber + let nextModelConstructionInputWithThreshold = modelConstruction {iterationNumber = currentIterationNumber + 1,upperTaxonomyLimit = currentUpperTaxonomyLimit, taxRecords = newTaxEntries,taxonomicContext = currentTaxonomicContext} logMessage (iterationSummaryLog nextModelConstructionInputWithThreshold) (tempDirPath staticOptions) logVerboseMessage (verbositySwitch staticOptions) (show nextModelConstructionInputWithThreshold) (tempDirPath staticOptions) ---- writeFile (iterationDirectory ++ "done") "" - nextModelConstruction <- modelConstructer staticOptions nextModelConstructionInputWithThreshold - return nextModelConstruction - else do + modelConstructer staticOptions nextModelConstructionInputWithThreshold + else if (alignmentModeInfernal modelConstruction) then do logVerboseMessage (verbositySwitch staticOptions) ("Alignment construction with candidates - infernal mode\n") (tempDirPath staticOptions) --prepare next iteration - let nextModelConstructionInput = constructNext currentIterationNumber modelConstruction alignmentResults currentUpperTaxonomyLimit currentTaxonomicContext [] currentPotentialMembers True - constructModel nextModelConstructionInput staticOptions + let nextModelConstructionInput = constructNext currentIterationNumber modelConstruction alignmentResults currentUpperTaxonomyLimit currentTaxonomicContext [] currentPotentialMembers True + constructModel nextModelConstructionInput staticOptions writeFile (iterationDirectory ++ "done") "" logMessage (iterationSummaryLog nextModelConstructionInput) (tempDirPath staticOptions) logVerboseMessage (verbositySwitch staticOptions) (show nextModelConstructionInput) (tempDirPath staticOptions) @@ -332,7 +335,7 @@ logVerboseMessage (verbositySwitch staticOptions) ("Alignment construction with candidates - initial mode\n") (tempDirPath staticOptions) --First round enough candidates are available for modelconstruction, alignmentModeInfernal is set to true after this iteration --prepare next iteration - let nextModelConstructionInput = constructNext currentIterationNumber modelConstruction alignmentResults currentUpperTaxonomyLimit currentTaxonomicContext [] currentPotentialMembers False + let nextModelConstructionInput = constructNext currentIterationNumber modelConstruction alignmentResults currentUpperTaxonomyLimit currentTaxonomicContext [] currentPotentialMembers False constructModel nextModelConstructionInput staticOptions currentSelectedQueries <- selectQueries staticOptions modelConstruction alignmentResults --select queries @@ -340,104 +343,101 @@ logMessage (iterationSummaryLog nextModelConstructionInputWithInfernalMode) (tempDirPath staticOptions) logVerboseMessage (verbositySwitch staticOptions) (show nextModelConstructionInputWithInfernalMode) (tempDirPath staticOptions) writeFile (iterationDirectory ++ "done") "" - nextModelConstruction <- modelConstructer staticOptions nextModelConstructionInputWithInfernalMode + nextModelConstruction <- modelConstructer staticOptions nextModelConstructionInputWithInfernalMode return nextModelConstruction - + alignmentConstructionWithoutCandidates :: Maybe Taxon -> Maybe Int -> StaticOptions -> ModelConstruction -> IO ModelConstruction alignmentConstructionWithoutCandidates currentTaxonomicContext upperTaxLimit staticOptions modelConstruction = do - let currentIterationNumber = (iterationNumber modelConstruction) - let iterationDirectory = (tempDirPath staticOptions) ++ (show currentIterationNumber) ++ "/" + let currentIterationNumber = iterationNumber modelConstruction + let iterationDirectory = tempDirPath staticOptions ++ show currentIterationNumber ++ "/" --Found no new candidates in this iteration, reusing previous modelconstruction with increased upperTaxonomyLimit - let nextModelConstructionInputWithThreshold = modelConstruction {iterationNumber = (currentIterationNumber + 1),upperTaxonomyLimit = upperTaxLimit,taxonomicContext = currentTaxonomicContext} + let nextModelConstructionInputWithThreshold = modelConstruction {iterationNumber = currentIterationNumber + 1,upperTaxonomyLimit = upperTaxLimit,taxonomicContext = currentTaxonomicContext} --copy model and alignment from last iteration in place if present - let previousIterationCMPath = (tempDirPath staticOptions) ++ (show (currentIterationNumber - 1)) ++ "/model.cm" + let previousIterationCMPath = tempDirPath staticOptions ++ show (currentIterationNumber - 1) ++ "/model.cm" previousIterationCMexists <- doesFileExist previousIterationCMPath if previousIterationCMexists then do - logVerboseMessage (verbositySwitch staticOptions) ("Alignment construction no candidates - previous cm\n") (tempDirPath staticOptions) - let previousIterationFastaPath = (tempDirPath staticOptions) ++ (show (currentIterationNumber - 1)) ++ "/model.fa" - let previousIterationAlignmentPath = (tempDirPath staticOptions) ++ (show (currentIterationNumber - 1)) ++ "/model.stockholm" - let thisIterationFastaPath = (tempDirPath staticOptions) ++ (show (currentIterationNumber)) ++ "/model.fa" - let thisIterationAlignmentPath = (tempDirPath staticOptions) ++ (show (currentIterationNumber)) ++ "/model.stockholm" - let thisIterationCMPath = (tempDirPath staticOptions) ++ (show (currentIterationNumber)) ++ "/model.cm" + logVerboseMessage (verbositySwitch staticOptions) "Alignment construction no candidates - previous cm\n" (tempDirPath staticOptions) + let previousIterationFastaPath = tempDirPath staticOptions ++ show (currentIterationNumber - 1) ++ "/model.fa" + let previousIterationAlignmentPath = tempDirPath staticOptions ++ show (currentIterationNumber - 1) ++ "/model.stockholm" + let thisIterationFastaPath = tempDirPath staticOptions ++ show (currentIterationNumber) ++ "/model.fa" + let thisIterationAlignmentPath = tempDirPath staticOptions ++ show (currentIterationNumber) ++ "/model.stockholm" + let thisIterationCMPath = tempDirPath staticOptions ++ show (currentIterationNumber) ++ "/model.cm" copyFile previousIterationFastaPath thisIterationFastaPath copyFile previousIterationAlignmentPath thisIterationAlignmentPath copyFile previousIterationCMPath thisIterationCMPath logMessage (iterationSummaryLog nextModelConstructionInputWithThreshold) (tempDirPath staticOptions) logVerboseMessage (verbositySwitch staticOptions) (show nextModelConstructionInputWithThreshold) (tempDirPath staticOptions) writeFile (iterationDirectory ++ "done") "" - nextModelConstruction <- modelConstructer staticOptions nextModelConstructionInputWithThreshold - return nextModelConstruction + modelConstructer staticOptions nextModelConstructionInputWithThreshold else do - logVerboseMessage (verbositySwitch staticOptions) ("Alignment construction no candidates - no previous iteration cm\n") (tempDirPath staticOptions) + logVerboseMessage (verbositySwitch staticOptions) "Alignment construction no candidates - no previous iteration cm\n" (tempDirPath staticOptions) logMessage (iterationSummaryLog nextModelConstructionInputWithThreshold) (tempDirPath staticOptions) logVerboseMessage (verbositySwitch staticOptions) (show nextModelConstructionInputWithThreshold) (tempDirPath staticOptions) ---- writeFile (iterationDirectory ++ "done") "" - nextModelConstruction <- modelConstructer staticOptions nextModelConstructionInputWithThreshold - return nextModelConstruction - + modelConstructer staticOptions nextModelConstructionInputWithThreshold + findTaxonomyStart :: Maybe String -> String -> Sequence -> IO Int findTaxonomyStart inputBlastDatabase temporaryDirectory querySequence = do let queryIndexString = "1" - let hitNumberQuery = buildHitNumberQuery "&HITLIST_SIZE=10" + let hitNumberQuery = buildHitNumberQuery "&HITLIST_SIZE=10" let registrationInfo = buildRegistration "RNAlien" "florian.eggenhofer@univie.ac.at" let blastQuery = BlastHTTPQuery (Just "ncbi") (Just "blastn") inputBlastDatabase [querySequence] (Just (hitNumberQuery ++ registrationInfo)) (Just (5400000000 :: Int)) - logMessage ("No tax id provided - Sending find taxonomy start blast query \n") temporaryDirectory + logMessage "No tax id provided - Sending find taxonomy start blast query \n" temporaryDirectory blastOutput <- CE.catch (blastHTTP blastQuery) (\e -> do let err = show (e :: CE.IOException) logWarning ("Warning: Blast attempt failed:" ++ " " ++ err) temporaryDirectory error "findTaxonomyStart: Blast attempt failed" return (Left "")) - let logFileDirectoryPath = temporaryDirectory ++ "taxonomystart" ++ "/" + let logFileDirectoryPath = temporaryDirectory ++ "taxonomystart" ++ "/" createDirectory logFileDirectoryPath writeFile (logFileDirectoryPath ++ "/" ++ queryIndexString ++ "_1blastOutput") (show blastOutput) - logEither blastOutput temporaryDirectory - let blastHitsArePresent = either (\_ -> False) blastMatchesPresent blastOutput - if (blastHitsArePresent) + logEither blastOutput temporaryDirectory + let blastHitsArePresent = either (const False) blastMatchesPresent blastOutput + if blastHitsArePresent then do let rightBlast = fromRight blastOutput let bestHit = getBestHit rightBlast bestBlastHitTaxIdOutput <- retrieveBlastHitTaxIdEntrez [bestHit] let taxIdFromEntrySummaries = extractTaxIdFromEntrySummaries (snd bestBlastHitTaxIdOutput) - if (null taxIdFromEntrySummaries) then (error "findTaxonomyStart: - head: empty list of taxonomy entry summary for best hit") else return () + Control.Monad.when (null taxIdFromEntrySummaries) (error "findTaxonomyStart: - head: empty list of taxonomy entry summary for best hit") let rightBestTaxIdResult = head taxIdFromEntrySummaries - logMessage ("Initial TaxId: " ++ (show rightBestTaxIdResult) ++ "\n") temporaryDirectory + logMessage ("Initial TaxId: " ++ show rightBestTaxIdResult ++ "\n") temporaryDirectory CE.evaluate rightBestTaxIdResult else error "Find taxonomy start: Could not find blast hits to use as a taxonomic starting point" searchCandidates :: StaticOptions -> Maybe String -> Int -> Maybe Int -> Maybe Int -> Double -> [Sequence] -> IO SearchResult searchCandidates staticOptions finaliterationprefix iterationnumber upperTaxLimit lowerTaxLimit expectThreshold inputQuerySequences = do --let fastaSeqData = seqdata _querySequence - if (null inputQuerySequences) then error "searchCandidates: - head: empty list of query sequences" else return () + Control.Monad.when (null inputQuerySequences) $ error "searchCandidates: - head: empty list of query sequences" let queryLength = fromIntegral (seqlength (head inputQuerySequences)) let queryIndexString = "1" - let entrezTaxFilter = buildTaxFilterQuery upperTaxLimit lowerTaxLimit + let entrezTaxFilter = buildTaxFilterQuery upperTaxLimit lowerTaxLimit logVerboseMessage (verbositySwitch staticOptions) ("entrezTaxFilter" ++ show entrezTaxFilter ++ "\n") (tempDirPath staticOptions) let hitNumberQuery = buildHitNumberQuery "&HITLIST_SIZE=5000&EXPECT=" ++ show expectThreshold let registrationInfo = buildRegistration "RNAlien" "florian.eggenhofer@univie.ac.at" - let softmaskFilter = if (blastSoftmaskingToggle staticOptions) then "&FILTER=True&FILTER=m" else "" + let softmaskFilter = if blastSoftmaskingToggle staticOptions then "&FILTER=True&FILTER=m" else "" let blastQuery = BlastHTTPQuery (Just "ncbi") (Just "blastn") (blastDatabase staticOptions) inputQuerySequences (Just (hitNumberQuery ++ entrezTaxFilter ++ softmaskFilter ++ registrationInfo)) (Just (5400000000 :: Int)) --appendFile "/scratch/egg/blasttest/queries" ("\nBlast query:\n" ++ show blastQuery ++ "\n") - logVerboseMessage (verbositySwitch staticOptions) ("Sending blast query " ++ (show iterationnumber) ++ "\n") (tempDirPath staticOptions) + logVerboseMessage (verbositySwitch staticOptions) ("Sending blast query " ++ show iterationnumber ++ "\n") (tempDirPath staticOptions) blastOutput <- CE.catch (blastHTTP blastQuery) (\e -> do let err = show (e :: CE.IOException) logWarning ("Warning: Blast attempt failed:" ++ " " ++ err) (tempDirPath staticOptions) return (Left "")) - let logFileDirectoryPath = (tempDirPath staticOptions) ++ (show iterationnumber) ++ "/" ++ (fromMaybe "" finaliterationprefix) ++ "log" - logDirectoryPresent <- doesDirectoryExist logFileDirectoryPath - if (not logDirectoryPresent) - then createDirectory (logFileDirectoryPath) else return () + let logFileDirectoryPath = tempDirPath staticOptions ++ show iterationnumber ++ "/" ++ fromMaybe "" finaliterationprefix ++ "log" + logDirectoryPresent <- doesDirectoryExist logFileDirectoryPath + Control.Monad.when (not logDirectoryPresent) $ createDirectory (logFileDirectoryPath) writeFile (logFileDirectoryPath ++ "/" ++ queryIndexString ++ "_1blastOutput") (show blastOutput) - logEither blastOutput (tempDirPath staticOptions) + logEither blastOutput (tempDirPath staticOptions) let blastHitsArePresent = either (const False) blastMatchesPresent blastOutput if blastHitsArePresent then do let rightBlast = fromRight blastOutput - -- bestBlastHitTaxIdOutput <- retrieveBlastHitTaxIdEntrez [bestHit] - -- let taxIdFromEntrySummaries = extractTaxIdFromEntrySummaries bestBlastHitTaxIdOutput - -- if (null taxIdFromEntrySummaries) then (error "searchCandidates: - head: empty list of taxonomy entry summary for best hit") else return () - -- let rightBestTaxIdResult = head taxIdFromEntrySummaries - -- logVerboseMessage (verbositySwitch staticOptions) ("rightbestTaxIdResult: " ++ (show rightBestTaxIdResult) ++ "\n") (tempDirPath staticOptions) + -- bestBlastHitTaxIdOutput <- retrieveBlastHitTaxIdEntrez [bestHit] + -- let taxIdFromEntrySummaries = extractTaxIdFromEntrySummaries bestBlastHitTaxIdOutput + -- if (null taxIdFromEntrySummaries) then (error "searchCandidates: - head: empty list of taxonomy entry summary for best hit") else return () + -- let rightBestTaxIdResult = head taxIdFromEntrySummaries + -- logVerboseMessage (verbositySwitch staticOptions) ("rightbestTaxIdResult: " ++ (show rightBestTaxIdResult) ++ "\n") (tempDirPath staticOptions) let blastHits = concatMap hits (results rightBlast) writeFile (logFileDirectoryPath ++ "/" ++ queryIndexString ++ "_2blastHits") (showlines blastHits) --filter by length @@ -447,7 +447,7 @@ writeFile (logFileDirectoryPath ++ "/" ++ queryIndexString ++ "_3ablastHitsFilteredByLength") (showlines blastHitsFilteredByCoverage) --tag BlastHits with TaxId blastHitsWithTaxIdOutput <- retrieveBlastHitsTaxIdEntrez blastHitsFilteredByCoverage - let uncheckedBlastHitsWithTaxIdList = map (\(blasthits,taxIdout) -> (blasthits,extractTaxIdFromEntrySummaries taxIdout)) blastHitsWithTaxIdOutput + let uncheckedBlastHitsWithTaxIdList = map (Control.Arrow.second extractTaxIdFromEntrySummaries) blastHitsWithTaxIdOutput let checkedBlastHitsWithTaxId = filter (\(_,taxids) -> not (null taxids)) uncheckedBlastHitsWithTaxIdList --todo checked blasthittaxidswithblasthits need to be merged as taxid blasthit pairs let blastHitsWithTaxId = zip (concatMap fst checkedBlastHitsWithTaxId) (concatMap snd checkedBlastHitsWithTaxId) @@ -479,15 +479,15 @@ writeFile (logFileDirectoryPath ++ "/" ++ queryIndexString ++ "_10fullSequences") (showlines fullSequences) let maybeFractionEvalueMatch = getHitWithFractionEvalue rightBlast if isNothing maybeFractionEvalueMatch - then CE.evaluate (SearchResult [] Nothing) + then CE.evaluate (SearchResult [] Nothing) else do let fractionEvalueMatch = fromJust maybeFractionEvalueMatch let dbSize = computeDataBaseSize (e_val fractionEvalueMatch) (bits fractionEvalueMatch) (fromIntegral queryLength ::Double) CE.evaluate (SearchResult fullSequences (Just dbSize)) - else CE.evaluate (SearchResult [] Nothing) + else CE.evaluate (SearchResult [] Nothing) -- |Computes size of blast db in Mb -computeDataBaseSize :: Double -> Double -> Double -> Double +computeDataBaseSize :: Double -> Double -> Double -> Double computeDataBaseSize evalue bitscore querylength = ((evalue * 2 ** bitscore) / querylength)/10^(6 :: Integer) alignCandidates :: StaticOptions -> ModelConstruction -> String -> SearchResult -> IO ([(Sequence,Int,L.ByteString)],[(Sequence,Int,L.ByteString)]) @@ -496,7 +496,7 @@ createDirectoryIfMissing False (iterationDirectory ++ "log") if null (candidates searchResults) then do - writeFile (iterationDirectory ++ "log" ++ "/11candidates") ("No candidates to align") + writeFile (iterationDirectory ++ "log" ++ "/11candidates") "No candidates to align" return ([],[]) else do --refilter for similarity @@ -511,26 +511,26 @@ alignCandidatesInfernalMode :: StaticOptions -> ModelConstruction -> String -> Maybe Double -> [(Sequence,Int,L.ByteString)] -> IO ([(Sequence,Int,L.ByteString)],[(Sequence,Int,L.ByteString)]) alignCandidatesInfernalMode staticOptions modelConstruction multipleSearchResultPrefix blastDbSize filteredCandidates = do let iterationDirectory = tempDirPath staticOptions ++ show (iterationNumber modelConstruction) ++ "/" ++ multipleSearchResultPrefix - let candidateSequences = extractCandidateSequences filteredCandidates + let candidateSequences = extractCandidateSequences filteredCandidates logVerboseMessage (verbositySwitch staticOptions) "Alignment Mode Infernal\n" (tempDirPath staticOptions) - let indexedCandidateSequenceList = (V.toList candidateSequences) + let indexedCandidateSequenceList = V.toList candidateSequences let cmSearchFastaFilePaths = map (constructFastaFilePaths iterationDirectory) indexedCandidateSequenceList let cmSearchFilePaths = map (constructCMsearchFilePaths iterationDirectory) indexedCandidateSequenceList let covarianceModelPath = tempDirPath staticOptions ++ show (iterationNumber modelConstruction - 1) ++ "/" ++ "model.cm" mapM_ (\(number,_nucleotideSequence) -> writeFasta (iterationDirectory ++ show number ++ ".fa") [_nucleotideSequence]) indexedCandidateSequenceList - let zippedFastaCMSearchResultPaths = zip cmSearchFastaFilePaths cmSearchFilePaths + let zippedFastaCMSearchResultPaths = zip cmSearchFastaFilePaths cmSearchFilePaths --check with cmSearch mapM_ (uncurry (systemCMsearch (cpuThreads staticOptions) ("-Z " ++ show (fromJust blastDbSize)) covarianceModelPath)) zippedFastaCMSearchResultPaths - cmSearchResults <- mapM readCMSearch cmSearchFilePaths + cmSearchResults <- mapM readCMSearch cmSearchFilePaths writeFile (iterationDirectory ++ "cm_error") (concatMap show (lefts cmSearchResults)) - let rightCMSearchResults = rights cmSearchResults - let cmSearchCandidatesWithSequences = zip rightCMSearchResults filteredCandidates + let rightCMSearchResults = rights cmSearchResults + let cmSearchCandidatesWithSequences = zip rightCMSearchResults filteredCandidates let (trimmedSelectedCandidates,potentialCandidates,rejectedCandidates) = evaluePartitionTrimCMsearchHits (evalueThreshold modelConstruction) cmSearchCandidatesWithSequences createDirectoryIfMissing False (iterationDirectory ++ "log") writeFile (iterationDirectory ++ "log" ++ "/13selectedCandidates'") (showlines trimmedSelectedCandidates) writeFile (iterationDirectory ++ "log" ++ "/14rejectedCandidates'") (showlines rejectedCandidates) - writeFile (iterationDirectory ++ "log" ++ "/15potentialCandidates'") (showlines potentialCandidates) - CE.evaluate (map snd trimmedSelectedCandidates,map snd potentialCandidates) + writeFile (iterationDirectory ++ "log" ++ "/15potentialCandidates'") (showlines potentialCandidates) + return (map snd trimmedSelectedCandidates,map snd potentialCandidates) alignCandidatesInitialMode :: StaticOptions -> ModelConstruction -> String -> [(Sequence,Int,L.ByteString)] -> IO ([(Sequence,Int,L.ByteString)],[(Sequence,Int,L.ByteString)]) alignCandidatesInitialMode staticOptions modelConstruction multipleSearchResultPrefix filteredCandidates = do @@ -543,7 +543,7 @@ --write Fasta sequences let inputFastaFilepath = iterationDirectory ++ "input.fa" let inputFoldFilepath = iterationDirectory ++ "input.fold" - writeFasta (iterationDirectory ++ "input.fa") ([inputFasta modelConstruction]) + writeFasta (iterationDirectory ++ "input.fa") [inputFasta modelConstruction] logMessage (showlines (V.toList candidateSequences)) (tempDirPath staticOptions) V.mapM_ (\(number,nucleotideSequence') -> writeFasta (iterationDirectory ++ show number ++ ".fa") [nucleotideSequence']) candidateSequences let candidateFastaFilepath = V.toList (V.map (\(number,_) -> iterationDirectory ++ show number ++ ".fa") candidateSequences) @@ -567,21 +567,21 @@ let consensusMFE = map (alignmentConsensusMinimumFreeEnergy . fromRight) alifoldResults let sciidfraction = map (\(consMFE,averMFE,seqId) -> (consMFE/averMFE)/seqId) (zip3 consensusMFE averageMFEs (V.toList sequenceIdentities)) let idlog = concatMap (\(sciidfraction',consMFE,averMFE,seqId) -> show sciidfraction' ++ "," ++ show consMFE ++ "," ++ show averMFE ++ "," ++ show seqId ++ "\n")(zip4 sciidfraction consensusMFE averageMFEs (V.toList sequenceIdentities)) - writeFile (iterationDirectory ++ "log" ++ "/idlog") (idlog) + writeFile (iterationDirectory ++ "log" ++ "/idlog") idlog let alignedCandidates = zip sciidfraction filteredCandidates writeFile (iterationDirectory ++ "log" ++ "/zscores") (showlines alignedCandidates) let (selectedCandidates,rejectedCandidates) = partition (\(sciidfraction',_) -> sciidfraction' > nSCICutoff staticOptions) alignedCandidates mapM_ print (zip3 consensusMFE averageMFEs (V.toList sequenceIdentities)) writeFile (iterationDirectory ++ "log" ++ "/13selectedCandidates") (showlines selectedCandidates) writeFile (iterationDirectory ++ "log" ++ "/14rejectedCandidates") (showlines rejectedCandidates) - CE.evaluate (map snd selectedCandidates,[]) + return (map snd selectedCandidates,[]) setClusterNumber :: Int -> Int setClusterNumber x - | x <= 5 = x - | otherwise = 5 + | x <= 5 = x + | otherwise = 5 -findCutoffforClusterNumber :: Dendrogram a -> Int -> Distance -> Distance +findCutoffforClusterNumber :: Dendrogram a -> Int -> Distance -> Distance findCutoffforClusterNumber clustaloDendrogram numberOfClusters currentCutoff | currentClusterNumber >= numberOfClusters = currentCutoff | otherwise = findCutoffforClusterNumber clustaloDendrogram numberOfClusters (currentCutoff-0.01) @@ -592,20 +592,20 @@ selectQueries staticOptions modelConstruction selectedCandidates = do logVerboseMessage (verbositySwitch staticOptions) "SelectQueries\n" (tempDirPath staticOptions) --Extract sequences from modelconstruction - let alignedSequences = extractAlignedSequences (iterationNumber modelConstruction) modelConstruction + let alignedSequences = extractAlignedSequences (iterationNumber modelConstruction) modelConstruction let candidateSequences = extractQueryCandidates selectedCandidates let iterationDirectory = tempDirPath staticOptions ++ show (iterationNumber modelConstruction) ++ "/" let stockholmFilepath = iterationDirectory ++ "model" ++ ".stockholm" let alignmentSequences = map snd (V.toList (V.concat [candidateSequences,alignedSequences])) if length alignmentSequences > 3 - then do + then if (querySelectionMethod staticOptions) == "clustering" then do --write Fasta sequences writeFasta (iterationDirectory ++ "query" ++ ".fa") alignmentSequences let fastaFilepath = iterationDirectory ++ "query" ++ ".fa" let clustaloFilepath = iterationDirectory ++ "query" ++ ".clustalo" - let clustaloDistMatrixPath = iterationDirectory ++ "query" ++ ".matrix" + let clustaloDistMatrixPath = iterationDirectory ++ "query" ++ ".matrix" alignSequences "clustalo" ("--full --distmat-out=" ++ clustaloDistMatrixPath ++ " ") [fastaFilepath] [] [clustaloFilepath] [] idsDistancematrix <- readClustaloDistMatrix clustaloDistMatrixPath logEither idsDistancematrix (tempDirPath staticOptions) @@ -630,7 +630,7 @@ --Stockholm sequnces contain conservation annotation from cmalign in infernal mode let currentSelectedSequences = if (blastSoftmaskingToggle staticOptions) then stockholmSelectedSequences else fastaSelectedSequences --let currentSelectedQueries = concatMap (\querySeqId -> filter (\alignedSeq -> L.unpack (unSL (seqid alignedSeq)) == querySeqId) alignmentSequences) querySeqIds - logVerboseMessage (verbositySwitch staticOptions) ("SelectedQueries: " ++ show currentSelectedSequences ++ "\n") (tempDirPath staticOptions) + logVerboseMessage (verbositySwitch staticOptions) ("SelectedQueries: " ++ show currentSelectedSequences ++ "\n") (tempDirPath staticOptions) writeFile (tempDirPath staticOptions ++ show (iterationNumber modelConstruction) ++ "/log" ++ "/13selectedQueries") (showlines currentSelectedSequences) CE.evaluate currentSelectedSequences else do @@ -663,7 +663,7 @@ let stockholmFilepath = outputDirectory ++ "model" ++ ".stockholm" --- let reformatedClustalFilepath = outputDirectory ++ "model" ++ ".clustal.reformated" let updatedStructureStockholmFilepath = outputDirectory ++ "newstructuremodel" ++ ".stockholm" - let cmalignCMFilepath = (tempDirPath staticOptions) ++ show (iterationNumber modelConstruction - 2) ++ "/" ++ "model" ++ ".cm" + let cmalignCMFilepath = tempDirPath staticOptions ++ show (iterationNumber modelConstruction - 2) ++ "/" ++ "model" ++ ".cm" let cmFilepath = outputDirectory ++ "model" ++ ".cm" let cmCalibrateFilepath = outputDirectory ++ "model" ++ ".cmcalibrate" let cmBuildFilepath = outputDirectory ++ "model" ++ ".cmbuild" @@ -703,7 +703,7 @@ inputAln <- readFile stockholmFilepath inputRNAalifold <- readRNAalifold alifoldFilepath if isLeft inputRNAalifold - then do + then return (show (fromLeft inputRNAalifold)) else do let alifoldstructure = alignmentConsensusDotBracket (fromRight inputRNAalifold) @@ -711,7 +711,7 @@ let structureIndices = V.toList (V.findIndices isStructureLine seedLinesVector) let updatedStructureElements = updateStructureElements seedLinesVector alifoldstructure structureIndices let newVector = seedLinesVector V.// updatedStructureElements - let newVectorString = concatMap (++ "\n") (V.toList newVector) + let newVectorString = unlines (V.toList newVector) writeFile updatedStructureStockholmFilepath newVectorString return [] @@ -736,14 +736,14 @@ iterationSummaryLog mC = output where upperTaxonomyLimitOutput = maybe "not set" show (upperTaxonomyLimit mC) output = "Upper taxonomy id limit: " ++ upperTaxonomyLimitOutput ++ ", Collected members: " ++ show (length (concatMap sequenceRecords (taxRecords mC))) ++ "\n" - + -- | Used for passing progress to Alien server iterationSummary :: ModelConstruction -> StaticOptions -> IO() iterationSummary mC sO = do --iteration -- tax limit -- bitscore cutoff -- blastresult -- aligned seqs --queries --fa link --aln link --cm link let upperTaxonomyLimitOutput = maybe "not set" show (upperTaxonomyLimit mC) let output = show (iterationNumber mC) ++ "," ++ upperTaxonomyLimitOutput ++ "," ++ show (length (concatMap sequenceRecords (taxRecords mC))) - writeFile ((tempDirPath sO) ++ "/log/" ++ show (iterationNumber mC) ++ ".log") output + writeFile (tempDirPath sO ++ "/log/" ++ show (iterationNumber mC) ++ ".log") output -- | Used for passing progress to Alien server resultSummary :: ModelConstruction -> StaticOptions -> IO() @@ -751,16 +751,16 @@ --iteration -- tax limit -- bitscore cutoff -- blastresult -- aligned seqs --queries --fa link --aln link --cm link let upperTaxonomyLimitOutput = maybe "not set" show (upperTaxonomyLimit mC) let output = show (iterationNumber mC) ++ "," ++ upperTaxonomyLimitOutput ++ "," ++ show (length (concatMap sequenceRecords (taxRecords mC))) - writeFile (tempDirPath sO ++ "/log/result" ++ ".log") output - -readClustaloDistMatrix :: String -> IO (Either ParseError ([String],Matrix Double)) -readClustaloDistMatrix = parseFromFile genParserClustaloDistMatrix - + writeFile (tempDirPath sO ++ "/log/result" ++ ".log") output + +readClustaloDistMatrix :: String -> IO (Either ParseError ([String],Matrix Double)) +readClustaloDistMatrix = parseFromFile genParserClustaloDistMatrix + genParserClustaloDistMatrix :: GenParser Char st ([String],Matrix Double) genParserClustaloDistMatrix = do _ <- many1 digit newline - clustaloDistRow <- many1 (try genParserClustaloDistRow) + clustaloDistRow <- many1 (try genParserClustaloDistRow) eof return (map fst clustaloDistRow,fromLists (map snd clustaloDistRow)) @@ -801,21 +801,21 @@ --filterIdenticalSequencesWithOrigin [] _ = [] -- | Filter a list of similar extended blast hits -filterIdenticalSequences :: [(Sequence,Int,L.ByteString)] -> Double -> [(Sequence,Int,L.ByteString)] +filterIdenticalSequences :: [(Sequence,Int,L.ByteString)] -> Double -> [(Sequence,Int,L.ByteString)] filterIdenticalSequences (headSequence:rest) identitycutoff = result - where filteredSequences = filter (\x -> sequenceIdentity (firstOfTriple headSequence) (firstOfTriple x) < identitycutoff) rest + where filteredSequences = filter (\x -> sequenceIdentity (firstOfTriple headSequence) (firstOfTriple x) < identitycutoff) rest result = headSequence:filterIdenticalSequences filteredSequences identitycutoff filterIdenticalSequences [] _ = [] -- | Filter sequences too similar to already aligned sequences -filterWithCollectedSequences :: [(Sequence,Int,L.ByteString)] -> [Sequence] -> Double -> [(Sequence,Int,L.ByteString)] -filterWithCollectedSequences inputCandidates collectedSequences identitycutoff = filter (isUnSimilarSequence collectedSequences identitycutoff . firstOfTriple) inputCandidates +filterWithCollectedSequences :: [(Sequence,Int,L.ByteString)] -> [Sequence] -> Double -> [(Sequence,Int,L.ByteString)] +filterWithCollectedSequences inputCandidates collectedSequences identitycutoff = filter (isUnSimilarSequence collectedSequences identitycutoff . firstOfTriple) inputCandidates --filterWithCollectedSequences [] [] _ = [] -- | Filter alignment entries by similiarity filterIdenticalSequences' :: [Sequence] -> Double -> [Sequence] filterIdenticalSequences' (headEntry:rest) identitycutoff = result - where filteredEntries = filter (\ x -> (sequenceIdentity headEntry x) < identitycutoff) rest + where filteredEntries = filter (\ x -> sequenceIdentity headEntry x < identitycutoff) rest result = headEntry:filterIdenticalSequences' filteredEntries identitycutoff filterIdenticalSequences' [] _ = [] @@ -827,14 +827,14 @@ --filterIdenticalAlignmentEntry [] _ = [] isUnSimilarSequence :: [Sequence] -> Double -> Sequence -> Bool -isUnSimilarSequence collectedSequences identitycutoff checkSequence = any (\ x -> (sequenceIdentity checkSequence x) < identitycutoff) collectedSequences - +isUnSimilarSequence collectedSequences identitycutoff checkSequence = any (\ x -> sequenceIdentity checkSequence x < identitycutoff) collectedSequences + firstOfTriple :: (t, t1, t2) -> t -firstOfTriple (a,_,_) = a +firstOfTriple (a,_,_) = a -- | Check if the result field of BlastResult is filled and if hits are present blastMatchesPresent :: BlastResult -> Bool -blastMatchesPresent blastResult +blastMatchesPresent blastResult | null resultList = False | otherwise = True where resultList = concatMap matches (concatMap hits (results blastResult)) @@ -848,8 +848,8 @@ maximumDistance = maximum [T.length text1, T.length text2] distanceDouble = toInteger ( fromJust distance ) identityPercent = 1 - (fromIntegral distanceDouble/fromIntegral maximumDistance) - - + + -- | Compute identity of sequences -- stringIdentity :: String -> String -> Double -- stringIdentity string1 string2 = identityPercent @@ -869,7 +869,7 @@ identityPercent = 100 - ((fromIntegral distance/fromIntegral maximumDistance) * (read "100" ::Double)) getTaxonomicContextEntrez :: Maybe Int -> Maybe Taxon -> IO (Maybe Taxon) -getTaxonomicContextEntrez upperTaxLimit currentTaxonomicContext = do +getTaxonomicContextEntrez upperTaxLimit currentTaxonomicContext = if isJust upperTaxLimit then if isJust currentTaxonomicContext then return currentTaxonomicContext @@ -878,12 +878,12 @@ else return Nothing setTaxonomicContextEntrez :: Int -> Maybe Taxon -> Maybe Int -> (Maybe Int, Maybe Int) -setTaxonomicContextEntrez currentIterationNumber currentTaxonomicContext subTreeTaxId +setTaxonomicContextEntrez currentIterationNumber currentTaxonomicContext subTreeTaxId | currentIterationNumber == 0 = (subTreeTaxId, Nothing) | otherwise = setUpperLowerTaxLimitEntrez (fromJust subTreeTaxId) (fromJust currentTaxonomicContext) - + -- setTaxonomic Context for next candidate search, the upper bound of the last search become the lower bound of the next -setUpperLowerTaxLimitEntrez :: Int -> Taxon -> (Maybe Int, Maybe Int) +setUpperLowerTaxLimitEntrez :: Int -> Taxon -> (Maybe Int, Maybe Int) setUpperLowerTaxLimitEntrez subTreeTaxId currentTaxonomicContext = (upperLimit,lowerLimit) where upperLimit = raiseTaxIdLimitEntrez subTreeTaxId currentTaxonomicContext lowerLimit = Just subTreeTaxId @@ -895,7 +895,7 @@ lineageExVector = V.fromList (lineageEx taxon) --the input taxid is not part of the lineage, therefor we look for further taxids in the lineage after we used the parent tax id of the input node parentNodeTaxId = if subTreeTaxId == taxonTaxId taxon then Just (taxonParentTaxId taxon) else linageNodeTaxId - + constructNext :: Int -> ModelConstruction -> [(Sequence,Int,L.ByteString)] -> Maybe Int -> Maybe Taxon -> [Sequence] -> [SearchResult] -> Bool -> ModelConstruction constructNext currentIterationNumber modelconstruction alignmentResults upperTaxLimit inputTaxonomicContext inputSelectedQueries inputPotentialMembers toggleInfernalAlignmentModeTrue = nextModelConstruction where newIterationNumber = currentIterationNumber + 1 @@ -903,11 +903,11 @@ potMembers = potentialMembers modelconstruction ++ inputPotentialMembers currentAlignmentMode = toggleInfernalAlignmentModeTrue || alignmentModeInfernal modelconstruction nextModelConstruction = ModelConstruction newIterationNumber (inputFasta modelconstruction) taxEntries upperTaxLimit inputTaxonomicContext (evalueThreshold modelconstruction) currentAlignmentMode inputSelectedQueries potMembers - + buildTaxRecords :: [(Sequence,Int,L.ByteString)] -> Int -> [TaxonomyRecord] buildTaxRecords alignmentResults currentIterationNumber = taxonomyRecords where taxIdGroups = groupBy sameTaxIdAlignmentResult alignmentResults - taxonomyRecords = map (buildTaxRecord currentIterationNumber) taxIdGroups + taxonomyRecords = map (buildTaxRecord currentIterationNumber) taxIdGroups sameTaxIdAlignmentResult :: (Sequence,Int,L.ByteString) -> (Sequence,Int,L.ByteString) -> Bool sameTaxIdAlignmentResult (_,taxId1,_) (_,taxId2,_) = taxId1 == taxId2 @@ -918,8 +918,8 @@ seqRecords = map (buildSeqRecord currentIterationNumber) entries taxRecord = TaxonomyRecord recordTaxId seqRecords -buildSeqRecord :: Int -> (Sequence,Int,L.ByteString) -> SequenceRecord -buildSeqRecord currentIterationNumber (parsedFasta,_,seqSubject) = SequenceRecord parsedFasta currentIterationNumber seqSubject +buildSeqRecord :: Int -> (Sequence,Int,L.ByteString) -> SequenceRecord +buildSeqRecord currentIterationNumber (parsedFasta,_,seqSubject) = SequenceRecord parsedFasta currentIterationNumber seqSubject -- | Partitions sequences by containing a cmsearch hit and extracts the hit region as new sequence evaluePartitionTrimCMsearchHits :: Double -> [(CMsearch,(Sequence, Int, L.ByteString))] -> ([(CMsearch,(Sequence, Int, L.ByteString))],[(CMsearch,(Sequence, Int, L.ByteString))],[(CMsearch,(Sequence, Int, L.ByteString))]) @@ -928,8 +928,8 @@ (prefilteredCandidates,rejectedCandidates) = partition (\(cmSearchResult,_) -> any (\hitScore' -> potentialMemberseValueThreshold >= hitEvalue hitScore') (cmsearchHits cmSearchResult)) cmSearchCandidatesWithSequences (selectedCandidates,potentialCandidates) = partition (\(cmSearchResult,_) -> any (\hitScore' -> eValueThreshold >= hitEvalue hitScore') (cmsearchHits cmSearchResult)) prefilteredCandidates trimmedSelectedCandidates = map (\(cmSearchResult,inputSequence) -> (cmSearchResult,trimCMsearchHit cmSearchResult inputSequence)) selectedCandidates - - + + trimCMsearchHit :: CMsearch -> (Sequence, Int, L.ByteString) -> (Sequence, Int, L.ByteString) trimCMsearchHit cmSearchResult (inputSequence,b,c) = (subSequence,b,c) where hitScoreEntry = head (cmsearchHits cmSearchResult) @@ -941,21 +941,21 @@ -- | Extract a substring with coordinates from cmsearch, first nucleotide has index 1 cmSearchsubString :: Int -> Int -> String -> String -cmSearchsubString startSubString endSubString inputString +cmSearchsubString startSubString endSubString inputString | startSubString < endSubString = take (endSubString - (startSubString -1))(drop (startSubString - 1) inputString) | startSubString > endSubString = take (reverseEnd - (reverseStart - 1))(drop (reverseStart - 1 ) (reverse inputString)) | otherwise = take (endSubString - (startSubString -1))(drop (startSubString - 1) inputString) where stringLength = length inputString reverseStart = stringLength - (startSubString + 1) reverseEnd = stringLength - (endSubString - 1) - -extractQueries :: Int -> ModelConstruction -> [Sequence] + +extractQueries :: Int -> ModelConstruction -> [Sequence] extractQueries foundSequenceNumber modelconstruction - | foundSequenceNumber < 3 = [fastaSeqData] - | otherwise = querySequences' + | foundSequenceNumber < 3 = [fastaSeqData] + | otherwise = querySequences' where fastaSeqData = inputFasta modelconstruction querySequences' = selectedQueries modelconstruction - + extractQueryCandidates :: [(Sequence,Int,L.ByteString)] -> V.Vector (Int,Sequence) extractQueryCandidates querycandidates = indexedSeqences where sequences = map (\(candidateSequence,_,_) -> candidateSequence) querycandidates @@ -966,7 +966,7 @@ | isNothing upperTaxLimit = "" | isNothing lowerTaxLimit = "&ENTREZ_QUERY=" ++ encodedTaxIDQuery (fromJust upperTaxLimit) | otherwise = "&ENTREZ_QUERY=" ++ "%28txid" ++ show (fromJust upperTaxLimit) ++ "%5BORGN%5D%29" ++ "NOT" ++ "%28txid" ++ show (fromJust lowerTaxLimit) ++ "%5BORGN%5D&EQ_OP%29" - + buildHitNumberQuery :: String -> String buildHitNumberQuery hitNumber | hitNumber == "" = "" @@ -981,14 +981,14 @@ -- | Create session id for RNAlien createSessionID :: Maybe String -> IO String -createSessionID sessionIdentificator = do +createSessionID sessionIdentificator = if isJust sessionIdentificator then return (fromJust sessionIdentificator) else do randomNumber <- randomIO :: IO Int16 let sessionId = randomid (abs (randomNumber)) return sessionId - + -- | Run external locarna command and read the output into the corresponding datatype systemlocarna :: String -> (String,String,String,String) -> IO ExitCode systemlocarna options (inputFilePath1, inputFilePath2, clustalformatoutputFilePath, outputFilePath) = system ("locarna " ++ options ++ " --clustal=" ++ clustalformatoutputFilePath ++ " " ++ inputFilePath1 ++ " " ++ inputFilePath2 ++ " > " ++ outputFilePath) @@ -996,11 +996,11 @@ -- | Run external mlocarna command and read the output into the corresponding datatype, there is also a folder created at the location of the input fasta file systemMlocarna :: String -> (String,String) -> IO ExitCode systemMlocarna options (inputFilePath, outputFilePath) = system ("mlocarna " ++ options ++ " " ++ inputFilePath ++ " > " ++ outputFilePath) - + -- | Run external mlocarna command and read the output into the corresponding datatype, there is also a folder created at the location of the input fasta file, the job is terminated after the timeout provided in seconds systemMlocarnaWithTimeout :: String -> String -> (String,String) -> IO ExitCode systemMlocarnaWithTimeout timeout options (inputFilePath, outputFilePath) = system ("timeout " ++ timeout ++"s "++ "mlocarna " ++ options ++ " " ++ inputFilePath ++ " > " ++ outputFilePath) - + -- | Run external clustalo command and return the Exitcode systemClustalw2 :: String -> (String,String,String) -> IO ExitCode systemClustalw2 options (inputFilePath, outputFilePath, summaryFilePath) = system ("clustalw2 " ++ options ++ "-INFILE=" ++ inputFilePath ++ " -OUTFILE=" ++ outputFilePath ++ ">" ++ summaryFilePath) @@ -1011,8 +1011,8 @@ -- | Run external CMbuild command and read the output into the corresponding datatype systemCMbuild :: String -> String -> String -> String -> IO ExitCode -systemCMbuild options alignmentFilepath modelFilepath outputFilePath = system ("cmbuild " ++ options ++ " " ++ modelFilepath ++ " " ++ alignmentFilepath ++ " > " ++ outputFilePath) - +systemCMbuild options alignmentFilepath modelFilepath outputFilePath = system ("cmbuild " ++ options ++ " " ++ modelFilepath ++ " " ++ alignmentFilepath ++ " > " ++ outputFilePath) + -- | Run CMCompare and read the output into the corresponding datatype systemCMcompare :: String -> String -> String -> IO ExitCode systemCMcompare model1path model2path outputFilePath = system ("CMCompare -q " ++ model1path ++ " " ++ model2path ++ " >" ++ outputFilePath) @@ -1026,14 +1026,14 @@ systemCMstat covarianceModelPath outputPath = system ("cmstat " ++ covarianceModelPath ++ " > " ++ outputPath) -- | Run CMcalibrate and return exitcode -systemCMcalibrate :: String -> Int -> String -> String -> IO ExitCode -systemCMcalibrate mode cpus covarianceModelPath outputPath +systemCMcalibrate :: String -> Int -> String -> String -> IO ExitCode +systemCMcalibrate mode cpus covarianceModelPath outputPath | mode == "fast" = system ("cmcalibrate --beta 1E-4 --cpu " ++ show cpus ++ " " ++ covarianceModelPath ++ "> " ++ outputPath) | otherwise = system ("cmcalibrate --cpu " ++ show cpus ++ " " ++ covarianceModelPath ++ "> " ++ outputPath) -- | Run CMcalibrate and return exitcode -systemCMalign :: String -> String -> String -> String -> IO ExitCode +systemCMalign :: String -> String -> String -> String -> IO ExitCode systemCMalign options filePathCovarianceModel filePathSequence filePathAlignment = system ("cmalign " ++ options ++ " " ++ filePathCovarianceModel ++ " " ++ filePathSequence ++ "> " ++ filePathAlignment) compareCM :: String -> String -> String -> IO (Either String Double) @@ -1053,18 +1053,18 @@ let bitscore2 = read (decodedCmCompareOutput !! 3) :: Double let minmax = minimum [bitscore1,bitscore2] return (Right minmax) - + readInt :: String -> Int readInt = read readDouble :: String -> Double readDouble = read - + extractCandidateSequences :: [(Sequence,Int,L.ByteString)] -> V.Vector (Int,Sequence) extractCandidateSequences candidates' = indexedSeqences where sequences = map (\(inputSequence,_,_) -> inputSequence) candidates' indexedSeqences = V.map (\(number,inputSequence) -> (number + 1,inputSequence))(V.indexed (V.fromList sequences)) - + extractAlignedSequences :: Int -> ModelConstruction -> V.Vector (Int,Sequence) extractAlignedSequences iterationnumber modelconstruction | iterationnumber == 0 = V.map (\(number,seq') -> (number + 1,seq')) (V.indexed (V.fromList [inputSequence])) @@ -1076,7 +1076,7 @@ indexedSeqRecords = V.map (\(number,seq') -> (number + 1,seq')) (V.indexed (V.fromList (inputSequence : map nucleotideSequence seqRecords))) filterByParentTaxId :: [(BlastHit,Int)] -> Bool -> [(BlastHit,Int)] -filterByParentTaxId blastHitsWithParentTaxId singleHitPerParentTaxId +filterByParentTaxId blastHitsWithParentTaxId singleHitPerParentTaxId | singleHitPerParentTaxId = singleBlastHitperParentTaxId | otherwise = blastHitsWithParentTaxId where blastHitsWithParentTaxIdSortedByParentTaxId = sortBy compareTaxId blastHitsWithParentTaxId @@ -1084,7 +1084,7 @@ singleBlastHitperParentTaxId = map (maximumBy compareHitEValue) blastHitsWithParentTaxIdGroupedByParentTaxId filterByHitLength :: [BlastHit] -> Int -> Bool -> [BlastHit] -filterByHitLength blastHits queryLength filterOn +filterByHitLength blastHits queryLength filterOn | filterOn = filteredBlastHits | otherwise = blastHits where filteredBlastHits = filter (hitLengthCheck queryLength) blastHits @@ -1099,13 +1099,13 @@ maxHtoHSP = fromJust (find (\hsp -> maxHto == h_to hsp) blastMatches) minHonQuery = q_from minHfromHSP maxHonQuery = q_to maxHtoHSP - startCoordinate = minHfrom - minHonQuery - endCoordinate = maxHto + (queryLength - maxHonQuery) + startCoordinate = minHfrom - minHonQuery + endCoordinate = maxHto + (queryLength - maxHonQuery) fullSeqLength = endCoordinate - startCoordinate lengthStatus = fullSeqLength < (queryLength * 3) filterByCoverage :: [BlastHit] -> Int -> Bool -> [BlastHit] -filterByCoverage blastHits queryLength filterOn +filterByCoverage blastHits queryLength filterOn | filterOn = filteredBlastHits | otherwise = blastHits where filteredBlastHits = filter (coverageCheck queryLength) blastHits @@ -1115,8 +1115,8 @@ coverageCheck queryLength blastHit = coverageStatus where blastMatches = matches blastHit maxIdentity = fromIntegral (maximum (map (snd . Bio.BlastXML.identity) blastMatches)) - coverageStatus = (maxIdentity/(fromIntegral queryLength))* (100 :: Double) >= (80 :: Double) - + coverageStatus = (maxIdentity/fromIntegral queryLength)* (100 :: Double) >= (80 :: Double) + -- | Wrapper for retrieveFullSequence that rerequests incomplete return sequees retrieveFullSequences :: StaticOptions -> [(String,Int,Int,String,T.Text,Int,L.ByteString)] -> IO [(Sequence,Int,L.ByteString)] retrieveFullSequences staticOptions requestedSequences = do @@ -1133,21 +1133,21 @@ let unwrappedRetrievals = map (\(x,y,z) -> (fromJust x,y,z)) (map fst successfulRetrievals ++ reRetrievedSequences) CE.evaluate unwrappedRetrievals else CE.evaluate (map (\(x,y,z) -> (fromJust x,y,z)) fullSequences) - + retrieveFullSequence :: String -> (String,Int,Int,String,T.Text,Int,L.ByteString) -> IO (Maybe Sequence,Int,L.ByteString) retrieveFullSequence temporaryDirectoryPath (nucleotideId,seqStart,seqStop,strand,_,taxid,subject') = do let program' = Just "efetch" let database' = Just "nucleotide" let registrationInfo = buildRegistration "RNAlien" "florian.eggenhofer@univie.ac.at" let queryString = "id=" ++ nucleotideId ++ "&seq_start=" ++ show seqStart ++ "&seq_stop=" ++ show seqStop ++ "&rettype=fasta" ++ "&strand=" ++ strand ++ registrationInfo - let entrezQuery = EntrezHTTPQuery program' database' queryString + let entrezQuery = EntrezHTTPQuery program' database' queryString result <- CE.catch (entrezHTTP entrezQuery) (\e -> do let err = show (e :: CE.IOException) logWarning ("Warning: Full sequence retrieval failed:" ++ " " ++ err) temporaryDirectoryPath return []) if null result then return (Nothing,taxid,subject') - else do + else if null ((mkSeqs . L.lines) (L.pack result)) then return (Nothing,taxid,subject') else do @@ -1155,9 +1155,9 @@ if L.null (unSD (seqdata parsedFasta)) then return (Nothing,taxid,subject') else CE.evaluate (Just parsedFasta,taxid,subject') - + getRequestedSequenceElement :: Int -> (BlastHit,Int) -> (String,Int,Int,String,T.Text,Int,L.ByteString) -getRequestedSequenceElement queryLength (blastHit,taxid) +getRequestedSequenceElement queryLength (blastHit,taxid) | blastHitIsReverseComplement (blastHit,taxid) = getReverseRequestedSequenceElement queryLength (blastHit,taxid) | otherwise = getForwardRequestedSequenceElement queryLength (blastHit,taxid) @@ -1180,10 +1180,10 @@ minHonQuery = q_from blastMatch maxHonQuery = q_to blastMatch --unsafe coordinates may exceed length of available sequence - unsafestartcoordinate = minHfrom - minHonQuery - unsafeendcoordinate = maxHto + (queryLength - maxHonQuery) + unsafestartcoordinate = minHfrom - minHonQuery + unsafeendcoordinate = maxHto + (queryLength - maxHonQuery) startcoordinate = lowerBoundryCoordinateSetter 0 unsafestartcoordinate - endcoordinate = upperBoundryCoordinateSetter blastHitOriginSequenceLength unsafeendcoordinate + endcoordinate = upperBoundryCoordinateSetter blastHitOriginSequenceLength unsafeendcoordinate strand = "1" ---- --blastMatches = matches blastHit @@ -1214,19 +1214,19 @@ getReverseRequestedSequenceElement :: Int -> (BlastHit,Int) -> (String,Int,Int,String,T.Text,Int,L.ByteString) getReverseRequestedSequenceElement queryLength (blastHit,taxid) = (geneIdentifier',startcoordinate,endcoordinate,strand,accession',taxid,subjectBlast) where accession' = extractAccession blastHit - subjectBlast = unSL (subject blastHit) + subjectBlast = unSL (subject blastHit) geneIdentifier' = extractGeneId blastHit blastMatch = head (matches blastHit) - blastHitOriginSequenceLength = slength blastHit - maxHfrom = h_from blastMatch + blastHitOriginSequenceLength = slength blastHit + maxHfrom = h_from blastMatch minHto = h_to blastMatch minHonQuery = q_from blastMatch maxHonQuery = q_to blastMatch --unsafe coordinates may exceed length of avialable sequence - unsafestartcoordinate = maxHfrom + minHonQuery - unsafeendcoordinate = minHto - (queryLength - maxHonQuery) - startcoordinate = lowerBoundryCoordinateSetter 0 unsafeendcoordinate - endcoordinate = upperBoundryCoordinateSetter blastHitOriginSequenceLength unsafestartcoordinate + unsafestartcoordinate = maxHfrom + minHonQuery + unsafeendcoordinate = minHto - (queryLength - maxHonQuery) + startcoordinate = lowerBoundryCoordinateSetter 0 unsafeendcoordinate + endcoordinate = upperBoundryCoordinateSetter blastHitOriginSequenceLength unsafestartcoordinate strand = "2" -- --blastMatches = matches blastHit @@ -1252,7 +1252,7 @@ alignSequences :: String -> String -> [String] -> [String] -> [String] -> [String] -> IO () alignSequences program' options fastaFilepaths fastaFilepaths2 alignmentFilepaths summaryFilepaths = do let zipped4Filepaths = zip4 fastaFilepaths fastaFilepaths2 alignmentFilepaths summaryFilepaths - let zipped3Filepaths = zip3 fastaFilepaths alignmentFilepaths summaryFilepaths + let zipped3Filepaths = zip3 fastaFilepaths alignmentFilepaths summaryFilepaths let zippedFilepaths = zip fastaFilepaths alignmentFilepaths let timeout = "3600" case program' of @@ -1267,32 +1267,32 @@ constructCMsearchFilePaths :: String -> (Int, Sequence) -> String constructCMsearchFilePaths currentDirectory (fastaIdentifier, _) = currentDirectory ++ show fastaIdentifier ++".cmsearch" - + -- Smaller e-Values are greater, the maximum function is applied -compareHitEValue :: (BlastHit,Int) -> (BlastHit,Int) -> Ordering +compareHitEValue :: (BlastHit,Int) -> (BlastHit,Int) -> Ordering compareHitEValue (hit1,_) (hit2,_) | hitEValue hit1 > hitEValue hit2 = LT | hitEValue hit1 < hitEValue hit2 = GT -- in case of equal evalues the first hit is selected - | hitEValue hit1 == hitEValue hit2 = GT + | hitEValue hit1 == hitEValue hit2 = GT -- comparing (hitEValue . Down . fst) -compareHitEValue (_,_) (_,_) = EQ +compareHitEValue (_,_) (_,_) = EQ -compareTaxId :: (BlastHit,Int) -> (BlastHit,Int) -> Ordering +compareTaxId :: (BlastHit,Int) -> (BlastHit,Int) -> Ordering compareTaxId (_,taxId1) (_,taxId2) | taxId1 > taxId2 = LT | taxId1 < taxId2 = GT -- in case of equal evalues the first hit is selected | taxId1 == taxId2 = EQ compareTaxId (_,_) (_,_) = EQ - + sameTaxId :: (BlastHit,Int) -> (BlastHit,Int) -> Bool sameTaxId (_,taxId1) (_,taxId2) = taxId1 == taxId2 -- | NCBI uses the e-Value of the best HSP as the Hits e-Value hitEValue :: BlastHit -> Double hitEValue hit = minimum (map e_val (matches hit)) - + convertFastaFoldStockholm :: Sequence -> String -> String convertFastaFoldStockholm fastasequence foldedStructure = stockholmOutput where alnHeader = "# STOCKHOLM 1.0\n\n" @@ -1306,7 +1306,7 @@ structureString = "#=GC SS_cons" ++ replicate (spacerLength' - 12) ' ' ++ foldedStructure ++ "\n" bottom = "//" stockholmOutput = alnHeader ++ entrystring ++ structureString ++ bottom - + convertClustaltoStockholm :: StructuralClustalAlignment -> T.Text convertClustaltoStockholm parsedMlocarnaAlignment = stockholmOutput where alnHeader = T.pack "# STOCKHOLM 1.0\n\n" @@ -1316,21 +1316,21 @@ maxIdentifierLenght = maximum (map (T.length . entrySequenceIdentifier) clustalAlignment) spacerLength' = maxIdentifierLenght + 2 stockholmEntries = T.concat (map (buildStockholmAlignmentEntries spacerLength') mergedEntries) - structureString = (T.pack "#=GC SS_cons") `T.append` T.replicate (spacerLength' - 12) (T.pack " ") `T.append` secondaryStructureTrack parsedMlocarnaAlignment `T.append` (T.pack "\n") + structureString = T.pack "#=GC SS_cons" `T.append` T.replicate (spacerLength' - 12) (T.pack " ") `T.append` secondaryStructureTrack parsedMlocarnaAlignment `T.append` T.pack "\n" bottom = T.pack "//" stockholmOutput = alnHeader `T.append` stockholmEntries `T.append` structureString `T.append` bottom mergeEntry :: [ClustalAlignmentEntry] -> T.Text -> ClustalAlignmentEntry mergeEntry clustalAlignment uniqueId = mergedEntry where idEntries = filter (\entry -> entrySequenceIdentifier entry==uniqueId) clustalAlignment - mergedSeq = foldr ((T.append) . entryAlignedSequence) (T.pack "") idEntries + mergedSeq = foldr (T.append . entryAlignedSequence) (T.pack "") idEntries mergedEntry = ClustalAlignmentEntry uniqueId mergedSeq buildStockholmAlignmentEntries :: Int -> ClustalAlignmentEntry -> T.Text buildStockholmAlignmentEntries inputSpacerLength entry = entrystring where idLength = T.length (T.filter (/= '\n') (entrySequenceIdentifier entry)) spacer = T.replicate (inputSpacerLength - idLength) (T.pack " ") - entrystring = entrySequenceIdentifier entry `T.append` spacer `T.append` entryAlignedSequence entry `T.append` (T.pack "\n") + entrystring = entrySequenceIdentifier entry `T.append` spacer `T.append` entryAlignedSequence entry `T.append` T.pack "\n" retrieveTaxonomicContextEntrez :: Int -> IO (Maybe Taxon) retrieveTaxonomicContextEntrez inputTaxId = do @@ -1339,7 +1339,7 @@ let taxIdString = show inputTaxId let registrationInfo = buildRegistration "RNAlien" "florian.eggenhofer@univie.ac.at" let queryString = "id=" ++ taxIdString ++ registrationInfo - let entrezQuery = EntrezHTTPQuery program' database' queryString + let entrezQuery = EntrezHTTPQuery program' database' queryString result <- entrezHTTP entrezQuery if null result then do @@ -1350,10 +1350,10 @@ --print taxon if null (lineageEx taxon) then error "Retrieved taxonomic context taxon from NCBI Entrez with empty lineage, cannot proceed." - else CE.evaluate (Just taxon) + else return (Just taxon) retrieveParentTaxIdEntrez :: [(BlastHit,Int)] -> IO [(BlastHit,Int)] -retrieveParentTaxIdEntrez blastHitsWithHitTaxids = do +retrieveParentTaxIdEntrez blastHitsWithHitTaxids = if not (null blastHitsWithHitTaxids) then do let program' = Just "efetch" @@ -1364,7 +1364,7 @@ let taxIdQuery = intercalate "," taxIdStrings let registrationInfo = buildRegistration "RNAlien" "florian.eggenhofer@univie.ac.at" let queryString = "id=" ++ taxIdQuery ++ registrationInfo - let entrezQuery = EntrezHTTPQuery program' database' queryString + let entrezQuery = EntrezHTTPQuery program' database' queryString result <- entrezHTTP entrezQuery let parentTaxIds = readEntrezParentIds result if null parentTaxIds @@ -1384,10 +1384,10 @@ retrieveBlastHitsTaxIdEntrez blastHits = do let splits = portionListElements blastHits 20 mapM retrieveBlastHitTaxIdEntrez splits - + retrieveBlastHitTaxIdEntrez :: [BlastHit] -> IO ([BlastHit],String) -retrieveBlastHitTaxIdEntrez blastHits = do +retrieveBlastHitTaxIdEntrez blastHits = if not (null blastHits) then do let geneIds = map extractGeneId blastHits @@ -1395,7 +1395,7 @@ let registrationInfo = buildRegistration "RNAlien" "florian.eggenhofer@univie.ac.at" let query' = "id=" ++ idList ++ registrationInfo let entrezQuery = EntrezHTTPQuery (Just "esummary") (Just "nucleotide") query' - threadDelay 10000000 + threadDelay 10000000 result <- entrezHTTP entrezQuery CE.evaluate (blastHits,result) else return (blastHits,"") @@ -1415,7 +1415,7 @@ extractAccession currentBlastHit = accession' where splitedFields = T.splitOn (T.pack "|") (DTE.decodeUtf8 (L.toStrict (hitId currentBlastHit))) accession' = splitedFields !! 3 - + extractGeneId :: BlastHit -> String extractGeneId currentBlastHit = nucleotideId where truncatedId = drop 3 (L.unpack (hitId currentBlastHit)) @@ -1426,18 +1426,18 @@ extractTaxIdfromDocumentSummary documentSummary = itemContent (fromJust (find (\item -> "TaxId" == itemName item) (summaryItems documentSummary))) getBestHit :: BlastResult -> BlastHit -getBestHit blastResult +getBestHit blastResult | null (concatMap hits (results blastResult)) = error "getBestHit - head: empty list" | otherwise = head (hits (head (results blastResult))) -- Blast returns low evalues with zero instead of the exact number getHitWithFractionEvalue :: BlastResult -> Maybe BlastMatch -getHitWithFractionEvalue blastResult +getHitWithFractionEvalue blastResult | null (concatMap hits (results blastResult)) = Nothing | otherwise = find (\match -> e_val match /= (0 ::Double)) (concatMap matches (concatMap hits (results blastResult))) showlines :: Show a => [a] -> String -showlines = concatMap (\x -> show x ++ "\n") +showlines = concatMap (\x -> show x ++ "\n") logMessage :: String -> String -> IO () logMessage logoutput temporaryDirectoryPath = appendFile (temporaryDirectoryPath ++ "Log") logoutput @@ -1446,10 +1446,10 @@ logWarning logoutput temporaryDirectoryPath = appendFile (temporaryDirectoryPath ++ "log/warnings") logoutput logVerboseMessage :: Bool -> String -> String -> IO () -logVerboseMessage verboseTrue logoutput temporaryDirectoryPath +logVerboseMessage verboseTrue logoutput temporaryDirectoryPath | verboseTrue = appendFile (temporaryDirectoryPath ++ "Log") (show logoutput) | otherwise = return () - + logEither :: (Show a) => Either a b -> String -> IO () logEither (Left logoutput) temporaryDirectoryPath = appendFile (temporaryDirectoryPath ++ "Log") (show logoutput) logEither _ _ = return () @@ -1462,7 +1462,7 @@ checks <- mapM checkTool additionaltools if not (null (lefts checks)) then return (Left (concat (lefts checks))) - else do + else do logMessage ("Tools : " ++ intercalate "," tools ++ "\n") temporaryDirectoryPath return (Right "Tools ok") @@ -1475,11 +1475,11 @@ --_ <- system ("clustalo --version >" ++ clustaloversionpath) _ <- system ("mlocarna --version >" ++ mlocarnaversionpath) _ <- system ("RNAfold --version >" ++ rnafoldversionpath) - _ <- system ("cmcalibrate -h >" ++ infernalversionpath) + _ <- system ("cmcalibrate -h >" ++ infernalversionpath) -- _ <- system ("RNAz" ++ rnazversionpath) -- _ <- system ("CMCompare >" ++ infernalversionpath) mlocarnaversion <- readFile mlocarnaversionpath - rnafoldversion <- readFile rnafoldversionpath + rnafoldversion <- readFile rnafoldversionpath infernalversionOutput <- readFile infernalversionpath let infernalversion = lines infernalversionOutput !! 1 if inputQuerySelectionMethod == "clustering" @@ -1499,7 +1499,7 @@ if isJust toolcheck then return (Right (fromJust toolcheck)) else return (Left ("RNAlien could not find "++ tool ++ " in your $PATH and has to abort.\n")) - + constructTaxonomyRecordsCSVTable :: ModelConstruction -> String constructTaxonomyRecordsCSVTable modelconstruction = csvtable where tableheader = "Taxonomy Id;Added in Iteration Step;Entry Header\n" @@ -1528,36 +1528,33 @@ evaluateConstructionResult staticOptions mCResult = do let evaluationDirectoryFilepath = tempDirPath staticOptions ++ "evaluation/" createDirectoryIfMissing False evaluationDirectoryFilepath - let fastaFilepath = tempDirPath staticOptions ++ "result.fa" - let clustalFilepath = evaluationDirectoryFilepath ++ "result.clustal" let reformatedClustalPath = evaluationDirectoryFilepath ++ "result.clustal.reformated" let cmFilepath = tempDirPath staticOptions ++ "result.cm" - let resultSequences = (inputFasta mCResult):map nucleotideSequence (concatMap sequenceRecords (taxRecords mCResult)) - let resultNumber = length resultSequences - let rnaCentralQueries = map RCH.buildSequenceViaMD5Query resultSequences + let resultSequences = inputFasta mCResult:map nucleotideSequence (concatMap sequenceRecords (taxRecords mCResult)) + let resultNumber = length resultSequences + let rnaCentralQueries = map RCH.buildSequenceViaMD5Query resultSequences rnaCentralEntries <- RCH.getRNACentralEntries rnaCentralQueries let rnaCentralEvaluationResult = RCH.showRNAcentralAlienEvaluation rnaCentralEntries writeFile (tempDirPath staticOptions ++ "result.rnacentral") rnaCentralEvaluationResult - systemCMalign ("--outformat=Clustal --cpu " ++ show (cpuThreads staticOptions)) cmFilepath fastaFilepath clustalFilepath let resultModelStatistics = tempDirPath staticOptions ++ "result.cmstat" systemCMstat cmFilepath resultModelStatistics inputcmStat <- readCMstat resultModelStatistics let cmstatString = cmstatEvalOutput inputcmStat if resultNumber > 1 - then do + then do let resultRNAz = tempDirPath staticOptions ++ "result.rnaz" let resultRNAcode = tempDirPath staticOptions ++ "result.rnacode" + let resultClustalFilepath = tempDirPath staticOptions ++ "result.clustal" let seqNumber = 6 :: Int let optimalIdentity = 80 :: Double let maximalIdentity = 99 :: Double let referenceSequence = True - --rnazClustalpath <- preprocessClustalForRNAcodeExternal clustalFilepath reformatedClustalPath seqNumber (truncate optimalIdentity) (truncate maximalIdentity) referenceSequence - preprocessingOutput <- preprocessClustalForRNAz clustalFilepath reformatedClustalPath seqNumber optimalIdentity maximalIdentity referenceSequence + preprocessingOutput <- preprocessClustalForRNAz resultClustalFilepath reformatedClustalPath seqNumber optimalIdentity maximalIdentity referenceSequence if isRight preprocessingOutput then do let rightPreprocessingOutput = fromRight preprocessingOutput let rnazClustalpath = snd rightPreprocessingOutput - systemRNAz "-l" rnazClustalpath resultRNAz + systemRNAz "-l" rnazClustalpath resultRNAz inputRNAz <- readRNAz resultRNAz let rnaZString = rnaZEvalOutput inputRNAz RC.systemRNAcode " -t " rnazClustalpath resultRNAcode @@ -1565,29 +1562,29 @@ let rnaCodeString = rnaCodeEvalOutput inputRNAcode return ("\nEvaluation of RNAlien result :\nCMstat statistics for result.cm\n" ++ cmstatString ++ "\nRNAz statistics for result alignment: " ++ rnaZString ++ "\nRNAcode output for result alignment:\n" ++ rnaCodeString ++ "\nSequences found by RNAlien with RNAcentral entry:\n" ++ rnaCentralEvaluationResult) else do - logWarning ("Running RNAz for result evalution encountered a problem:" ++ fromLeft preprocessingOutput) (tempDirPath staticOptions) + logWarning ("Running RNAz for result evalution encountered a problem:" ++ fromLeft preprocessingOutput) (tempDirPath staticOptions) return ("\nEvaluation of RNAlien result :\nCMstat statistics for result.cm\n" ++ cmstatString ++ "\nRNAz statistics for result alignment: Running RNAz for result evalution encountered a problem\n" ++ fromLeft preprocessingOutput ++ "\n" ++ "Sequences found by RNAlien with RNAcentral entry:\n" ++ rnaCentralEvaluationResult) else do - logWarning "Message: RNAlien could not find additional covariant sequences\n Could not run RNAz statistics. Could not run RNAz statistics with a single sequence.\n" (tempDirPath staticOptions) + logWarning "Message: RNAlien could not find additional covariant sequences\n Could not run RNAz statistics. Could not run RNAz statistics with a single sequence.\n" (tempDirPath staticOptions) return ("\nEvaluation of RNAlien result :\nCMstat statistics for result.cm\n" ++ cmstatString ++ "\nRNAlien could not find additional covariant sequences. Could not run RNAz statistics with a single sequence.\n\nSequences found by RNAlien with RNAcentral entry:\n" ++ rnaCentralEvaluationResult) -cmstatEvalOutput :: Either ParseError CMstat -> String +cmstatEvalOutput :: Either ParseError CMstat -> String cmstatEvalOutput inputcmstat | isRight inputcmstat = cmstatString | otherwise = show (fromLeft inputcmstat) - where cmStat = fromRight inputcmstat + where cmStat = fromRight inputcmstat cmstatString = " Sequence Number: " ++ show (statSequenceNumber cmStat)++ "\n" ++ " Effective Sequences: " ++ show (statEffectiveSequences cmStat)++ "\n" ++ " Consensus length: " ++ show (statConsensusLength cmStat) ++ "\n" ++ " Expected maximum hit-length: " ++ show (statW cmStat) ++ "\n" ++ " Basepairs: " ++ show (statBasepairs cmStat)++ "\n" ++ " Bifurcations: " ++ show (statBifurcations cmStat) ++ "\n" ++ " Modeltype: " ++ show (statModel cmStat) ++ "\n" ++ " Relative Entropy CM: " ++ show (relativeEntropyCM cmStat) ++ "\n" ++ " Relative Entropy HMM: " ++ show (relativeEntropyHMM cmStat) ++ "\n" -rnaZEvalOutput :: Either ParseError RNAz -> String -rnaZEvalOutput inputRNAz +rnaZEvalOutput :: Either ParseError RNAz -> String +rnaZEvalOutput inputRNAz | isRight inputRNAz = rnazString | otherwise = show (fromLeft inputRNAz) where rnaZ = fromRight inputRNAz rnazString = " Mean pairwise identity: " ++ show (meanPairwiseIdentity rnaZ) ++ "\n Shannon entropy: " ++ show (shannonEntropy rnaZ) ++ "\n GC content: " ++ show (gcContent rnaZ) ++ "\n Mean single sequence minimum free energy: " ++ show (meanSingleSequenceMinimumFreeEnergy rnaZ) ++ "\n Consensus minimum free energy: " ++ show (consensusMinimumFreeEnergy rnaZ) ++ "\n Energy contribution: " ++ show (energyContribution rnaZ) ++ "\n Covariance contribution: " ++ show (covarianceContribution rnaZ) ++ "\n Combinations pair: " ++ show (combinationsPair rnaZ) ++ "\n Mean z-score: " ++ show (meanZScore rnaZ) ++ "\n Structure conservation index: " ++ show (structureConservationIndex rnaZ) ++ "\n Background model: " ++ backgroundModel rnaZ ++ "\n Decision model: " ++ decisionModel rnaZ ++ "\n SVM decision value: " ++ show (svmDecisionValue rnaZ) ++ "\n SVM class propability: " ++ show (svmRNAClassProbability rnaZ) ++ "\n Prediction: " ++ prediction rnaZ -rnaCodeEvalOutput :: Either ParseError RC.RNAcode -> String -rnaCodeEvalOutput inputRNAcode +rnaCodeEvalOutput :: Either ParseError RC.RNAcode -> String +rnaCodeEvalOutput inputRNAcode | isRight inputRNAcode = rnaCodeString | otherwise = show (fromLeft inputRNAcode) where rnaCode = fromRight inputRNAcode @@ -1595,7 +1592,7 @@ rnaCodeEntries = concatMap showRNACodeHits (RC.rnacodeHits rnaCode) showRNACodeHits :: RC.RNAcodeHit -> String -showRNACodeHits rnacodeHit = show (RC.hss rnacodeHit) ++ "\t" ++ show (RC.frame rnacodeHit) ++ "\t" ++ show (RC.hitLength rnacodeHit) ++ "\t"++ show (RC.from rnacodeHit) ++ "\t" ++ show (RC.to rnacodeHit) ++ "\t" ++ (RC.name rnacodeHit) ++ "\t" ++ show (RC.start rnacodeHit) ++ "\t" ++ show (RC.end rnacodeHit) ++ "\t" ++ show (RC.score rnacodeHit) ++ show (RC.pvalue rnacodeHit) ++ "\n" +showRNACodeHits rnacodeHit = show (RC.hss rnacodeHit) ++ "\t" ++ show (RC.frame rnacodeHit) ++ "\t" ++ show (RC.hitLength rnacodeHit) ++ "\t"++ show (RC.from rnacodeHit) ++ "\t" ++ show (RC.to rnacodeHit) ++ "\t" ++ RC.name rnacodeHit ++ "\t" ++ show (RC.start rnacodeHit) ++ "\t" ++ show (RC.end rnacodeHit) ++ "\t" ++ show (RC.score rnacodeHit) ++ show (RC.pvalue rnacodeHit) ++ "\n" -- | Call for external preprocessClustalForRNAz preprocessClustalForRNAzExternal :: String -> String -> Int -> Int -> Int -> Bool -> IO (Either String (String,String)) @@ -1610,7 +1607,7 @@ let optimalIdentityOption = " -i " ++ show optimalIdentity ++ " " let maximalIdentityOption = " --max-id=" ++ show maximalIdenity ++ " " let referenceSequenceOption = if referenceSequence then " " else " -x " - let syscall = ("rnazSelectSeqs.pl " ++ reformatedClustalPath ++ " " ++ sequenceNumberOption ++ optimalIdentityOption ++ maximalIdentityOption ++ referenceSequenceOption ++ " >" ++ selectedClustalpath) + let syscall = "rnazSelectSeqs.pl " ++ reformatedClustalPath ++ " " ++ sequenceNumberOption ++ optimalIdentityOption ++ maximalIdentityOption ++ referenceSequenceOption ++ " >" ++ selectedClustalpath --putStrLn syscall system syscall selectedClustalText <- readFile selectedClustalpath @@ -1625,14 +1622,14 @@ let headerClustalTextLines = T.unlines (take 2 clustalTextLines) let headerlessClustalTextLines = T.unlines (drop 2 clustalTextLines) let reformatedClustalText = T.map reformatRNACodeAln headerlessClustalTextLines - TI.writeFile reformatedClustalPath (headerClustalTextLines `T.append` (T.singleton '\n') `T.append` reformatedClustalText) + TI.writeFile reformatedClustalPath (headerClustalTextLines `T.append` T.singleton '\n' `T.append` reformatedClustalText) --select representative entries from result.Clustal with select_sequences let selectedClustalpath = clustalFilepath ++ ".selected" let sequenceNumberOption = " -n " ++ show seqenceNumber ++ " " let optimalIdentityOption = " -i " ++ show optimalIdentity ++ " " let maximalIdentityOption = " --max-id=" ++ show maximalIdenity ++ " " let referenceSequenceOption = if referenceSequence then " " else " -x " - let syscall = ("rnazSelectSeqs.pl " ++ reformatedClustalPath ++ " " ++ sequenceNumberOption ++ optimalIdentityOption ++ maximalIdentityOption ++ referenceSequenceOption ++ " >" ++ selectedClustalpath) + let syscall = "rnazSelectSeqs.pl " ++ reformatedClustalPath ++ " " ++ sequenceNumberOption ++ optimalIdentityOption ++ maximalIdentityOption ++ referenceSequenceOption ++ " >" ++ selectedClustalpath --putStrLn syscall system syscall selectedClustalText <- readFile selectedClustalpath @@ -1641,14 +1638,14 @@ preprocessClustalForRNAz :: String -> String -> Int -> Double -> Double -> Bool -> IO (Either String (String,String)) preprocessClustalForRNAz clustalFilepath _ seqenceNumber optimalIdentity maximalIdenity referenceSequence = do clustalText <- TI.readFile clustalFilepath - let clustalTextLines = T.lines clustalText + let clustalTextLines = T.lines clustalText parsedClustalInput <- readClustalAlignment clustalFilepath - let selectedClustalpath = clustalFilepath ++ ".selected" + let selectedClustalpath = clustalFilepath ++ ".selected" if length clustalTextLines > 5 - then do + then if isRight parsedClustalInput then do - let (idMatrix,filteredClustalInput) = rnaCodeSelectSeqs2 (fromRight parsedClustalInput) seqenceNumber optimalIdentity maximalIdenity referenceSequence + let (idMatrix,filteredClustalInput) = rnaCodeSelectSeqs2 (fromRight parsedClustalInput) seqenceNumber optimalIdentity maximalIdenity referenceSequence writeFile selectedClustalpath (show filteredClustalInput) let formatedIdMatrix = show (fmap formatIdMatrix idMatrix) return (Right (formatedIdMatrix,selectedClustalpath)) @@ -1658,7 +1655,7 @@ let headerClustalTextLines = T.unlines (take 2 clustalLines) let headerlessClustalTextLines = T.unlines (drop 2 clustalLines) let reformatedClustalText = T.map reformatRNACodeAln headerlessClustalTextLines - TI.writeFile selectedClustalpath (headerClustalTextLines `T.append` (T.singleton '\n') `T.append` reformatedClustalText) + TI.writeFile selectedClustalpath (headerClustalTextLines `T.append` T.singleton '\n' `T.append` reformatedClustalText) return (Right ([],clustalFilepath)) formatIdMatrix :: Maybe (Int,Int,Double) -> String @@ -1672,7 +1669,7 @@ where entryVector = V.fromList (alignmentEntries currentClustalAlignment) entrySequences = V.map entryAlignedSequence entryVector entryReformatedSequences = V.map (T.map reformatRNACodeAln) entrySequences - totalSeqNumber = (V.length entryVector) + totalSeqNumber = V.length entryVector identityMatrix = computeSequenceIdentityMatrix entryReformatedSequences entryIdentityVector = V.map fromJust (V.filter isJust (getMatrixAsVector identityMatrix)) entryIdentities = V.toList entryIdentityVector @@ -1682,36 +1679,36 @@ prefilteredEntries = allEntries \\ entriesToDiscard --Optimize mean pairwise similarity (greedily) - remove worst sequence until desired number is reached costList = map (computeEntryCost optimalIdentity entryIdentityVector) prefilteredEntries - sortedCostList = (sortBy compareEntryCost2 costList) - sortedIndices = map fst sortedCostList + sortedCostList = sortBy compareEntryCost2 costList + sortedIndices = map fst sortedCostList --selectedEntryIndices = [1] ++ map fst (take (targetSeqNumber -1) sortedCostList) - selectedEntryIndices = selectEntryIndices referenceSequence targetSeqNumber sortedIndices + selectedEntryIndices = selectEntryIndices referenceSequence targetSeqNumber sortedIndices selectedEntries = map (\ind -> entryVector V.! (ind-1)) selectedEntryIndices selectedEntryHeader = map entrySequenceIdentifier selectedEntries reformatedSelectedEntryHeader = map (T.map reformatRNACodeId) selectedEntryHeader selectedEntrySequences = map (\ind -> entryReformatedSequences V.! (ind-1)) selectedEntryIndices --gapfreeEntrySequences = T.transpose (T.filter (\a -> not (T.all isGap a)) (T.transpose selectedEntrySequences)) - gapfreeEntrySequences = T.transpose (filter (\a -> not (T.all isGap a)) (T.transpose selectedEntrySequences)) - gapfreeEntries = map (\(a,b) -> ClustalAlignmentEntry a b)(zip reformatedSelectedEntryHeader gapfreeEntrySequences) + gapfreeEntrySequences = T.transpose (filter (not . T.all isGap) (T.transpose selectedEntrySequences)) + gapfreeEntries = map (uncurry ClustalAlignmentEntry)(zip reformatedSelectedEntryHeader gapfreeEntrySequences) emptyConservationTrack = setEmptyConservationTrack gapfreeEntries (conservationTrack currentClustalAlignment) newClustalAlignment = currentClustalAlignment {alignmentEntries = gapfreeEntries, conservationTrack = emptyConservationTrack} selectEntryIndices :: Bool -> Int -> [Int] -> [Int] selectEntryIndices referenceSequence targetSeqNumber sortedIndices - | referenceSequence = if elem (1 :: Int) firstX then firstX else 1:firstXm1 + | referenceSequence = if (1 :: Int) `elem` firstX then firstX else 1:firstXm1 | otherwise = firstX where firstXm1 = take (targetSeqNumber - 1) sortedIndices firstX = take targetSeqNumber sortedIndices - + setEmptyConservationTrack :: [ClustalAlignmentEntry] -> T.Text -> T.Text setEmptyConservationTrack alnentries currentConservationTrack - | null alnentries = currentConservationTrack - | otherwise = newConservationTrack + | null alnentries = currentConservationTrack + | otherwise = newConservationTrack where trackLength = T.length (entryAlignedSequence (head alnentries)) newConservationTrack = T.replicate (trackLength + 0) (T.pack " ") - + isGap :: Char -> Bool -isGap a +isGap a | a == '-' = True | otherwise = False @@ -1721,7 +1718,7 @@ entryIdentities = getEntryIdentities currentIndex allIdentities getEntryIdentities :: Int -> V.Vector (Int,Int,Double) -> V.Vector (Int,Int,Double) -getEntryIdentities currentIndex allIdentities = (V.filter (isIIdx currentIndex) allIdentities) V.++ (V.filter (isJIdx currentIndex) allIdentities) +getEntryIdentities currentIndex allIdentities = V.filter (isIIdx currentIndex) allIdentities V.++ V.filter (isJIdx currentIndex) allIdentities isIIdx :: Int -> (Int,Int,Double) -> Bool isIIdx currentIdx (i,_,_) = currentIdx == i @@ -1731,13 +1728,13 @@ computeCost :: Double -> (Int,Int,Double) -> Double computeCost optimalIdentity (_,_,c) = (c - optimalIdentity) * (c - optimalIdentity) -compareEntryCost2 :: (Int, Double) -> (Int, Double) -> Ordering +compareEntryCost2 :: (Int, Double) -> (Int, Double) -> Ordering compareEntryCost2 (_,costA) (_,costB) = compare costA costB -- TODO change to vector preFilterIdentityMatrix :: Double -> Int -> Int-> [Int] -> [(Int,Int,Double)] -> [Int] preFilterIdentityMatrix identityCutoff minSeqNumber totalSeqNumber filteredIds entryIdentities - | (totalSeqNumber - (length filteredIds)) <= minSeqNumber = [] + | (totalSeqNumber - length filteredIds) <= minSeqNumber = [] | identityCutoff == (100 :: Double) = [] | Prelude.null entryIdentities = [] | otherwise = entryresult ++ preFilterIdentityMatrix identityCutoff minSeqNumber totalSeqNumber (filteredIds ++ entryresult) (tail entryIdentities) @@ -1746,46 +1743,46 @@ checkIdentityEntry :: Double -> [Int] -> (Int,Int,Double) -> [Int] checkIdentityEntry identityCutoff filteredIds (i,j,ident) - | elem i filteredIds = [] - | elem j filteredIds = [] + | i `elem` filteredIds = [] + | j `elem` filteredIds = [] | ident > identityCutoff = [j] | otherwise = [] - + computeSequenceIdentityMatrix :: V.Vector T.Text -> Matrix (Maybe (Int,Int,Double)) computeSequenceIdentityMatrix entryVector = matrix (V.length entryVector) (V.length entryVector) (computeSequenceIdentityEntry entryVector) -- Computes Sequence identity once for each pair and not vs itself computeSequenceIdentityEntry :: V.Vector T.Text -> (Int,Int) -> Maybe (Int,Int,Double) computeSequenceIdentityEntry entryVector (row,col) - | i < j = Just $ (row,col,ident) + | i < j = Just (row,col,ident) | otherwise = Nothing where i=row-1 j=col-1 --gaps in both sequences need to be removed, because they count as match - ientry = (entryVector V.! i) - jentry = (entryVector V.! j) + ientry = entryVector V.! i + jentry = entryVector V.! j (gfi,gfj) = unzip (filter notDoubleGap (T.zip ientry jentry)) gfitext = T.pack gfi - gfjtext = T.pack gfj + gfjtext = T.pack gfj --ident=stringIdentity gfi gfj ident=textIdentity gfitext gfjtext - + notDoubleGap :: (Char,Char) -> Bool notDoubleGap (a,b) | a == '-' && b == '-' = False | otherwise = True -reformatRNACodeId :: Char -> Char +reformatRNACodeId :: Char -> Char reformatRNACodeId c | c == ':' = '-' | c == '|' = '-' | c == '.' = '-' | c == '~' = '-' | c == '_' = '-' - | c == '/' = '-' + | c == '/' = '-' | otherwise = c - -reformatRNACodeAln :: Char -> Char + +reformatRNACodeAln :: Char -> Char reformatRNACodeAln c | c == ':' = '-' | c == '|' = '-' @@ -1799,7 +1796,7 @@ | c == 'a' = 'A' | otherwise = c -reformatAln :: Char -> Char +reformatAln :: Char -> Char reformatAln c | c == '.' = '-' | c == '~' = '-' @@ -1819,7 +1816,7 @@ response <- N.httpLbs req manager let sta = N.responseStatus response if statusIsSuccessful sta - then return (Right ("Network connection with NCBI server was successful")) + then return (Right "Network connection with NCBI server was successful") else return (Left ("Could not connect to NCBI server \"https://www.ncbi.nlm.nih.gov\". Response Code: " ++ show (statusCode sta) ++ " \n" ++ B.unpack (statusMessage sta))) -- | Blast evalue is set stricter in inital alignment mode @@ -1832,7 +1829,7 @@ reformatFasta input = Seq (seqheader input) updatedSequence Nothing where updatedSequence = SeqData (L.pack (map reformatFastaSequence (L.unpack (unSD (seqdata input))))) -reformatFastaSequence :: Char -> Char +reformatFastaSequence :: Char -> Char reformatFastaSequence c | c == '.' = '-' | c == '~' = '-' @@ -1846,7 +1843,7 @@ | otherwise = c setRestrictedTaxonomyLimits :: String -> (Maybe Int,Maybe Int) -setRestrictedTaxonomyLimits trestriction +setRestrictedTaxonomyLimits trestriction | trestriction == "bacteria" = (Just (2 :: Int), Nothing) | trestriction == "archea" = (Just (2157 :: Int), Nothing) | trestriction == "eukaryia" = (Just (2759 :: Int), Nothing) @@ -1871,22 +1868,22 @@ --let splitIds = map E.encodeUtf8 (TL.splitOn (TL.pack ",") (TL.pack sequenceIds)) let filteredEntries = concatMap (filterSequencesById alnEntries) sequenceIds return filteredEntries - + extractAlignmentSequences :: TL.Text -> [Sequence] extractAlignmentSequences seedFamilyAln = rfamIDAndseedFamilySequences where seedFamilyAlnLines = TL.lines seedFamilyAln -- remove empty lines from splitting - seedFamilyNonEmpty = filter (\alnline -> not (TL.empty == alnline)) seedFamilyAlnLines + seedFamilyNonEmpty = filter (\alnline -> TL.empty /= alnline) seedFamilyAlnLines -- remove annotation and spacer lines - seedFamilyIdSeqLines = filter (\alnline -> ((not ((TL.head alnline) == '#'))) && (not ((TL.head alnline) == ' ')) && (not ((TL.head alnline) == '/'))) seedFamilyNonEmpty + seedFamilyIdSeqLines = filter (\alnline -> not ((TL.head alnline) == '#') && not ((TL.head alnline) == ' ') && not ((TL.head alnline) == '/')) seedFamilyNonEmpty -- put id and corresponding seq of each line into a list and remove whitspaces seedFamilyIdandSeqLines = map TL.words seedFamilyIdSeqLines -- linewise tuples with id and seq without alinment characters - . - seedFamilyIdandSeqLineTuples = map (\alnline -> ((head alnline),(filterAlnChars (last alnline)))) seedFamilyIdandSeqLines + seedFamilyIdandSeqLineTuples = map (\alnline -> (head alnline,filterAlnChars (last alnline))) seedFamilyIdandSeqLines -- line tuples sorted by id seedFamilyIdandSeqTupleSorted = sortBy (\tuple1 tuple2 -> compare (fst tuple1) (fst tuple2)) seedFamilyIdandSeqLineTuples -- line tuples grouped by id - seedFamilyIdandSeqTupleGroups = groupBy (\tuple1 tuple2 -> (fst tuple1) == (fst tuple2)) seedFamilyIdandSeqTupleSorted + seedFamilyIdandSeqTupleGroups = groupBy (\tuple1 tuple2 -> fst tuple1 == fst tuple2) seedFamilyIdandSeqTupleSorted seedFamilySequences = map mergeIdSeqTuplestoSequence seedFamilyIdandSeqTupleGroups rfamIDAndseedFamilySequences = seedFamilySequences @@ -1894,10 +1891,10 @@ filterSequencesById alignmentSequences sequenceId = filter (sequenceHasId sequenceId) alignmentSequences sequenceHasId :: L.ByteString -> Sequence -> Bool -sequenceHasId sequenceId currentSequence = sequenceId == (unSL (seqid currentSequence)) +sequenceHasId sequenceId currentSequence = sequenceId == unSL (seqid currentSequence) filterAlnChars :: TL.Text -> TL.Text -filterAlnChars cs = TL.filter (\c -> (not (c == '-')) && (not (c == '.'))) cs +filterAlnChars cs = TL.filter (\c -> not (c == '-') && not (c == '.')) cs mergeIdSeqTuplestoSequence :: [(TL.Text,TL.Text)] -> Sequence mergeIdSeqTuplestoSequence tuplelist = currentSequence diff --git a/src/Bio/RNAlienStatistics.hs b/src/Bio/RNAlienStatistics.hs --- a/src/Bio/RNAlienStatistics.hs +++ b/src/Bio/RNAlienStatistics.hs @@ -4,14 +4,14 @@ -- | Statistics for RNAlien Results -- dist/build/RNAlienStatistics/RNAlienStatistics -s bitscore -i /scratch/egg/temp/cm13676/1/model.cm -r /home/mescalin/egg/current/Data/AlienTest/cms/BsrG.cm -g /scratch/egg/temp/AlienSearch/genomes/ -o /scratch/egg/temp/AlienStatistics module Main where - -import System.Console.CmdArgs + +import System.Console.CmdArgs import Data.Either.Unwrap import System.Process import qualified Data.ByteString.Lazy.Char8 as L import Bio.RNAlienLibrary import System.Directory -import Bio.Core.Sequence +import Bio.Core.Sequence import Bio.Sequence.Fasta import Data.List import qualified System.FilePath as FP @@ -20,7 +20,7 @@ import Bio.RNAzParser import qualified Bio.RNAcodeParser as RC -data Options = Options +data Options = Options { alienCovarianceModelPath :: String, alienrnazPath :: String, alienrnacodePath :: String, @@ -29,7 +29,7 @@ rfamFastaFilePath :: String, alienFastaFilePath :: String, rfamModelName :: String, - rfamModelId :: String, + rfamModelId :: String, rfamThreshold :: Double, alienThreshold :: Double, databaseSize :: Maybe Double, @@ -49,7 +49,7 @@ rfamCovarianceModelPath = def &= name "r" &= help "Path to rfamCovarianceModelPath", rfamFastaFilePath = def &= name "g" &= help "Path to rfamFastaFile", rfamModelName = def &= name "n" &= help "Rfam model name", - rfamModelId = def &= name "d" &= help "Rfam model id", + rfamModelId = def &= name "d" &= help "Rfam model id", alienFastaFilePath = def &= name "a" &= help "Path to alienFastaFile", outputDirectoryPath = def &= name "o" &= help "Path to output directory", alienThreshold = 20 &= name "t" &= help "Bitscore threshold for RNAlien model hits on Rfam fasta, default 20", @@ -59,16 +59,16 @@ thresholdSelection = "bitscore" &= name "s" &= help "Selection method, (bitscore, evalue), default bitscore", linkScores = False &= name "l" &= help "Triggers computation of linkscores via CMCompare", threads = 1 &= name "c" &= help "Number of available cpu slots/cores, default 1" - } &= summary "RNAlienStatistics" &= help "Florian Eggenhofer - >2013" &= verbosity + } &= summary "RNAlienStatistics" &= help "Florian Eggenhofer - >2013" &= verbosity --cmSearchFasta threads rfamCovarianceModelPath outputDirectoryPath "Rfam" False genomesDirectoryPath cmSearchFasta :: Int -> String -> Double -> Maybe Double -> Int -> String -> String -> String -> String -> IO [CMsearchHit] cmSearchFasta benchmarkIndex thresholdSelection thresholdScore databaseSize cpuThreads covarianceModelPath outputDirectory modelType fastapath = do createDirectoryIfMissing False (outputDirectory ++ "/" ++ modelType) - _ <- systemCMsearch cpuThreads (maybe "" (\dbs -> " -Z " ++ show dbs ++ " ") databaseSize) covarianceModelPath fastapath (outputDirectory ++ "/" ++ modelType ++ "/" ++ (show benchmarkIndex) ++ ".cmsearch") + _ <- systemCMsearch cpuThreads (maybe "" (\dbs -> " -Z " ++ show dbs ++ " ") databaseSize) covarianceModelPath fastapath (outputDirectory ++ "/" ++ modelType ++ "/" ++ show benchmarkIndex ++ ".cmsearch") --_ <- systemCMsearch cpuThreads " " covarianceModelPath fastapath (outputDirectory ++ "/" ++ modelType ++ "/" ++ (show benchmarkIndex) ++ ".cmsearch") - result <- readCMSearch (outputDirectory ++ "/" ++ modelType ++ "/" ++ (show benchmarkIndex) ++ ".cmsearch") - if (isLeft result) + result <- readCMSearch (outputDirectory ++ "/" ++ modelType ++ "/" ++ show benchmarkIndex ++ ".cmsearch") + if isLeft result then do print (fromLeft result) return [] @@ -82,10 +82,10 @@ cmSearchesFasta :: Int -> String -> Double -> Maybe Double -> Int -> String -> String -> String -> String -> IO [CMsearchHit] cmSearchesFasta benchmarkIndex thresholdSelection thresholdScore databaseSize cpuThreads covarianceModelPath outputDirectory modelType fastapath = do createDirectoryIfMissing False (outputDirectory ++ "/" ++ modelType) - _ <- systemCMsearch cpuThreads (maybe "" (\dbs -> " -Z " ++ show dbs ++ " ") databaseSize) covarianceModelPath fastapath (outputDirectory ++ "/" ++ modelType ++ "/" ++ (show benchmarkIndex) ++ ".cmsearch") + _ <- systemCMsearch cpuThreads (maybe "" (\dbs -> " -Z " ++ show dbs ++ " ") databaseSize) covarianceModelPath fastapath (outputDirectory ++ "/" ++ modelType ++ "/" ++ show benchmarkIndex ++ ".cmsearch") --_ <- systemCMsearch cpuThreads " " covarianceModelPath fastapath (outputDirectory ++ "/" ++ modelType ++ "/" ++ (show benchmarkIndex) ++ ".cmsearch") - result <- readCMSearches (outputDirectory ++ "/" ++ modelType ++ "/" ++ (show benchmarkIndex) ++ ".cmsearch") - if (isLeft result) + result <- readCMSearches (outputDirectory ++ "/" ++ modelType ++ "/" ++ show benchmarkIndex ++ ".cmsearch") + if isLeft result then do print (fromLeft result) return [] @@ -119,37 +119,37 @@ fastaSequences <- readFasta fastaInputPath let trimmedSequence = trimCMsearchSequence cmsearch (head fastaSequences) writeFasta fastaOutputPath [trimmedSequence] - + trimCMsearchSequence :: CMsearch -> Sequence -> Sequence trimCMsearchSequence cmSearchResult inputSequence = subSequence where hitScoreEntry = head (cmsearchHits cmSearchResult) sequenceString = L.unpack (unSD (seqdata inputSequence)) sequenceSubstring = cmSearchsubString (hitStart hitScoreEntry) (hitEnd hitScoreEntry) sequenceString - newSequenceHeader = L.pack ((L.unpack (unSL (seqheader inputSequence))) ++ "cmS_" ++ (show (hitStart hitScoreEntry)) ++ "_" ++ (show (hitEnd hitScoreEntry)) ++ "_" ++ (show (hitStrand hitScoreEntry))) - subSequence = Seq (SeqLabel newSequenceHeader) (SeqData (L.pack sequenceSubstring)) Nothing + newSequenceHeader = L.pack (L.unpack (unSL (seqheader inputSequence)) ++ "cmS_" ++ show (hitStart hitScoreEntry) ++ "_" ++ show (hitEnd hitScoreEntry) ++ "_" ++ show (hitStrand hitScoreEntry)) + subSequence = Seq (SeqLabel newSequenceHeader) (SeqData (L.pack sequenceSubstring)) Nothing --With paralogs allowed cmSearchSameHit :: CMsearchHit -> CMsearchHit -> Bool cmSearchSameHit hitscore1 hitscore2 | unpackedSeqHeader1 == unpackedSeqHeader2 = True | otherwise = False - where unpackedSeqHeader1 = (L.unpack (hitSequenceHeader hitscore1)) - unpackedSeqHeader2 = (L.unpack (hitSequenceHeader hitscore2)) + where unpackedSeqHeader1 = L.unpack (hitSequenceHeader hitscore1) + unpackedSeqHeader2 = L.unpack (hitSequenceHeader hitscore2) cmSearchSameOrganism :: CMsearchHit -> CMsearchHit -> Bool cmSearchSameOrganism hitscore1 hitscore2 | hitOrganism1 == hitOrganism2 = True | otherwise = False - where unpackedSeqHeader1 = (L.unpack (hitSequenceHeader hitscore1)) - unpackedSeqHeader2 = (L.unpack (hitSequenceHeader hitscore2)) + where unpackedSeqHeader1 = L.unpack (hitSequenceHeader hitscore1) + unpackedSeqHeader2 = L.unpack (hitSequenceHeader hitscore2) separationcharacter1 = selectSeparationChar unpackedSeqHeader1 separationcharacter2 = selectSeparationChar unpackedSeqHeader2 - hitOrganism1 = (DS.splitOn separationcharacter1 unpackedSeqHeader1) !! 0 - hitOrganism2 = (DS.splitOn separationcharacter2 unpackedSeqHeader2) !! 0 + hitOrganism1 = head (DS.splitOn separationcharacter1 unpackedSeqHeader1) + hitOrganism2 = head (DS.splitOn separationcharacter2 unpackedSeqHeader2) selectSeparationChar :: String -> String selectSeparationChar inputString - | any (\a -> a == ':') inputString = ":" + | any ((== ':')) inputString = ":" | otherwise = "/" main :: IO () @@ -171,61 +171,61 @@ _ <- system ("cat " ++ rfamFastaFilePath ++ " | grep '>' | wc -l >" ++ outputDirectoryPath ++ FP.takeFileName rfamFastaFilePath ++ ".entries") _ <- system ("cat " ++ alienFastaFilePath ++ " | grep '>' | wc -l >" ++ outputDirectoryPath ++ FP.takeFileName alienFastaFilePath ++ ".entries") rfamFastaEntries <- readFile (outputDirectoryPath ++ FP.takeFileName rfamFastaFilePath ++ ".entries") - alienFastaEntries <- readFile (outputDirectoryPath ++ FP.takeFileName alienFastaFilePath ++ ".entries") + alienFastaEntries <- readFile (outputDirectoryPath ++ FP.takeFileName alienFastaFilePath ++ ".entries") let rfamFastaEntriesNumber = read rfamFastaEntries :: Int let alienFastaEntriesNumber = read alienFastaEntries :: Int - rfamonAlienResults <- cmSearchesFasta benchmarkIndex thresholdSelection rfamThreshold databaseSize threads rfamCovarianceModelPath outputDirectoryPath "rfamOnAlien" alienFastaFilePath - alienonRfamResults <- cmSearchFasta benchmarkIndex thresholdSelection alienThreshold databaseSize threads alienCovarianceModelPath outputDirectoryPath "alienOnRfam" rfamFastaFilePath + rfamonAlienResults <- cmSearchesFasta benchmarkIndex thresholdSelection rfamThreshold databaseSize threads rfamCovarianceModelPath outputDirectoryPath "rfamOnAlien" alienFastaFilePath + alienonRfamResults <- cmSearchFasta benchmarkIndex thresholdSelection alienThreshold databaseSize threads alienCovarianceModelPath outputDirectoryPath "alienOnRfam" rfamFastaFilePath let rfamonAlienResultsNumber = length rfamonAlienResults let alienonRfamResultsNumber = length alienonRfamResults let rfamonAlienRecovery = (fromIntegral rfamonAlienResultsNumber :: Double) / (fromIntegral alienFastaEntriesNumber :: Double) let alienonRfamRecovery = (fromIntegral alienonRfamResultsNumber :: Double) / (fromIntegral rfamFastaEntriesNumber :: Double) - if (verbose == Loud) + if verbose == Loud then do putStrLn ("BenchmarkIndex: " ++ show benchmarkIndex) putStrLn ("RfamModelName: " ++ rfamModelName) putStrLn ("RfamModelId: " ++ rfamModelId) - putStrLn ("Linkscore: " ++ (either id show linkscore)) - putStrLn ("rfamMaxLinkScore: " ++ (either id show rfamMaxLinkScore)) - putStrLn ("alienMaxLinkscore: " ++ (either id show alienMaxLinkscore)) + putStrLn ("Linkscore: " ++ either id show linkscore) + putStrLn ("rfamMaxLinkScore: " ++ either id show rfamMaxLinkScore) + putStrLn ("alienMaxLinkscore: " ++ either id show alienMaxLinkscore) putStrLn ("rfamGatheringThreshold: " ++ show rfamThreshold) - putStrLn ("alienGatheringThreshold: " ++ show alienThreshold) + putStrLn ("alienGatheringThreshold: " ++ show alienThreshold) putStrLn ("rfamFastaEntriesNumber: " ++ show rfamFastaEntriesNumber) - putStrLn ("alienFastaEntriesNumber: " ++ show alienFastaEntriesNumber) + putStrLn ("alienFastaEntriesNumber: " ++ show alienFastaEntriesNumber) putStrLn ("rfamonAlienResultsNumber: " ++ show rfamonAlienResultsNumber) putStrLn ("alienonRfamResultsNumber: " ++ show alienonRfamResultsNumber) - putStrLn ("RfamonAlienRecovery: " ++ show rfamonAlienRecovery) + putStrLn ("RfamonAlienRecovery: " ++ show rfamonAlienRecovery) putStrLn ("AlienonRfamRecovery: " ++ show alienonRfamRecovery) print rnazString print rnacodeString print cmStatString - else do + else putStrLn (show benchmarkIndex ++ "\t" ++ rfamModelName ++ "\t" ++ rfamModelId ++ "\t" ++ (either id show linkscore) ++ "\t" ++ (either id show rfamMaxLinkScore) ++ "\t" ++ (either id show alienMaxLinkscore) ++ "\t" ++ show rfamThreshold ++ "\t" ++ show alienThreshold ++ "\t" ++ show rfamFastaEntriesNumber ++ "\t" ++ show alienFastaEntriesNumber ++ "\t" ++ show rfamonAlienResultsNumber ++ "\t" ++ show alienonRfamResultsNumber ++ "\t" ++ printf "%.2f" rfamonAlienRecovery ++ "\t" ++ printf "%.2f" alienonRfamRecovery ++ "\t" ++ rnazString ++ "\t" ++ rnacodeString ++ "\t" ++ cmStatString) else do --compute linkscore alienMaxLinkscore <- if linkScores then compareCM alienCovarianceModelPath alienCovarianceModelPath outputDirectoryPath else return ( Left "-") _ <- system ("cat " ++ alienFastaFilePath ++ " | grep '>' | wc -l >" ++ outputDirectoryPath ++ FP.takeFileName alienFastaFilePath ++ ".entries") - alienFastaEntries <- readFile (outputDirectoryPath ++ FP.takeFileName alienFastaFilePath ++ ".entries") - let alienFastaEntriesNumber = read alienFastaEntries :: Int - if (verbose == Loud) + alienFastaEntries <- readFile (outputDirectoryPath ++ FP.takeFileName alienFastaFilePath ++ ".entries") + let alienFastaEntriesNumber = read alienFastaEntries :: Int + if verbose == Loud then do - putStrLn ("BenchmarkIndex:") - putStrLn ("RfamModelName: -") - putStrLn ("RfamModelId: -") - putStrLn ("Linkscore: -") - putStrLn ("rfamMaxLinkScore: -") - putStrLn ("alienMaxLinkscore: " ++ (either id show alienMaxLinkscore)) - putStrLn ("rfamGatheringThreshold: -") - putStrLn ("alienGatheringThreshold: -") - putStrLn ("rfamFastaEntriesNumber: -") - putStrLn ("alienFastaEntriesNumber: " ++ show alienFastaEntriesNumber) - putStrLn ("rfamonAlienResultsNumber: -") - putStrLn ("alienonRfamResultsNumber: -") - putStrLn ("RfamonAlienRecovery: -") - putStrLn ("AlienonRfamRecovery: -") + putStrLn "BenchmarkIndex:" + putStrLn "RfamModelName: -" + putStrLn "RfamModelId: -" + putStrLn "Linkscore: -" + putStrLn "rfamMaxLinkScore: -" + putStrLn ("alienMaxLinkscore: " ++ either id show alienMaxLinkscore) + putStrLn "rfamGatheringThreshold: -" + putStrLn "alienGatheringThreshold: -" + putStrLn "rfamFastaEntriesNumber: -" + putStrLn ("alienFastaEntriesNumber: " ++ show alienFastaEntriesNumber) + putStrLn "rfamonAlienResultsNumber: -" + putStrLn "alienonRfamResultsNumber: -" + putStrLn "RfamonAlienRecovery: -" + putStrLn "AlienonRfamRecovery: -" print rnazString print cmStatString - else do + else putStrLn (show benchmarkIndex ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ (either id show alienMaxLinkscore) ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ show alienFastaEntriesNumber ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ rnazString ++ "\t" ++ rnacodeString ++ "\t" ++ cmStatString) rnazOutput :: Verbosity -> String -> IO String @@ -237,14 +237,14 @@ if isRight inputRNAz then do let rnaZ = fromRight inputRNAz - if (verbose == Loud) + if verbose == Loud then do - let output = "Mean pairwise identity: " ++ show (meanPairwiseIdentity rnaZ) ++ "\n Shannon entropy: " ++ show (shannonEntropy rnaZ) ++ "\n GC content: " ++ show (gcContent rnaZ) ++ "\n Mean single sequence minimum free energy: " ++ show (meanSingleSequenceMinimumFreeEnergy rnaZ) ++ "\n Consensus minimum free energy: " ++ show (consensusMinimumFreeEnergy rnaZ) ++ "\n Energy contribution: " ++ show (energyContribution rnaZ) ++ "\n Covariance contribution: " ++ show (covarianceContribution rnaZ) ++ "\n Combinations pair: " ++ show (combinationsPair rnaZ) ++ "\n Mean z-score: " ++ show (meanZScore rnaZ) ++ "\n Structure conservation index: " ++ show (structureConservationIndex rnaZ) ++ "\n Background model: " ++ backgroundModel rnaZ ++ "\n Decision model: " ++ decisionModel rnaZ ++ "\n SVM decision value: " ++ show (svmDecisionValue rnaZ) ++ "\n SVM class propability: " ++ show (svmRNAClassProbability rnaZ) ++ "\n Prediction: " ++ (prediction rnaZ) + let output = "Mean pairwise identity: " ++ show (meanPairwiseIdentity rnaZ) ++ "\n Shannon entropy: " ++ show (shannonEntropy rnaZ) ++ "\n GC content: " ++ show (gcContent rnaZ) ++ "\n Mean single sequence minimum free energy: " ++ show (meanSingleSequenceMinimumFreeEnergy rnaZ) ++ "\n Consensus minimum free energy: " ++ show (consensusMinimumFreeEnergy rnaZ) ++ "\n Energy contribution: " ++ show (energyContribution rnaZ) ++ "\n Covariance contribution: " ++ show (covarianceContribution rnaZ) ++ "\n Combinations pair: " ++ show (combinationsPair rnaZ) ++ "\n Mean z-score: " ++ show (meanZScore rnaZ) ++ "\n Structure conservation index: " ++ show (structureConservationIndex rnaZ) ++ "\n Background model: " ++ backgroundModel rnaZ ++ "\n Decision model: " ++ decisionModel rnaZ ++ "\n SVM decision value: " ++ show (svmDecisionValue rnaZ) ++ "\n SVM class propability: " ++ show (svmRNAClassProbability rnaZ) ++ "\n Prediction: " ++ prediction rnaZ return output else do - let output = show (meanPairwiseIdentity rnaZ) ++ "\t" ++ show (shannonEntropy rnaZ) ++ "\t" ++ show (gcContent rnaZ) ++ "\t" ++ show (meanSingleSequenceMinimumFreeEnergy rnaZ) ++ "\t" ++ show (consensusMinimumFreeEnergy rnaZ) ++ "\t" ++ show (energyContribution rnaZ) ++ "\t" ++ show (covarianceContribution rnaZ) ++ "\t" ++ show (combinationsPair rnaZ) ++ "\t" ++ show (meanZScore rnaZ) ++ "\t" ++ show (structureConservationIndex rnaZ) ++ "\t" ++ show (svmDecisionValue rnaZ) ++ "\t" ++ show (svmRNAClassProbability rnaZ) ++ "\t" ++ (prediction rnaZ) + let output = show (meanPairwiseIdentity rnaZ) ++ "\t" ++ show (shannonEntropy rnaZ) ++ "\t" ++ show (gcContent rnaZ) ++ "\t" ++ show (meanSingleSequenceMinimumFreeEnergy rnaZ) ++ "\t" ++ show (consensusMinimumFreeEnergy rnaZ) ++ "\t" ++ show (energyContribution rnaZ) ++ "\t" ++ show (covarianceContribution rnaZ) ++ "\t" ++ show (combinationsPair rnaZ) ++ "\t" ++ show (meanZScore rnaZ) ++ "\t" ++ show (structureConservationIndex rnaZ) ++ "\t" ++ show (svmDecisionValue rnaZ) ++ "\t" ++ show (svmRNAClassProbability rnaZ) ++ "\t" ++ prediction rnaZ return output - else do + else if (verbose == Loud) then do let output = "Mean pairwise identity: " ++ " - \n Shannon entropy: " ++ " - \n GC content: " ++ " - \n Mean single sequence minimum free energy: " ++ " - \n Consensus minimum free energy: " ++ " - \n Energy contribution: " ++ " - \n Covariance contribution: " ++ " - \n Combinations pair: " ++ " - \n Mean z-score: " ++ " - \n Structure conservation index: " ++ " - \n Background model: " ++ " - \n Decision model: " ++ " - \n SVM decision value: " ++ " - \n SVM class propability: " ++ " - \n Prediction: " ++ " - \n" @@ -252,7 +252,7 @@ else do let output = "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" return output - else do + else if (verbose == Loud) then do let output = "Mean pairwise identity: " ++ " - \n Shannon entropy: " ++ " - \n GC content: " ++ " - \n Mean single sequence minimum free energy: " ++ " - \n Consensus minimum free energy: " ++ " - \n Energy contribution: " ++ " - \n Covariance contribution: " ++ " - \n Combinations pair: " ++ " - \n Mean z-score: " ++ " - \n Structure conservation index: " ++ " - \n Background model: " ++ " - \n Decision model: " ++ " - \n SVM decision value: " ++ " - \n SVM class propability: " ++ " - \n Prediction: " ++ " - \n" @@ -270,14 +270,14 @@ if isRight inputCMstat then do let cmStat = fromRight inputCMstat - if (verbose == Loud) + if verbose == Loud then do - let output = "statSequenceNumber: " ++ (show (statSequenceNumber cmStat)) ++ "\nstatEffectiveSequences: " ++ (show (statEffectiveSequences cmStat)) ++ "\nstatConsensusLength: " ++ (show (statConsensusLength cmStat)) ++ "\nstatW: " ++ show (statW cmStat) ++ "\nstatBasepairs: " ++ show (statBasepairs cmStat) ++ "\nstatBifurcations: " ++ (show (statBifurcations cmStat)) ++ "\nstatModel: " ++ (statModel cmStat) ++ "\nrelativeEntropyCM: " ++ show (relativeEntropyCM cmStat) ++ "\nrelativeEntropyHMM: " ++ show (relativeEntropyHMM cmStat) + let output = "statSequenceNumber: " ++ show (statSequenceNumber cmStat) ++ "\nstatEffectiveSequences: " ++ show (statEffectiveSequences cmStat) ++ "\nstatConsensusLength: " ++ show (statConsensusLength cmStat) ++ "\nstatW: " ++ show (statW cmStat) ++ "\nstatBasepairs: " ++ show (statBasepairs cmStat) ++ "\nstatBifurcations: " ++ show (statBifurcations cmStat) ++ "\nstatModel: " ++ statModel cmStat ++ "\nrelativeEntropyCM: " ++ show (relativeEntropyCM cmStat) ++ "\nrelativeEntropyHMM: " ++ show (relativeEntropyHMM cmStat) return output else do - let output = (show (statSequenceNumber cmStat)) ++ "\t" ++ (show (statEffectiveSequences cmStat)) ++ "\t" ++ (show (statConsensusLength cmStat)) ++ "\t" ++ show (statW cmStat) ++ "\t" ++ show (statBasepairs cmStat) ++ "\t" ++ (show (statBifurcations cmStat)) ++ "\t" ++ (statModel cmStat) ++ "\t" ++ show (relativeEntropyCM cmStat) ++ "\t" ++ show (relativeEntropyHMM cmStat) + let output = show (statSequenceNumber cmStat) ++ "\t" ++ show (statEffectiveSequences cmStat) ++ "\t" ++ show (statConsensusLength cmStat) ++ "\t" ++ show (statW cmStat) ++ "\t" ++ show (statBasepairs cmStat) ++ "\t" ++ show (statBifurcations cmStat) ++ "\t" ++ statModel cmStat ++ "\t" ++ show (relativeEntropyCM cmStat) ++ "\t" ++ show (relativeEntropyHMM cmStat) return output - else do + else if (verbose == Loud) then do let output = "statSequenceNumber: -" ++ "\nstatEffectiveSequences: -" ++ "\nstatConsensusLength: -" ++ "\nstatW: -" ++ "\nstatBasepairs: -" ++ "\nstatBifurcations: -" ++ "\nstatModel: -" ++ "\nrelativeEntropyCM: -" ++ "\nrelativeEntropyHMM: -" @@ -285,7 +285,7 @@ else do let output = "-\t" ++ "-\t" ++ "-\t" ++ "-\t" ++ "-\t" ++ "-\t" ++ "-\t" ++ "-\t" ++ "-" return output - else do + else if (verbose == Loud) then do let output = "statSequenceNumber: -" ++ "\nstatEffectiveSequences: -" ++ "\nstatConsensusLength: -" ++ "\nstatW: -" ++ "\nstatBasepairs: -" ++ "\nstatBifurcations: -" ++ "\nstatModel: -" ++ "\nrelativeEntropyCM: -" ++ "\nrelativeEntropyHMM: -" @@ -305,14 +305,14 @@ let rnaCode = fromRight inputRNACode let lowestPvalue = minimum (map RC.pvalue (RC.rnacodeHits rnaCode)) let rnaCodeClassification = if lowestPvalue < 0.05 then "PROTEIN" else "OTHER" - if (verbose == Loud) - then do - let output = "RNAcode lowest p-value: " ++ (show lowestPvalue) ++ "\nrnaCodeClassification: " ++ rnaCodeClassification + if verbose == Loud + then do + let output = "RNAcode lowest p-value: " ++ show lowestPvalue ++ "\nrnaCodeClassification: " ++ rnaCodeClassification return output else do - let output = (show lowestPvalue) ++ "\t" ++ rnaCodeClassification + let output = show lowestPvalue ++ "\t" ++ rnaCodeClassification return output - else do + else if (verbose == Loud) then do let output = "RNAcode lowest p-value: " ++ "-" ++ "\nrnaCodeClassification: " ++ "-" @@ -321,11 +321,11 @@ let output = "-\t" ++ "-" --let output = show (fromLeft inputRNACode) return output - else do + else if (verbose == Loud) then do let output = "RNAcode lowest p-value: " ++ "-" ++ "\nrnaCodeClassification: " ++ "-" return output else do let output = "-\t" ++ "-" - return output + return output diff --git a/src/Bio/SelectSequences.hs b/src/Bio/SelectSequences.hs deleted file mode 100644 --- a/src/Bio/SelectSequences.hs +++ /dev/null @@ -1,54 +0,0 @@ -{-# LANGUAGE RecordWildCards #-} -{-# LANGUAGE DeriveDataTypeable #-} - --- | Select Sequences --- Testcommand: SelectSequences -i /path/to/test.clustal -module Main where - -import System.Console.CmdArgs -import Bio.RNAlienLibrary -import Data.Either.Unwrap - -data Options = Options - { inputClustalPath :: String, - toogleExternalSelectSequences :: Bool, - seqenceNumber :: Int, - optimalIdentity :: Double, - maximalIdenity :: Double, - referenceSequence :: Bool, - distanceMatrixPath :: String - } deriving (Show,Data,Typeable) - -options :: Options -options = Options - { inputClustalPath = def &= name "c" &= help "Path to input clustal file", - toogleExternalSelectSequences = False &= name "e" &= help "Use only replacement of alignment characters and external 'selectSequence.pl'. Default: False", - seqenceNumber = (6 :: Int) &= name "n" &= help "Number of sequences in the output alignment. (Default: 6)", - optimalIdentity = (80 :: Double) &= name "i" &= help "Optimize for this percentage of mean pairwise identity (Default: 80)", - maximalIdenity = (95 :: Double) &= name "m" &= help "Sequences with a higher percentage of pairwise Identity will be removed. (Default: 95)", - referenceSequence = True &= name "x" &= help "The first sequence (=reference sequence) is always present in the output alignment per default. Default: True", - distanceMatrixPath = "" &= name "d" &= help "Path to distance matrix output file, only internal for interal sequence selection, e.g. /home/user/distmat (Default: )" - } &= summary "SelectSequences" &= help "Florian Eggenhofer 2016" &= verbosity - -main :: IO () -main = do - Options{..} <- cmdArgs options - let reformatedClustalPath = inputClustalPath ++ ".reformated" - if toogleExternalSelectSequences - then do - resultStatus <- preprocessClustalForRNAzExternal inputClustalPath reformatedClustalPath seqenceNumber (truncate optimalIdentity) (truncate maximalIdenity) referenceSequence - if (isRight resultStatus) - then do - let (idMatrix,resultAln) = fromRight resultStatus - putStr resultAln - if null distanceMatrixPath - then return () - else (writeFile distanceMatrixPath idMatrix) - else (print ("A problem occured selecting sequences: " ++ fromLeft resultStatus)) - else do - resultStatus <- preprocessClustalForRNAz inputClustalPath reformatedClustalPath seqenceNumber optimalIdentity maximalIdenity referenceSequence - if (isRight resultStatus) - then do - let (_,resultAln) = fromRight resultStatus - putStr resultAln - else (print ("A problem occured selecting sequences: " ++ fromLeft resultStatus)) diff --git a/src/Bio/cmsearchToBED.hs b/src/Bio/cmsearchToBED.hs --- a/src/Bio/cmsearchToBED.hs +++ b/src/Bio/cmsearchToBED.hs @@ -4,10 +4,10 @@ -- | Convert cmsearch output to Browser Extensible Data (BED) format -- Testcommand: cmsearchToBED -i /path/to/test.clustal module Main where -import Prelude -import System.Console.CmdArgs +import Prelude +import System.Console.CmdArgs import Bio.RNAlienLibrary -import Data.Either.Unwrap +import Data.Either.Unwrap import qualified Data.ByteString.Lazy.Char8 as L import qualified Data.Text as T import Data.List @@ -24,17 +24,17 @@ instance Show Bed where show (Bed _browserPostition _browserSettings _bedName _bedDescription _bedVisibility _bedItemRgb _bedEntries) = a ++ b ++ c ++ d ++ e ++ f ++ g - where a = "browser position " ++ (T.unpack _browserPostition) ++ "\n" - b = (T.unpack _browserSettings) ++ "\n" + where a = "browser position " ++ T.unpack _browserPostition ++ "\n" + b = T.unpack _browserSettings ++ "\n" c = "track name=\"" ++ T.unpack _bedName ++ "\" " d = "description=\"" ++ T.unpack _bedDescription ++ "\" " e = "visibility=" ++ show _bedVisibility ++ " " f = "itemRgb=\"" ++ itemRbg ++ "\"\n" itemRbg = if _bedItemRgb then "On" else "Off" g = concatMap show _bedEntries - - -data BedEntry = BedEntry + + +data BedEntry = BedEntry { chrom :: T.Text, chromStart :: Int, chromEnd :: Int, @@ -47,24 +47,24 @@ blockCount :: Maybe Int, blockSizes :: Maybe [Int], blockStarts :: Maybe [Int] - } deriving (Eq, Read) + } deriving (Eq, Read) instance Show BedEntry where show (BedEntry _chrom _chromStart _chromEnd _chromName _score _strand _thickStart _thickEnd _color _blockCount _blockSizes _blockStarts) = a ++ b ++ c ++ d ++ e ++ f ++ g ++ h ++ i ++ j ++ k ++ l - where a = T.unpack _chrom ++ "\t" - b = show _chromStart ++ "\t" + where a = T.unpack _chrom ++ "\t" + b = show _chromStart ++ "\t" c = show _chromEnd ++ "\t" d = maybe "" T.unpack _chromName ++ "\t" e = maybe "" show _score ++ "\t" - f = maybe "" (\s -> [s]) _strand ++ "\t" + f = maybe "" ((: [])) _strand ++ "\t" g = maybe "" show _thickStart ++ "\t" h = maybe "" show _thickEnd ++ "\t" i = maybe "" T.unpack _color ++ "\t" - j = maybe "" show _blockCount ++ "\t" - k = maybe "" (\sizes -> intercalate "," (map show sizes)) _blockSizes ++ "\t" - l = maybe "" (\starts -> intercalate "," (map show starts)) _blockStarts ++ "\n" - -data Options = Options + j = maybe "" show _blockCount ++ "\t" + k = maybe "" (intercalate "," . map show) _blockSizes ++ "\t" + l = maybe "" (intercalate "," . map show) _blockStarts ++ "\n" + +data Options = Options { cmsearchPath :: String, inputBrowserSettings :: String, inputBedVisibility :: Int, @@ -87,16 +87,16 @@ inputTrackColor = "255,0,0" &= name "c" &= help "RGB Color of the track. Default: 255,0,0", sortBed = True &= name "s" &= help "Sort entries of Bed file by start end end cooridinates. Default: True", withHeader = True &= name "w" &= help "Output contains bed header. Default: True" - } &= summary "cmsearchToBED - Converts cmsearch file hits to BED file entries" &= help "Florian Eggenhofer 2016" &= verbosity - + } &= summary "cmsearchToBED - Converts cmsearch file hits to BED file entries" &= help "Florian Eggenhofer 2016" &= verbosity + main :: IO () main = do Options{..} <- cmdArgs options parsedCmsearch <- readCMSearch cmsearchPath - if (isRight parsedCmsearch) + if isRight parsedCmsearch then do let outputBED = convertcmSearchToBED (fromRight parsedCmsearch) inputBrowserSettings inputTrackName inputTrackDescription inputTrackColor inputBedVisibility inputItemRgb sortBed - if (isRight outputBED) + if isRight outputBED then if withHeader then print (fromRight outputBED) @@ -104,7 +104,7 @@ let output = concatMap show (bedEntries (fromRight outputBED)) putStr output else putStr (fromLeft outputBED) - else (putStr ("A problem occured converting from cmsearch to BED format:\n " ++ show (fromLeft parsedCmsearch))) + else putStr ("A problem occured converting from cmsearch to BED format:\n " ++ show (fromLeft parsedCmsearch)) --convertcmSearchToBED :: CMsearch -> String -> String -> Either String String --convertcmSearchToBED inputcmsearch trackName trackColor @@ -125,13 +125,13 @@ bedEntries = map (cmsearchHitToBEDentry trackName trackColor) cmHits sortedBedEntries = if sortBed then sortBy orderBedEntry bedEntries else bedEntries currentBrowserPosition = T.unpack (chrom firstEntry) ++ ":" ++ show (chromStart firstEntry) ++ "-" ++ show (chromEnd firstEntry) - firstEntry = (head sortedBedEntries) + firstEntry = head sortedBedEntries bed = Bed (T.pack currentBrowserPosition) (T.pack inputBrowserSettings) (T.pack trackName) (T.pack trackDescription) inputBedVisibility inputItemRgb sortedBedEntries cmsearchHitToBEDentry :: String -> String -> CMsearchHit -> BedEntry cmsearchHitToBEDentry hitName hitColor cmHit = entry where entry = BedEntry chromosome entrystart entryend (Just (T.pack hitName)) entryscore entrystrand thickstart thickend entrycolor blocks blockSize blockStart - chromosome = T.pack (L.unpack (hitSequenceHeader cmHit)) + chromosome = T.pack (L.unpack (hitSequenceHeader cmHit)) --entryline = L.unpack (hitSequenceHeader cmHit) ++ "\t" ++ entryStart cmHit ++ "\t" ++ entryEnd cmHit++ "\t" ++ (hitName) ++ "\t" ++ "0" ++ "\t" ++ [(hitStrand cmHit)] ++ "\t" ++ show (hitStart cmHit) ++ "\t" ++ show (hitEnd cmHit) ++ "\t" ++ hitColor ++ "\n" entrystart = if hitStrand cmHit == '+' then hitStart cmHit else hitEnd cmHit entryend = if hitStrand cmHit == '+' then hitEnd cmHit else hitStart cmHit @@ -141,10 +141,10 @@ thickend = Just entryend entrycolor = Just (T.pack hitColor) blocks = Just (1 :: Int) - blockSize = Just [(entryend - entrystart)] - blockStart = Just [(0 :: Int)] - + blockSize = Just [entryend - entrystart] + blockStart = Just [0 :: Int] + --cmsearchHitToBEDentry :: String -> String -> CMsearchHit -> String --cmsearchHitToBEDentry hitName hitColor cmHit = entryline -- where entryline = L.unpack (hitSequenceHeader cmHit) ++ "\t" ++ entryStart cmHit ++ "\t" ++ entryEnd cmHit++ "\t" ++ (hitName) ++ "\t" ++ "0" ++ "\t" ++ [(hitStrand cmHit)] ++ "\t" ++ show (hitStart cmHit) ++ "\t" ++ show (hitEnd cmHit) ++ "\t" ++ hitColor ++ "\n" @@ -153,13 +153,13 @@ entryStart :: CMsearchHit -> String entryStart cmHit - | (hitStrand cmHit) == '+' = show (hitStart cmHit) + | hitStrand cmHit == '+' = show (hitStart cmHit) | otherwise = show (hitEnd cmHit) entryEnd :: CMsearchHit -> String entryEnd cmHit - | (hitStrand cmHit) == '+' = show (hitEnd cmHit) - | otherwise = show (hitStart cmHit) + | hitStrand cmHit == '+' = show (hitEnd cmHit) + | otherwise = show (hitStart cmHit) orderBedEntry :: BedEntry -> BedEntry -> Ordering orderBedEntry firstHit secondHit