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.0 +version: 1.1.1 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.0 - tag: 1.1.0 + location: https://github.com/eggzilla/RNAlien/tree/1.1.1 + tag: 1.1.1 executable RNAlien Hs-Source-Dirs: ./src/Bio/ @@ -77,7 +77,7 @@ Library Hs-Source-Dirs: ./src/ ghc-options: -Wall -O2 -fno-warn-unused-do-bind - build-depends: base >=4.5 && <5, cmdargs, ViennaRNAParser>=1.2.8, process, directory, blastxml>=0.3.2, biofasta, parsec, random, BlastHTTP, biocore, bytestring, Taxonomy >= 1.0.1, either-unwrap, containers, ClustalParser>=1.1.0, EntrezHTTP>=1.0.1, vector, edit-distance, cassava, matrix, hierarchical-clustering, filepath, HTTP, http-conduit, hxt, network, aeson, text, transformers >= 0.4 && <0.5, pureMD5 + build-depends: base >=4.5 && <5, cmdargs, ViennaRNAParser>=1.2.9, process, directory, blastxml>=0.3.2, biofasta, parsec, random, BlastHTTP, biocore, bytestring, Taxonomy >= 1.0.1, either-unwrap, containers, ClustalParser>=1.1.0, EntrezHTTP>=1.0.1, vector, edit-distance, cassava, matrix, hierarchical-clustering, filepath, HTTP, http-conduit, hxt, network, aeson, text, transformers >= 0.4 && <0.5, pureMD5 Exposed-Modules: Bio.RNAlienData Bio.RNAlienLibrary Bio.RNAcentralHTTP diff --git a/src/Bio/RNAlien.hs b/src/Bio/RNAlien.hs --- a/src/Bio/RNAlien.hs +++ b/src/Bio/RNAlien.hs @@ -26,6 +26,7 @@ lengthFilter :: Bool, coverageFilter :: Bool, singleHitperTax :: Bool, + blastSoftmasking :: Bool, inputQuerySelectionMethod :: String, inputQueryNumber :: Int, threads :: Int, @@ -44,12 +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", inputQuerySelectionMethod = "filtering" &= name "m" &= help "Method for selection of queries (filtering,clustering). Default: filtering", inputQueryNumber = (5 :: Int) &= name "n" &= help "Number of queries used for candidate search. Default: 5", threads = 1 &= name "c" &= help "Number of available cpu slots/cores. Default: 1", taxonomyRestriction = Nothing &= name "r" &= help "Restrict search space to taxonomic kingdom (bacteria,archea,eukaryia). Default: not set", sessionIdentificator = Nothing &= name "d" &= help "Optional session id that is used instead of automatically generated one." - } &= summary "RNAlien version 1.0.0" &= help "Florian Eggenhofer, Ivo L. Hofacker, Christian Höner zu Siederdissen - 2013 - 2016" &= verbosity + } &= summary "RNAlien version 1.1.1" &= help "Florian Eggenhofer, Ivo L. Hofacker, Christian Höner zu Siederdissen - 2013 - 2016" &= verbosity main :: IO () main = do @@ -62,7 +64,7 @@ createDirectoryIfMissing False temporaryDirectoryPath createDirectoryIfMissing False (temporaryDirectoryPath ++ "log") -- Create Log files - writeFile (temporaryDirectoryPath ++ "Log") ("RNAlien 1.0.0" ++ "\n") + writeFile (temporaryDirectoryPath ++ "Log") ("RNAlien 1.1.1" ++ "\n") writeFile (temporaryDirectoryPath ++ "log/warnings") ("") logMessage ("Timestamp: " ++ (show timestamp) ++ "\n") temporaryDirectoryPath logMessage ("Temporary Directory: " ++ temporaryDirectoryPath ++ "\n") temporaryDirectoryPath @@ -91,7 +93,7 @@ let inputSequence = reformatFasta (head inputFasta) initialTaxId <- setInitialTaxId inputBlastDatabase temporaryDirectoryPath inputTaxId inputSequence let checkedTaxonomyRestriction = checkTaxonomyRestriction taxonomyRestriction - let staticOptions = StaticOptions temporaryDirectoryPath sessionId (fromJust inputnSCICutoff) inputTaxId singleHitperTax inputQuerySelectionMethod inputQueryNumber lengthFilter coverageFilter threads inputBlastDatabase checkedTaxonomyRestriction (setVerbose verboseLevel) + let staticOptions = StaticOptions temporaryDirectoryPath sessionId (fromJust inputnSCICutoff) inputTaxId singleHitperTax inputQuerySelectionMethod inputQueryNumber lengthFilter coverageFilter blastSoftmasking threads inputBlastDatabase checkedTaxonomyRestriction (setVerbose verboseLevel) let initialization = ModelConstruction iterationNumber inputSequence [] initialTaxId Nothing (fromJust inputEvalueCutoff) False [] [] logMessage (show initialization) temporaryDirectoryPath modelConstructionResults <- modelConstructer staticOptions initialization diff --git a/src/Bio/RNAlienData.hs b/src/Bio/RNAlienData.hs --- a/src/Bio/RNAlienData.hs +++ b/src/Bio/RNAlienData.hs @@ -18,6 +18,7 @@ queryNumber :: Int, lengthFilterToggle :: Bool, coverageFilterToggle :: Bool, + blastSoftmaskingToggle :: Bool, cpuThreads :: Int, blastDatabase :: Maybe String, taxRestriction :: Maybe String, diff --git a/src/Bio/RNAlienLibrary.hs b/src/Bio/RNAlienLibrary.hs --- a/src/Bio/RNAlienLibrary.hs +++ b/src/Bio/RNAlienLibrary.hs @@ -404,7 +404,7 @@ 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 = "&FILTER=True&FILTER=m" + 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)) --appendFile "/scratch/egg/blasttest/queries" ("\nBlast query:\n" ++ show blastQuery ++ "\n") logVerboseMessage (verbositySwitch staticOptions) ("Sending blast query " ++ (show iterationnumber) ++ "\n") (tempDirPath staticOptions) diff --git a/src/Bio/RNAlienStatistics.hs b/src/Bio/RNAlienStatistics.hs --- a/src/Bio/RNAlienStatistics.hs +++ b/src/Bio/RNAlienStatistics.hs @@ -18,10 +18,12 @@ import qualified Data.List.Split as DS import Text.Printf import Bio.RNAzParser +import qualified Bio.RNAcodeParser as RC data Options = Options { alienCovarianceModelPath :: String, alienrnazPath :: String, + alienrnacodePath :: String, aliencmstatPath :: String, rfamCovarianceModelPath :: String, rfamFastaFilePath :: String, @@ -30,6 +32,7 @@ rfamModelId :: String, rfamThreshold :: Double, alienThreshold :: Double, + databaseSize :: Maybe Double, outputDirectoryPath :: String, benchmarkIndex :: Int, thresholdSelection :: String, @@ -41,6 +44,7 @@ options = Options { alienCovarianceModelPath = def &= name "i" &= help "Path to alienCovarianceModelPath", alienrnazPath = def &= name "z" &= help "Path to alienRNAzResult", + alienrnacodePath = def &= name "w" &= help "Path to alienRNAcodeResult", aliencmstatPath = def &= name "m" &= help "Path to aliencmstatResult", rfamCovarianceModelPath = def &= name "r" &= help "Path to rfamCovarianceModelPath", rfamFastaFilePath = def &= name "g" &= help "Path to rfamFastaFile", @@ -50,6 +54,7 @@ outputDirectoryPath = def &= name "o" &= help "Path to output directory", alienThreshold = 20 &= name "t" &= help "Bitscore threshold for RNAlien model hits on Rfam fasta, default 20", rfamThreshold = 20 &= name "x" &= help "Bitscore threshold for Rfam model hits on Alien fasta, default 20", + databaseSize = Nothing &= name "k" &= help "Cmsearch database size in mega bases. default not set", benchmarkIndex = 1 &= name "b" &= help "Index used to identify sRNA tagged RNA families", thresholdSelection = "bitscore" &= name "s" &= help "Selection method, (bitscore, evalue), default bitscore", linkScores = False &= name "l" &= help "Triggers computation of linkscores via CMCompare", @@ -57,10 +62,11 @@ } &= summary "RNAlienStatistics" &= help "Florian Eggenhofer - >2013" &= verbosity --cmSearchFasta threads rfamCovarianceModelPath outputDirectoryPath "Rfam" False genomesDirectoryPath -cmSearchFasta :: Int -> String -> Double -> Int -> String -> String -> String -> String -> IO [CMsearchHit] -cmSearchFasta benchmarkIndex thresholdSelection thresholdScore cpuThreads covarianceModelPath outputDirectory modelType fastapath = do +cmSearchFasta :: Int -> String -> Double -> Maybe Double -> Int -> String -> String -> String -> String -> IO [CMsearchHit] +cmSearchFasta benchmarkIndex thresholdSelection thresholdScore databaseSize cpuThreads covarianceModelPath outputDirectory modelType fastapath = do createDirectoryIfMissing False (outputDirectory ++ "/" ++ modelType) - _ <- systemCMsearch cpuThreads " -Z 1000 " covarianceModelPath fastapath (outputDirectory ++ "/" ++ modelType ++ "/" ++ (show benchmarkIndex) ++ ".cmsearch") + _ <- systemCMsearch cpuThreads (maybe "" (\dbs -> " -Z " ++ show dbs ++ " ") databaseSize) covarianceModelPath fastapath (outputDirectory ++ "/" ++ modelType ++ "/" ++ (show benchmarkIndex) ++ ".cmsearch") + --_ <- systemCMsearch cpuThreads " " covarianceModelPath fastapath (outputDirectory ++ "/" ++ modelType ++ "/" ++ (show benchmarkIndex) ++ ".cmsearch") result <- readCMSearch (outputDirectory ++ "/" ++ modelType ++ "/" ++ (show benchmarkIndex) ++ ".cmsearch") if (isLeft result) then do @@ -69,14 +75,15 @@ else do let rightResults = fromRight result let significantHits = filterCMsearchHits thresholdSelection thresholdScore rightResults - let organismUniquesignificantHits = nubBy cmSearchSameOrganism significantHits - return organismUniquesignificantHits + let uniquesignificantHits = nubBy cmSearchSameHit significantHits + return uniquesignificantHits --cmSearchFasta threads rfamCovarianceModelPath outputDirectoryPath "Rfam" False genomesDirectoryPath -cmSearchesFasta :: Int -> String -> Double -> Int -> String -> String -> String -> String -> IO [CMsearchHit] -cmSearchesFasta benchmarkIndex thresholdSelection thresholdScore cpuThreads covarianceModelPath outputDirectory modelType fastapath = do +cmSearchesFasta :: Int -> String -> Double -> Maybe Double -> Int -> String -> String -> String -> String -> IO [CMsearchHit] +cmSearchesFasta benchmarkIndex thresholdSelection thresholdScore databaseSize cpuThreads covarianceModelPath outputDirectory modelType fastapath = do createDirectoryIfMissing False (outputDirectory ++ "/" ++ modelType) - _ <- systemCMsearch cpuThreads " -Z 1000 " covarianceModelPath fastapath (outputDirectory ++ "/" ++ modelType ++ "/" ++ (show benchmarkIndex) ++ ".cmsearch") + _ <- systemCMsearch cpuThreads (maybe "" (\dbs -> " -Z " ++ show dbs ++ " ") databaseSize) covarianceModelPath fastapath (outputDirectory ++ "/" ++ modelType ++ "/" ++ (show benchmarkIndex) ++ ".cmsearch") + --_ <- systemCMsearch cpuThreads " " covarianceModelPath fastapath (outputDirectory ++ "/" ++ modelType ++ "/" ++ (show benchmarkIndex) ++ ".cmsearch") result <- readCMSearches (outputDirectory ++ "/" ++ modelType ++ "/" ++ (show benchmarkIndex) ++ ".cmsearch") if (isLeft result) then do @@ -85,8 +92,11 @@ else do let rightResults = fromRight result let significantHits = filterCMsearchHits thresholdSelection thresholdScore rightResults - let organismUniquesignificantHits = nubBy cmSearchSameOrganism significantHits - return organismUniquesignificantHits + --putStrLn ("significant Hits " ++ show (length significantHits)) + let uniquesignificantHits = nubBy cmSearchSameHit significantHits + --putStrLn ("unique significant Hits " ++ show (length uniquesignificantHits)) + --let organismUniquesignificantHits = nubBy cmSearchSameOrganism significantHits + return uniquesignificantHits filterCMsearchHits :: String -> Double -> CMsearch -> [CMsearchHit] filterCMsearchHits thresholdSelection thresholdScore cmSearchResult @@ -118,6 +128,14 @@ newSequenceHeader = L.pack ((L.unpack (unSL (seqheader inputSequence))) ++ "cmS_" ++ (show (hitStart hitScoreEntry)) ++ "_" ++ (show (hitEnd hitScoreEntry)) ++ "_" ++ (show (hitStrand hitScoreEntry))) subSequence = Seq (SeqLabel newSequenceHeader) (SeqData (L.pack sequenceSubstring)) Nothing +--With paralogs allowed +cmSearchSameHit :: CMsearchHit -> CMsearchHit -> Bool +cmSearchSameHit hitscore1 hitscore2 + | unpackedSeqHeader1 == unpackedSeqHeader2 = True + | otherwise = False + where unpackedSeqHeader1 = (L.unpack (hitSequenceHeader hitscore1)) + unpackedSeqHeader2 = (L.unpack (hitSequenceHeader hitscore2)) + cmSearchSameOrganism :: CMsearchHit -> CMsearchHit -> Bool cmSearchSameOrganism hitscore1 hitscore2 | hitOrganism1 == hitOrganism2 = True @@ -139,7 +157,8 @@ Options{..} <- cmdArgs options rfamModelExists <- doesFileExist rfamCovarianceModelPath verbose <- getVerbosity - rnazString <- rnazOutput verbose alienrnazPath + rnazString <- rnazOutput verbose alienrnazPath + rnacodeString <- rnaCodeOutput verbose alienrnacodePath cmStatString <- cmStatOutput verbose aliencmstatPath if rfamModelExists then do @@ -155,8 +174,8 @@ alienFastaEntries <- readFile (outputDirectoryPath ++ FP.takeFileName alienFastaFilePath ++ ".entries") let rfamFastaEntriesNumber = read rfamFastaEntries :: Int let alienFastaEntriesNumber = read alienFastaEntries :: Int - rfamonAlienResults <- cmSearchesFasta benchmarkIndex thresholdSelection rfamThreshold threads rfamCovarianceModelPath outputDirectoryPath "rfamOnAlien" alienFastaFilePath - alienonRfamResults <- cmSearchFasta benchmarkIndex thresholdSelection alienThreshold threads alienCovarianceModelPath outputDirectoryPath "alienOnRfam" rfamFastaFilePath + rfamonAlienResults <- cmSearchesFasta benchmarkIndex thresholdSelection rfamThreshold databaseSize threads rfamCovarianceModelPath outputDirectoryPath "rfamOnAlien" alienFastaFilePath + alienonRfamResults <- cmSearchFasta benchmarkIndex thresholdSelection alienThreshold databaseSize threads alienCovarianceModelPath outputDirectoryPath "alienOnRfam" rfamFastaFilePath let rfamonAlienResultsNumber = length rfamonAlienResults let alienonRfamResultsNumber = length alienonRfamResults let rfamonAlienRecovery = (fromIntegral rfamonAlienResultsNumber :: Double) / (fromIntegral alienFastaEntriesNumber :: Double) @@ -178,9 +197,10 @@ putStrLn ("RfamonAlienRecovery: " ++ show rfamonAlienRecovery) putStrLn ("AlienonRfamRecovery: " ++ show alienonRfamRecovery) print rnazString + print rnacodeString print cmStatString else do - putStrLn (show benchmarkIndex ++ "\t" ++ rfamModelName ++ "\t" ++ rfamModelId ++ "\t" ++ (either id show linkscore) ++ "\t" ++ (either id show rfamMaxLinkScore) ++ "\t" ++ (either id show alienMaxLinkscore) ++ "\t" ++ show rfamThreshold ++ "\t" ++ show alienThreshold ++ "\t" ++ show rfamFastaEntriesNumber ++ "\t" ++ show alienFastaEntriesNumber ++ "\t" ++ show rfamonAlienResultsNumber ++ "\t" ++ show alienonRfamResultsNumber ++ "\t" ++ printf "%.2f" rfamonAlienRecovery ++ "\t" ++ printf "%.2f" alienonRfamRecovery ++ "\t" ++ rnazString ++ "\t" ++ cmStatString) + putStrLn (show benchmarkIndex ++ "\t" ++ rfamModelName ++ "\t" ++ rfamModelId ++ "\t" ++ (either id show linkscore) ++ "\t" ++ (either id show rfamMaxLinkScore) ++ "\t" ++ (either id show alienMaxLinkscore) ++ "\t" ++ show rfamThreshold ++ "\t" ++ show alienThreshold ++ "\t" ++ show rfamFastaEntriesNumber ++ "\t" ++ show alienFastaEntriesNumber ++ "\t" ++ show rfamonAlienResultsNumber ++ "\t" ++ show alienonRfamResultsNumber ++ "\t" ++ printf "%.2f" rfamonAlienRecovery ++ "\t" ++ printf "%.2f" alienonRfamRecovery ++ "\t" ++ rnazString ++ "\t" ++ rnacodeString ++ "\t" ++ cmStatString) else do --compute linkscore alienMaxLinkscore <- if linkScores then compareCM alienCovarianceModelPath alienCovarianceModelPath outputDirectoryPath else return ( Left "-") @@ -206,7 +226,7 @@ print rnazString print cmStatString else do - putStrLn (show benchmarkIndex ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ (either id show alienMaxLinkscore) ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ show alienFastaEntriesNumber ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ rnazString ++ "\t" ++ cmStatString) + putStrLn (show benchmarkIndex ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ (either id show alienMaxLinkscore) ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ show alienFastaEntriesNumber ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ rnazString ++ "\t" ++ rnacodeString ++ "\t" ++ cmStatString) rnazOutput :: Verbosity -> String -> IO String rnazOutput verbose rnazPath = do @@ -274,4 +294,38 @@ let output = "-\t" ++ "-\t" ++ "-\t" ++ "-\t" ++ "-\t" ++ "-\t" ++ "-\t" ++ "-\t" ++ "-" return output - +rnaCodeOutput :: Verbosity -> String -> IO String +rnaCodeOutput verbose rnaCodePath = do + rnacodePresent <- doesFileExist rnaCodePath + if rnacodePresent + then do + inputRNACode <- RC.readRNAcodeTabular rnaCodePath + if isRight inputRNACode + then do + let rnaCode = fromRight inputRNACode + let lowestPvalue = minimum (map RC.pvalue (RC.rnacodeHits rnaCode)) + let rnaCodeClassification = if lowestPvalue < 0.05 then "PROTEIN" else "OTHER" + if (verbose == Loud) + then do + let output = "RNAcode lowest p-value: " ++ (show lowestPvalue) ++ "\nrnaCodeClassification: " ++ rnaCodeClassification + return output + else do + let output = (show lowestPvalue) ++ "\t" ++ rnaCodeClassification + return output + else do + if (verbose == Loud) + then do + let output = "RNAcode lowest p-value: " ++ "-" ++ "\nrnaCodeClassification: " ++ "-" + return output + else do + let output = "-\t" ++ "-" + --let output = show (fromLeft inputRNACode) + return output + else do + if (verbose == Loud) + then do + let output = "RNAcode lowest p-value: " ++ "-" ++ "\nrnaCodeClassification: " ++ "-" + return output + else do + let output = "-\t" ++ "-" + return output