RNAlien 1.1.1 → 1.1.3
raw patch · 4 files changed
+107/−35 lines, 4 files
Files
- RNAlien.cabal +3/−3
- src/Bio/RNAlien.hs +3/−3
- src/Bio/RNAlienData.hs +2/−2
- src/Bio/RNAlienLibrary.hs +99/−27
RNAlien.cabal view
@@ -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/
src/Bio/RNAlien.hs view
@@ -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
src/Bio/RNAlienData.hs view
@@ -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
src/Bio/RNAlienLibrary.hs view
@@ -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+