diff --git a/RNAlien.cabal b/RNAlien.cabal --- a/RNAlien.cabal +++ b/RNAlien.cabal @@ -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 @@ . * . - * - . * . * + . + * . * . @@ -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 diff --git a/src/Bio/InfernalParser.hs b/src/Bio/InfernalParser.hs new file mode 100644 --- /dev/null +++ b/src/Bio/InfernalParser.hs @@ -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 diff --git a/src/Bio/RNAcentralHTTP.hs b/src/Bio/RNAcentralHTTP.hs new file mode 100644 --- /dev/null +++ b/src/Bio/RNAcentralHTTP.hs @@ -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 diff --git a/src/Bio/RNAlien.hs b/src/Bio/RNAlien.hs --- a/src/Bio/RNAlien.hs +++ b/src/Bio/RNAlien.hs @@ -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 diff --git a/src/Bio/RNAlienData.hs b/src/Bio/RNAlienData.hs --- a/src/Bio/RNAlienData.hs +++ b/src/Bio/RNAlienData.hs @@ -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 } diff --git a/src/Bio/RNAlienLibrary.hs b/src/Bio/RNAlienLibrary.hs --- a/src/Bio/RNAlienLibrary.hs +++ b/src/Bio/RNAlienLibrary.hs @@ -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 == '.' = '-' diff --git a/src/Bio/RNAlienStatistics.hs b/src/Bio/RNAlienStatistics.hs --- a/src/Bio/RNAlienStatistics.hs +++ b/src/Bio/RNAlienStatistics.hs @@ -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