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.1.1 +version: 1.1.3 synopsis: Unsupervized construction of RNA family models description: RNAlien is a tool for automatic construction of RNAfamily models from a single sequence. . @@ -53,8 +53,8 @@ source-repository this type: git - location: https://github.com/eggzilla/RNAlien/tree/1.1.1 - tag: 1.1.1 + location: https://github.com/eggzilla/RNAlien/tree/1.1.3 + tag: 1.1.3 executable RNAlien Hs-Source-Dirs: ./src/Bio/ diff --git a/src/Bio/RNAlien.hs b/src/Bio/RNAlien.hs --- a/src/Bio/RNAlien.hs +++ b/src/Bio/RNAlien.hs @@ -45,13 +45,13 @@ lengthFilter = True &= name "l" &= help "Filter blast hits per genomic length. Default: True", coverageFilter = True &= name "a" &= help "Filter blast hits by coverage of at least 80%. Default: True", singleHitperTax = False &= name "s" &= help "Only the best blast hit per taxonomic entry is considered. Default: False", - blastSoftmasking = True &= name "f" &= help "Toggles blast softmasking, meaning exclusion of low complexity (repetative) regions in lookup table. Default: True", + blastSoftmasking = False &= name "f" &= help "Toggles blast query softmasking, meaning masking of non-conserved regions on the query. 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.1.1" &= help "Florian Eggenhofer, Ivo L. Hofacker, Christian Höner zu Siederdissen - 2013 - 2016" &= verbosity + } &= summary "RNAlien version 1.1.2" &= help "Florian Eggenhofer, Ivo L. Hofacker, Christian Höner zu Siederdissen - 2013 - 2016" &= verbosity main :: IO () main = do @@ -64,7 +64,7 @@ createDirectoryIfMissing False temporaryDirectoryPath createDirectoryIfMissing False (temporaryDirectoryPath ++ "log") -- Create Log files - writeFile (temporaryDirectoryPath ++ "Log") ("RNAlien 1.1.1" ++ "\n") + writeFile (temporaryDirectoryPath ++ "Log") ("RNAlien 1.1.2" ++ "\n") writeFile (temporaryDirectoryPath ++ "log/warnings") ("") logMessage ("Timestamp: " ++ (show timestamp) ++ "\n") temporaryDirectoryPath logMessage ("Temporary Directory: " ++ temporaryDirectoryPath ++ "\n") temporaryDirectoryPath diff --git a/src/Bio/RNAlienData.hs b/src/Bio/RNAlienData.hs --- a/src/Bio/RNAlienData.hs +++ b/src/Bio/RNAlienData.hs @@ -35,7 +35,7 @@ taxonomicContext :: Maybe Taxon, evalueThreshold :: Double, alignmentModeInfernal :: Bool, - selectedQueries :: [String], + selectedQueries :: [Sequence], potentialMembers :: [SearchResult] } @@ -47,7 +47,7 @@ d = "Upper taxonomy limit: " ++ maybe "not set" show _upperTaxonomyLimit ++ "\n" e = "Taxonomic Context: " ++ maybe "not set" show _taxonomicContext ++ "\n" g = "Evalue cutoff: " ++ show _evalueThreshold ++ "\n" - h = "Selected queries: \n" ++ concatMap (\x -> x ++ "\n") _selectedQueries + h = "Selected queries: \n" ++ concatMap show _selectedQueries i = "Potential Members: \n" ++ concatMap show _potentialMembers data TaxonomyRecord = TaxonomyRecord diff --git a/src/Bio/RNAlienLibrary.hs b/src/Bio/RNAlienLibrary.hs --- a/src/Bio/RNAlienLibrary.hs +++ b/src/Bio/RNAlienLibrary.hs @@ -70,6 +70,9 @@ import qualified Data.Text as T import qualified Data.Text.IO as TI import qualified Data.Text.Encoding as DTE +import qualified Data.Text.Lazy.Encoding as E +import qualified Data.Text.Lazy as TL +import qualified Data.Text.Lazy.IO as TIO -- | Initial RNA family model construction - generates iteration number, seed alignment and model modelConstructer :: StaticOptions -> ModelConstruction -> IO ModelConstruction @@ -306,28 +309,31 @@ nextModelConstruction <- modelConstructer staticOptions nextModelConstructionInputWithThreshold return nextModelConstruction else do - --select queries - currentSelectedQueries <- selectQueries staticOptions modelConstruction alignmentResults if (alignmentModeInfernal modelConstruction) then do logVerboseMessage (verbositySwitch staticOptions) ("Alignment construction with candidates - infernal mode\n") (tempDirPath staticOptions) --prepare next iteration - let nextModelConstructionInput = constructNext currentIterationNumber modelConstruction alignmentResults currentUpperTaxonomyLimit currentTaxonomicContext currentSelectedQueries currentPotentialMembers True + let nextModelConstructionInput = constructNext currentIterationNumber modelConstruction alignmentResults currentUpperTaxonomyLimit currentTaxonomicContext [] currentPotentialMembers True constructModel nextModelConstructionInput staticOptions writeFile (iterationDirectory ++ "done") "" logMessage (iterationSummaryLog nextModelConstructionInput) (tempDirPath staticOptions) - logVerboseMessage (verbositySwitch staticOptions) (show nextModelConstructionInput) (tempDirPath staticOptions) ---- - nextModelConstruction <- modelConstructer staticOptions nextModelConstructionInput + logVerboseMessage (verbositySwitch staticOptions) (show nextModelConstructionInput) (tempDirPath staticOptions) + --select queries + currentSelectedQueries <- selectQueries staticOptions modelConstruction alignmentResults + let nextModelConstructionInputWithQueries = nextModelConstructionInput {selectedQueries = currentSelectedQueries} + nextModelConstruction <- modelConstructer staticOptions nextModelConstructionInputWithQueries return nextModelConstruction else do logVerboseMessage (verbositySwitch staticOptions) ("Alignment construction with candidates - initial mode\n") (tempDirPath staticOptions) --First round enough candidates are available for modelconstruction, alignmentModeInfernal is set to true after this iteration --prepare next iteration - let nextModelConstructionInput = constructNext currentIterationNumber modelConstruction alignmentResults currentUpperTaxonomyLimit currentTaxonomicContext currentSelectedQueries currentPotentialMembers False - constructModel nextModelConstructionInput staticOptions - let nextModelConstructionInputWithInfernalMode = nextModelConstructionInput {alignmentModeInfernal = True} + let nextModelConstructionInput = constructNext currentIterationNumber modelConstruction alignmentResults currentUpperTaxonomyLimit currentTaxonomicContext [] currentPotentialMembers False + constructModel nextModelConstructionInput staticOptions + currentSelectedQueries <- selectQueries staticOptions modelConstruction alignmentResults + --select queries + let nextModelConstructionInputWithInfernalMode = nextModelConstructionInput {alignmentModeInfernal = True, selectedQueries = currentSelectedQueries} logMessage (iterationSummaryLog nextModelConstructionInputWithInfernalMode) (tempDirPath staticOptions) - logVerboseMessage (verbositySwitch staticOptions) (show nextModelConstructionInputWithInfernalMode) (tempDirPath staticOptions) ---- + logVerboseMessage (verbositySwitch staticOptions) (show nextModelConstructionInputWithInfernalMode) (tempDirPath staticOptions) writeFile (iterationDirectory ++ "done") "" nextModelConstruction <- modelConstructer staticOptions nextModelConstructionInputWithInfernalMode return nextModelConstruction @@ -395,17 +401,17 @@ else error "Find taxonomy start: Could not find blast hits to use as a taxonomic starting point" searchCandidates :: StaticOptions -> Maybe String -> Int -> Maybe Int -> Maybe Int -> Double -> [Sequence] -> IO SearchResult -searchCandidates staticOptions finaliterationprefix iterationnumber upperTaxLimit lowerTaxLimit expectThreshold querySequences' = do +searchCandidates staticOptions finaliterationprefix iterationnumber upperTaxLimit lowerTaxLimit expectThreshold inputQuerySequences = do --let fastaSeqData = seqdata _querySequence - if (null querySequences') then error "searchCandidates: - head: empty list of query sequences" else return () - let queryLength = fromIntegral (seqlength (head querySequences')) + if (null inputQuerySequences) then error "searchCandidates: - head: empty list of query sequences" else return () + let queryLength = fromIntegral (seqlength (head inputQuerySequences)) let queryIndexString = "1" let entrezTaxFilter = buildTaxFilterQuery upperTaxLimit lowerTaxLimit logVerboseMessage (verbositySwitch staticOptions) ("entrezTaxFilter" ++ show entrezTaxFilter ++ "\n") (tempDirPath staticOptions) let hitNumberQuery = buildHitNumberQuery "&HITLIST_SIZE=5000&EXPECT=" ++ show expectThreshold let registrationInfo = buildRegistration "RNAlien" "florian.eggenhofer@univie.ac.at" let softmaskFilter = if (blastSoftmaskingToggle staticOptions) then "&FILTER=True&FILTER=m" else "" - let blastQuery = BlastHTTPQuery (Just "ncbi") (Just "blastn") (blastDatabase staticOptions) querySequences' (Just (hitNumberQuery ++ entrezTaxFilter ++ softmaskFilter ++ registrationInfo)) (Just (5400000000 :: Int)) + let blastQuery = BlastHTTPQuery (Just "ncbi") (Just "blastn") (blastDatabase staticOptions) inputQuerySequences (Just (hitNumberQuery ++ entrezTaxFilter ++ softmaskFilter ++ registrationInfo)) (Just (5400000000 :: Int)) --appendFile "/scratch/egg/blasttest/queries" ("\nBlast query:\n" ++ show blastQuery ++ "\n") logVerboseMessage (verbositySwitch staticOptions) ("Sending blast query " ++ (show iterationnumber) ++ "\n") (tempDirPath staticOptions) blastOutput <- CE.catch (blastHTTP blastQuery) @@ -577,13 +583,14 @@ where currentClusterNumber = length (cutAt clustaloDendrogram currentCutoff) -- 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 -> [(Sequence,Int,L.ByteString)] -> IO [Sequence] selectQueries staticOptions modelConstruction selectedCandidates = do logVerboseMessage (verbositySwitch staticOptions) "SelectQueries\n" (tempDirPath staticOptions) --Extract sequences from modelconstruction let alignedSequences = extractAlignedSequences (iterationNumber modelConstruction) modelConstruction let candidateSequences = extractQueryCandidates selectedCandidates let iterationDirectory = tempDirPath staticOptions ++ show (iterationNumber modelConstruction) ++ "/" + let stockholmFilepath = iterationDirectory ++ "model" ++ ".stockholm" let alignmentSequences = map snd (V.toList (V.concat [candidateSequences,alignedSequences])) if length alignmentSequences > 3 then do @@ -611,17 +618,32 @@ 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 + let currentSelectedSequenceIds = map L.pack (take (queryNumber staticOptions) (concatMap (take 1 . elements) cutDendrogram)) + --let alignedSequences = fastaSeqData:map nucleotideSequence (concatMap sequenceRecords (taxRecords modelconstruction)) + let fastaSelectedSequences = concatMap (filterSequenceById alignmentSequences) currentSelectedSequenceIds + stockholmSelectedSequences <- extractAlignmentSequencesByIds stockholmFilepath currentSelectedSequenceIds + --Stockholm sequnces contain conservation annotation from cmalign in infernal mode + let currentSelectedSequences = if (blastSoftmaskingToggle staticOptions) then stockholmSelectedSequences else fastaSelectedSequences + --let currentSelectedQueries = concatMap (\querySeqId -> filter (\alignedSeq -> L.unpack (unSL (seqid alignedSeq)) == querySeqId) alignmentSequences) querySeqIds + logVerboseMessage (verbositySwitch staticOptions) ("SelectedQueries: " ++ show currentSelectedSequences ++ "\n") (tempDirPath staticOptions) + writeFile (tempDirPath staticOptions ++ show (iterationNumber modelConstruction) ++ "/log" ++ "/13selectedQueries") (showlines currentSelectedSequences) + CE.evaluate currentSelectedSequences else do - let currentSelectedSequences = filterIdenticalSequences' alignmentSequences (95 :: Double) - let currentSelectedQueries = map (L.unpack . unSL . seqid) (take (queryNumber staticOptions) currentSelectedSequences) - CE.evaluate currentSelectedQueries - + let fastaSelectedSequences = filterIdenticalSequences' alignmentSequences (95 :: Double) + let currentSelectedSequenceIds = map (unSL . seqid) (take (queryNumber staticOptions) fastaSelectedSequences) + stockholmSelectedSequences <- extractAlignmentSequencesByIds stockholmFilepath currentSelectedSequenceIds + let currentSelectedSequences = if (blastSoftmaskingToggle staticOptions) then stockholmSelectedSequences else fastaSelectedSequences + writeFile (tempDirPath staticOptions ++ show (iterationNumber modelConstruction) ++ "/log" ++ "/13selectedQueries") (showlines currentSelectedSequences) + CE.evaluate currentSelectedSequences else return [] + +filterSequenceById :: [Sequence] -> L.ByteString-> [Sequence] +filterSequenceById alignmentSequences querySequenceId = filter (seqenceHasId querySequenceId) alignmentSequences + +seqenceHasId :: L.ByteString -> Sequence -> Bool +seqenceHasId querySequenceId alignmentSequence = unSL (seqid alignmentSequence) == querySequenceId + constructModel :: ModelConstruction -> StaticOptions -> IO String constructModel modelConstruction staticOptions = do --Extract sequences from modelconstruction @@ -856,7 +878,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,L.ByteString)] -> Maybe Int -> Maybe Taxon -> [String] -> [SearchResult] -> Bool -> ModelConstruction +constructNext :: Int -> ModelConstruction -> [(Sequence,Int,L.ByteString)] -> Maybe Int -> Maybe Taxon -> [Sequence] -> [SearchResult] -> Bool -> ModelConstruction constructNext currentIterationNumber modelconstruction alignmentResults upperTaxLimit inputTaxonomicContext inputSelectedQueries inputPotentialMembers toggleInfernalAlignmentModeTrue = nextModelConstruction where newIterationNumber = currentIterationNumber + 1 taxEntries = taxRecords modelconstruction ++ buildTaxRecords alignmentResults currentIterationNumber @@ -914,9 +936,7 @@ | foundSequenceNumber < 3 = [fastaSeqData] | otherwise = querySequences' where fastaSeqData = inputFasta modelconstruction - querySeqIds = selectedQueries modelconstruction - alignedSequences = fastaSeqData:map nucleotideSequence (concatMap sequenceRecords (taxRecords modelconstruction)) - querySequences' = concatMap (\querySeqId -> filter (\alignedSeq -> L.unpack (unSL (seqid alignedSeq)) == querySeqId) alignedSequences) querySeqIds + querySequences' = selectedQueries modelconstruction extractQueryCandidates :: [(Sequence,Int,L.ByteString)] -> V.Vector (Int,Sequence) extractQueryCandidates querycandidates = indexedSeqences @@ -1462,8 +1482,18 @@ csvtable = tableheader ++ tablebody constructTaxonomyRecordCSVEntries :: TaxonomyRecord -> String -constructTaxonomyRecordCSVEntries taxRecord = concatMap (\seqrec -> show (recordTaxonomyId taxRecord) ++ ";" ++ show (aligned seqrec) ++ ";" ++ filter (/= ';') (L.unpack (unSL (seqheader (nucleotideSequence seqrec)))) ++ "\n") (sequenceRecords taxRecord) +constructTaxonomyRecordCSVEntries taxRecord = concatMap (constructTaxonomyRecordCSVEntry taxIdString) (sequenceRecords taxRecord) + where taxIdString = show (recordTaxonomyId taxRecord) +constructTaxonomyRecordCSVEntry :: String -> SequenceRecord -> String +constructTaxonomyRecordCSVEntry taxIdString seqrec = taxIdString ++ ";" ++ show (aligned seqrec) ++ ";" ++ filter checkTaxonomyRecordCSVChar (L.unpack (unSL (seqheader (nucleotideSequence seqrec)))) ++ "\n" + +checkTaxonomyRecordCSVChar :: Char -> Bool +checkTaxonomyRecordCSVChar c + | c == '"' = False + | c == ';' = False + | otherwise = True + setVerbose :: Verbosity -> Bool setVerbose verbosityLevel | verbosityLevel == Loud = True @@ -1672,3 +1702,45 @@ | restrictionString == "bacteria" = Just "bacteria" | restrictionString == "eukaryia" = Just "eukaryia" | otherwise = Nothing + +extractAlignmentSequencesByIds :: String -> [L.ByteString] -> IO [Sequence] +extractAlignmentSequencesByIds stockholmFilePath sequenceIds = do + inputSeedAln <- TIO.readFile stockholmFilePath + let alnEntries = extractAlignmentSequences inputSeedAln + --let splitIds = map E.encodeUtf8 (TL.splitOn (TL.pack ",") (TL.pack sequenceIds)) + let filteredEntries = concatMap (filterSequencesById alnEntries) sequenceIds + return filteredEntries + +extractAlignmentSequences :: TL.Text -> [Sequence] +extractAlignmentSequences seedFamilyAln = rfamIDAndseedFamilySequences + where seedFamilyAlnLines = TL.lines seedFamilyAln + -- remove empty lines from splitting + seedFamilyNonEmpty = filter (\alnline -> not (TL.empty == alnline)) seedFamilyAlnLines + -- remove annotation and spacer lines + seedFamilyIdSeqLines = filter (\alnline -> ((not ((TL.head alnline) == '#'))) && (not ((TL.head alnline) == ' ')) && (not ((TL.head alnline) == '/'))) seedFamilyNonEmpty + -- put id and corresponding seq of each line into a list and remove whitspaces + seedFamilyIdandSeqLines = map TL.words seedFamilyIdSeqLines + -- linewise tuples with id and seq without alinment characters - . + seedFamilyIdandSeqLineTuples = map (\alnline -> ((head alnline),(filterAlnChars (last alnline)))) seedFamilyIdandSeqLines + -- line tuples sorted by id + seedFamilyIdandSeqTupleSorted = sortBy (\tuple1 tuple2 -> compare (fst tuple1) (fst tuple2)) seedFamilyIdandSeqLineTuples + -- line tuples grouped by id + seedFamilyIdandSeqTupleGroups = groupBy (\tuple1 tuple2 -> (fst tuple1) == (fst tuple2)) seedFamilyIdandSeqTupleSorted + seedFamilySequences = map mergeIdSeqTuplestoSequence seedFamilyIdandSeqTupleGroups + rfamIDAndseedFamilySequences = seedFamilySequences + +filterSequencesById :: [Sequence] -> L.ByteString -> [Sequence] +filterSequencesById alignmentSequences sequenceId = filter (sequenceHasId sequenceId) alignmentSequences + +sequenceHasId :: L.ByteString -> Sequence -> Bool +sequenceHasId sequenceId currentSequence = sequenceId == (unSL (seqid currentSequence)) + +filterAlnChars :: TL.Text -> TL.Text +filterAlnChars cs = TL.filter (\c -> (not (c == '-')) && (not (c == '.'))) cs + +mergeIdSeqTuplestoSequence :: [(TL.Text,TL.Text)] -> Sequence +mergeIdSeqTuplestoSequence tuplelist = currentSequence + where seqId = fst (head tuplelist) + seqData = TL.concat (map snd tuplelist) + currentSequence = Seq (SeqLabel (E.encodeUtf8 seqId)) (SeqData (E.encodeUtf8 seqData)) Nothing +