RNAlien 1.2.8 → 1.2.9
raw patch · 3 files changed
+242/−72 lines, 3 filesdep +text-metricsdep ~ClustalParser
Dependencies added: text-metrics
Dependency ranges changed: ClustalParser
Files
- RNAlien.cabal +4/−4
- src/Bio/RNAlienLibrary.hs +213/−60
- src/Bio/SelectSequences.hs +25/−8
RNAlien.cabal view
@@ -1,5 +1,5 @@ name: RNAlien-version: 1.2.8+version: 1.2.9 synopsis: Unsupervized construction of RNA family models description: RNAlien is a tool for automatic construction of RNAfamily models from a single sequence. .@@ -50,8 +50,8 @@ source-repository this type: git- location: https://github.com/eggzilla/RNAlien/tree/1.2.8- tag: 1.2.8+ location: https://github.com/eggzilla/RNAlien/tree/1.2.9+ tag: 1.2.9 executable RNAlien Hs-Source-Dirs: ./src/Bio/@@ -86,7 +86,7 @@ Library Hs-Source-Dirs: ./src/ ghc-options: -Wall -O2 -fno-warn-unused-do-bind- build-depends: base >=4.5 && <5, cmdargs, ViennaRNAParser>=1.3.1, process, directory, blastxml>=0.3.2, biofasta, parsec, random, BlastHTTP>=1.2.1, biocore, bytestring, Taxonomy >= 1.0.2, either-unwrap, containers, ClustalParser>=1.1.0, EntrezHTTP>=1.0.3, vector, edit-distance, cassava, matrix, hierarchical-clustering, filepath, HTTP, http-conduit, hxt, network, aeson, text, transformers, pureMD5, http-types+ build-depends: base >=4.5 && <5, cmdargs, ViennaRNAParser>=1.3.1, process, directory, blastxml>=0.3.2, biofasta, parsec, random, BlastHTTP>=1.2.1, biocore, bytestring, Taxonomy >= 1.0.2, either-unwrap, containers, ClustalParser>=1.2.0, EntrezHTTP>=1.0.3, vector, edit-distance, cassava, matrix, hierarchical-clustering, filepath, HTTP, http-conduit, hxt, network, aeson, text, transformers, pureMD5, http-types, text-metrics Exposed-Modules: Bio.RNAlienData Bio.RNAlienLibrary Bio.RNAcentralHTTP
src/Bio/RNAlienLibrary.hs view
@@ -1,5 +1,5 @@ -- | This module contains functions for RNAlien-+{-# LANGUAGE RankNTypes #-} module Bio.RNAlienLibrary ( module Bio.RNAlienData, createSessionID,@@ -24,6 +24,7 @@ checkNCBIConnection, preprocessClustalForRNAz, preprocessClustalForRNAzExternal,+ preprocessClustalForRNAcodeExternal, rnaZEvalOutput, reformatFasta, checkTaxonomyRestriction,@@ -75,6 +76,8 @@ import qualified Data.Text.Lazy.Encoding as E import qualified Data.Text.Lazy as TL import qualified Data.Text.Lazy.IO as TIO+import Text.Printf+import qualified Data.Text.Metrics as TM -- | Initial RNA family model construction - generates iteration number, seed alignment and model modelConstructer :: StaticOptions -> ModelConstruction -> IO ModelConstruction@@ -689,7 +692,7 @@ mlocarnaAlignment <- readStructuralClustalAlignment locarnaFilepath logEither mlocarnaAlignment (tempDirPath staticOptions) let stockholAlignment = convertClustaltoStockholm (fromRight mlocarnaAlignment)- writeFile stockholmFilepath stockholAlignment+ TI.writeFile stockholmFilepath stockholAlignment _ <- systemCMbuild cmBuildOptions stockholmFilepath cmFilepath cmBuildFilepath _ <- systemCMcalibrate "fast" (cpuThreads staticOptions) cmFilepath cmCalibrateFilepath return cmFilepath@@ -816,12 +819,12 @@ result = headEntry:filterIdenticalSequences' filteredEntries identitycutoff filterIdenticalSequences' [] _ = [] --- | Filter alignment entries by similiarity -filterIdenticalAlignmentEntry :: [ClustalAlignmentEntry] -> Double -> [ClustalAlignmentEntry]-filterIdenticalAlignmentEntry (headEntry:rest) identitycutoff = result- where filteredEntries = filter (\x -> (stringIdentity (entryAlignedSequence headEntry) (entryAlignedSequence x)) < identitycutoff) rest- result = headEntry:filterIdenticalAlignmentEntry filteredEntries identitycutoff-filterIdenticalAlignmentEntry [] _ = []+---- | Filter alignment entries by similiarity +--filterIdenticalAlignmentEntry :: [ClustalAlignmentEntry] -> Double -> [ClustalAlignmentEntry]+--filterIdenticalAlignmentEntry (headEntry:rest) identitycutoff = result+-- where filteredEntries = filter (\x -> (stringIdentity (entryAlignedSequence headEntry) (entryAlignedSequence x)) < identitycutoff) rest+-- result = headEntry:filterIdenticalAlignmentEntry filteredEntries identitycutoff+--filterIdenticalAlignmentEntry [] _ = [] isUnSimilarSequence :: [Sequence] -> Double -> Sequence -> Bool isUnSimilarSequence collectedSequences identitycutoff checkSequence = any (\ x -> (sequenceIdentity checkSequence x) < identitycutoff) collectedSequences@@ -835,13 +838,26 @@ | null resultList = False | otherwise = True where resultList = concatMap matches (concatMap hits (results blastResult))- + -- | Compute identity of sequences+textIdentity :: T.Text -> T.Text -> Double+textIdentity text1 text2 = identityPercent+ where distance = TM.hamming text1 text2+ --Replication of RNAz select sequences requires only allowing substitutions+ --costs = ED.defaultEditCosts {ED.deletionCosts = ED.ConstantCost 100,ED.insertionCosts = ED.ConstantCost 100,ED.transpositionCosts = ED.ConstantCost 100}+ maximumDistance = maximum [T.length text1, T.length text2]+ distanceDouble = toInteger ( fromJust distance )+ identityPercent = 1 - (fromIntegral distanceDouble/fromIntegral maximumDistance)+ + +-- | Compute identity of sequences stringIdentity :: String -> String -> Double stringIdentity string1 string2 = identityPercent- where distance = ED.levenshteinDistance ED.defaultEditCosts string1 string2+ where distance = ED.levenshteinDistance costs string1 string2+ --Replication of RNAz select sequences requires only allowing substitutions+ costs = ED.defaultEditCosts {ED.deletionCosts = ED.ConstantCost 100,ED.insertionCosts = ED.ConstantCost 100,ED.transpositionCosts = ED.ConstantCost 100} maximumDistance = maximum [length string1,length string2]- identityPercent = 100 - ((fromIntegral distance/fromIntegral maximumDistance) * (read "100" ::Double))+ identityPercent = 1 - (fromIntegral distance/fromIntegral maximumDistance) -- | Compute identity of sequences sequenceIdentity :: Sequence -> Sequence -> Double@@ -1291,30 +1307,30 @@ bottom = "//" stockholmOutput = alnHeader ++ entrystring ++ structureString ++ bottom -convertClustaltoStockholm :: StructuralClustalAlignment -> String+convertClustaltoStockholm :: StructuralClustalAlignment -> T.Text convertClustaltoStockholm parsedMlocarnaAlignment = stockholmOutput- where alnHeader = "# STOCKHOLM 1.0\n\n"+ where alnHeader = T.pack "# STOCKHOLM 1.0\n\n" clustalAlignment = structuralAlignmentEntries parsedMlocarnaAlignment uniqueIds = nub (map entrySequenceIdentifier clustalAlignment) mergedEntries = map (mergeEntry clustalAlignment) uniqueIds- maxIdentifierLenght = maximum (map (length . entrySequenceIdentifier) clustalAlignment)+ maxIdentifierLenght = maximum (map (T.length . entrySequenceIdentifier) clustalAlignment) spacerLength' = maxIdentifierLenght + 2- stockholmEntries = concatMap (buildStockholmAlignmentEntries spacerLength') mergedEntries- structureString = "#=GC SS_cons" ++ replicate (spacerLength' - 12) ' ' ++ secondaryStructureTrack parsedMlocarnaAlignment ++ "\n"- bottom = "//"- stockholmOutput = alnHeader ++ stockholmEntries ++ structureString ++ bottom+ stockholmEntries = T.concat (map (buildStockholmAlignmentEntries spacerLength') mergedEntries)+ structureString = (T.pack "#=GC SS_cons") `T.append` T.replicate (spacerLength' - 12) (T.pack " ") `T.append` secondaryStructureTrack parsedMlocarnaAlignment `T.append` (T.pack "\n")+ bottom = T.pack "//"+ stockholmOutput = alnHeader `T.append` stockholmEntries `T.append` structureString `T.append` bottom -mergeEntry :: [ClustalAlignmentEntry] -> String -> ClustalAlignmentEntry+mergeEntry :: [ClustalAlignmentEntry] -> T.Text -> ClustalAlignmentEntry mergeEntry clustalAlignment uniqueId = mergedEntry where idEntries = filter (\entry -> entrySequenceIdentifier entry==uniqueId) clustalAlignment- mergedSeq = foldr ((++) . entryAlignedSequence) "" idEntries+ mergedSeq = foldr ((T.append) . entryAlignedSequence) (T.pack "") idEntries mergedEntry = ClustalAlignmentEntry uniqueId mergedSeq -buildStockholmAlignmentEntries :: Int -> ClustalAlignmentEntry -> String+buildStockholmAlignmentEntries :: Int -> ClustalAlignmentEntry -> T.Text buildStockholmAlignmentEntries inputSpacerLength entry = entrystring- where idLength = length (filter (/= '\n') (entrySequenceIdentifier entry))- spacer = replicate (inputSpacerLength - idLength) ' '- entrystring = entrySequenceIdentifier entry ++ spacer ++ entryAlignedSequence entry ++ "\n"+ where idLength = T.length (T.filter (/= '\n') (entrySequenceIdentifier entry))+ spacer = T.replicate (inputSpacerLength - idLength) (T.pack " ")+ entrystring = entrySequenceIdentifier entry `T.append` spacer `T.append` entryAlignedSequence entry `T.append` (T.pack "\n") retrieveTaxonomicContextEntrez :: Int -> IO (Maybe Taxon) retrieveTaxonomicContextEntrez inputTaxId = do@@ -1524,19 +1540,26 @@ then do let resultRNAz = tempDirPath staticOptions ++ "result.rnaz" let resultRNAcode = tempDirPath staticOptions ++ "result.rnacode"- rnazClustalpath <- preprocessClustalForRNAcodeExternal clustalFilepath reformatedClustalPath- if isRight rnazClustalpath+ let seqNumber = 6 :: Int+ let optimalIdentity = 80 :: Double+ let maximalIdentity = 99 :: Double+ let referenceSequence = True+ --rnazClustalpath <- preprocessClustalForRNAcodeExternal clustalFilepath reformatedClustalPath seqNumber (truncate optimalIdentity) (truncate maximalIdentity) referenceSequence+ preprocessingOutput <- preprocessClustalForRNAz clustalFilepath reformatedClustalPath seqNumber optimalIdentity maximalIdentity referenceSequence + if isRight preprocessingOutput then do- systemRNAz "-l" (fromRight rnazClustalpath) resultRNAz + let rightPreprocessingOutput = fromRight preprocessingOutput+ let rnazClustalpath = snd rightPreprocessingOutput+ systemRNAz "-l" rnazClustalpath resultRNAz inputRNAz <- readRNAz resultRNAz let rnaZString = rnaZEvalOutput inputRNAz- RC.systemRNAcode " -t " (fromRight rnazClustalpath) resultRNAcode+ RC.systemRNAcode " -t " rnazClustalpath resultRNAcode inputRNAcode <- RC.readRNAcodeTabular resultRNAcode let rnaCodeString = rnaCodeEvalOutput inputRNAcode return ("\nEvaluation of RNAlien result :\nCMstat statistics for result.cm\n" ++ cmstatString ++ "\nRNAz statistics for result alignment: " ++ rnaZString ++ "\nRNAcode output for result alignment:\n" ++ rnaCodeString ++ "\nSequences found by RNAlien with RNAcentral entry:\n" ++ rnaCentralEvaluationResult) else do- logWarning ("Running RNAz for result evalution encountered a problem:" ++ fromLeft rnazClustalpath) (tempDirPath staticOptions) - return ("\nEvaluation of RNAlien result :\nCMstat statistics for result.cm\n" ++ cmstatString ++ "\nRNAz statistics for result alignment: Running RNAz for result evalution encountered a problem\n" ++ fromLeft rnazClustalpath ++ "\n" ++ "Sequences found by RNAlien with RNAcentral entry:\n" ++ rnaCentralEvaluationResult)+ logWarning ("Running RNAz for result evalution encountered a problem:" ++ fromLeft preprocessingOutput) (tempDirPath staticOptions) + return ("\nEvaluation of RNAlien result :\nCMstat statistics for result.cm\n" ++ cmstatString ++ "\nRNAz statistics for result alignment: Running RNAz for result evalution encountered a problem\n" ++ fromLeft preprocessingOutput ++ "\n" ++ "Sequences found by RNAlien with RNAcentral entry:\n" ++ rnaCentralEvaluationResult) else do logWarning "Message: RNAlien could not find additional covariant sequences\n Could not run RNAz statistics. Could not run RNAz statistics with a single sequence.\n" (tempDirPath staticOptions) return ("\nEvaluation of RNAlien result :\nCMstat statistics for result.cm\n" ++ cmstatString ++ "\nRNAlien could not find additional covariant sequences. Could not run RNAz statistics with a single sequence.\n\nSequences found by RNAlien with RNAcentral entry:\n" ++ rnaCentralEvaluationResult)@@ -1568,20 +1591,27 @@ showRNACodeHits rnacodeHit = show (RC.hss rnacodeHit) ++ "\t" ++ show (RC.frame rnacodeHit) ++ "\t" ++ show (RC.hitLength rnacodeHit) ++ "\t"++ show (RC.from rnacodeHit) ++ "\t" ++ show (RC.to rnacodeHit) ++ "\t" ++ (RC.name rnacodeHit) ++ "\t" ++ show (RC.start rnacodeHit) ++ "\t" ++ show (RC.end rnacodeHit) ++ "\t" ++ show (RC.score rnacodeHit) ++ show (RC.pvalue rnacodeHit) ++ "\n" -- | Call for external preprocessClustalForRNAz-preprocessClustalForRNAzExternal :: String -> String -> IO (Either String String)-preprocessClustalForRNAzExternal clustalFilepath reformatedClustalPath = do+preprocessClustalForRNAzExternal :: String -> String -> Int -> Int -> Int -> Bool -> IO (Either String (String,String))+preprocessClustalForRNAzExternal clustalFilepath reformatedClustalPath seqenceNumber optimalIdentity maximalIdenity referenceSequence = do clustalText <- TI.readFile clustalFilepath --change clustal format for rnazSelectSeqs.pl let reformatedClustalText = T.map reformatAln clustalText TI.writeFile reformatedClustalPath reformatedClustalText --select representative entries from result.Clustal with select_sequences let selectedClustalpath = clustalFilepath ++ ".selected"- system ("rnazSelectSeqs.pl " ++ reformatedClustalPath ++ " >" ++ selectedClustalpath)- return (Right selectedClustalpath)+ let sequenceNumberOption = " -n " ++ show seqenceNumber ++ " "+ let optimalIdentityOption = " -i " ++ show optimalIdentity ++ " "+ let maximalIdentityOption = " --max-id=" ++ show maximalIdenity ++ " "+ let referenceSequenceOption = if referenceSequence then " " else " -x "+ let syscall = ("rnazSelectSeqs.pl " ++ reformatedClustalPath ++ " " ++ sequenceNumberOption ++ optimalIdentityOption ++ maximalIdentityOption ++ referenceSequenceOption ++ " >" ++ selectedClustalpath)+ --putStrLn syscall+ system syscall+ selectedClustalText <- readFile selectedClustalpath+ return (Right ([],selectedClustalText)) -- | Call for external preprocessClustalForRNAcode - RNAcode additionally to RNAz requirements does not accept pipe,underscore, doublepoint symbols-preprocessClustalForRNAcodeExternal :: String -> String -> IO (Either String String)-preprocessClustalForRNAcodeExternal clustalFilepath reformatedClustalPath = do+preprocessClustalForRNAcodeExternal :: String -> String -> Int -> Int -> Int -> Bool -> IO (Either String (String,String))+preprocessClustalForRNAcodeExternal clustalFilepath reformatedClustalPath seqenceNumber optimalIdentity maximalIdenity referenceSequence = do clustalText <- TI.readFile clustalFilepath --change clustal format for rnazSelectSeqs.pl let clustalTextLines = T.lines clustalText@@ -1591,40 +1621,163 @@ TI.writeFile reformatedClustalPath (headerClustalTextLines `T.append` (T.singleton '\n') `T.append` reformatedClustalText) --select representative entries from result.Clustal with select_sequences let selectedClustalpath = clustalFilepath ++ ".selected"- system ("rnazSelectSeqs.pl " ++ reformatedClustalPath ++ " >" ++ selectedClustalpath)- return (Right selectedClustalpath)+ let sequenceNumberOption = " -n " ++ show seqenceNumber ++ " "+ let optimalIdentityOption = " -i " ++ show optimalIdentity ++ " "+ let maximalIdentityOption = " --max-id=" ++ show maximalIdenity ++ " "+ let referenceSequenceOption = if referenceSequence then " " else " -x "+ let syscall = ("rnazSelectSeqs.pl " ++ reformatedClustalPath ++ " " ++ sequenceNumberOption ++ optimalIdentityOption ++ maximalIdentityOption ++ referenceSequenceOption ++ " >" ++ selectedClustalpath)+ --putStrLn syscall+ system syscall+ selectedClustalText <- readFile selectedClustalpath+ return (Right ([],selectedClustalText)) --- | RNAz can process 500 sequences at max. Using rnazSelectSeqs to isolate representative sample. rnazSelectSeqs only accepts - gap characters, alignment is reformatted accordingly.-preprocessClustalForRNAz :: String -> String -> IO (Either String String)-preprocessClustalForRNAz clustalFilepath reformatedClustalPath = do+preprocessClustalForRNAz :: String -> String -> Int -> Double -> Double -> Bool -> IO (Either String (String,String))+preprocessClustalForRNAz clustalFilepath _ seqenceNumber optimalIdentity maximalIdenity referenceSequence = do clustalText <- TI.readFile clustalFilepath- let clustalTextLines = T.lines clustalText- if length clustalTextLines > 500+ let clustalTextLines = T.lines clustalText + parsedClustalInput <- readClustalAlignment clustalFilepath+ let selectedClustalpath = clustalFilepath ++ ".selected" + if length clustalTextLines > 5 then do - --change clustal format for rnazSelectSeqs.pl- let reformatedClustalString = T.map reformatAln clustalText- TI.writeFile reformatedClustalPath reformatedClustalString- --select representative entries from result.Clustal with select_sequences- let selectedClustalpath = clustalFilepath ++ ".selected"- parsedClustalInput <- readClustalAlignment clustalFilepath if isRight parsedClustalInput then do- let filteredClustalInput = rnaZSelectSeqs (fromRight parsedClustalInput) 500 99+ let (idMatrix,filteredClustalInput) = rnaCodeSelectSeqs2 (fromRight parsedClustalInput) seqenceNumber optimalIdentity maximalIdenity referenceSequence writeFile selectedClustalpath (show filteredClustalInput)- return (Right selectedClustalpath)+ let formatedIdMatrix = show (fmap formatIdMatrix idMatrix)+ return (Right (formatedIdMatrix,selectedClustalpath)) else return (Left (show (fromLeft parsedClustalInput)))- else return (Right clustalFilepath)+ else do+ let clustalTextLines = T.lines clustalText+ let headerClustalTextLines = T.unlines (take 2 clustalTextLines)+ let headerlessClustalTextLines = T.unlines (drop 2 clustalTextLines)+ let reformatedClustalText = T.map reformatRNACodeAln headerlessClustalTextLines+ TI.writeFile selectedClustalpath (headerClustalTextLines `T.append` (T.singleton '\n') `T.append` reformatedClustalText)+ return (Right ([],clustalFilepath)) +formatIdMatrix :: Maybe (Int,Int,Double) -> String+formatIdMatrix (Just (_,_,c)) = printf "%.2f" c+formatIdMatrix _ = "-" --- Iteratively removes sequences with decreasing similarity until target number of alignment entries is reached.-rnaZSelectSeqs :: ClustalAlignment -> Int -> Double -> ClustalAlignment-rnaZSelectSeqs currentClustalAlignment targetEntries identityCutoff- | targetEntries < numberOfEntries = rnaZSelectSeqs filteredAlignment targetEntries (identityCutoff - 1)- | otherwise = currentClustalAlignment- where numberOfEntries = length (alignmentEntries currentClustalAlignment) - filteredEntries = filterIdenticalAlignmentEntry (alignmentEntries currentClustalAlignment) identityCutoff - filteredAlignment = ClustalAlignment filteredEntries (conservationTrack currentClustalAlignment) +-- | Sequence preselection for RNAz and RNAcode +rnaCodeSelectSeqs2 :: ClustalAlignment -> Int -> Double -> Double -> Bool -> (Matrix (Maybe (Int,Int,Double)),ClustalAlignment)+rnaCodeSelectSeqs2 currentClustalAlignment targetSeqNumber optimalIdentity maximalIdentity referenceSequence = (identityMatrix,newClustalAlignment)+ where entryVector = V.fromList (alignmentEntries currentClustalAlignment)+ entrySequences = V.map entryAlignedSequence entryVector+ entryReformatedSequences = V.map (T.map reformatRNACodeAln) entrySequences+ totalSeqNumber = (V.length entryVector)+ identityMatrix = computeSequenceIdentityMatrix entryReformatedSequences+ entryIdentityVector = V.map fromJust (V.filter isJust (getMatrixAsVector identityMatrix))+ entryIdentities = V.toList entryIdentityVector+ --Similarity filter - filter too similar sequences until alive seqs are less then minSeqs+ entriesToDiscard = preFilterIdentityMatrix maximalIdentity targetSeqNumber totalSeqNumber [] entryIdentities+ allEntries = [1..totalSeqNumber]+ prefilteredEntries = allEntries \\ entriesToDiscard+ --Optimize mean pairwise similarity (greedily) - remove worst sequence until desired number is reached+ costList = map (computeEntryCost optimalIdentity entryIdentityVector) prefilteredEntries+ sortedCostList = (sortBy compareEntryCost2 costList)+ sortedIndices = map fst sortedCostList + --selectedEntryIndices = [1] ++ map fst (take (targetSeqNumber -1) sortedCostList)+ selectedEntryIndices = selectEntryIndices referenceSequence targetSeqNumber sortedIndices + selectedEntries = map (\ind -> entryVector V.! (ind-1)) selectedEntryIndices+ selectedEntryHeader = map entrySequenceIdentifier selectedEntries+ reformatedSelectedEntryHeader = map (T.map reformatRNACodeId) selectedEntryHeader+ selectedEntrySequences = map (\ind -> entryReformatedSequences V.! (ind-1)) selectedEntryIndices+ --gapfreeEntrySequences = T.transpose (T.filter (\a -> not (T.all isGap a)) (T.transpose selectedEntrySequences))+ gapfreeEntrySequences = T.transpose (filter (\a -> not (T.all isGap a)) (T.transpose selectedEntrySequences)) + gapfreeEntries = map (\(a,b) -> ClustalAlignmentEntry a b)(zip reformatedSelectedEntryHeader gapfreeEntrySequences)+ emptyConservationTrack = setEmptyConservationTrack gapfreeEntries (conservationTrack currentClustalAlignment)+ newClustalAlignment = currentClustalAlignment {alignmentEntries = gapfreeEntries, conservationTrack = emptyConservationTrack}++selectEntryIndices :: Bool -> Int -> [Int] -> [Int]+selectEntryIndices referenceSequence targetSeqNumber sortedIndices+ | referenceSequence = if elem (1 :: Int) firstX then firstX else 1:firstXm1+ | otherwise = firstX+ where firstXm1 = take (targetSeqNumber - 1) sortedIndices+ firstX = take targetSeqNumber sortedIndices+ +setEmptyConservationTrack :: [ClustalAlignmentEntry] -> T.Text -> T.Text+setEmptyConservationTrack alnentries currentConservationTrack+ | null alnentries = currentConservationTrack + | otherwise = newConservationTrack + where trackLength = T.length (entryAlignedSequence (head alnentries))+ newConservationTrack = T.replicate (trackLength + 0) (T.pack " ")+ +isGap :: Char -> Bool+isGap a + | a == '-' = True+ | otherwise = False++computeEntryCost :: Double -> V.Vector (Int,Int,Double) -> Int -> (Int,Double)+computeEntryCost optimalIdentity allIdentities currentIndex = (currentIndex,entryCost)+ where entryCost = V.sum (V.map (computeCost optimalIdentity) entryIdentities)+ entryIdentities = getEntryIdentities currentIndex allIdentities++getEntryIdentities :: Int -> V.Vector (Int,Int,Double) -> V.Vector (Int,Int,Double)+getEntryIdentities currentIndex allIdentities = (V.filter (isIIdx currentIndex) allIdentities) V.++ (V.filter (isJIdx currentIndex) allIdentities)++isIIdx :: Int -> (Int,Int,Double) -> Bool+isIIdx currentIdx (i,_,_) = currentIdx == i+isJIdx :: Int -> (Int,Int,Double) -> Bool+isJIdx currentIdx (_,j,_) = currentIdx == j++computeCost :: Double -> (Int,Int,Double) -> Double+computeCost optimalIdentity (_,_,c) = (c - optimalIdentity) * (c - optimalIdentity)++compareEntryCost2 :: (Int, Double) -> (Int, Double) -> Ordering +compareEntryCost2 (_,costA) (_,costB) = compare costA costB++-- TODO change to vector+preFilterIdentityMatrix :: Double -> Int -> Int-> [Int] -> [(Int,Int,Double)] -> [Int]+preFilterIdentityMatrix identityCutoff minSeqNumber totalSeqNumber filteredIds entryIdentities+ | (totalSeqNumber - (length filteredIds)) <= minSeqNumber = []+ | identityCutoff == (100 :: Double) = []+ | Prelude.null entryIdentities = []+ | otherwise = entryresult ++ preFilterIdentityMatrix identityCutoff minSeqNumber totalSeqNumber (filteredIds ++ entryresult) (tail entryIdentities)+ where currentEntry = head entryIdentities+ entryresult = checkIdentityEntry identityCutoff filteredIds currentEntry++checkIdentityEntry :: Double -> [Int] -> (Int,Int,Double) -> [Int]+checkIdentityEntry identityCutoff filteredIds (i,j,ident)+ | elem i filteredIds = []+ | elem j filteredIds = []+ | ident > identityCutoff = [j]+ | otherwise = []+ +computeSequenceIdentityMatrix :: V.Vector T.Text -> Matrix (Maybe (Int,Int,Double))+computeSequenceIdentityMatrix entryVector = matrix (V.length entryVector) (V.length entryVector) (computeSequenceIdentityEntry entryVector)++-- Computes Sequence identity once for each pair and not vs itself+computeSequenceIdentityEntry :: V.Vector T.Text -> (Int,Int) -> Maybe (Int,Int,Double)+computeSequenceIdentityEntry entryVector (row,col)+ | i < j = Just $ (row,col,ident)+ | otherwise = Nothing+ where i=row-1+ j=col-1+ --gaps in both sequences need to be removed, because they count as match+ ientry = (entryVector V.! i)+ jentry = (entryVector V.! j)+ (gfi,gfj) = unzip (filter notDoubleGap (T.zip ientry jentry))+ gfitext = T.pack gfi+ gfjtext = T.pack gfj + --ident=stringIdentity gfi gfj+ ident=textIdentity gfitext gfjtext+ +notDoubleGap :: (Char,Char) -> Bool+notDoubleGap (a,b)+ | a == '-' && b == '-' = False+ | otherwise = True++reformatRNACodeId :: Char -> Char +reformatRNACodeId c+ | c == ':' = '-'+ | c == '|' = '-'+ | c == '.' = '-'+ | c == '~' = '-'+ | c == '_' = '-'+ | c == '/' = '-' + | otherwise = c+ reformatRNACodeAln :: Char -> Char reformatRNACodeAln c | c == ':' = '-'
src/Bio/SelectSequences.hs view
@@ -11,14 +11,24 @@ data Options = Options { inputClustalPath :: String,- toogleExternalSelectSequences :: Bool+ toogleExternalSelectSequences :: Bool,+ seqenceNumber :: Int,+ optimalIdentity :: Double,+ maximalIdenity :: Double,+ referenceSequence :: Bool,+ distanceMatrixPath :: String } deriving (Show,Data,Typeable) options :: Options options = Options- { inputClustalPath = def &= name "i" &= help "Path to input clustal file",- toogleExternalSelectSequences = False &= name "e" &= help "Use only replacement of alignment characters and external 'selectSequence.pl'. Default: False" - } &= summary "SelectSequences" &= help "Florian Eggenhofer 2015" &= verbosity + { inputClustalPath = def &= name "c" &= help "Path to input clustal file",+ toogleExternalSelectSequences = False &= name "e" &= help "Use only replacement of alignment characters and external 'selectSequence.pl'. Default: False",+ seqenceNumber = (6 :: Int) &= name "n" &= help "Number of sequences in the output alignment. (Default: 6)",+ optimalIdentity = (80 :: Double) &= name "i" &= help "Optimize for this percentage of mean pairwise identity (Default: 80)",+ maximalIdenity = (95 :: Double) &= name "m" &= help "Sequences with a higher percentage of pairwise Identity will be removed. (Default: 95)",+ referenceSequence = True &= name "x" &= help "The first sequence (=reference sequence) is always present in the output alignment per default. Default: True",+ distanceMatrixPath = "" &= name "d" &= help "Path to distance matrix output file, only internal for interal sequence selection, e.g. /home/user/distmat (Default: )" + } &= summary "SelectSequences" &= help "Florian Eggenhofer 2016" &= verbosity main :: IO () main = do@@ -26,12 +36,19 @@ let reformatedClustalPath = inputClustalPath ++ ".reformated" if toogleExternalSelectSequences then do- resultStatus <- preprocessClustalForRNAzExternal inputClustalPath reformatedClustalPath+ resultStatus <- preprocessClustalForRNAzExternal inputClustalPath reformatedClustalPath seqenceNumber (truncate optimalIdentity) (truncate maximalIdenity) referenceSequence if (isRight resultStatus)- then (return ())+ then do+ let (idMatrix,resultAln) = fromRight resultStatus+ putStr resultAln+ if null distanceMatrixPath+ then return ()+ else (writeFile distanceMatrixPath idMatrix) else (print ("A problem occured selecting sequences: " ++ fromLeft resultStatus)) else do- resultStatus <- preprocessClustalForRNAz inputClustalPath reformatedClustalPath+ resultStatus <- preprocessClustalForRNAz inputClustalPath reformatedClustalPath seqenceNumber optimalIdentity maximalIdenity referenceSequence if (isRight resultStatus)- then (return ())+ then do+ let (_,resultAln) = fromRight resultStatus+ putStr resultAln else (print ("A problem occured selecting sequences: " ++ fromLeft resultStatus))