RNAlien 1.0.0 → 1.1.0

raw patch · 7 files changed

+720/−380 lines, 7 filesdep +aesondep +http-conduitdep +hxtdep ~EntrezHTTPdep ~ViennaRNAParserPVP ok

version bump matches the API change (PVP)

Dependencies added: aeson, http-conduit, hxt, network, pureMD5, text, transformers

Dependency ranges changed: EntrezHTTP, ViennaRNAParser

API changes (from Hackage documentation)

Files

• RNAlien.cabal +15/−4
• src/Bio/InfernalParser.hs +340/−0
• src/Bio/RNAcentralHTTP.hs +111/−0
• src/Bio/RNAlien.hs +10/−7
• src/Bio/RNAlienData.hs +8/−9
• src/Bio/RNAlienLibrary.hs +182/−327
• src/Bio/RNAlienStatistics.hs +54/−33

@@ -2,7 +2,7 @@ -- see http://haskell.org/cabal/users-guide/  name:                RNAlien-version:             1.0.0+version:             1.1.0 synopsis:            Unsupervized construction of RNA family models description:         RNAlien is a tool for automatic construction of RNAfamily models from a single sequence.                      .@@ -26,11 +26,11 @@ 		     .                      * <http://infernal.janelia.org/ Infernal> 		     .-                     * <http://www.clustal.org/omega/#Download clustal-omega>-		     . 		     * <http://www.bioinf.uni-freiburg.de/Software/LocARNA/#download Locarna> 		     . 		     * <https://www.tbi.univie.ac.at/~wash/RNAz/ RNAz>+                     .+                     * <http://wash.github.io/rnacode/ RNAcode>                         		     . 		     * <http://www.tbi.univie.ac.at/RNA/index.html#download ViennaRNA package>                      .@@ -47,6 +47,15 @@ build-type:          Simple cabal-version:       >=1.8 +source-repository head+  type:     git+  location: https://github.com/eggzilla/RNAlien++source-repository this+  type:     git+  location: https://github.com/eggzilla/RNAlien/tree/1.1.0+  tag:      1.1.0+                      executable RNAlien   Hs-Source-Dirs:      ./src/Bio/   main-is:	       RNAlien.hs   @@ -68,8 +77,10 @@ Library   Hs-Source-Dirs:      ./src/   ghc-options:         -Wall -O2 -fno-warn-unused-do-bind-  build-depends:       base >=4.5 && <5, cmdargs, ViennaRNAParser>=1.2.6, process, directory, blastxml>=0.3.2, biofasta, parsec, random, BlastHTTP, biocore, bytestring, Taxonomy >= 1.0.1, either-unwrap, containers, ClustalParser>=1.1.0, EntrezHTTP>=1.0.0, split, vector, edit-distance, cassava, matrix, hierarchical-clustering, filepath, HTTP+  build-depends:       base >=4.5 && <5, cmdargs, ViennaRNAParser>=1.2.8, process, directory, blastxml>=0.3.2, biofasta, parsec, random, BlastHTTP, biocore, bytestring, Taxonomy >= 1.0.1, either-unwrap, containers, ClustalParser>=1.1.0, EntrezHTTP>=1.0.1, vector, edit-distance, cassava, matrix, hierarchical-clustering, filepath, HTTP, http-conduit, hxt, network, aeson, text, transformers >= 0.4 && <0.5, pureMD5   Exposed-Modules:     Bio.RNAlienData                        Bio.RNAlienLibrary+                       Bio.RNAcentralHTTP+                       Bio.InfernalParser  

@@ -0,0 +1,340 @@+-- | This module contains parsing functions for Infernal programs++module Bio.InfernalParser (+                           module Bio.RNAlienData,                          +                           readCMSearch,+                           readCMSearches,+                           parseCMSearch,+                           parseCMSearches,+                           parseCMstat,+                           readCMstat+                           )+where++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" ++-- | 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 +  parsedFile <- parseFromFile genParserCMSearch filePath+  CE.evaluate parsedFile ++-- | parse from input filePath                      +readCMSearches :: String -> IO (Either ParseError CMsearch)             +readCMSearches filePath = do +  parsedFile <- parseFromFile genParserCMSearches filePath+  CE.evaluate parsedFile ++genParserCMSearches :: GenParser Char st CMsearch+genParserCMSearches = do+  string "# cmsearch :: search CM(s) against a sequence database"+  newline+  string "# INFERNAL "+  many1 (noneOf "\n")+  newline       +  string "# Copyright (C) 201"+  many1 (noneOf "\n")+  newline       +  string "# Freely distributed under the GNU General Public License (GPLv3)."+  newline       +  string "# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -"+  newline+  string "# query CM file:"+  many1 space+  queryCMfile' <- many1 (noneOf "\n")+  newline+  string "# target sequence database:"+  many1 space      +  targetSequenceDatabase' <- many1 (noneOf "\n")+  newline+  optional (try (genParserCMsearchHeaderField "# CM configuration"))+  optional (try (genParserCMsearchHeaderField "# database size is set to"))+  optional (try (genParserCMsearchHeaderField "# truncated sequence detection"))+  string "# number of worker threads:"+  many1 space+  numberOfWorkerThreads' <- many1 (noneOf "\n")+  newline+  string "# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -"+  newline+  optional newline+  cmSearchesHits <- many1 (try genParserMultipleCMSearch)+  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       +  string "# Copyright (C) 201"+  many1 (noneOf "\n")+  newline       +  string "# Freely distributed under the GNU General Public License (GPLv3)."+  newline       +  string "# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -"+  newline+  string "# query CM file:"+  many1 space+  queryCMfile' <- many1 (noneOf "\n")+  newline+  string "# target sequence database:"+  many1 space      +  targetSequenceDatabase' <- many1 (noneOf "\n")+  newline+  optional (try (genParserCMsearchHeaderField "# CM configuration"))+  optional (try (genParserCMsearchHeaderField "# database size is set to"))+  optional (try (genParserCMsearchHeaderField "# truncated sequence detection"))+  string "# number of worker threads:"+  many1 space+  numberOfWorkerThreads' <- many1 (noneOf "\n")+  newline+  string "# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -"+  newline+  optional newline+  string "Query:"+  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 +  string "E-value"+  many1 space        +  string "score"+  many1 space +  string "bias"+  many1 space +  string "sequence"+  many1 space  +  string "start"+  many1 space +  string "end"+  many1 space +  string "mdl"+  many1 space +  string "trunc"+  many1 space +  string "gc"+  many1 space +  string "description"+  newline+  string " -"+  many1 (try (oneOf " -"))+  newline+  optional (try (string " ------ inclusion threshold ------"))+  many newline+  hitScores' <- many (try genParserCMsearchHit) --`endBy` (try (string "Hit alignments:"))+  optional (try genParserCMsearchEmptyHit)+  -- this is followed by hit alignments and internal cmsearch statistics which are not parsed+  many anyChar+  eof+  return $ CMsearch queryCMfile' targetSequenceDatabase' numberOfWorkerThreads' hitScores'++-- | Parsing function for CMSearches with multiple querymodels in one modelfile, e.g. clans+genParserMultipleCMSearch :: GenParser Char st [CMsearchHit]+genParserMultipleCMSearch = do+  --optional newline+  --optional string "//"+  string "Query:"+  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 +  string "E-value"+  many1 space        +  string "score"+  many1 space +  string "bias"+  many1 space +  string "sequence"+  many1 space  +  string "start"+  many1 space +  string "end"+  many1 space +  string "mdl"+  many1 space +  string "trunc"+  many1 space +  string "gc"+  many1 space +  string "description"+  newline+  string " -"+  many1 (try (oneOf " -"))+  newline+  optional (try (string " ------ inclusion threshold ------"))+  many newline+  hitScores' <- many (try genParserCMsearchHit) --`endBy` (try (string "Hit alignments:"))+  optional (try genParserCMsearchEmptyHit)+  -- this is followed by hit alignments and internal cmsearch statistics which are not parsed+  --many anyChar+  manyTill anyChar (try (string "//\n"))+  return hitScores'++genParserCMsearchHeaderField :: String -> GenParser Char st String+genParserCMsearchHeaderField fieldname = do+  string (fieldname ++ ":")+  many1 space+  many1 (noneOf "\n")+  newline+  return []++genParserCMsearchEmptyHit :: GenParser Char st [CMsearchHit]+genParserCMsearchEmptyHit = do+  string "   [No hits detected that satisfy reporting thresholds]"+  newline+  optional (try newline)+  return []++genParserCMsearchHit :: GenParser Char st CMsearchHit+genParserCMsearchHit = do+  many1 space+  string "("     +  hitRank' <- many1 digit+  string ")"+  many1 space+  hitSignificant' <- choice [char '!', char '?']+  many1 space                  +  hitEValue' <- many1 (oneOf "0123456789.e-")+  many1 space             +  hitScore'  <- many1 (oneOf "0123456789.e-")+  many1 space   +  hitBias' <- many1 (oneOf "0123456789.e-")+  many1 space+  hitSequenceHeader' <- many1 (noneOf " ")+  many1 space                +  hitStart' <- many1 digit+  many1 space+  hitEnd' <- many1 digit+  many1 space            +  hitStrand' <- choice [char '+', char '-', char '.']+  many1 space              +  hitModel' <- many1 letter+  many1 space          +  hitTruncation' <- many1 (choice [alphaNum, char '\''])+  many1 space                   +  hitGCcontent' <- many1 (oneOf "0123456789.e-")+  many1 space                +  hitDescription' <- many1 (noneOf "\n")     +  newline+  optional (try (string " ------ inclusion threshold ------"))+  optional (try newline)+  return $ CMsearchHit (readInt hitRank') hitSignificant' (readDouble hitEValue') (readDouble hitScore') (readDouble hitBias') (L.pack hitSequenceHeader') (readInt hitStart') (readInt hitEnd') hitStrand' (L.pack hitModel') (L.pack hitTruncation') (readDouble hitGCcontent') (L.pack hitDescription')++-- | parse from input filePath              +parseCMstat :: String -> Either ParseError CMstat+parseCMstat = parse genParserCMstat "parseCMstat"++-- | parse from input filePath                      +readCMstat :: String -> IO (Either ParseError CMstat)             +readCMstat filePath = do +  parsedFile <- parseFromFile genParserCMstat filePath+  CE.evaluate parsedFile +                      +genParserCMstat :: GenParser Char st CMstat+genParserCMstat = do+  string "# cmstat :: display summary statistics for CMs"+  newline+  string "# INFERNAL "+  many1 (noneOf "\n")+  newline       +  string "# Copyright (C) 201"+  many1 (noneOf "\n")+  newline       +  string "# Freely distributed under the GNU General Public License (GPLv3)."+  newline       +  string "# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -"+  newline+  char '#'+  many1 (char ' ')+  string "rel entropy"+  newline+  char '#'+  many1 (char ' ')+  many1 (char '-')+  newline+  char '#'+  many1 space +  string "idx"+  many1 space        +  string "name"+  many1 space +  string "accession"+  many1 space +  string "nseq"+  many1 space  +  string "eff_nseq"+  many1 space +  string "clen"+  many1 space +  string "W"+  many1 space +  string "bps"+  many1 space +  string "bifs"+  many1 space +  string "model"+  many1 space +  string "cm"+  many1 space+  string "hmm"+  newline+  string "#"+  many1 (try (oneOf " -"))+  newline+  many1 space     +  _statIndex <- many1 digit+  many1 space+  _statName <- many1 letter+  many1 space                  +  _statAccession <- many1 (noneOf " ")+  many1 space             +  _statSequenceNumber <- many1 digit+  many1 space   +  _statEffectiveSequences <- many1 (oneOf "0123456789.e-")+  many1 space+  _statConsensusLength <- many digit+  many1 space                +  _statW <- many1 digit+  many1 space+  _statBasepaires <- many1 digit+  many1 space            +  _statBifurcations <- many1 digit+  many1 space              +  _statModel <- many1 letter+  many1 space          +  _relativeEntropyCM <- many1 (oneOf "0123456789.e-")+  many1 space                   +  _relativeEntropyHMM <- many1 (oneOf "0123456789.e-")+  newline+  char '#'+  newline+  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+readInt = read++readDouble :: String -> Double+readDouble = read

@@ -0,0 +1,111 @@+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE Arrows #-}+{-# LANGUAGE DeriveGeneric #-}++-- | Interface for the RNAcentral REST webservice.+--   +module Bio.RNAcentralHTTP (rnaCentralHTTP,+                      buildSequenceViaMD5Query,+                      getRNACentralEntries,+                      showRNAcentralAlienEvaluation,+                      RNAcentralEntryResponse,+                      RNAcentralEntry+                      ) where++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.Sequence.Fasta+import Data.Either+import Data.Aeson.Types++--Datatypes+-- | Data structure for RNAcentral entry response+data RNAcentralEntryResponse = RNAcentralEntryResponse+  {+    _count :: Int,+    _next :: Maybe Text,+    _previous :: Maybe Text,+    results :: [RNAcentralEntry]+  }+  deriving (Show, Eq, Generic)++instance ToJSON RNAcentralEntryResponse where+  toJSON = genericToJSON defaultOptions+  --toEncoding = genericToEncoding defaultOptions++instance FromJSON RNAcentralEntryResponse ++data RNAcentralEntry = RNAcentralEntry+  {+    _url :: Text,+    rnacentral_id :: Text,+    md5 :: Text,+    _sequence :: Text,+    length :: Int,+    _xrefs :: Text,+    _publications :: Text+  }+  deriving (Show, Eq, Generic)++instance ToJSON RNAcentralEntry where+  toJSON = genericToJSON defaultOptions+  --toEncoding = genericToEncoding defaultOptions++instance FromJSON RNAcentralEntry ++-- | Send query and parse return XML +startSession :: String -> IO (Either String RNAcentralEntryResponse)+startSession query' = do+  requestXml <- withSocketsDo+      $ sendQuery query'+  let eitherErrorResponse = eitherDecode requestXml :: Either String RNAcentralEntryResponse+  return eitherErrorResponse+  +-- | Send query and return response XML+sendQuery :: String -> IO L8.ByteString+sendQuery query' = simpleHttp ("http://rnacentral.org/api/v1/rna/" ++ query')++-- | Function for querying the RNAcentral REST interface.+rnaCentralHTTP :: String -> IO (Either String RNAcentralEntryResponse)+rnaCentralHTTP query' = do+  startSession query'++-- | Function for delayed queries to the RNAcentral REST interface. Enforces the maximum 20 requests per second policy.+delayedRNACentralHTTP :: String -> IO (Either String RNAcentralEntryResponse)+delayedRNACentralHTTP query' = do+  threadDelay 55000+  startSession query'++getRNACentralEntries :: [String] -> IO [(Either String RNAcentralEntryResponse)]+getRNACentralEntries queries = do+  responses <- mapM delayedRNACentralHTTP queries+  return responses++buildSequenceViaMD5Query :: Sequence -> String+buildSequenceViaMD5Query s = qString+  where querySequence = unSD (seqdata s)+        querySequenceUreplacedwithT = L8.map bsreplaceUT querySequence+        md5Sequence = M.md5 querySequenceUreplacedwithT+        qString = "?md5=" ++ (show md5Sequence)++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 resulthead ++ "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"++bsreplaceUT :: Char -> Char+bsreplaceUT a+  | a == 'U' = 'T'+  | otherwise = a

@@ -26,6 +26,8 @@     lengthFilter :: Bool,     coverageFilter :: Bool,     singleHitperTax :: Bool,+    inputQuerySelectionMethod :: String,+    inputQueryNumber :: Int,     threads :: Int,     taxonomyRestriction :: Maybe String,     sessionIdentificator :: Maybe String@@ -41,11 +43,13 @@     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 = True &= name "s" &= help "Only the best blast hit per taxonomic entry is considered. Default: True",+    singleHitperTax = False &= name "s" &= help "Only the best blast hit per taxonomic entry is considered. Default: False",+    inputQuerySelectionMethod = "filtering" &= name "m" &= help "Method for selection of queries (filtering,clustering). Default: filtering",+    inputQueryNumber = (5 :: Int) &= name "n" &= help "Number of queries used for candidate search. Default: 5",     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."-  } &= summary "RNAlien version 1.0.0" &= help "Florian Eggenhofer, Ivo L. Hofacker, Christian Höner zu Siederdissen - 2013 - 2015" &= verbosity       +  } &= summary "RNAlien version 1.0.0" &= help "Florian Eggenhofer, Ivo L. Hofacker, Christian Höner zu Siederdissen - 2013 - 2016" &= verbosity                         main :: IO () main = do@@ -75,8 +79,8 @@           logMessage "Error: Input fasta file is empty.\n" temporaryDirectoryPath         else do           let iterationNumber = 0-          let tools = ["clustalo","mlocarna","RNAfold","RNAalifold","cmcalibrate","cmstat","cmbuild","RNAz"]-          toolsCheck <- checkTools tools temporaryDirectoryPath+          let tools = ["mlocarna","RNAfold","RNAalifold","cmcalibrate","cmstat","cmbuild","RNAz","RNAcode"]+          toolsCheck <- checkTools tools inputQuerySelectionMethod temporaryDirectoryPath           -- Check required commandline tools           if isLeft toolsCheck             then do @@ -87,13 +91,12 @@               let inputSequence = reformatFasta (head inputFasta)               initialTaxId <- setInitialTaxId inputBlastDatabase temporaryDirectoryPath inputTaxId inputSequence               let checkedTaxonomyRestriction = checkTaxonomyRestriction taxonomyRestriction-              let staticOptions = StaticOptions temporaryDirectoryPath sessionId (fromJust inputnSCICutoff) inputTaxId singleHitperTax lengthFilter coverageFilter threads inputBlastDatabase checkedTaxonomyRestriction (setVerbose verboseLevel)+              let staticOptions = StaticOptions temporaryDirectoryPath sessionId (fromJust inputnSCICutoff) inputTaxId singleHitperTax inputQuerySelectionMethod inputQueryNumber lengthFilter coverageFilter 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-              let resultNumber = (length (concatMap sequenceRecords (taxRecords modelConstructionResults))) + 1-              resultEvaluation <- evaluateConstructionResult staticOptions resultNumber+              resultEvaluation <- evaluateConstructionResult staticOptions modelConstructionResults               appendFile (temporaryDirectoryPath ++ "Log") resultEvaluation               writeFile (temporaryDirectoryPath ++ "result.csv") resultTaxonomyRecordsCSVTable               resultSummary modelConstructionResults staticOptions

@@ -14,6 +14,8 @@     nSCICutoff :: Double,     userTaxId :: Maybe Int,     singleHitperTaxToggle :: Bool,+    querySelectionMethod :: String,+    queryNumber :: Int,     lengthFilterToggle :: Bool,     coverageFilterToggle :: Bool,     cpuThreads :: Int,@@ -62,17 +64,14 @@     nucleotideSequence :: Sequence,     -- 0 is unaligned, number is the iteration the sequence has been included into the alignment     aligned  :: Int,-    recordDescription :: String,-    -- Is the sequence derived from the blast hit coordinates (B) or from a corresponding genbank feature (G)-    sequenceOrigin :: Char    +    recordDescription :: L.ByteString      }   instance Show SequenceRecord where-  show (SequenceRecord _nucleotideSequence _aligned _recordDescription _sequenceOrigin) = a ++ b ++ c ++ d -    where a = "SequenceRecord TaxonomyId: " ++ show _recordDescription ++ "\n" -          b = "Sequence Origin: " ++ _recordDescription ++ "\n" -          c = "Aligned in iteration: " ++ show _aligned ++ "\n" -          d = "Sequence Origin: " ++ 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,@@ -101,7 +100,7 @@   } deriving (Show, Eq, Read)   data SearchResult = SearchResult-  { candidates :: [(Sequence,Int,String,Char)],+  { candidates :: [(Sequence,Int,L.ByteString)],     blastDatabaseSize :: Maybe Double   } 

@@ -13,6 +13,7 @@                            checkTools,                            systemCMsearch,                            readCMSearch,+                           readCMSearches,                            compareCM,                            parseCMSearch,                            cmSearchsubString,@@ -46,9 +47,8 @@ import Data.Either.Unwrap import Data.Maybe import Bio.EntrezHTTP -import qualified Data.List.Split as DS import System.Exit-import Data.Either (lefts,rights)+import Data.Either (lefts,rights,Either) import qualified Text.EditDistance as ED    import qualified Data.Vector as V import Control.Concurrent @@ -64,6 +64,12 @@ import Bio.RNAalifoldParser import Bio.RNAzParser import Network.HTTP+import qualified Bio.RNAcodeParser as RC+import Bio.RNAcentralHTTP+import Bio.InfernalParser+import qualified Data.Text as T+import qualified Data.Text.IO as TI+import qualified Data.Text.Encoding as DTE  -- | Initial RNA family model construction - generates iteration number, seed alignment and model modelConstructer :: StaticOptions -> ModelConstruction -> IO ModelConstruction@@ -192,13 +198,15 @@       let alignmentSequences = map snd (V.toList (V.concat [alignedSequences]))       writeFasta preliminaryFastaPath alignmentSequences       let cmBuildFilepath = iterationDirectory ++ "model" ++ ".cmbuild"+      let refinedAlignmentFilepath = iterationDirectory ++ "modelrefined" ++ ".stockholm"+      let cmBuildOptions ="--refine " ++ refinedAlignmentFilepath       let foldFilepath = iterationDirectory ++ "model" ++ ".fold"       _ <- systemRNAfold preliminaryFastaPath foldFilepath       foldoutput <- readRNAfold foldFilepath       let seqStructure = foldSecondaryStructure (fromRight foldoutput)       let stockholAlignment = convertFastaFoldStockholm (head alignmentSequences) seqStructure       writeFile preliminaryAlignmentPath stockholAlignment-      _ <- systemCMbuild preliminaryAlignmentPath preliminaryCMPath cmBuildFilepath+      _ <- systemCMbuild cmBuildOptions preliminaryAlignmentPath preliminaryCMPath cmBuildFilepath       _ <- systemCMcalibrate "fast" (cpuThreads staticOptions) preliminaryCMPath preliminaryCMLogPath       resultModelConstruction <- reevaluatePotentialMembers staticOptions nextModelConstructionInput       return resultModelConstruction@@ -247,6 +255,7 @@       let lastIterationAlignmentPath = outputDirectory ++ show (currentIterationNumber - 1)  ++ "/model.stockholm"       let lastIterationCMPath = outputDirectory ++ show (currentIterationNumber - 1)++ "/model.cm"       copyFile lastIterationCMPath resultCMPath+      --copyFile lastIterationCMPath (resultCMPath ++ ".bak1")       copyFile lastIterationFastaPath resultFastaPath       copyFile lastIterationAlignmentPath resultAlignmentPath       _ <- systemCMcalibrate "standard" (cpuThreads staticOptions) resultCMPath resultCMLogPath@@ -260,8 +269,10 @@       let nextModelConstructionInput = constructNext currentIterationNumber modelConstruction alignmentResults Nothing Nothing [] [] (alignmentModeInfernal modelConstruction)       constructModel nextModelConstructionInput staticOptions       copyFile lastIterationCMPath resultCMPath+      --debug+      --copyFile lastIterationCMPath (resultCMPath ++ ".bak2")       copyFile lastIterationFastaPath resultFastaPath-      copyFile lastIterationAlignmentPath resultAlignmentPath +      copyFile lastIterationAlignmentPath resultAlignmentPath       logMessage (iterationSummaryLog nextModelConstructionInput) outputDirectory       logVerboseMessage (verbositySwitch staticOptions) (show nextModelConstructionInput) outputDirectory       _ <- systemCMcalibrate "standard" (cpuThreads staticOptions) resultCMPath resultCMLogPath@@ -393,7 +404,9 @@   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 blastQuery = BlastHTTPQuery (Just "ncbi") (Just "blastn") (blastDatabase staticOptions) querySequences'  (Just (hitNumberQuery ++ entrezTaxFilter ++ registrationInfo)) (Just (5400000000 :: Int))+  let softmaskFilter = "&FILTER=True&FILTER=m"+  let blastQuery = BlastHTTPQuery (Just "ncbi") (Just "blastn") (blastDatabase staticOptions) querySequences'  (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)   blastOutput <- CE.catch (blastHTTP blastQuery)                        (\e -> do let err = show (e :: CE.IOException)@@ -422,7 +435,7 @@        let blastHitsFilteredByCoverage = filterByCoverage blastHitsFilteredByLength queryLength (coverageFilterToggle staticOptions)        writeFile (logFileDirectoryPath ++ "/" ++ queryIndexString  ++ "_3ablastHitsFilteredByLength") (showlines blastHitsFilteredByCoverage)        --tag BlastHits with TaxId-       blastHitsWithTaxIdOutput <- retrieveBlastHitsTaxIdEntrez blastHitsFilteredByLength+       blastHitsWithTaxIdOutput <- retrieveBlastHitsTaxIdEntrez blastHitsFilteredByCoverage        let uncheckedBlastHitsWithTaxIdList = map (\(blasthits,taxIdout) -> (blasthits,extractTaxIdFromEntrySummaries taxIdout)) blastHitsWithTaxIdOutput        let checkedBlastHitsWithTaxId = filter (\(_,taxids) -> not (null taxids)) uncheckedBlastHitsWithTaxIdList        --todo checked blasthittaxidswithblasthits need to be merged as taxid blasthit pairs@@ -451,7 +464,7 @@          else do            writeFile (logFileDirectoryPath ++ "/" ++ queryIndexString ++ "_10afullSequencesWithSimilars") (showlines fullSequencesWithSimilars)            let fullSequences = filterIdenticalSequences fullSequencesWithSimilars 100-           let fullSequencesWithOrigin = map (\(parsedFasta,taxid,seqSubject) -> (parsedFasta,taxid,seqSubject,'B')) fullSequences+           --let fullSequencesWithOrigin = map (\(parsedFasta,taxid,seqSubject) -> (parsedFasta,taxid,seqSubject,'B')) fullSequences            writeFile (logFileDirectoryPath ++ "/" ++ queryIndexString ++ "_10fullSequences") (showlines fullSequences)            let maybeFractionEvalueMatch = getHitWithFractionEvalue rightBlast            if isNothing maybeFractionEvalueMatch@@ -459,14 +472,14 @@              else do                let fractionEvalueMatch = fromJust maybeFractionEvalueMatch                let dbSize = computeDataBaseSize (e_val fractionEvalueMatch) (bits fractionEvalueMatch) (fromIntegral queryLength ::Double)-               CE.evaluate (SearchResult fullSequencesWithOrigin (Just dbSize))+               CE.evaluate (SearchResult fullSequences (Just dbSize))      else CE.evaluate (SearchResult [] Nothing)    -- |Computes size of blast db in Mb  computeDataBaseSize :: Double -> Double -> Double -> Double  computeDataBaseSize evalue bitscore querylength = ((evalue * 2 ** bitscore) / querylength)/10^(6 :: Integer) -alignCandidates :: StaticOptions -> ModelConstruction -> String -> SearchResult -> IO ([(Sequence,Int,String,Char)],[(Sequence,Int,String,Char)])+alignCandidates :: StaticOptions -> ModelConstruction -> String -> SearchResult -> IO ([(Sequence,Int,L.ByteString)],[(Sequence,Int,L.ByteString)]) alignCandidates staticOptions modelConstruction multipleSearchResultPrefix searchResults = do   let iterationDirectory = tempDirPath staticOptions ++ show (iterationNumber modelConstruction) ++ "/" ++ multipleSearchResultPrefix   createDirectoryIfMissing False (iterationDirectory ++ "log")@@ -484,7 +497,7 @@         then alignCandidatesInfernalMode staticOptions modelConstruction multipleSearchResultPrefix (blastDatabaseSize searchResults) filteredCandidates         else alignCandidatesInitialMode staticOptions modelConstruction multipleSearchResultPrefix filteredCandidates -alignCandidatesInfernalMode :: StaticOptions -> ModelConstruction -> String -> Maybe Double -> [(Sequence,Int,String,Char)] -> IO ([(Sequence,Int,String,Char)],[(Sequence,Int,String,Char)])+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 @@ -508,7 +521,7 @@   writeFile (iterationDirectory ++ "log" ++ "/15potentialCandidates'") (showlines potentialCandidates)                                            CE.evaluate (map snd trimmedSelectedCandidates,map snd potentialCandidates) -alignCandidatesInitialMode :: StaticOptions -> ModelConstruction -> String -> [(Sequence,Int,String,Char)] -> IO ([(Sequence,Int,String,Char)],[(Sequence,Int,String,Char)])+alignCandidatesInitialMode :: StaticOptions -> ModelConstruction -> String -> [(Sequence,Int,L.ByteString)] -> IO ([(Sequence,Int,L.ByteString)],[(Sequence,Int,L.ByteString)]) alignCandidatesInitialMode staticOptions modelConstruction multipleSearchResultPrefix filteredCandidates = do   let iterationDirectory = tempDirPath staticOptions ++ show (iterationNumber modelConstruction) ++ "/" ++ multipleSearchResultPrefix   --writeFile (iterationDirectory ++ "log" ++ "/11bcandidates") (showlines filteredCandidates)@@ -562,8 +575,9 @@   | currentClusterNumber >= numberOfClusters = currentCutoff   | otherwise = findCutoffforClusterNumber clustaloDendrogram numberOfClusters (currentCutoff-0.01)     where currentClusterNumber = length (cutAt clustaloDendrogram currentCutoff)-                -selectQueries :: StaticOptions -> ModelConstruction -> [(Sequence,Int,String,Char)] -> IO [String]++-- Selects Query sequence ids from all collected seqeuences. Queries are then fetched by extractQueries function.+selectQueries :: StaticOptions -> ModelConstruction -> [(Sequence,Int,L.ByteString)] -> IO [String] selectQueries staticOptions modelConstruction selectedCandidates = do   logVerboseMessage (verbositySwitch staticOptions) "SelectQueries\n" (tempDirPath staticOptions)   --Extract sequences from modelconstruction@@ -573,32 +587,39 @@   let alignmentSequences = map snd (V.toList (V.concat [candidateSequences,alignedSequences]))   if length alignmentSequences > 3     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" -      alignSequences "clustalo" ("--full --distmat-out=" ++ clustaloDistMatrixPath ++ " ") [fastaFilepath] [] [clustaloFilepath] []-      idsDistancematrix <- readClustaloDistMatrix clustaloDistMatrixPath-      logEither idsDistancematrix (tempDirPath staticOptions)-      let (clustaloIds,clustaloDistMatrix) = fromRight idsDistancematrix-      logVerboseMessage (verbositySwitch staticOptions) ("Clustalid: " ++ intercalate "," clustaloIds ++ "\n") (tempDirPath staticOptions)-      logVerboseMessage (verbositySwitch staticOptions) ("Distmatrix: " ++ show clustaloDistMatrix ++ "\n") (tempDirPath staticOptions)-      let clustaloDendrogram = dendrogram UPGMA clustaloIds (getDistanceMatrixElements clustaloIds clustaloDistMatrix)-      logVerboseMessage (verbositySwitch staticOptions) ("ClustaloDendrogram: " ++ show  clustaloDendrogram ++ "\n") (tempDirPath staticOptions)-      logVerboseMessage (verbositySwitch staticOptions) ("ClustaloDendrogram: " ++ show clustaloDistMatrix ++ "\n") (tempDirPath staticOptions)-      let numberOfClusters = setClusterNumber (length alignmentSequences)-      logVerboseMessage (verbositySwitch staticOptions) ("numberOfClusters: " ++ show numberOfClusters ++ "\n") (tempDirPath staticOptions)-      let dendrogramStartCutDistance = 1 :: Double-      let dendrogramCutDistance' = findCutoffforClusterNumber clustaloDendrogram numberOfClusters dendrogramStartCutDistance-      logVerboseMessage (verbositySwitch staticOptions) ("dendrogramCutDistance': " ++ show dendrogramCutDistance' ++ "\n") (tempDirPath staticOptions)-      let cutDendrogram = cutAt clustaloDendrogram dendrogramCutDistance'-      --putStrLn "cutDendrogram: "-      --print cutDendrogram-      let currentSelectedQueries = take 5 (concatMap (take 1 . elements) cutDendrogram)-      logVerboseMessage (verbositySwitch staticOptions) ("SelectedQueries: " ++ show currentSelectedQueries ++ "\n") (tempDirPath staticOptions)                       -      writeFile (tempDirPath staticOptions ++ show (iterationNumber modelConstruction) ++ "/log" ++ "/13selectedQueries") (showlines currentSelectedQueries)-      CE.evaluate currentSelectedQueries+      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" +          alignSequences "clustalo" ("--full --distmat-out=" ++ clustaloDistMatrixPath ++ " ") [fastaFilepath] [] [clustaloFilepath] []+          idsDistancematrix <- readClustaloDistMatrix clustaloDistMatrixPath+          logEither idsDistancematrix (tempDirPath staticOptions)+          let (clustaloIds,clustaloDistMatrix) = fromRight idsDistancematrix+          logVerboseMessage (verbositySwitch staticOptions) ("Clustalid: " ++ intercalate "," clustaloIds ++ "\n") (tempDirPath staticOptions)+          logVerboseMessage (verbositySwitch staticOptions) ("Distmatrix: " ++ show clustaloDistMatrix ++ "\n") (tempDirPath staticOptions)+          let clustaloDendrogram = dendrogram UPGMA clustaloIds (getDistanceMatrixElements clustaloIds clustaloDistMatrix)+          logVerboseMessage (verbositySwitch staticOptions) ("ClustaloDendrogram: " ++ show  clustaloDendrogram ++ "\n") (tempDirPath staticOptions)+          logVerboseMessage (verbositySwitch staticOptions) ("ClustaloDendrogram: " ++ show clustaloDistMatrix ++ "\n") (tempDirPath staticOptions)+          let numberOfClusters = setClusterNumber (length alignmentSequences)+          logVerboseMessage (verbositySwitch staticOptions) ("numberOfClusters: " ++ show numberOfClusters ++ "\n") (tempDirPath staticOptions)+          let dendrogramStartCutDistance = 1 :: Double+          let dendrogramCutDistance' = findCutoffforClusterNumber clustaloDendrogram numberOfClusters dendrogramStartCutDistance+          logVerboseMessage (verbositySwitch staticOptions) ("dendrogramCutDistance': " ++ show dendrogramCutDistance' ++ "\n") (tempDirPath staticOptions)+          let cutDendrogram = cutAt clustaloDendrogram dendrogramCutDistance'+          --putStrLn "cutDendrogram: "+          --print cutDendrogram+          let currentSelectedQueries = take (queryNumber staticOptions) (concatMap (take 1 . elements) cutDendrogram)+          logVerboseMessage (verbositySwitch staticOptions) ("SelectedQueries: " ++ show currentSelectedQueries ++ "\n") (tempDirPath staticOptions)                       +          writeFile (tempDirPath staticOptions ++ show (iterationNumber modelConstruction) ++ "/log" ++ "/13selectedQueries") (showlines currentSelectedQueries)+          CE.evaluate currentSelectedQueries+        else do+          let currentSelectedSequences = filterIdenticalSequences' alignmentSequences (95 :: Double)+          let currentSelectedQueries = map (L.unpack . unSL . seqid) (take (queryNumber staticOptions) currentSelectedSequences)+          CE.evaluate currentSelectedQueries+               else return []  constructModel :: ModelConstruction -> StaticOptions -> IO String@@ -620,6 +641,8 @@   let cmCalibrateFilepath = outputDirectory ++ "model" ++ ".cmcalibrate"   let cmBuildFilepath = outputDirectory ++ "model" ++ ".cmbuild"   let alifoldFilepath = outputDirectory ++ "model" ++ ".alifold"+  let refinedAlignmentFilepath = outputDirectory ++ "modelrefined" ++ ".stockholm"+  let cmBuildOptions ="--refine " ++ refinedAlignmentFilepath   if alignmentModeInfernal modelConstruction      then do        logVerboseMessage (verbositySwitch staticOptions) "Construct Model - infernal mode\n" (tempDirPath staticOptions)@@ -628,12 +651,12 @@        replaceStatus <- replaceStockholmStructure stockholmFilepath alifoldFilepath updatedStructureStockholmFilepath        if null replaceStatus          then do-           systemCMbuild updatedStructureStockholmFilepath cmFilepath cmBuildFilepath+           systemCMbuild cmBuildOptions updatedStructureStockholmFilepath cmFilepath cmBuildFilepath            systemCMcalibrate "fast" (cpuThreads staticOptions) cmFilepath cmCalibrateFilepath            return cmFilepath          else do            logWarning ("Warning: A problem occured updating the secondary structure of iteration " ++ show (iterationNumber modelConstruction)  ++ " stockholm alignment: " ++ replaceStatus) (tempDirPath staticOptions)-           systemCMbuild stockholmFilepath cmFilepath cmBuildFilepath+           systemCMbuild cmBuildOptions stockholmFilepath cmFilepath cmBuildFilepath            systemCMcalibrate "fast" (cpuThreads staticOptions) cmFilepath cmCalibrateFilepath            return cmFilepath      else do@@ -643,7 +666,7 @@        logEither mlocarnaAlignment (tempDirPath staticOptions)        let stockholAlignment = convertClustaltoStockholm (fromRight mlocarnaAlignment)        writeFile stockholmFilepath stockholAlignment-       _ <- systemCMbuild stockholmFilepath cmFilepath cmBuildFilepath+       _ <- systemCMbuild cmBuildOptions stockholmFilepath cmFilepath cmBuildFilepath        _ <- systemCMcalibrate "fast" (cpuThreads staticOptions) cmFilepath cmCalibrateFilepath        return cmFilepath @@ -740,7 +763,7 @@ filterDuplicates modelConstruction inputSearchResult = uniqueSearchResult   where alignedSequences = map snd (V.toList (extractAlignedSequences (iterationNumber modelConstruction) modelConstruction))         collectedIdentifiers = map seqid alignedSequences-        uniques = filter (\(s,_,_,_) -> notElem (seqid s) collectedIdentifiers) (candidates inputSearchResult)+        uniques = filter (\(s,_,_) -> notElem (seqid s) collectedIdentifiers) (candidates inputSearchResult)         uniqueSearchResult = SearchResult uniques (blastDatabaseSize inputSearchResult)  -- | Filter a list of similar extended blast hits   @@ -751,35 +774,37 @@ --filterIdenticalSequencesWithOrigin [] _ = []  -- | Filter a list of similar extended blast hits   -filterIdenticalSequences :: [(Sequence,Int,String)] -> Double -> [(Sequence,Int,String)]                            +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          result = headSequence:filterIdenticalSequences filteredSequences identitycutoff filterIdenticalSequences [] _ = []  -- | Filter sequences too similar to already aligned sequences-filterWithCollectedSequences :: [(Sequence,Int,String,Char)] -> [Sequence] -> Double -> [(Sequence,Int,String,Char)]                            -filterWithCollectedSequences inputCandidates collectedSequences identitycutoff = filter (isUnSimilarSequence collectedSequences identitycutoff . firstOfQuadruple) 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+        result = headEntry:filterIdenticalSequences' filteredEntries identitycutoff+filterIdenticalSequences' [] _ = []++-- | Filter alignment entries by similiarity   filterIdenticalAlignmentEntry :: [ClustalAlignmentEntry] -> Double -> [ClustalAlignmentEntry] filterIdenticalAlignmentEntry (headEntry:rest) identitycutoff = result   where filteredEntries = filter (\x -> (stringIdentity (entryAlignedSequence headEntry) (entryAlignedSequence x)) < identitycutoff) rest         result = headEntry:filterIdenticalAlignmentEntry filteredEntries identitycutoff filterIdenticalAlignmentEntry [] _ = [] - isUnSimilarSequence :: [Sequence] -> Double -> Sequence -> Bool isUnSimilarSequence collectedSequences identitycutoff checkSequence = any (\ x -> (sequenceIdentity checkSequence x) < identitycutoff) collectedSequences--                 +                 firstOfTriple :: (t, t1, t2) -> t firstOfTriple (a,_,_) = a  -firstOfQuadruple :: (t, t1, t2, t3) -> t-firstOfQuadruple (a,_,_,_) = a - -- | Check if the result field of BlastResult is filled and if hits are present blastMatchesPresent :: BlastResult -> Bool blastMatchesPresent blastResult @@ -831,7 +856,7 @@         --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, String, Char)] -> Maybe Int -> Maybe Taxon  -> [String] -> [SearchResult] -> Bool -> ModelConstruction+constructNext :: Int -> ModelConstruction -> [(Sequence,Int,L.ByteString)] -> Maybe Int -> Maybe Taxon  -> [String] -> [SearchResult] -> Bool -> ModelConstruction constructNext currentIterationNumber modelconstruction alignmentResults upperTaxLimit inputTaxonomicContext inputSelectedQueries inputPotentialMembers toggleInfernalAlignmentModeTrue = nextModelConstruction   where newIterationNumber = currentIterationNumber + 1         taxEntries = taxRecords modelconstruction ++ buildTaxRecords alignmentResults currentIterationNumber@@ -839,25 +864,25 @@         currentAlignmentMode = toggleInfernalAlignmentModeTrue || alignmentModeInfernal modelconstruction         nextModelConstruction = ModelConstruction newIterationNumber (inputFasta modelconstruction) taxEntries upperTaxLimit inputTaxonomicContext (evalueThreshold modelconstruction) currentAlignmentMode inputSelectedQueries potMembers          -buildTaxRecords :: [(Sequence,Int,String,Char)] -> Int -> [TaxonomyRecord]+buildTaxRecords :: [(Sequence,Int,L.ByteString)] -> Int -> [TaxonomyRecord] buildTaxRecords alignmentResults currentIterationNumber = taxonomyRecords   where taxIdGroups = groupBy sameTaxIdAlignmentResult alignmentResults         taxonomyRecords = map (buildTaxRecord currentIterationNumber) taxIdGroups     -sameTaxIdAlignmentResult :: (Sequence,Int,String,Char) -> (Sequence,Int,String,Char) -> Bool-sameTaxIdAlignmentResult (_,taxId1,_,_) (_,taxId2,_,_) = taxId1 == taxId2+sameTaxIdAlignmentResult :: (Sequence,Int,L.ByteString) -> (Sequence,Int,L.ByteString) -> Bool+sameTaxIdAlignmentResult (_,taxId1,_) (_,taxId2,_) = taxId1 == taxId2 -buildTaxRecord :: Int -> [(Sequence,Int,String,Char)] -> TaxonomyRecord+buildTaxRecord :: Int -> [(Sequence,Int,L.ByteString)] -> TaxonomyRecord buildTaxRecord currentIterationNumber entries = taxRecord-  where recordTaxId = (\(_,currentTaxonomyId,_,_) -> currentTaxonomyId) (head entries)+  where recordTaxId = (\(_,currentTaxonomyId,_) -> currentTaxonomyId) (head entries)         seqRecords = map (buildSeqRecord currentIterationNumber)  entries         taxRecord = TaxonomyRecord recordTaxId seqRecords -buildSeqRecord :: Int -> (Sequence,Int,String,Char) -> SequenceRecord -buildSeqRecord currentIterationNumber (parsedFasta,_,seqSubject,seqOrigin) = SequenceRecord parsedFasta currentIterationNumber seqSubject seqOrigin   +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, String, Char))] -> ([(CMsearch,(Sequence, Int, String, Char))],[(CMsearch,(Sequence, Int, String, Char))],[(CMsearch,(Sequence, Int, String, Char))])+evaluePartitionTrimCMsearchHits :: Double -> [(CMsearch,(Sequence, Int, L.ByteString))] -> ([(CMsearch,(Sequence, Int, L.ByteString))],[(CMsearch,(Sequence, Int, L.ByteString))],[(CMsearch,(Sequence, Int, L.ByteString))]) evaluePartitionTrimCMsearchHits eValueThreshold cmSearchCandidatesWithSequences = (trimmedSelectedCandidates,potentialCandidates,rejectedCandidates)   where potentialMemberseValueThreshold = eValueThreshold * 1000         (prefilteredCandidates,rejectedCandidates) = partition (\(cmSearchResult,_) -> any (\hitScore' -> potentialMemberseValueThreshold >= hitEvalue hitScore') (cmsearchHits cmSearchResult)) cmSearchCandidatesWithSequences@@ -865,8 +890,8 @@         trimmedSelectedCandidates = map (\(cmSearchResult,inputSequence) -> (cmSearchResult,trimCMsearchHit cmSearchResult inputSequence)) selectedCandidates                  -trimCMsearchHit :: CMsearch -> (Sequence, Int, String, Char) -> (Sequence, Int, String, Char)-trimCMsearchHit cmSearchResult (inputSequence,b,c,d) = (subSequence,b,c,d)+trimCMsearchHit :: CMsearch -> (Sequence, Int, L.ByteString) -> (Sequence, Int, L.ByteString)+trimCMsearchHit cmSearchResult (inputSequence,b,c) = (subSequence,b,c)   where hitScoreEntry = head (cmsearchHits cmSearchResult)         sequenceString = L.unpack (unSD (seqdata inputSequence))         sequenceSubstring = cmSearchsubString (hitStart hitScoreEntry) (hitEnd hitScoreEntry) sequenceString@@ -893,9 +918,9 @@         alignedSequences = fastaSeqData:map nucleotideSequence (concatMap sequenceRecords (taxRecords modelconstruction))          querySequences' = concatMap (\querySeqId -> filter (\alignedSeq -> L.unpack (unSL (seqid alignedSeq)) == querySeqId) alignedSequences) querySeqIds         -extractQueryCandidates :: [(Sequence,Int,String,Char)] -> V.Vector (Int,Sequence)+extractQueryCandidates :: [(Sequence,Int,L.ByteString)] -> V.Vector (Int,Sequence) extractQueryCandidates querycandidates = indexedSeqences-  where sequences = map (\(candidateSequence,_,_,_) -> candidateSequence) querycandidates+  where sequences = map (\(candidateSequence,_,_) -> candidateSequence) querycandidates         indexedSeqences = V.map (\(number,candidateSequence) -> (number + 1,candidateSequence))(V.indexed (V.fromList sequences))  buildTaxFilterQuery :: Maybe Int -> Maybe Int -> String@@ -916,6 +941,7 @@ randomid :: Int16 -> String randomid number = "cm" ++ show number +-- | Create session id for RNAlien createSessionID :: Maybe String -> IO String createSessionID sessionIdentificator = do   if isJust sessionIdentificator@@ -946,8 +972,8 @@ systemClustalo options (inputFilePath, outputFilePath) = system ("clustalo " ++ options ++ "--infile=" ++ inputFilePath ++ " >" ++ outputFilePath)  -- | Run external CMbuild command and read the output into the corresponding datatype -systemCMbuild ::  String -> String -> String -> IO ExitCode-systemCMbuild alignmentFilepath modelFilepath outputFilePath = system ("cmbuild " ++ modelFilepath ++ " " ++ alignmentFilepath  ++ " > " ++ outputFilePath)  +systemCMbuild ::  String -> String -> String -> String -> IO ExitCode+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@@ -972,7 +998,7 @@ systemCMalign :: String -> String -> String -> String -> IO ExitCode  systemCMalign options filePathCovarianceModel filePathSequence filePathAlignment = system ("cmalign " ++ options ++ " " ++ filePathCovarianceModel ++ " " ++ filePathSequence ++ "> " ++ filePathAlignment) -compareCM :: String -> String -> String -> IO Double+compareCM :: String -> String -> String -> IO (Either String Double) compareCM rfamCMPath resultCMpath outputDirectory = do   let myOptions = defaultDecodeOptions {       decDelimiter = fromIntegral (ord ' ')@@ -988,7 +1014,7 @@   let bitscore1 = read (decodedCmCompareOutput !! 2) :: Double   let bitscore2 = read (decodedCmCompareOutput !! 3) :: Double   let minmax = minimum [bitscore1,bitscore2]-  return minmax+  return (Right minmax)                                                                   readInt :: String -> Int readInt = read@@ -996,233 +1022,9 @@ readDouble :: String -> Double readDouble = read  --- | parse from input filePath              -parseCMSearch :: String -> Either ParseError CMsearch-parseCMSearch = parse genParserCMsearch "parseCMsearch" ---- | parse from input filePath                      -readCMSearch :: String -> IO (Either ParseError CMsearch)             -readCMSearch filePath = do -  parsedFile <- parseFromFile genParserCMsearch filePath-  CE.evaluate parsedFile -                      -genParserCMsearch :: GenParser Char st CMsearch-genParserCMsearch = do-  string "# cmsearch :: search CM(s) against a sequence database"-  newline-  string "# INFERNAL "-  many1 (noneOf "\n")-  newline       -  string "# Copyright (C) 201"-  many1 (noneOf "\n")-  newline       -  string "# Freely distributed under the GNU General Public License (GPLv3)."-  newline       -  string "# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -"-  newline-  string "# query CM file:"-  many1 space-  queryCMfile' <- many1 (noneOf "\n")-  newline-  string "# target sequence database:"-  many1 space      -  targetSequenceDatabase' <- many1 (noneOf "\n")-  newline-  optional (try (genParserCMsearchHeaderField "# CM configuration"))-  optional (try (genParserCMsearchHeaderField "# database size is set to"))-  optional (try (genParserCMsearchHeaderField "# truncated sequence detection"))-  string "# number of worker threads:"-  many1 space-  numberOfWorkerThreads' <- many1 (noneOf "\n")-  newline-  string "# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -"-  newline-  optional newline-  string "Query:"-  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 -  string "E-value"-  many1 space        -  string "score"-  many1 space -  string "bias"-  many1 space -  string "sequence"-  many1 space  -  string "start"-  many1 space -  string "end"-  many1 space -  string "mdl"-  many1 space -  string "trunc"-  many1 space -  string "gc"-  many1 space -  string "description"-  newline-  string " -"-  many1 (try (oneOf " -"))-  newline-  optional (try (string " ------ inclusion threshold ------"))-  many newline-  hitScores' <- many (try genParserCMsearchHit) --`endBy` (try (string "Hit alignments:"))-  optional (try genParserCMsearchEmptyHit)-  -- this is followed by hit alignments and internal cmsearch statistics which are not parsed-  many anyChar-  eof-  return $ CMsearch queryCMfile' targetSequenceDatabase' numberOfWorkerThreads' hitScores'--genParserCMsearchHeaderField :: String -> GenParser Char st String-genParserCMsearchHeaderField fieldname = do-  string (fieldname ++ ":")-  many1 space-  many1 (noneOf "\n")-  newline-  return []--genParserCMsearchEmptyHit :: GenParser Char st [CMsearchHit]-genParserCMsearchEmptyHit = do-  string "   [No hits detected that satisfy reporting thresholds]"-  newline-  optional (try newline)-  return []--genParserCMsearchHit :: GenParser Char st CMsearchHit-genParserCMsearchHit = do-  many1 space-  string "("     -  hitRank' <- many1 digit-  string ")"-  many1 space-  hitSignificant' <- choice [char '!', char '?']-  many1 space                  -  hitEValue' <- many1 (oneOf "0123456789.e-")-  many1 space             -  hitScore'  <- many1 (oneOf "0123456789.e-")-  many1 space   -  hitBias' <- many1 (oneOf "0123456789.e-")-  many1 space-  hitSequenceHeader' <- many1 (noneOf " ")-  many1 space                -  hitStart' <- many1 digit-  many1 space-  hitEnd' <- many1 digit-  many1 space            -  hitStrand' <- choice [char '+', char '-', char '.']-  many1 space              -  hitModel' <- many1 letter-  many1 space          -  hitTruncation' <- many1 (choice [alphaNum, char '\''])-  many1 space                   -  hitGCcontent' <- many1 (oneOf "0123456789.e-")-  many1 space                -  hitDescription' <- many1 (noneOf "\n")     -  newline-  optional (try (string " ------ inclusion threshold ------"))-  optional (try newline)-  return $ CMsearchHit (readInt hitRank') hitSignificant' (readDouble hitEValue') (readDouble hitScore') (readDouble hitBias') (L.pack hitSequenceHeader') (readInt hitStart') (readInt hitEnd') hitStrand' (L.pack hitModel') (L.pack hitTruncation') (readDouble hitGCcontent') (L.pack hitDescription')---- | parse from input filePath              -parseCMstat :: String -> Either ParseError CMstat-parseCMstat = parse genParserCMstat "parseCMstat"---- | parse from input filePath                      -readCMstat :: String -> IO (Either ParseError CMstat)             -readCMstat filePath = do -  parsedFile <- parseFromFile genParserCMstat filePath-  CE.evaluate parsedFile -                      -genParserCMstat :: GenParser Char st CMstat-genParserCMstat = do-  string "# cmstat :: display summary statistics for CMs"-  newline-  string "# INFERNAL "-  many1 (noneOf "\n")-  newline       -  string "# Copyright (C) 201"-  many1 (noneOf "\n")-  newline       -  string "# Freely distributed under the GNU General Public License (GPLv3)."-  newline       -  string "# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -"-  newline-  char '#'-  many1 (char ' ')-  string "rel entropy"-  newline-  char '#'-  many1 (char ' ')-  many1 (char '-')-  newline-  char '#'-  many1 space -  string "idx"-  many1 space        -  string "name"-  many1 space -  string "accession"-  many1 space -  string "nseq"-  many1 space  -  string "eff_nseq"-  many1 space -  string "clen"-  many1 space -  string "W"-  many1 space -  string "bps"-  many1 space -  string "bifs"-  many1 space -  string "model"-  many1 space -  string "cm"-  many1 space-  string "hmm"-  newline-  string "#"-  many1 (try (oneOf " -"))-  newline-  many1 space     -  _statIndex <- many1 digit-  many1 space-  _statName <- many1 letter-  many1 space                  -  _statAccession <- many1 (noneOf " ")-  many1 space             -  _statSequenceNumber <- many1 digit-  many1 space   -  _statEffectiveSequences <- many1 (oneOf "0123456789.e-")-  many1 space-  _statConsensusLength <- many digit-  many1 space                -  _statW <- many1 digit-  many1 space-  _statBasepaires <- many1 digit-  many1 space            -  _statBifurcations <- many1 digit-  many1 space              -  _statModel <- many1 letter-  many1 space          -  _relativeEntropyCM <- many1 (oneOf "0123456789.e-")-  many1 space                   -  _relativeEntropyHMM <- many1 (oneOf "0123456789.e-")-  newline-  char '#'-  newline-  eof  -  return $ CMstat (readInt _statIndex) _statName _statAccession (readInt _statSequenceNumber) (readDouble _statEffectiveSequences) (readInt _statConsensusLength) (readInt _statW) (readInt _statBasepaires) (readInt _statBifurcations) _statModel (readDouble _relativeEntropyCM) (readDouble _relativeEntropyHMM)-   -extractCandidateSequences :: [(Sequence,Int,String,Char)] -> V.Vector (Int,Sequence)+extractCandidateSequences :: [(Sequence,Int,L.ByteString)] -> V.Vector (Int,Sequence) extractCandidateSequences candidates' = indexedSeqences-  where sequences = map (\(inputSequence,_,_,_) -> inputSequence) candidates'+  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)@@ -1277,9 +1079,8 @@          maxIdentity = fromIntegral (maximum (map (snd . Bio.BlastXML.identity) blastMatches))          coverageStatus = (maxIdentity/(fromIntegral queryLength))* (100 :: Double) >= (80 :: Double)          -   -- | Wrapper for retrieveFullSequence that rerequests incomplete return sequees-retrieveFullSequences :: StaticOptions -> [(String,Int,Int,String,String,Int,String)] -> IO [(Sequence,Int,String)]+retrieveFullSequences :: StaticOptions -> [(String,Int,Int,String,T.Text,Int,L.ByteString)] -> IO [(Sequence,Int,L.ByteString)] retrieveFullSequences staticOptions requestedSequences = do   fullSequences <- mapM (retrieveFullSequence (tempDirPath staticOptions)) requestedSequences   if any (isNothing . firstOfTriple) fullSequences@@ -1295,7 +1096,7 @@       CE.evaluate unwrappedRetrievals     else CE.evaluate (map (\(x,y,z) -> (fromJust x,y,z)) fullSequences)         -retrieveFullSequence :: String -> (String,Int,Int,String,String,Int,String) -> IO (Maybe Sequence,Int,String)+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"@@ -1317,7 +1118,7 @@             then return (Nothing,taxid,subject')             else CE.evaluate (Just parsedFasta,taxid,subject')  -getRequestedSequenceElement :: Int -> (BlastHit,Int) -> (String,Int,Int,String,String,Int,String)+getRequestedSequenceElement :: Int -> (BlastHit,Int) -> (String,Int,Int,String,T.Text,Int,L.ByteString) getRequestedSequenceElement queryLength (blastHit,taxid)    | blastHitIsReverseComplement (blastHit,taxid) = getReverseRequestedSequenceElement queryLength (blastHit,taxid)   | otherwise = getForwardRequestedSequenceElement queryLength (blastHit,taxid)@@ -1329,11 +1130,11 @@         firstHSPto = h_to blastMatch         isReverse = firstHSPfrom > firstHSPto -getForwardRequestedSequenceElement :: Int -> (BlastHit,Int) -> (String,Int,Int,String,String,Int,String)+getForwardRequestedSequenceElement :: Int -> (BlastHit,Int) -> (String,Int,Int,String,T.Text,Int,L.ByteString) getForwardRequestedSequenceElement queryLength (blastHit,taxid) = (geneIdentifier',startcoordinate,endcoordinate,strand,accession',taxid,subjectBlast)-   where    accession' = L.unpack (extractAccession blastHit)-            subjectBlast = L.unpack (unSL (subject blastHit))-            geneIdentifier' = extractGeneId blastHit--+   where    accession' = extractAccession blastHit+            subjectBlast = unSL (subject blastHit)+            geneIdentifier' = extractGeneId blastHit             blastMatch = head (matches blastHit)             blastHitOriginSequenceLength = slength blastHit             minHfrom = h_from blastMatch@@ -1372,10 +1173,10 @@   | currentValue > upperBoundry = upperBoundry   | otherwise = currentValue -getReverseRequestedSequenceElement :: Int -> (BlastHit,Int) -> (String,Int,Int,String,String,Int,String)+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' = L.unpack (extractAccession blastHit)-           subjectBlast = L.unpack (unSL (subject blastHit))           +   where   accession' = extractAccession blastHit+           subjectBlast = unSL (subject blastHit)                    geneIdentifier' = extractGeneId blastHit            blastMatch = head (matches blastHit)            blastHitOriginSequenceLength = slength blastHit               @@ -1572,10 +1373,10 @@         hitTaxIdStrings = map extractTaxIdfromDocumentSummary blastHitSummaries         hitTaxIds = map readInt hitTaxIdStrings -extractAccession :: BlastHit -> L.ByteString+extractAccession :: BlastHit -> T.Text extractAccession currentBlastHit = accession'-  where splitedFields = DS.splitOn "|" (L.unpack (hitId currentBlastHit))-        accession' =  L.pack (splitedFields !! 3) +  where splitedFields = T.splitOn (T.pack "|") (DTE.decodeUtf8 (L.toStrict (hitId currentBlastHit)))+        accession' =  splitedFields !! 3          extractGeneId :: BlastHit -> String extractGeneId currentBlastHit = nucleotideId@@ -1615,10 +1416,12 @@ logEither (Left logoutput) temporaryDirectoryPath = appendFile (temporaryDirectoryPath ++ "Log") (show logoutput) logEither  _ _ = return () -checkTools :: [String] -> String -> IO (Either String String)-checkTools tools temporaryDirectoryPath = do+checkTools :: [String] -> String -> String -> IO (Either String String)+checkTools tools temporaryDirectoryPath selectedQuerySelectionMethod = do+  -- if queryselectionmethod is set to clustering then also check for clustal omega+  let additionaltools = if selectedQuerySelectionMethod == "clustering" then tools ++ ["clustalo"] else tools   -- check if all tools are available via PATH or Left-  checks <- mapM checkTool tools+  checks <- mapM checkTool additionaltools   if not (null (lefts checks))     then return (Left (concat (lefts checks)))     else do  @@ -1666,35 +1469,45 @@   | verbosityLevel == Loud = True   | otherwise = False -evaluateConstructionResult :: StaticOptions -> Int -> IO String-evaluateConstructionResult staticOptions entryNumber = do+evaluateConstructionResult :: StaticOptions -> ModelConstruction -> IO String+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 = map nucleotideSequence (concatMap sequenceRecords (taxRecords mCResult))+  let resultNumber = length resultSequences + 1 +  let rnaCentralQueries = map buildSequenceViaMD5Query resultSequences    +  rnaCentralEntries <- getRNACentralEntries rnaCentralQueries+  let rnaCentralEvaluationResult = 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 entryNumber > 1+  if resultNumber > 1     then do        let resultRNAz = tempDirPath staticOptions ++ "result.rnaz"-      rnazClustalpath <- preprocessClustalForRNAzExternal clustalFilepath reformatedClustalPath+      let resultRNAcode = tempDirPath staticOptions ++ "result.rnacode"+      rnazClustalpath <- preprocessClustalForRNAcodeExternal clustalFilepath reformatedClustalPath       if isRight rnazClustalpath         then do           systemRNAz "-l" (fromRight rnazClustalpath) resultRNAz            inputRNAz <- readRNAz resultRNAz           let rnaZString = rnaZEvalOutput inputRNAz-          return ("\nEvaluation of RNAlien result :\nCMstat statistics for result.cm\n" ++ cmstatString ++ "\nRNAz statistics for result alignment: " ++ rnaZString)+          RC.systemRNAcode " -t " (fromRight rnazClustalpath) resultRNAcode+          inputRNAcode <- RC.readRNAcodeTabular resultRNAcode+          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 rnazClustalpath) (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 rnazClustalpath)+          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 rnazClustalpath ++ "\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) -      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")+      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 @@ -1711,27 +1524,54 @@     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 +  | isRight inputRNAcode = rnaCodeString+  | otherwise = show (fromLeft inputRNAcode)+    where rnaCode = fromRight inputRNAcode+          rnaCodeString = "HSS\tFrame\tLength\tFrom\tTo\tName\tStart\tEnd\tScore\tP\n" ++ rnaCodeEntries+          rnaCodeEntries = concatMap showRNACodeHits (RC.rnacodeHits rnaCode)++showRNACodeHits :: RC.RNAcodeHit -> String+showRNACodeHits rnacodeHit = show (RC.hss rnacodeHit) ++ "\t" ++ show (RC.frame rnacodeHit) ++ "\t" ++ show (RC.length 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 -> IO (Either String String) preprocessClustalForRNAzExternal clustalFilepath reformatedClustalPath = do-  clustalString <- readFile clustalFilepath+  clustalText <- TI.readFile clustalFilepath   --change clustal format for rnazSelectSeqs.pl-  let reformatedClustalString = map reformatAln clustalString-  writeFile reformatedClustalPath reformatedClustalString+  let reformatedClustalText = T.map reformatAln clustalText+  TI.writeFile reformatedClustalPath reformatedClustalText   --select representative entries from result.Clustal with select_sequences   let selectedClustalpath = clustalFilepath ++ ".selected"   system ("rnazSelectSeqs.pl " ++ reformatedClustalPath ++ " >" ++ selectedClustalpath)   return (Right selectedClustalpath) +-- | Call for external preprocessClustalForRNAcode - RNAcode additionally to RNAz requirements does not accept pipe,underscore, doublepoint symbols+preprocessClustalForRNAcodeExternal :: String -> String -> IO (Either String String)+preprocessClustalForRNAcodeExternal clustalFilepath reformatedClustalPath = do+  clustalText <- TI.readFile clustalFilepath+  --change clustal format for rnazSelectSeqs.pl+  let clustalTextLines = T.lines clustalText+  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)+  --select representative entries from result.Clustal with select_sequences+  let selectedClustalpath = clustalFilepath ++ ".selected"+  system ("rnazSelectSeqs.pl " ++ reformatedClustalPath ++ " >" ++ selectedClustalpath)+  return (Right selectedClustalpath)+ -- | RNAz can process 500 sequences at max. Using rnazSelectSeqs to isolate representative sample. rnazSelectSeqs only accepts - gap characters, alignment is reformatted accordingly. preprocessClustalForRNAz :: String -> String -> IO (Either String String) preprocessClustalForRNAz clustalFilepath reformatedClustalPath = do-  clustalString <- readFile clustalFilepath-  if length (lines clustalString) > 500+  clustalText <- TI.readFile clustalFilepath+  let clustalTextLines = T.lines clustalText+  if length clustalTextLines > 500     then do        --change clustal format for rnazSelectSeqs.pl-      let reformatedClustalString = map reformatAln clustalString-      writeFile reformatedClustalPath reformatedClustalString+      let reformatedClustalString = T.map reformatAln clustalText+      TI.writeFile reformatedClustalPath reformatedClustalString       --select representative entries from result.Clustal with select_sequences       let selectedClustalpath = clustalFilepath ++ ".selected"       parsedClustalInput <- readClustalAlignment clustalFilepath@@ -1743,6 +1583,7 @@         else return (Left (show (fromLeft parsedClustalInput)))     else return (Right clustalFilepath) + -- Iteratively removes sequences with decreasing similarity until target number of alignment entries is reached. rnaZSelectSeqs :: ClustalAlignment -> Int -> Double -> ClustalAlignment rnaZSelectSeqs currentClustalAlignment targetEntries identityCutoff@@ -1751,7 +1592,21 @@   where numberOfEntries =  length (alignmentEntries currentClustalAlignment)          filteredEntries = filterIdenticalAlignmentEntry (alignmentEntries currentClustalAlignment) identityCutoff          filteredAlignment = ClustalAlignment filteredEntries (conservationTrack currentClustalAlignment)- ++reformatRNACodeAln :: Char -> Char +reformatRNACodeAln c+  | c == ':' = '-'+  | c == '|' = '-'+  | c == '.' = '-'+  | c == '~' = '-'+  | c == '_' = '-'+  | c == 'u' = 'U'+  | c == 't' = 'T'+  | c == 'g' = 'G'+  | c == 'c' = 'C'+  | c == 'a' = 'A'+  | otherwise = c+ reformatAln :: Char -> Char  reformatAln c   | c == '.' = '-'

@@ -2,7 +2,7 @@ {-# LANGUAGE DeriveDataTypeable #-}  -- | Statistics for RNAlien Results--- dist/build/RNAlienStatistics/RNAlienStatistics -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+-- 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      @@ -32,6 +32,8 @@     alienThreshold :: Double,     outputDirectoryPath :: String,     benchmarkIndex :: Int,+    thresholdSelection :: String,+    linkScores :: Bool,     threads :: Int   } deriving (Show,Data,Typeable) @@ -49,14 +51,16 @@     alienThreshold = 20 &= name "t" &= help "Bitscore threshold for RNAlien model hits on Rfam fasta, default 20",     rfamThreshold = 20 &= name "x" &= help "Bitscore threshold for Rfam model hits on Alien fasta, default 20",     benchmarkIndex = 1 &= name "b" &= help "Index used to identify sRNA tagged RNA families",+    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 devel version" &= help "Florian Eggenhofer - >2013" &= verbosity       +  } &= summary "RNAlienStatistics" &= help "Florian Eggenhofer - >2013" &= verbosity         --cmSearchFasta threads rfamCovarianceModelPath outputDirectoryPath "Rfam" False genomesDirectoryPath-cmSearchFasta :: Int -> Double -> Int -> String -> String -> String -> String -> IO [CMsearchHit]-cmSearchFasta benchmarkIndex thresholdScore cpuThreads covarianceModelPath outputDirectory modelType fastapath = do+cmSearchFasta :: Int -> String -> Double -> Int -> String -> String -> String -> String -> IO [CMsearchHit]+cmSearchFasta benchmarkIndex thresholdSelection thresholdScore cpuThreads covarianceModelPath outputDirectory modelType fastapath = do   createDirectoryIfMissing False (outputDirectory ++ "/" ++ modelType)-  _ <- systemCMsearch cpuThreads "" covarianceModelPath fastapath (outputDirectory ++ "/" ++ modelType ++ "/" ++ (show benchmarkIndex) ++ ".cmsearch")+  _ <- systemCMsearch cpuThreads " -Z 1000 " covarianceModelPath fastapath (outputDirectory ++ "/" ++ modelType ++ "/" ++ (show benchmarkIndex) ++ ".cmsearch")   result <- readCMSearch (outputDirectory ++ "/" ++ modelType ++ "/" ++ (show benchmarkIndex) ++ ".cmsearch")   if (isLeft result)      then do@@ -64,14 +68,40 @@        return []      else do        let rightResults = fromRight result-       let (significantHits,_) = partitionCMsearchHitsByScore thresholdScore rightResults+       let significantHits = filterCMsearchHits thresholdSelection thresholdScore rightResults        let organismUniquesignificantHits = nubBy cmSearchSameOrganism significantHits        return organismUniquesignificantHits -partitionCMsearchHitsByScore :: Double -> CMsearch -> ([CMsearchHit],[CMsearchHit])-partitionCMsearchHitsByScore thresholdScore cmSearchResult = (selected,rejected)-  where (selected,rejected) = partition (\hit -> hitScore hit >= thresholdScore) (cmsearchHits cmSearchResult)+--cmSearchFasta threads rfamCovarianceModelPath outputDirectoryPath "Rfam" False genomesDirectoryPath+cmSearchesFasta :: Int -> String -> Double -> Int -> String -> String -> String -> String -> IO [CMsearchHit]+cmSearchesFasta benchmarkIndex thresholdSelection thresholdScore cpuThreads covarianceModelPath outputDirectory modelType fastapath = do+  createDirectoryIfMissing False (outputDirectory ++ "/" ++ modelType)+  _ <- systemCMsearch cpuThreads " -Z 1000 " covarianceModelPath fastapath (outputDirectory ++ "/" ++ modelType ++ "/" ++ (show benchmarkIndex) ++ ".cmsearch")+  result <- readCMSearches (outputDirectory ++ "/" ++ modelType ++ "/" ++ (show benchmarkIndex) ++ ".cmsearch")+  if (isLeft result)+     then do+       print (fromLeft result)+       return []+     else do+       let rightResults = fromRight result+       let significantHits = filterCMsearchHits thresholdSelection thresholdScore rightResults+       let organismUniquesignificantHits = nubBy cmSearchSameOrganism significantHits+       return organismUniquesignificantHits +filterCMsearchHits :: String -> Double -> CMsearch -> [CMsearchHit]+filterCMsearchHits thresholdSelection thresholdScore cmSearchResult+  | thresholdSelection == "bitscore" = bitscorefiltered+  | otherwise =  evaluefiltered+  where bitscorefiltered = filter (\hit -> hitScore hit >= thresholdScore) (cmsearchHits cmSearchResult)+        evaluefiltered = filter (\hit -> hitEvalue hit <= thresholdScore) (cmsearchHits cmSearchResult)++partitionCMsearchHits :: String -> Double -> CMsearch -> ([CMsearchHit],[CMsearchHit])+partitionCMsearchHits thresholdSelection thresholdScore cmSearchResult+  | thresholdSelection == "bitscore" = (bitscoreselected,bitscorerejected)+  | otherwise =  (evalueselected,evaluerejected)+  where (bitscoreselected,bitscorerejected) = partition (\hit -> hitScore hit >= thresholdScore) (cmsearchHits cmSearchResult)+        (evalueselected,evaluerejected) = partition (\hit -> hitEvalue hit <= thresholdScore) (cmsearchHits cmSearchResult)+ trimCMsearchFastaFile :: String -> String -> String -> CMsearch -> String -> IO () trimCMsearchFastaFile genomesDirectory outputFolder modelType cmsearch fastafile  = do   let fastaInputPath = genomesDirectory ++ "/" ++ fastafile@@ -114,17 +144,19 @@   if rfamModelExists     then do       --compute linkscore-      linkscore <- compareCM rfamCovarianceModelPath alienCovarianceModelPath outputDirectoryPath-      rfamMaxLinkScore <- compareCM rfamCovarianceModelPath rfamCovarianceModelPath outputDirectoryPath-      alienMaxLinkscore <- compareCM alienCovarianceModelPath alienCovarianceModelPath outputDirectoryPath+      linkscore <- if linkScores+        then compareCM rfamCovarianceModelPath alienCovarianceModelPath outputDirectoryPath+        else return (Left "-")+      rfamMaxLinkScore <- if linkScores then compareCM rfamCovarianceModelPath rfamCovarianceModelPath outputDirectoryPath else return (Left "-")+      alienMaxLinkscore <- if linkScores then compareCM alienCovarianceModelPath alienCovarianceModelPath outputDirectoryPath else return (Left "-")       _ <- 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")                           let rfamFastaEntriesNumber = read rfamFastaEntries :: Int       let alienFastaEntriesNumber = read alienFastaEntries :: Int-      rfamonAlienResults <- cmSearchFasta benchmarkIndex rfamThreshold threads rfamCovarianceModelPath outputDirectoryPath "rfamOnAlien" alienFastaFilePath -      alienonRfamResults <- cmSearchFasta benchmarkIndex alienThreshold threads alienCovarianceModelPath outputDirectoryPath "alienOnRfam" rfamFastaFilePath  +      rfamonAlienResults <- cmSearchesFasta benchmarkIndex thresholdSelection rfamThreshold threads rfamCovarianceModelPath outputDirectoryPath "rfamOnAlien" alienFastaFilePath +      alienonRfamResults <- cmSearchFasta benchmarkIndex thresholdSelection alienThreshold threads alienCovarianceModelPath outputDirectoryPath "alienOnRfam" rfamFastaFilePath         let rfamonAlienResultsNumber = length rfamonAlienResults       let alienonRfamResultsNumber = length alienonRfamResults       let rfamonAlienRecovery = (fromIntegral rfamonAlienResultsNumber :: Double) / (fromIntegral alienFastaEntriesNumber :: Double)@@ -134,9 +166,9 @@           putStrLn ("BenchmarkIndex: " ++ show benchmarkIndex)           putStrLn ("RfamModelName: " ++ rfamModelName)           putStrLn ("RfamModelId: " ++ rfamModelId)-          putStrLn ("Linkscore: " ++ show linkscore)-          putStrLn ("rfamMaxLinkScore: " ++ show rfamMaxLinkScore)-          putStrLn ("alienMaxLinkscore: " ++ 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 ("rfamFastaEntriesNumber: " ++ show rfamFastaEntriesNumber)@@ -148,24 +180,13 @@           print rnazString           print cmStatString         else do-          putStrLn (show benchmarkIndex ++ "\t" ++ rfamModelName ++ "\t" ++ rfamModelId ++ "\t" ++ show linkscore ++ "\t" ++ show rfamMaxLinkScore ++ "\t" ++ 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" ++ cmStatString)+          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" ++ cmStatString)     else do       --compute linkscore-      --linkscore <- compareCM rfamCovarianceModelPath alienCovarianceModelPath outputDirectoryPath-      --rfamMaxLinkScore <- compareCM rfamCovarianceModelPath rfamCovarianceModelPath outputDirectoryPath-      alienMaxLinkscore <- compareCM alienCovarianceModelPath alienCovarianceModelPath outputDirectoryPath-      --_ <- system ("cat " ++ rfamFastaFilePath ++ " | grep '>' | wc -l >" ++ outputDirectoryPath ++ FP.takeFileName rfamFastaFilePath ++ ".entries")+      alienMaxLinkscore <- if linkScores then compareCM alienCovarianceModelPath alienCovarianceModelPath outputDirectoryPath else return ( Left "-")       _ <- 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")                    -      --let rfamFastaEntriesNumber = read rfamFastaEntries :: Int-      let alienFastaEntriesNumber = read alienFastaEntries :: Int-      --rfamonAlienResults <- cmSearchFasta benchmarkIndex rfamThreshold threads rfamCovarianceModelPath outputDirectoryPath "rfamOnAlien" alienFastaFilePath -      --alienonRfamResults <- cmSearchFasta benchmarkIndex alienThreshold 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)  +      let alienFastaEntriesNumber = read alienFastaEntries :: Int         if (verbose == Loud)         then do           putStrLn ("BenchmarkIndex:")@@ -173,7 +194,7 @@           putStrLn ("RfamModelId: -")           putStrLn ("Linkscore: -")           putStrLn ("rfamMaxLinkScore: -")-          putStrLn ("alienMaxLinkscore: " ++ show alienMaxLinkscore)    +          putStrLn ("alienMaxLinkscore: " ++  (either id show alienMaxLinkscore))               putStrLn ("rfamGatheringThreshold: -")           putStrLn ("alienGatheringThreshold: -")            putStrLn ("rfamFastaEntriesNumber: -")@@ -185,7 +206,7 @@           print rnazString           print cmStatString         else do-          putStrLn (show benchmarkIndex ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ show alienMaxLinkscore ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ show alienFastaEntriesNumber ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-"  ++ "\t" ++ "-" ++ "\t" ++ rnazString ++ "\t" ++ cmStatString)+          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" ++ cmStatString)  rnazOutput :: Verbosity -> String -> IO String rnazOutput verbose rnazPath = do