RNAlien 1.1.0 → 1.1.1
raw patch · 5 files changed
+81/−24 lines, 5 filesdep ~ViennaRNAParser
Dependency ranges changed: ViennaRNAParser
Files
- RNAlien.cabal +4/−4
- src/Bio/RNAlien.hs +5/−3
- src/Bio/RNAlienData.hs +1/−0
- src/Bio/RNAlienLibrary.hs +1/−1
- src/Bio/RNAlienStatistics.hs +70/−16
RNAlien.cabal view
@@ -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
src/Bio/RNAlien.hs view
@@ -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
src/Bio/RNAlienData.hs view
@@ -18,6 +18,7 @@ queryNumber :: Int, lengthFilterToggle :: Bool, coverageFilterToggle :: Bool,+ blastSoftmaskingToggle :: Bool, cpuThreads :: Int, blastDatabase :: Maybe String, taxRestriction :: Maybe String,
src/Bio/RNAlienLibrary.hs view
@@ -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)
src/Bio/RNAlienStatistics.hs view
@@ -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