RNAlien 1.6.0 → 1.7.0
raw patch · 5 files changed
+474/−112 lines, 5 filesdep +silentlydep ~BiobaseBlastdep ~BlastHTTPdep ~aesonnew-component:exe:RNAlienScanPVP ok
version bump matches the API change (PVP)
Dependencies added: silently
Dependency ranges changed: BiobaseBlast, BlastHTTP, aeson
API changes (from Hackage documentation)
+ Biobase.RNAlien.Library: scanModelConstructer :: StaticOptions -> ModelConstruction -> IO ModelConstruction
+ Biobase.RNAlien.Types: [genomeFastas] :: ModelConstruction -> [Fasta () ()]
- Biobase.RNAlien.Types: ModelConstruction :: Int -> [Fasta () ()] -> [TaxonomyRecord] -> Maybe Int -> Maybe Taxon -> Double -> Bool -> [Fasta () ()] -> [SearchResult] -> ModelConstruction
+ Biobase.RNAlien.Types: ModelConstruction :: Int -> [Fasta () ()] -> [TaxonomyRecord] -> Maybe Int -> Maybe Taxon -> Double -> Bool -> [Fasta () ()] -> [SearchResult] -> [Fasta () ()] -> ModelConstruction
Files
- Biobase/RNAlien.hs +6/−5
- Biobase/RNAlien/Library.hs +322/−99
- Biobase/RNAlien/Types.hs +4/−2
- Biobase/RNAlienScan.hs +126/−0
- RNAlien.cabal +16/−6
Biobase/RNAlien.hs view
@@ -4,7 +4,7 @@ -- | Unsupervized construction of RNA family models -- For more information on RNA family models consult <http://> -- Usage example: RNAlien -i /path/input.fa -c 5 -o /outdir/--- Usage example offline mode: RNAlien -i /path/input.fa -b /backup/blast/nt_v5 -o /outdir/ -c 5 -t 1396 -j +-- Usage example offline mode: RNAlien -i /path/input.fa -b /backup/blast/nt_v5 -o /outdir/ -c 5 -t 1396 -j module Main where import System.Console.CmdArgs@@ -94,8 +94,8 @@ writeFile (temporaryDirectoryPath ++ "log/warnings") ("") logMessage ("Timestamp: " ++ (show timestamp) ++ "\n") temporaryDirectoryPath logMessage ("Temporary Directory: " ++ temporaryDirectoryPath ++ "\n") temporaryDirectoryPath- inputFasta <- readFastaFile inputFastaFilePath- if null inputFasta+ fastaInput <- readFastaFile inputFastaFilePath+ if null fastaInput then do putStrLn "Error: Input fasta file is empty." logMessage "Error: Input fasta file is empty.\n" temporaryDirectoryPath@@ -108,11 +108,12 @@ logMessage ("Error - Not all required tools could be found in $PATH: " ++ fromLeft toolsCheck ++ "\n") temporaryDirectoryPath else do logToolVersions inputQuerySelectionMethod temporaryDirectoryPath- let inputSequence = reformatFasta (head inputFasta)+ let reformatedFastaInput = map reformatFasta fastaInput+ let inputSequence = head reformatedFastaInput initialTaxId <- setInitialTaxId offlineMode threads inputBlastDatabase temporaryDirectoryPath inputTaxId inputSequence let checkedTaxonomyRestriction = checkTaxonomyRestriction taxonomyRestriction let staticOptions = StaticOptions temporaryDirectoryPath sessionId (fromJust inputnSCICutoff) inputTaxId singleHitperTax inputQuerySelectionMethod inputQueryNumber lengthFilter coverageFilter blastSoftmasking threads inputBlastDatabase checkedTaxonomyRestriction (setVerbose verboseLevel) offlineMode- let initialization = ModelConstruction iterationNumber inputFasta [] initialTaxId Nothing (fromJust inputEvalueCutoff) False [] []+ let initialization = ModelConstruction iterationNumber reformatedFastaInput [] initialTaxId Nothing (fromJust inputEvalueCutoff) False [] [] [] logMessage (show initialization) temporaryDirectoryPath modelConstructionResults <- modelConstructer staticOptions initialization let resultTaxonomyRecordsCSVTable = constructTaxonomyRecordsCSVTable modelConstructionResults
Biobase/RNAlien/Library.hs view
@@ -6,6 +6,7 @@ logMessage, logEither, modelConstructer,+ scanModelConstructer, constructTaxonomyRecordsCSVTable, resultSummary, setVerbose,@@ -83,6 +84,7 @@ import Data.Foldable import Biobase.Types.BioSequence import qualified Biobase.BLAST.Import as BBI+import System.IO.Silently -- | Initial RNA family model construction - generates iteration number, seed alignment and model modelConstructer :: StaticOptions -> ModelConstruction -> IO ModelConstruction@@ -105,15 +107,15 @@ logVerboseMessage (verbositySwitch staticOptions) ("Upper taxonomy limit: " ++ show upperTaxLimit ++ "\n " ++ "Lower taxonomy limit: "++ show lowerTaxLimit ++ "\n") (tempDirPath staticOptions) --search queries let expectThreshold = setBlastExpectThreshold modelConstruction- searchResults <- catchAll (searchCandidates staticOptions Nothing currentIterationNumber upperTaxLimit lowerTaxLimit expectThreshold queries)+ searchResults <- catchAll (searchCandidates staticOptions Nothing currentIterationNumber upperTaxLimit lowerTaxLimit expectThreshold Nothing queries) (\e -> do logWarning ("Warning: Search results iteration" ++ show (iterationNumber modelConstruction) ++ " - exception: " ++ show e) (tempDirPath staticOptions) return (SearchResult [] Nothing)) currentTaxonomicContext <- getTaxonomicContextEntrez upperTaxLimit (taxonomicContext modelConstruction) if null (candidates searchResults) then- alignmentConstructionWithoutCandidates currentTaxonomicContext upperTaxLimit staticOptions modelConstruction+ alignmentConstructionWithoutCandidates "alien" currentTaxonomicContext upperTaxLimit staticOptions modelConstruction else- alignmentConstructionWithCandidates currentTaxonomicContext upperTaxLimit searchResults staticOptions modelConstruction+ alignmentConstructionWithCandidates "alien" currentTaxonomicContext upperTaxLimit searchResults staticOptions modelConstruction else do logMessage "Message: Modelconstruction complete: Out of queries or taxonomic tree exhausted\n" (tempDirPath staticOptions) modelConstructionResult staticOptions modelConstruction@@ -158,7 +160,7 @@ then do --taxonomic restriction let (upperTaxLimit,lowerTaxLimit) = setRestrictedTaxonomyLimits (fromJust (taxRestriction staticOptions))- restrictedCandidates <- catchAll (searchCandidates staticOptions (taxRestriction staticOptions) currentIterationNumber upperTaxLimit lowerTaxLimit expectThreshold queries)+ restrictedCandidates <- catchAll (searchCandidates staticOptions (taxRestriction staticOptions) currentIterationNumber upperTaxLimit lowerTaxLimit expectThreshold Nothing queries) (\e -> do logWarning ("Warning: Search results iteration" ++ show currentIterationNumber ++ " - exception: " ++ show e) outputDirectory return (SearchResult [] Nothing)) let uniqueCandidates = filterDuplicates modelConstruction restrictedCandidates@@ -170,7 +172,7 @@ else do --taxonomic context archea let (upperTaxLimit1,lowerTaxLimit1) = (Just (2157 :: Int), Nothing)- candidates1 <- catchAll (searchCandidates staticOptions (Just "archea") currentIterationNumber upperTaxLimit1 lowerTaxLimit1 expectThreshold queries)+ candidates1 <- catchAll (searchCandidates staticOptions (Just "archea") currentIterationNumber upperTaxLimit1 lowerTaxLimit1 expectThreshold Nothing queries) (\e -> do logWarning ("Warning: Search results iteration" ++ show currentIterationNumber ++ " - exception: " ++ show e) outputDirectory return (SearchResult [] Nothing)) let uniqueCandidates1 = filterDuplicates modelConstruction candidates1@@ -179,7 +181,7 @@ return ([],[])) --taxonomic context bacteria let (upperTaxLimit2,lowerTaxLimit2) = (Just (2 :: Int), Nothing)- candidates2 <- catchAll (searchCandidates staticOptions (Just "bacteria") currentIterationNumber upperTaxLimit2 lowerTaxLimit2 expectThreshold queries)+ candidates2 <- catchAll (searchCandidates staticOptions (Just "bacteria") currentIterationNumber upperTaxLimit2 lowerTaxLimit2 expectThreshold Nothing queries) (\e -> do logWarning ("Warning: Search results iteration" ++ show currentIterationNumber ++ " - exception: " ++ show e) outputDirectory return (SearchResult [] Nothing)) let uniqueCandidates2 = filterDuplicates modelConstruction candidates2@@ -188,7 +190,7 @@ return ([],[])) --taxonomic context eukaryia let (upperTaxLimit3,lowerTaxLimit3) = (Just (2759 :: Int), Nothing)- candidates3 <- catchAll (searchCandidates staticOptions (Just "eukaryia") currentIterationNumber upperTaxLimit3 lowerTaxLimit3 expectThreshold queries)+ candidates3 <- catchAll (searchCandidates staticOptions (Just "eukaryia") currentIterationNumber upperTaxLimit3 lowerTaxLimit3 expectThreshold Nothing queries) (\e -> do logWarning ("Warning: Search results iteration" ++ show currentIterationNumber ++ " - exception: " ++ show e) outputDirectory return (SearchResult [] Nothing)) let uniqueCandidates3 = filterDuplicates modelConstruction candidates3@@ -301,8 +303,8 @@ --------------------------------------------------------- -alignmentConstructionWithCandidates :: Maybe Taxon -> Maybe Int -> SearchResult -> StaticOptions -> ModelConstruction -> IO ModelConstruction-alignmentConstructionWithCandidates currentTaxonomicContext currentUpperTaxonomyLimit searchResults staticOptions modelConstruction = do+alignmentConstructionWithCandidates :: String -> Maybe Taxon -> Maybe Int -> SearchResult -> StaticOptions -> ModelConstruction -> IO ModelConstruction+alignmentConstructionWithCandidates alienType currentTaxonomicContext currentUpperTaxonomyLimit searchResults staticOptions modelConstruction = do --candidates usedUpperTaxonomyLimit blastDatabaseSize let currentIterationNumber = iterationNumber modelConstruction let iterationDirectory = tempDirPath staticOptions ++ show currentIterationNumber ++ "/"@@ -320,10 +322,18 @@ --prepare next iteration let newTaxEntries = taxRecords modelConstruction ++ buildTaxRecords alignmentResults currentIterationNumber let nextModelConstructionInputWithThreshold = modelConstruction {iterationNumber = currentIterationNumber + 1,upperTaxonomyLimit = currentUpperTaxonomyLimit, taxRecords = newTaxEntries,taxonomicContext = currentTaxonomicContext}- logMessage (iterationSummaryLog nextModelConstructionInputWithThreshold) (tempDirPath staticOptions)- logVerboseMessage (verbositySwitch staticOptions) (show nextModelConstructionInputWithThreshold) (tempDirPath staticOptions) ----+ let nextGenomeFastas = tail (genomeFastas modelConstruction)+ let nextScanModelConstructionInputWithThreshold = modelConstruction {iterationNumber = currentIterationNumber + 1,upperTaxonomyLimit = currentUpperTaxonomyLimit, taxRecords = newTaxEntries,taxonomicContext = currentTaxonomicContext, genomeFastas = nextGenomeFastas} writeFile (iterationDirectory ++ "done") ""- modelConstructer staticOptions nextModelConstructionInputWithThreshold+ if alienType == "alien"+ then do+ logMessage (iterationSummaryLog nextModelConstructionInputWithThreshold) (tempDirPath staticOptions)+ logVerboseMessage (verbositySwitch staticOptions) (show nextModelConstructionInputWithThreshold) (tempDirPath staticOptions)+ modelConstructer staticOptions nextModelConstructionInputWithThreshold+ else do+ logMessage (iterationSummaryLog nextScanModelConstructionInputWithThreshold) (tempDirPath staticOptions)+ logVerboseMessage (verbositySwitch staticOptions) (show nextScanModelConstructionInputWithThreshold) (tempDirPath staticOptions)+ scanModelConstructer staticOptions nextScanModelConstructionInputWithThreshold else if (alignmentModeInfernal modelConstruction) then do@@ -332,12 +342,21 @@ 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) --select queries currentSelectedQueries <- selectQueries staticOptions modelConstruction alignmentResults- let nextModelConstructionInputWithQueries = nextModelConstructionInput {selectedQueries = currentSelectedQueries}- nextModelConstruction <- modelConstructer staticOptions nextModelConstructionInputWithQueries+ nextModelConstruction <- if alienType == "alien"+ then do+ logMessage (iterationSummaryLog nextModelConstructionInput) (tempDirPath staticOptions)+ logVerboseMessage (verbositySwitch staticOptions) (show nextModelConstructionInput) (tempDirPath staticOptions)+ let nextModelConstructionInputWithQueries = nextModelConstructionInput {selectedQueries = currentSelectedQueries}+ modelConstructer staticOptions nextModelConstructionInputWithQueries+ else do+ let nextGenomeFastas = tail (genomeFastas modelConstruction)+ let nextScanModelConstructionInputWithQueries = nextModelConstructionInput {selectedQueries = currentSelectedQueries, genomeFastas = nextGenomeFastas}+ logMessage (iterationSummaryLog nextScanModelConstructionInputWithQueries) (tempDirPath staticOptions)+ logVerboseMessage (verbositySwitch staticOptions) (show nextScanModelConstructionInputWithQueries) (tempDirPath staticOptions)++ scanModelConstructer staticOptions nextScanModelConstructionInputWithQueries return nextModelConstruction else do logVerboseMessage (verbositySwitch staticOptions) ("Alignment construction with candidates - initial mode\n") (tempDirPath staticOptions)@@ -347,22 +366,32 @@ 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) writeFile (iterationDirectory ++ "done") ""- nextModelConstruction <- modelConstructer staticOptions nextModelConstructionInputWithInfernalMode+ nextModelConstruction <- if alienType == "alien"+ then do+ let nextModelConstructionInputWithInfernalMode = nextModelConstructionInput {alignmentModeInfernal = True, selectedQueries = currentSelectedQueries}+ logMessage (iterationSummaryLog nextModelConstructionInputWithInfernalMode) (tempDirPath staticOptions)+ logVerboseMessage (verbositySwitch staticOptions) (show nextModelConstructionInputWithInfernalMode) (tempDirPath staticOptions)+ modelConstructer staticOptions nextModelConstructionInputWithInfernalMode+ else do+ let nextGenomeFastas = tail (genomeFastas modelConstruction)+ let nextScanModelConstructionInputWithInfernalMode = nextModelConstructionInput {alignmentModeInfernal = True, selectedQueries = currentSelectedQueries, genomeFastas = nextGenomeFastas}+ logMessage (iterationSummaryLog nextScanModelConstructionInputWithInfernalMode) (tempDirPath staticOptions)+ logVerboseMessage (verbositySwitch staticOptions) (show nextScanModelConstructionInputWithInfernalMode) (tempDirPath staticOptions)+ scanModelConstructer staticOptions nextScanModelConstructionInputWithInfernalMode return nextModelConstruction -alignmentConstructionWithoutCandidates :: Maybe Taxon -> Maybe Int -> StaticOptions -> ModelConstruction -> IO ModelConstruction-alignmentConstructionWithoutCandidates currentTaxonomicContext upperTaxLimit staticOptions modelConstruction = do+alignmentConstructionWithoutCandidates :: String -> Maybe Taxon -> Maybe Int -> StaticOptions -> ModelConstruction -> IO ModelConstruction+alignmentConstructionWithoutCandidates alienType currentTaxonomicContext upperTaxLimit staticOptions modelConstruction = do let currentIterationNumber = iterationNumber modelConstruction let iterationDirectory = tempDirPath staticOptions ++ show currentIterationNumber ++ "/" --Found no new candidates in this iteration, reusing previous modelconstruction with increased upperTaxonomyLimit- let nextModelConstructionInputWithThreshold = modelConstruction {iterationNumber = currentIterationNumber + 1,upperTaxonomyLimit = upperTaxLimit,taxonomicContext = currentTaxonomicContext} --copy model and alignment from last iteration in place if present let previousIterationCMPath = tempDirPath staticOptions ++ show (currentIterationNumber - 1) ++ "/model.cm" previousIterationCMexists <- doesFileExist previousIterationCMPath+ let nextGenomeFastas = tail (genomeFastas modelConstruction)+ let nextModelConstructionInputWithThreshold = modelConstruction {iterationNumber = currentIterationNumber + 1,upperTaxonomyLimit = upperTaxLimit,taxonomicContext = currentTaxonomicContext}+ let nextScanModelConstructionInputWithThreshold = modelConstruction {iterationNumber = currentIterationNumber + 1,upperTaxonomyLimit = upperTaxLimit,taxonomicContext = currentTaxonomicContext, genomeFastas = nextGenomeFastas} if previousIterationCMexists then do logVerboseMessage (verbositySwitch staticOptions) "Alignment construction no candidates - previous cm\n" (tempDirPath staticOptions)@@ -374,17 +403,32 @@ copyFile previousIterationFastaPath thisIterationFastaPath copyFile previousIterationAlignmentPath thisIterationAlignmentPath copyFile previousIterationCMPath thisIterationCMPath- logMessage (iterationSummaryLog nextModelConstructionInputWithThreshold) (tempDirPath staticOptions)- logVerboseMessage (verbositySwitch staticOptions) (show nextModelConstructionInputWithThreshold) (tempDirPath staticOptions) writeFile (iterationDirectory ++ "done") ""- modelConstructer staticOptions nextModelConstructionInputWithThreshold+ if alienType == "alien"+ then do+ logMessage (iterationSummaryLog nextModelConstructionInputWithThreshold) (tempDirPath staticOptions)+ logVerboseMessage (verbositySwitch staticOptions) (show nextModelConstructionInputWithThreshold) (tempDirPath staticOptions)+ modelConstructer staticOptions nextModelConstructionInputWithThreshold+ else do+ logMessage (iterationSummaryLog nextScanModelConstructionInputWithThreshold) (tempDirPath staticOptions)+ logVerboseMessage (verbositySwitch staticOptions) (show nextScanModelConstructionInputWithThreshold) (tempDirPath staticOptions)+ scanModelConstructer staticOptions nextScanModelConstructionInputWithThreshold else do logVerboseMessage (verbositySwitch staticOptions) "Alignment construction no candidates - no previous iteration cm\n" (tempDirPath staticOptions) logMessage (iterationSummaryLog nextModelConstructionInputWithThreshold) (tempDirPath staticOptions) logVerboseMessage (verbositySwitch staticOptions) (show nextModelConstructionInputWithThreshold) (tempDirPath staticOptions) ---- writeFile (iterationDirectory ++ "done") ""- modelConstructer staticOptions nextModelConstructionInputWithThreshold+ if alienType == "alien"+ then do+ logMessage (iterationSummaryLog nextModelConstructionInputWithThreshold) (tempDirPath staticOptions)+ logVerboseMessage (verbositySwitch staticOptions) (show nextModelConstructionInputWithThreshold) (tempDirPath staticOptions)+ modelConstructer staticOptions nextModelConstructionInputWithThreshold+ else do+ logMessage (iterationSummaryLog nextScanModelConstructionInputWithThreshold) (tempDirPath staticOptions)+ logVerboseMessage (verbositySwitch staticOptions) (show nextScanModelConstructionInputWithThreshold) (tempDirPath staticOptions)+ scanModelConstructer staticOptions nextScanModelConstructionInputWithThreshold + findTaxonomyStart :: Bool -> Int -> Maybe String -> String -> Fasta () () -> IO Int findTaxonomyStart offlineMode threads inputBlastDatabase temporaryDirectory querySequence = do let queryIndexString = "1"@@ -418,8 +462,8 @@ CE.evaluate rightBestTaxIdResult else error "Find taxonomy start: Could not find blast hits to use as a taxonomic starting point" -searchCandidates :: StaticOptions -> Maybe String -> Int -> Maybe Int -> Maybe Int -> Double -> [Fasta () ()] -> IO SearchResult-searchCandidates staticOptions finaliterationprefix iterationnumber upperTaxLimit lowerTaxLimit expectThreshold inputQuerySequences = do+searchCandidates :: StaticOptions -> Maybe String -> Int -> Maybe Int -> Maybe Int -> Double -> Maybe String -> [Fasta () ()] -> IO SearchResult+searchCandidates staticOptions finaliterationprefix iterationnumber upperTaxLimit lowerTaxLimit expectThreshold maybeGenomeFasta inputQuerySequences = do Control.Monad.when (null inputQuerySequences) $ error "searchCandidates: - head: empty list of query sequences" let queryLength = fromIntegral (B.length (_bioSequence (_fasta (head inputQuerySequences)))) let queryIndexString = "1"@@ -428,12 +472,13 @@ 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) inputQuerySequences (Just (hitNumberQuery ++ entrezTaxFilter ++ softmaskFilter ++ registrationInfo)) (Just (5400000000 :: Int))+ let blastQuery = if (isJust maybeGenomeFasta) then BlastHTTPQuery (Just "ncbi") (Just "blastn") maybeGenomeFasta inputQuerySequences (Just (hitNumberQuery ++ entrezTaxFilter ++ softmaskFilter ++ registrationInfo)) (Just (5400000000 :: Int)) else 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) let logFileDirectoryPath = tempDirPath staticOptions ++ show iterationnumber ++ "/" ++ fromMaybe "" finaliterationprefix ++ "log" logDirectoryPresent <- doesDirectoryExist logFileDirectoryPath Control.Monad.when (not logDirectoryPresent) $ createDirectory (logFileDirectoryPath)+ -- print "Searching" --- blastOutput <- if (offline staticOptions) then CE.catch (blast logFileDirectoryPath (cpuThreads staticOptions) upperTaxLimit lowerTaxLimit (Just expectThreshold) (blastSoftmaskingToggle staticOptions) blastQuery) (\e -> do let err = show (e :: CE.IOException)@@ -461,30 +506,9 @@ writeFile (logFileDirectoryPath ++ "/" ++ queryIndexString ++ "_3blastHitsFilteredByLength") (showlines blastHitsFilteredByLength) let blastHitsFilteredByCoverage = filterByCoverage blastHitsFilteredByLength queryLength (coverageFilterToggle staticOptions) writeFile (logFileDirectoryPath ++ "/" ++ queryIndexString ++ "_3ablastHitsFilteredByLength") (showlines blastHitsFilteredByCoverage)- --tag BlastHits with TaxId- --blastHitsWithTaxIdOutput <- retrieveBlastHitsTaxIdEntrez blastHitsFilteredByCoverage- let blastHitsWithTaxId = extractBlastHitsTaxId blastHitsFilteredByCoverage- --let uncheckedBlastHitsWithTaxIdList = map (Control.Arrow.second extractTaxIdFromEntrySummaries) blastHitsWithTaxIdOutput- --let checkedBlastHitsWithTaxId = filter (\(_,taxids) -> not (null taxids)) uncheckedBlastHitsWithTaxIdList- --todo checked blasthittaxidswithblasthits need to be merged as taxid blasthit pairs- --let blastHitsWithTaxId = zip (concatMap fst checkedBlastHitsWithTaxId) (concatMap snd checkedBlastHitsWithTaxId)- blastHitsWithParentTaxIdOutput <- retrieveParentTaxIdsEntrez blastHitsWithTaxId- --let blastHitsWithParentTaxId = concat blastHitsWithParentTaxIdOutput- -- filter by ParentTaxId (only one hit per TaxId)- let blastHitsFilteredByParentTaxIdWithParentTaxId = filterByParentTaxId blastHitsWithParentTaxIdOutput True- --let blastHitsFilteredByParentTaxId = map fst blastHitsFilteredByParentTaxIdWithParentTaxId- writeFile (logFileDirectoryPath ++ "/" ++ queryIndexString ++ "_4blastHitsFilteredByParentTaxId") (showlines blastHitsFilteredByParentTaxIdWithParentTaxId)- -- Filtering with TaxTree (only hits from the same subtree as besthit)- --let blastHitsWithTaxId = zip blastHitsFilteredByParentTaxId blastHittaxIdList- --let (_, filteredBlastResults) = filterByNeighborhoodTreeConditional alignndmentModeInfernalToggle upperTaxLimit blastHitsWithTaxId (inputTaxNodes staticOptions) (fromJust upperTaxLimit) (singleHitperTaxToggle staticOptions)- --writeFile (logFileDirectoryPath ++ "/" ++ queryIndexString ++ "_5filteredBlastResults") (showlines filteredBlastResults)- -- Coordinate generation- let nonEmptyfilteredBlastResults = filter (\(blasthit,_) -> not (null (J._hsps blasthit))) blastHitsFilteredByParentTaxIdWithParentTaxId- let requestedSequenceElements = map (getRequestedSequenceElement queryLength) nonEmptyfilteredBlastResults- writeFile (logFileDirectoryPath ++ "/" ++ queryIndexString ++ "_6requestedSequenceElements") (showlines requestedSequenceElements)- -- Retrieval of full sequences from entrez- --fullSequencesWithSimilars <- retrieveFullSequences requestedSequenceElements- fullSequencesWithSimilars <- retrieveFullSequences staticOptions requestedSequenceElements+ fullSequencesWithSimilars <- if (isJust maybeGenomeFasta)+ then (scanFiltering blastHitsFilteredByCoverage logFileDirectoryPath queryIndexString queryLength maybeGenomeFasta staticOptions)+ else (alienFiltering blastHitsFilteredByCoverage logFileDirectoryPath queryIndexString queryLength staticOptions) if null fullSequencesWithSimilars then do writeFile (logFileDirectoryPath ++ "/" ++ queryIndexString ++ "_10afullSequencesWithSimilars") "No sequences retrieved"@@ -503,6 +527,35 @@ CE.evaluate (SearchResult fullSequences (Just dbSize)) else CE.evaluate (SearchResult [] Nothing) ++alienFiltering :: DS.Seq J.Hit -> [Char] -> [Char] -> Int -> StaticOptions -> IO [(Fasta () (), Int, B.ByteString)] +alienFiltering blastHitsFilteredByCoverage logFileDirectoryPath queryIndexString queryLength staticOptions = do+ --tag BlastHits with TaxId+ --blastHitsWithTaxIdOutput <- retrieveBlastHitsTaxIdEntrez blastHitsFilteredByCoverage+ let blastHitsWithTaxId = extractBlastHitsTaxId blastHitsFilteredByCoverage+ --let uncheckedBlastHitsWithTaxIdList = map (Control.Arrow.second extractTaxIdFromEntrySummaries) blastHitsWithTaxIdOutput+ --let checkedBlastHitsWithTaxId = filter (\(_,taxids) -> not (null taxids)) uncheckedBlastHitsWithTaxIdList+ --todo checked blasthittaxidswithblasthits need to be merged as taxid blasthit pairs+ --let blastHitsWithTaxId = zip (concatMap fst checkedBlastHitsWithTaxId) (concatMap snd checkedBlastHitsWithTaxId) + blastHitsWithParentTaxIdOutput <- retrieveParentTaxIdsEntrez blastHitsWithTaxId+ --let blastHitsWithParentTaxId = concat blastHitsWithParentTaxIdOutput+ -- filter by ParentTaxId (only one hit per TaxId)+ let blastHitsFilteredByParentTaxIdWithParentTaxId = filterByParentTaxId blastHitsWithParentTaxIdOutput True+ --let blastHitsFilteredByParentTaxId = map fst blastHitsFilteredByParentTaxIdWithParentTaxId+ writeFile (logFileDirectoryPath ++ "/" ++ queryIndexString ++ "_4blastHitsFilteredByParentTaxId") (showlines blastHitsFilteredByParentTaxIdWithParentTaxId)+ -- Filtering with TaxTree (only hits from the same subtree as besthit)+ --let blastHitsWithTaxId = zip blastHitsFilteredByParentTaxId blastHittaxIdList+ --let (_, filteredBlastResults) = filterByNeighborhoodTreeConditional alignmentModeInfernalToggle upperTaxLimit blastHitsWithTaxId (inputTaxNodes staticOptions) (fromJust upperTaxLimit) (singleHitperTaxToggle staticOptions)+ --writeFile (logFileDirectoryPath ++ "/" ++ queryIndexString ++ "_5filteredBlastResults") (showlines filteredBlastResults)+ -- Coordinate generation+ let nonEmptyfilteredBlastResults = filter (\(blasthit,_) -> not (null (J._hsps blasthit))) blastHitsFilteredByParentTaxIdWithParentTaxId+ let requestedSequenceElements = map (getRequestedSequenceElement queryLength) nonEmptyfilteredBlastResults+ writeFile (logFileDirectoryPath ++ "/" ++ queryIndexString ++ "_6requestedSequenceElements") (showlines requestedSequenceElements)+ -- Retrieval of full sequences from entrez+ --fullSequencesWithSimilars <- retrieveFullSequences requestedSequenceElements+ fullSequencesWithSimilars <- retrieveFullSequences staticOptions requestedSequenceElements+ return fullSequencesWithSimilars+ -- |Computes size of blast db in Mb computeDataBaseSize :: Double -> Double -> Double -> Double computeDataBaseSize evalue bitscore querylength = ((evalue * 2 ** bitscore) / querylength)/10^(6 :: Integer)@@ -922,7 +975,7 @@ taxEntries = taxRecords modelconstruction ++ buildTaxRecords alignmentResults currentIterationNumber potMembers = potentialMembers modelconstruction ++ inputPotentialMembers currentAlignmentMode = toggleInfernalAlignmentModeTrue || alignmentModeInfernal modelconstruction- nextModelConstruction = ModelConstruction newIterationNumber (inputFasta modelconstruction) taxEntries upperTaxLimit inputTaxonomicContext (evalueThreshold modelconstruction) currentAlignmentMode inputSelectedQueries potMembers+ nextModelConstruction = ModelConstruction newIterationNumber (inputFasta modelconstruction) taxEntries upperTaxLimit inputTaxonomicContext (evalueThreshold modelconstruction) currentAlignmentMode inputSelectedQueries potMembers [] buildTaxRecords :: [(Fasta () (),Int,B.ByteString)] -> Int -> [TaxonomyRecord] buildTaxRecords alignmentResults currentIterationNumber = taxonomyRecords@@ -1011,50 +1064,49 @@ -- | Run external locarna command and read the output into the corresponding datatype systemlocarna :: String -> (String,String,String,String) -> IO ExitCode-systemlocarna options (inputFilePath1, inputFilePath2, clustalformatoutputFilePath, outputFilePath) = system ("locarna " ++ options ++ " --clustal=" ++ clustalformatoutputFilePath ++ " " ++ inputFilePath1 ++ " " ++ inputFilePath2 ++ " > " ++ outputFilePath)+systemlocarna options (inputFilePath1, inputFilePath2, clustalformatoutputFilePath, outputFilePath) = silence (system ("locarna " ++ options ++ " --clustal=" ++ clustalformatoutputFilePath ++ " " ++ inputFilePath1 ++ " " ++ inputFilePath2 ++ " > " ++ outputFilePath)) -- | Run external mlocarna command and read the output into the corresponding datatype, there is also a folder created at the location of the input fasta file systemMlocarna :: String -> (String,String) -> IO ExitCode-systemMlocarna options (inputFilePath, outputFilePath) = system ("mlocarna " ++ options ++ " " ++ inputFilePath ++ " > " ++ outputFilePath)+systemMlocarna options (inputFilePath, outputFilePath) = silence $ system ("mlocarna " ++ options ++ " " ++ inputFilePath ++ " > " ++ outputFilePath) -- | Run external mlocarna command and read the output into the corresponding datatype, there is also a folder created at the location of the input fasta file, the job is terminated after the timeout provided in seconds systemMlocarnaWithTimeout :: String -> String -> (String,String) -> IO ExitCode-systemMlocarnaWithTimeout timeout options (inputFilePath, outputFilePath) = system ("timeout " ++ timeout ++"s "++ "mlocarna " ++ options ++ " " ++ inputFilePath ++ " > " ++ outputFilePath)+systemMlocarnaWithTimeout timeout options (inputFilePath, outputFilePath) = silence $ system ("timeout " ++ timeout ++"s "++ "mlocarna " ++ options ++ " " ++ inputFilePath ++ " > " ++ outputFilePath) -- | Run external clustalo command and return the Exitcode systemClustalw2 :: String -> (String,String,String) -> IO ExitCode-systemClustalw2 options (inputFilePath, outputFilePath, summaryFilePath) = system ("clustalw2 " ++ options ++ "-INFILE=" ++ inputFilePath ++ " -OUTFILE=" ++ outputFilePath ++ ">" ++ summaryFilePath)+systemClustalw2 options (inputFilePath, outputFilePath, summaryFilePath) = silence $ system ("clustalw2 " ++ options ++ "-INFILE=" ++ inputFilePath ++ " -OUTFILE=" ++ outputFilePath ++ ">" ++ summaryFilePath) -- | Run external clustalo command and return the Exitcode systemClustalo :: String -> (String,String) -> IO ExitCode-systemClustalo options (inputFilePath, outputFilePath) = system ("clustalo " ++ options ++ "--infile=" ++ inputFilePath ++ " >" ++ outputFilePath)+systemClustalo options (inputFilePath, outputFilePath) = silence $ system ("clustalo " ++ options ++ "--infile=" ++ inputFilePath ++ " >" ++ outputFilePath) -- | Run external CMbuild command and read the output into the corresponding datatype systemCMbuild :: String -> String -> String -> String -> IO ExitCode-systemCMbuild options alignmentFilepath modelFilepath outputFilePath = system ("cmbuild " ++ options ++ " " ++ modelFilepath ++ " " ++ alignmentFilepath ++ " > " ++ outputFilePath)+systemCMbuild options alignmentFilepath modelFilepath outputFilePath = silence $ system ("cmbuild " ++ options ++ " " ++ modelFilepath ++ " " ++ alignmentFilepath ++ " > " ++ outputFilePath) -- | Run CMCompare and read the output into the corresponding datatype systemCMcompare :: String -> String -> String -> IO ExitCode-systemCMcompare model1path model2path outputFilePath = system ("CMCompare -q " ++ model1path ++ " " ++ model2path ++ " >" ++ outputFilePath)+systemCMcompare model1path model2path outputFilePath = silence $ system ("CMCompare -q " ++ model1path ++ " " ++ model2path ++ " >" ++ outputFilePath) -- | Run CMsearch systemCMsearch :: Int -> String -> String -> String -> String -> IO ExitCode-systemCMsearch cpus options covarianceModelPath sequenceFilePath outputPath = system ("cmsearch --notrunc --cpu " ++ show cpus ++ " " ++ options ++ " -g " ++ covarianceModelPath ++ " " ++ sequenceFilePath ++ "> " ++ outputPath)+systemCMsearch cpus options covarianceModelPath sequenceFilePath outputPath = silence $ system ("cmsearch --notrunc --cpu " ++ show cpus ++ " " ++ options ++ " -g " ++ covarianceModelPath ++ " " ++ sequenceFilePath ++ "> " ++ outputPath) -- | Run CMstat systemCMstat :: String -> String -> IO ExitCode-systemCMstat covarianceModelPath outputPath = system ("cmstat " ++ covarianceModelPath ++ " > " ++ outputPath)+systemCMstat covarianceModelPath outputPath = silence $ system ("cmstat " ++ covarianceModelPath ++ " > " ++ outputPath) -- | Run CMcalibrate and return exitcode systemCMcalibrate :: String -> Int -> String -> String -> IO ExitCode systemCMcalibrate mode cpus covarianceModelPath outputPath- | mode == "fast" = system ("cmcalibrate --beta 1E-4 --cpu " ++ show cpus ++ " " ++ covarianceModelPath ++ "> " ++ outputPath)- | otherwise = system ("cmcalibrate --cpu " ++ show cpus ++ " " ++ covarianceModelPath ++ "> " ++ outputPath)-+ | mode == "fast" = silence $ system ("cmcalibrate --beta 1E-4 --cpu " ++ show cpus ++ " " ++ covarianceModelPath ++ "> " ++ outputPath)+ | otherwise = silence $ system ("cmcalibrate --cpu " ++ show cpus ++ " " ++ covarianceModelPath ++ "> " ++ outputPath) -- | Run CMcalibrate and return exitcode systemCMalign :: String -> String -> String -> String -> IO ExitCode-systemCMalign options filePathCovarianceModel filePathSequence filePathAlignment = system ("cmalign " ++ options ++ " " ++ filePathCovarianceModel ++ " " ++ filePathSequence ++ "> " ++ filePathAlignment)+systemCMalign options filePathCovarianceModel filePathSequence filePathAlignment = silence $ system ("cmalign " ++ options ++ " " ++ filePathCovarianceModel ++ " " ++ filePathSequence ++ "> " ++ filePathAlignment) compareCM :: String -> String -> String -> IO (Either String Double) compareCM rfamCMPath resultCMpath outputDirectory = do@@ -1175,7 +1227,7 @@ retrieveFullSequenceBlastDb blastDb temporaryDirectoryPath (nucleotideId,seqStart,seqStop,strand,_,taxid,subject') = do let sequencePath = temporaryDirectoryPath ++ "/" ++ nucleotideId ++ ".fa" let cmd = "blastdbcmd -db " ++ blastDb ++ " -range " ++ (show seqStart) ++ "-" ++ (show seqStop) ++ " -strand " ++ (setBlastDbStrand strand) ++ " -entry " ++ nucleotideId ++ " -outfmt %f -target_only -out " ++ sequencePath- print cmd+ --print cmd system(cmd) retrievedSequence <- readFastaFile sequencePath if null retrievedSequence@@ -1431,14 +1483,15 @@ else return [] -- | Wrapper functions that ensures that only 20 queries are sent per request-retrieveParentTaxIdsEntrez :: [(J.Hit,Int)] -> IO [(J.Hit,Int)]-retrieveParentTaxIdsEntrez taxIdwithBlastHits = do+retrieveParentTaxIdsEntrez :: [(J.Hit,Maybe Int)] -> IO [(J.Hit,Int)]+retrieveParentTaxIdsEntrez maybeTaxIdwithBlastHits = do+ let taxIdwithBlastHits = map (\(a,b) -> (a,fromJust b)) (filter (\(_,b) -> isJust b) maybeTaxIdwithBlastHits) let splits = portionListElements taxIdwithBlastHits 20 taxIdsOutput <- mapM retrieveParentTaxIdEntrez splits return (concat taxIdsOutput) -- | Extract taxids from JSON2 blasthit-extractBlastHitsTaxId :: DS.Seq J.Hit -> [(J.Hit,Int)]+extractBlastHitsTaxId :: DS.Seq J.Hit -> [(J.Hit,Maybe Int)] extractBlastHitsTaxId blastHits = do map (\a -> (a,J._taxid . head . J._description $ a)) (Data.Foldable.toList blastHits) @@ -1970,45 +2023,71 @@ blast :: String -> Int -> Maybe Int -> Maybe Int -> Maybe Double -> Bool -> BlastHTTPQuery -> IO (Either String J.BlastJSON2) blast _tempDirPath threads upperTaxIdLimit lowerTaxIdLimit expectThreshold _blastSoftmaskingToggle blastHTTPQuery = do+ --print "blast1" ---+ let selectedBlastDatabase = fromMaybe "" (Biobase.BLAST.HTTP.database blastHTTPQuery)+ --print ("blastdb: " ++ selectedBlastDatabase) --- --buildTaxonomyContext let upperTaxIdLimitPath = if isJust upperTaxIdLimit then _tempDirPath ++ "/upper.txids" else "" let lowerTaxIdLimitPath = if isJust lowerTaxIdLimit then _tempDirPath ++ "/lower.txids" else ""- when (isJust upperTaxIdLimit) $ systemGetSpeciesTaxId (fromJust upperTaxIdLimit) upperTaxIdLimitPath- when (isJust lowerTaxIdLimit) $ systemGetSpeciesTaxId (fromJust lowerTaxIdLimit) lowerTaxIdLimitPath- let positiveSetTaxIdLimitPath = _tempDirPath ++ "/postitiveset.txids"- if isJust lowerTaxIdLimit && isJust upperTaxIdLimit- then do- upperTaxIdsFile <- readFile upperTaxIdLimitPath- let upperTaxIds = lines upperTaxIdsFile- lowerTaxIdsFile <- readFile lowerTaxIdLimitPath- let lowerTaxIds = lines lowerTaxIdsFile- let positiveSetTaxIds = upperTaxIds \\ lowerTaxIds- let positiveSetTaxIdsFile = unlines positiveSetTaxIds- writeFile positiveSetTaxIdLimitPath positiveSetTaxIdsFile- else return ()- --sequenceSearch let fastaFilePath = _tempDirPath ++ "/blastQuery.fa" let blastResultFilePath = _tempDirPath ++ "/blastResult.json2"- let selectedBlastDatabase = fromMaybe "" (Biobase.BLAST.HTTP.database blastHTTPQuery) writeFastaFile fastaFilePath (querySequences blastHTTPQuery)- systemBlast threads selectedBlastDatabase upperTaxIdLimitPath lowerTaxIdLimitPath positiveSetTaxIdLimitPath expectThreshold _blastSoftmaskingToggle fastaFilePath blastResultFilePath- blastCmdResult <- BBI.blastCmdJSON2FromFile blastResultFilePath- --if isLeft blastResult then print (fromLeft blastResult) else print ""- if isRight blastCmdResult+ let isFastaDb = T.isSuffixOf (T.pack ".fa") (T.pack selectedBlastDatabase)+ if isFastaDb then do- let blastCmdOutput = J._blastcmdoutput2 (fromRight blastCmdResult)- when ((length blastCmdOutput) > 1) $ print "Blast output list with multiple elements"- if (not (null blastCmdOutput))- then (return (Right (J.BlastJSON2 (head blastCmdOutput)):: Either String J.BlastJSON2))- else (return (Left "Empty BlastOutput List" :: Either String J.BlastJSON2))- else (return (Left (fromLeft blastCmdResult) :: Either String J.BlastJSON2))+ systemBlast threads selectedBlastDatabase "" "" "" expectThreshold _blastSoftmaskingToggle fastaFilePath blastResultFilePath+ blastCmdResult <- BBI.blastCmdJSON2FromFile blastResultFilePath+ --if isLeft blastResult then print (fromLeft blastResult) else print ""+ if isRight blastCmdResult+ then do+ let blastCmdOutput = J._blastcmdoutput2 (fromRight blastCmdResult)+ --when ((length blastCmdOutput) > 1) $ print "Blast output list with multiple elements"+ if (not (null blastCmdOutput))+ then (return (Right (J.BlastJSON2 (head blastCmdOutput)):: Either String J.BlastJSON2))+ else (return (Left "Empty BlastOutput List" :: Either String J.BlastJSON2))+ else (return (Left (fromLeft blastCmdResult) :: Either String J.BlastJSON2))+ else do+ when (isJust upperTaxIdLimit) $ systemGetSpeciesTaxId (fromJust upperTaxIdLimit) upperTaxIdLimitPath+ when (isJust lowerTaxIdLimit) $ systemGetSpeciesTaxId (fromJust lowerTaxIdLimit) lowerTaxIdLimitPath+ let positiveSetTaxIdLimitPath = _tempDirPath ++ "/postitiveset.txids"+ if isJust lowerTaxIdLimit && isJust upperTaxIdLimit+ then do+ upperTaxIdsFile <- readFile upperTaxIdLimitPath+ let upperTaxIds = lines upperTaxIdsFile+ lowerTaxIdsFile <- readFile lowerTaxIdLimitPath+ let lowerTaxIds = lines lowerTaxIdsFile+ let positiveSetTaxIds = upperTaxIds \\ lowerTaxIds+ let positiveSetTaxIdsFile = unlines positiveSetTaxIds+ writeFile positiveSetTaxIdLimitPath positiveSetTaxIdsFile+ else return ()+ --sequenceSearch+ systemBlast threads selectedBlastDatabase upperTaxIdLimitPath lowerTaxIdLimitPath positiveSetTaxIdLimitPath expectThreshold _blastSoftmaskingToggle fastaFilePath blastResultFilePath+ blastCmdResult <- BBI.blastCmdJSON2FromFile blastResultFilePath+ --if isLeft blastResult then print (fromLeft blastResult) else print ""+ if isRight blastCmdResult+ then do+ let blastCmdOutput = J._blastcmdoutput2 (fromRight blastCmdResult)+ --when ((length blastCmdOutput) > 1) $ print "Blast output list with multiple elements"+ if (not (null blastCmdOutput))+ then (return (Right (J.BlastJSON2 (head blastCmdOutput)):: Either String J.BlastJSON2))+ else (return (Left "Empty BlastOutput List" :: Either String J.BlastJSON2))+ else (return (Left (fromLeft blastCmdResult) :: Either String J.BlastJSON2)) + + -- | Run external blast command systemBlast :: Int -> String -> String -> String -> String -> Maybe Double -> Bool -> String -> String -> IO ExitCode systemBlast threads _blastDatabase upperTaxLimitPath lowerTaxLimitPath positiveSetTaxIdLimitPath _evalueThreshold _blastSoftmaskingToggle queryFilepath outputFilePath = do let cmd = ("blastn " ++ threadedOption ++ expectThresholdOption ++ taxonomyOption ++ " " ++ softmaskOption ++ dbOption ++ " -query " ++ queryFilepath ++ " -outfmt 15 -out " ++ outputFilePath)- putStrLn cmd- system cmd+ if T.isSuffixOf (T.pack ".fa") (T.pack _blastDatabase)+ then do+ let makedbcmd = ("makeblastdb -in " ++ _blastDatabase ++ " -input_type fasta -dbtype nucl -parse_seqids ")+ (_,_) <- capture (system makedbcmd)+ (_,exitCode2) <- capture (system cmd)+ return exitCode2+ else do+ (_,exitCode3) <- capture (system cmd)+ return exitCode3 where threadedOption = " -num_threads " ++ show threads expectThresholdOption = if isJust _evalueThreshold then " -evalue " ++ show (fromJust _evalueThreshold) else "" dbOption = if null _blastDatabase then "" else " -db " ++ _blastDatabase ++ " "@@ -2029,4 +2108,148 @@ return () +------------------------------------------ RNAlienScan ------------------------------------ +-- | RNAlienScan RNA family model construction - generates iteration number, seed alignment and model+scanModelConstructer :: StaticOptions -> ModelConstruction -> IO ModelConstruction+scanModelConstructer staticOptions modelConstruction = do+ logMessage ("Iteration: " ++ show (iterationNumber modelConstruction) ++ "\n") (tempDirPath staticOptions)+ iterationSummary modelConstruction staticOptions+ let currentIterationNumber = iterationNumber modelConstruction+ let foundSequenceNumber = length (concatMap sequenceRecords (taxRecords modelConstruction))+ --extract queries+ let queries = extractQueries foundSequenceNumber modelConstruction+ logVerboseMessage (verbositySwitch staticOptions) ("Queries:" ++ show queries ++ "\n") (tempDirPath staticOptions)+ let iterationDirectory = tempDirPath staticOptions ++ show currentIterationNumber ++ "/"+ let lastGenome = null (genomeFastas modelConstruction)+ --Terminate on search space exhaustion+ if (not lastGenome)+ then do+ createDirectory iterationDirectory+ let (upperTaxLimit,lowerTaxLimit) = setTaxonomicContextEntrez currentIterationNumber (taxonomicContext modelConstruction) (upperTaxonomyLimit modelConstruction)+ logVerboseMessage (verbositySwitch staticOptions) ("Upper taxonomy limit: " ++ show upperTaxLimit ++ "\n " ++ "Lower taxonomy limit: "++ show lowerTaxLimit ++ "\n") (tempDirPath staticOptions)+ --search queries+ let expectThreshold = setBlastExpectThreshold modelConstruction+ let currentGenomeFasta = head (genomeFastas modelConstruction)+ let genomeFastaPath = (iterationDirectory ++ "genome.fa")+ writeFastaFile genomeFastaPath [currentGenomeFasta]+ searchResults <- catchAll (searchCandidates staticOptions Nothing currentIterationNumber upperTaxLimit lowerTaxLimit expectThreshold (Just genomeFastaPath) queries)+ (\e -> do logWarning ("Warning: Search results iteration" ++ show (iterationNumber modelConstruction) ++ " - exception: " ++ show e) (tempDirPath staticOptions)+ return (SearchResult [] Nothing))+ if null (candidates searchResults)+ then alignmentConstructionWithoutCandidates "scan" Nothing Nothing staticOptions modelConstruction+ else alignmentConstructionWithCandidates "scan" Nothing Nothing searchResults staticOptions modelConstruction+ else do+ logMessage "Message: Modelconstruction complete: Out of queries or taxonomic tree exhausted\n" (tempDirPath staticOptions)+ scanModelConstructionResult staticOptions modelConstruction++scanFiltering :: DS.Seq J.Hit -> [Char] -> [Char] -> Int -> Maybe String -> StaticOptions -> IO [(Fasta () (), Int, B.ByteString)] +scanFiltering blastHitsFilteredByCoverage logFileDirectoryPath queryIndexString queryLength maybeGenomeFasta staticOptions = do+ let nonEmptyfilteredBlastResults = filter (\(blasthit) -> not (null (J._hsps blasthit))) (Data.Foldable.toList blastHitsFilteredByCoverage)+ let dummyTaxId = replicate (length nonEmptyfilteredBlastResults) 0+ let blastResultsDummyTax = zip nonEmptyfilteredBlastResults dummyTaxId+ genomeFasta <- readFastaFile (fromJust maybeGenomeFasta)+ let sequenceByteString = _bioSequence . _fasta $ (head genomeFasta)+ let requestedSequenceElements = map (getRequestedSequenceElement queryLength) blastResultsDummyTax+ writeFile (logFileDirectoryPath ++ "/" ++ queryIndexString ++ "_6requestedSequenceElements") (showlines requestedSequenceElements)+ -- Retrieval of full sequences from genome+ let fullSequencesWithSimilars = retrieveGenomeFullSequences sequenceByteString staticOptions requestedSequenceElements+ return fullSequencesWithSimilars++-- | Wrapper for retrieveFullSequence that rerequests incomplete return sequees+retrieveGenomeFullSequences :: B.ByteString -> StaticOptions -> [(String,Int,Int,String,T.Text,Int,B.ByteString)] -> [(Fasta () (),Int,B.ByteString)]+retrieveGenomeFullSequences genomeSequence _ requestedSequences = map (retrieveGenomeFullSequence genomeSequence) requestedSequences++retrieveGenomeFullSequence :: B.ByteString -> (String,Int,Int,String,T.Text,Int,B.ByteString) -> ((Fasta () ()),Int,B.ByteString)+retrieveGenomeFullSequence sequenceByteString (nucleotideId,seqStart,seqStop,strand,_,_,subject') = (justFasta,0,subject')+ where retrievedSequence = byteStringSlice seqStart len sequenceByteString+ bioSequence = if strand == "1" then (BioSequence retrievedSequence) else (BioSequence rcretrievedSequence)+ currentFastaHeader= SequenceIdentifier (B.pack (nucleotideId ++ "_" ++ show seqStart ++ "_" ++ show seqStop ++ "_" ++ strand))+ justFasta = Fasta currentFastaHeader bioSequence+ len = if strand == "1" then seqStop - seqStart else seqStart - seqStop+ rcretrievedSequence = B.reverse (B.map complement' retrievedSequence)+ +complement' :: Char -> Char+complement' c+ | c == 'G' = 'C'+ | c == 'C' = 'G'+ | c == 'A' = 'T'+ | c == 'T' = 'A'+ | otherwise = 'N'++byteStringSlice :: Int -> Int -> B.ByteString -> B.ByteString+byteStringSlice start len = B.take len . B.drop start+++scanModelConstructionResult :: StaticOptions -> ModelConstruction -> IO ModelConstruction+scanModelConstructionResult staticOptions modelConstruction = do+ let currentIterationNumber = iterationNumber modelConstruction+ let outputDirectory = tempDirPath staticOptions+ logMessage ("Global search iteration: " ++ show currentIterationNumber ++ "\n") outputDirectory+ iterationSummary modelConstruction staticOptions+ let foundSequenceNumber = length (concatMap sequenceRecords (taxRecords modelConstruction))+ --extract queries+ --let querySeqIds = selectedQueries modelConstruction ---+ let queries = extractQueries foundSequenceNumber modelConstruction ---+ --let alignedSequences' = map nucleotideSequence (concatMap sequenceRecords (taxRecords modelConstruction)) ---+ logVerboseMessage (verbositySwitch staticOptions) ("Queries:" ++ show queries ++ "\n") outputDirectory+ let iterationDirectory = outputDirectory ++ show currentIterationNumber ++ "/"+ createDirectory iterationDirectory+ let logFileDirectoryPath = iterationDirectory ++ "log"+ createDirectoryIfMissing False logFileDirectoryPath+ let expectThreshold = setBlastExpectThreshold modelConstruction+ let (upperTaxLimit,lowerTaxLimit) = (Just (0 :: Int), Nothing)+ let currentGenomeFasta = genomeFastas modelConstruction+ let genomeFastaPath = (iterationDirectory ++ "genome.fa")+ writeFastaFile genomeFastaPath currentGenomeFasta+ candidates1 <- catchAll (searchCandidates staticOptions Nothing currentIterationNumber upperTaxLimit lowerTaxLimit expectThreshold (Just genomeFastaPath) queries)+ (\e -> do logWarning ("Warning: Search results iteration" ++ show currentIterationNumber ++ " - exception: " ++ show e) outputDirectory+ return (SearchResult [] Nothing))+ let uniqueCandidates = filterDuplicates modelConstruction candidates1+ (alignmentResults,currentPotentialMembers1) <- catchAll (alignCandidates staticOptions modelConstruction "" uniqueCandidates)+ (\e -> do logWarning ("Warning: Alignment results iteration" ++ show currentIterationNumber ++ " - exception: " ++ show e) outputDirectory+ return ([],[]))+ let currentPotentialMembers = [SearchResult currentPotentialMembers1 (blastDatabaseSize uniqueCandidates)] + let preliminaryFastaPath = iterationDirectory ++ "model.fa"+ let preliminaryCMPath = iterationDirectory ++ "model.cm"+ let preliminaryAlignmentPath = iterationDirectory ++ "model.stockholm"+ let preliminaryCMLogPath = iterationDirectory ++ "model.cm.log"+ let nextModelConstructionInput = constructNext currentIterationNumber modelConstruction alignmentResults Nothing Nothing [] currentPotentialMembers (alignmentModeInfernal modelConstruction)+ if (null alignmentResults) && not (alignmentModeInfernal modelConstruction)+ then do+ logVerboseMessage (verbositySwitch staticOptions) "Alignment result initial mode\n" outputDirectory+ logMessage "Message: No sequences found that statisfy filters. Try to reconstruct model with less strict cutoff parameters." outputDirectory+ let alignedSequences = extractAlignedSequences (iterationNumber modelConstruction) modelConstruction+ let alignmentSequences = map snd (V.toList (V.concat [alignedSequences]))+ writeFastaFile 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 cmBuildOptions preliminaryAlignmentPath preliminaryCMPath cmBuildFilepath+ _ <- systemCMcalibrate "fast" (cpuThreads staticOptions) preliminaryCMPath preliminaryCMLogPath+ reevaluatePotentialMembers staticOptions nextModelConstructionInput+ else+ if (alignmentModeInfernal modelConstruction)+ then do+ logVerboseMessage (verbositySwitch staticOptions) ("Alignment construction with candidates - infernal mode\n") outputDirectory+ constructModel nextModelConstructionInput staticOptions+ writeFile (iterationDirectory ++ "done") ""+ logMessage (iterationSummaryLog nextModelConstructionInput) outputDirectory+ logVerboseMessage (verbositySwitch staticOptions) (show nextModelConstructionInput) outputDirectory+ resultModelConstruction <- reevaluatePotentialMembers staticOptions nextModelConstructionInput+ return resultModelConstruction+ else do+ logVerboseMessage (verbositySwitch staticOptions) ("Alignment construction with candidates - initial mode\n") outputDirectory+ constructModel nextModelConstructionInput staticOptions+ let nextModelConstructionInputInfernalMode = nextModelConstructionInput {alignmentModeInfernal = True}+ logMessage (iterationSummaryLog nextModelConstructionInputInfernalMode) outputDirectory+ logVerboseMessage (verbositySwitch staticOptions) (show nextModelConstructionInputInfernalMode) outputDirectory+ writeFile (iterationDirectory ++ "done") ""+ resultModelConstruction <- reevaluatePotentialMembers staticOptions nextModelConstructionInputInfernalMode+ return resultModelConstruction
Biobase/RNAlien/Types.hs view
@@ -37,11 +37,12 @@ evalueThreshold :: Double, alignmentModeInfernal :: Bool, selectedQueries :: [Fasta () ()],- potentialMembers :: [SearchResult]+ potentialMembers :: [SearchResult],+ genomeFastas :: [Fasta () ()] } instance Show ModelConstruction where- show (ModelConstruction _iterationNumber _inputFasta _taxRecords _upperTaxonomyLimit _taxonomicContext _evalueThreshold _alignmentModeInfernal _selectedQueries _potentialMembers) = a ++ b ++ c ++ d ++ e ++ g ++ h ++ i+ show (ModelConstruction _iterationNumber _inputFasta _taxRecords _upperTaxonomyLimit _taxonomicContext _evalueThreshold _alignmentModeInfernal _selectedQueries _potentialMembers _genomeFastas) = a ++ b ++ c ++ d ++ e ++ g ++ h ++ i ++ j where a = "Modelconstruction iteration: " ++ show _iterationNumber ++ "\n" -- b = "Input fasta:\n" ++ concatMap (prettyPrintFasta 80) _inputFasta -- L.unpack (fastaHeader _inputFasta) ++ "\n" ++ L.unpack (fastaSequence _inputFasta) ++ "\n" b = "Input fasta:\n" ++ concatMap (convertString . fastaToByteString 80) _inputFasta@@ -51,6 +52,7 @@ g = "Evalue cutoff: " ++ show _evalueThreshold ++ "\n" h = "Selected queries: \n" ++ concatMap show _selectedQueries i = "Potential Members: \n" ++ concatMap show _potentialMembers+ j = "Number of genomes for RNAlienEgg: " ++ show (length _genomeFastas) data TaxonomyRecord = TaxonomyRecord { recordTaxonomyId :: Int,
+ Biobase/RNAlienScan.hs view
@@ -0,0 +1,126 @@+{-# LANGUAGE RecordWildCards #-}+{-# LANGUAGE DeriveDataTypeable #-}++-- | Unsupervized construction of RNA family models+-- For more information on RNA family models consult <http://>+-- Usage example: RNAlien -i /path/input.fa -c 5 -o /outdir/+-- Usage example offline mode: RNAlien -i /path/input.fa -b /backup/blast/nt_v5 -o /outdir/ -c 5 -t 1396 -j+module Main where++import System.Console.CmdArgs+import System.Directory+import Biobase.RNAlien.Types+import Biobase.RNAlien.Library+import Data.Maybe+import Data.Either.Unwrap+import Data.Time+import qualified System.FilePath as FP+import Paths_RNAlien (version)+import Data.Version (showVersion)+--import Biobase.Fasta.Streaming+import Control.Monad++data Options = Options+ { inputFastaFilePath :: String,+ inputGenomesFastaFilePath :: String,+ outputPath :: String,+ inputnSCICutoff :: Maybe Double,+ inputEvalueCutoff :: Maybe Double,+ lengthFilter :: Bool,+ coverageFilter :: Bool,+ singleHitperTax :: Bool,+ blastSoftmasking :: Bool,+ inputQuerySelectionMethod :: String,+ inputQueryNumber :: Int,+ threads :: Int,+ sessionIdentificator :: Maybe String,+ performEvaluation :: Bool,+ checkSetup :: Bool+ } deriving (Show,Data,Typeable)++options :: Options+options = Options+ { inputFastaFilePath = def &= name "i" &= help "Path to input fasta file",+ inputGenomesFastaFilePath = def &= name "b" &= help "Path to input genome fasta files",+ outputPath = def &= name "o" &= help "Path to output directory. Default: current working directory",+ inputnSCICutoff = Just (1 :: Double) &= name "z" &= help "Only candidate sequences with a normalized structure conservation index (nSCI) higher than this value are accepted. Default: 1",+ inputEvalueCutoff = Just (0.001 :: Double) &= name "e" &= help "Evalue cutoff for cmsearch filtering. Default: 0.001",+ 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 = 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",+ sessionIdentificator = Nothing &= name "d" &= help "Optional session id that is used instead of automatically generated one.",+ performEvaluation = True &= name "x" &= help "Perform evaluation step. Default: True",+ checkSetup = False &= name "g" &= help "Just prints installed tool versions and performs connection check. Default: False"+ } &= summary ("RNAlienScan " ++ alienVersion) &= help "Florian Eggenhofer - 2019" &= verbosity++main :: IO ()+main = do+ Options{..} <- cmdArgs options+ verboseLevel <- getVerbosity+ let tools = if inputQuerySelectionMethod == "clustering" then ["clustalo","mlocarna","RNAfold","RNAalifold","cmcalibrate","cmstat","cmbuild","RNAz","RNAcode"] else ["mlocarna","RNAfold","RNAalifold","cmcalibrate","cmstat","cmbuild","RNAz","RNAcode"]+ -- Generate SessionID+ sessionId <- createSessionID sessionIdentificator+ timestamp <- getCurrentTime+ currentWorkDirectory <- getCurrentDirectory+ let selectedOutputPath = if null outputPath then currentWorkDirectory else outputPath+ let temporaryDirectoryPath = FP.addTrailingPathSeparator selectedOutputPath ++ sessionId ++ "/"+ createDirectoryIfMissing False temporaryDirectoryPath+ networkCheck <- checkNCBIConnection+ if checkSetup+ then do+ toolsCheck <- checkTools tools inputQuerySelectionMethod temporaryDirectoryPath+ let setupCheckPath = temporaryDirectoryPath ++ "setupCheck"+ let toolCheckResult = either id id toolsCheck+ let networkCheckResult = either id id networkCheck+ writeFile setupCheckPath (toolCheckResult ++ "\n" ++ networkCheckResult ++ "\n")+ else+ if isLeft networkCheck+ then do+ putStrLn ("Error - Could not contact NCBI server: " ++ fromLeft networkCheck ++ "\n")+ logMessage ("Error - Could not contact NCBI server: " ++ fromLeft networkCheck ++ "\n") temporaryDirectoryPath+ else do+ createDirectoryIfMissing False (temporaryDirectoryPath ++ "log")+ -- Create Log files+ writeFile (temporaryDirectoryPath ++ "Log") ("RNAlienEgg " ++ alienVersion ++ "\n")+ writeFile (temporaryDirectoryPath ++ "log/warnings") ("")+ logMessage ("Timestamp: " ++ (show timestamp) ++ "\n") temporaryDirectoryPath+ logMessage ("Temporary Directory: " ++ temporaryDirectoryPath ++ "\n") temporaryDirectoryPath+ fastaInput <- readFastaFile inputFastaFilePath+ inputGenomesFasta <- readFastaFile inputGenomesFastaFilePath+ if null fastaInput+ then do+ putStrLn "Error: Input fasta file is empty."+ logMessage "Error: Input fasta file is empty.\n" temporaryDirectoryPath+ else do+ let iterationNumber = 0+ toolsCheck <- checkTools tools inputQuerySelectionMethod temporaryDirectoryPath+ if isLeft toolsCheck+ then do+ putStrLn ("Error - Not all required tools could be found in $PATH: " ++ fromLeft toolsCheck ++ "\n")+ logMessage ("Error - Not all required tools could be found in $PATH: " ++ fromLeft toolsCheck ++ "\n") temporaryDirectoryPath+ else do+ when (null inputGenomesFasta) (error "Please provide input genomes with the cmd line parameter -s")+ logToolVersions inputQuerySelectionMethod temporaryDirectoryPath+ let reformatedFastaInput = map reformatFasta fastaInput+ let staticOptions = StaticOptions temporaryDirectoryPath sessionId (fromJust inputnSCICutoff) Nothing singleHitperTax inputQuerySelectionMethod inputQueryNumber lengthFilter coverageFilter blastSoftmasking threads Nothing Nothing (setVerbose verboseLevel) True+ let initialization = ModelConstruction iterationNumber reformatedFastaInput [] Nothing Nothing (fromJust inputEvalueCutoff) False [] [] inputGenomesFasta+ logMessage (show initialization) temporaryDirectoryPath+ modelConstructionResults <- scanModelConstructer staticOptions initialization+ let resultTaxonomyRecordsCSVTable = constructTaxonomyRecordsCSVTable modelConstructionResults+ writeFile (temporaryDirectoryPath ++ "result.csv") resultTaxonomyRecordsCSVTable+ if performEvaluation+ then do+ resultEvaluation <- evaluateConstructionResult staticOptions modelConstructionResults+ appendFile (temporaryDirectoryPath ++ "Log") resultEvaluation+ resultSummary modelConstructionResults staticOptions+ writeFile (temporaryDirectoryPath ++ "done") ""+ else do+ resultSummary modelConstructionResults staticOptions+ writeFile (temporaryDirectoryPath ++ "done") ""++alienVersion :: String+alienVersion = showVersion version
RNAlien.cabal view
@@ -1,5 +1,5 @@ name: RNAlien-version: 1.6.0+version: 1.7.0 synopsis: Unsupervized construction of RNA family models description: RNAlien is a tool for automatic construction of RNAfamily models from a single sequence. .@@ -51,12 +51,12 @@ source-repository this type: git- location: https://github.com/eggzilla/RNAlien/tree/1.6.0- tag: 1.6.0+ location: https://github.com/eggzilla/RNAlien/tree/1.7.0+ tag: 1.7.0 executable RNAlien Hs-Source-Dirs: ./Biobase/- main-is: RNAlien.hs+ main-is: RNAlien.hs ghc-options: -Wall default-language: Haskell2010 other-modules: Paths_RNAlien@@ -64,6 +64,16 @@ random, containers, RNAlien, time, either-unwrap, filepath, BiobaseFasta == 0.3.0.* +executable RNAlienScan+ Hs-Source-Dirs: ./Biobase/+ main-is: RNAlienScan.hs+ ghc-options: -Wall+ default-language: Haskell2010+ other-modules: Paths_RNAlien+ build-depends: base >=4.5 && <5, cmdargs, directory,+ random, containers, RNAlien, time, either-unwrap, filepath,+ BiobaseFasta == 0.3.0.*+ executable RNAlienStatistics Hs-Source-Dirs: ./Biobase/ main-is: RNAlienStatistics.hs@@ -93,8 +103,8 @@ ghc-options: -Wall -fno-warn-unused-do-bind default-language: Haskell2010 build-depends: base >=4.5 && <5, cmdargs, ViennaRNAParser>=1.3.2, process, directory,- parsec, random, bytestring, Taxonomy >= 1.0.3, either-unwrap, containers, ClustalParser>=1.2.1, vector, edit-distance, cassava, matrix, hierarchical-clustering, filepath, HTTP, http-conduit, hxt, network<=2.8.0.0, aeson, text, transformers, pureMD5, http-types, text-metrics,- BiobaseTypes == 0.2.0.*, BiobaseFasta == 0.3.0.* , BiobaseBlast == 0.3.0.*, BlastHTTP >= 1.4.0, BiobaseHTTP == 1.1.0+ parsec, random, bytestring, Taxonomy >= 1.0.3, either-unwrap, containers, ClustalParser>=1.2.1, vector, edit-distance, cassava, matrix, hierarchical-clustering, filepath, HTTP, http-conduit, hxt, network<=2.8.0.0, aeson<=1.4.2.0, text, transformers, pureMD5, http-types, text-metrics,+ BiobaseTypes == 0.2.0.*, BiobaseFasta == 0.3.0.* , BiobaseBlast == 0.3.1.*, BlastHTTP >= 1.4.1, BiobaseHTTP == 1.1.0, silently Exposed-Modules: Biobase.RNAlien.Types Biobase.RNAlien.Library Biobase.RNAlien.RNAcentralHTTP