diff --git a/RNAlien.cabal b/RNAlien.cabal --- a/RNAlien.cabal +++ b/RNAlien.cabal @@ -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 diff --git a/src/Bio/RNAlienLibrary.hs b/src/Bio/RNAlienLibrary.hs --- a/src/Bio/RNAlienLibrary.hs +++ b/src/Bio/RNAlienLibrary.hs @@ -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 == ':' = '-' diff --git a/src/Bio/SelectSequences.hs b/src/Bio/SelectSequences.hs --- a/src/Bio/SelectSequences.hs +++ b/src/Bio/SelectSequences.hs @@ -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))