diff --git a/app/Main.hs b/app/Main.hs --- a/app/Main.hs +++ b/app/Main.hs @@ -39,6 +39,7 @@ , legacyFlag :: Bool , clipFastaFlag :: Bool , convertToAminoAcidsFlag :: Bool + , inputCodonTable :: CodonTable , inputFillIn :: FillInValue , inputStart :: Maybe Int , inputStop :: Maybe Int @@ -103,6 +104,15 @@ <> help "Whether to convert the filtered sequences to amino acids\ \ in the output. Applied last, even after add length." ) <*> option auto + ( long "codon-table" + <> metavar "[(CODON, AA)]" + <> value [] + <> help "The codon table used for converting to amino acids. This\ + \ is a list of codon AA pairs, so for example:\ + \ [(\"AAT\", '@'), (\"ATC\", '#')] for AATATCGGG would result\ + \ in @#G as defaults are used if the codon is not in\ + \ the custom table." ) + <*> option auto ( long "fill-in" <> short 'F' <> metavar "(FIELD, START, 'CHARACTER')" @@ -166,11 +176,13 @@ ( long "remove-highly-mutated" <> short 'h' <> help "Whether to remove highly mutated clone sequences (a third\ - \ of their sequence are different amino acids)." ) + \ of their sequence are different amino acids). Optionally\ + \ uses codon-table." ) <*> switch ( long "remove-stops" <> short 's' - <> help "Whether to remove sequences with stop codons" ) + <> help "Whether to remove sequences with stop codons (\"*\").\ + \ Optionally uses codon-table." ) <*> switch ( long "legacy-remove-duplicates" <> short 'd' @@ -344,7 +356,8 @@ x -- Remove clones with stops in the range - noStops x = not (removeStopsFlag opts) || hasNoStops genUnit stopRange x + noStops x = not (removeStopsFlag opts) + || hasNoStops genUnit (inputCodonTable opts) stopRange x -- Remove Ns from fasta list noNs = if removeTheNsFlag opts && (not . isAminoAcid $ genUnit) @@ -384,7 +397,7 @@ -- Convert to amino acids ntToaa = if convertToAminoAcidsFlag opts - then convertToAminoAcidsFastaSequence + then convertToAminoAcidsFastaSequence (inputCodonTable opts) else id -- Include sequence length in header at the end @@ -396,7 +409,9 @@ -- Remove highly mutated sequences removeHighMutations = if removeHighlyMutatedFlag opts - then filterHighlyMutatedEntry genUnit + then filterHighlyMutatedEntry + genUnit + (inputCodonTable opts) else id -- Extract mutations to a certain degree @@ -474,8 +489,9 @@ >> yield (T.pack "\n") ) -- want that newline at the end >-> PT.toHandle hOut - -- Finish up by closing file if written - unless (null . output $ opts) (IO.hClose hOut) + -- Finish up by closing handles just in case + IO.hClose hIn + IO.hClose hOut -- Legacy function modifyFastaCloneMap :: Options -> IO () @@ -528,6 +544,7 @@ (cloneMapNoStops, errorString) = if removeStopsFlag opts then removeStopsCloneMap genUnit + (inputCodonTable opts) stopRange cloneMapCustom else (cloneMapCustom, Nothing) @@ -549,6 +566,7 @@ $ genUnit ) then filterHighlyMutated genUnit + (inputCodonTable opts) cloneMapNoDuplicates else ( cloneMapNoDuplicates , Nothing ) @@ -572,6 +590,7 @@ (cloneMapAA, errorString3) = if (convertToAminoAcidsFlag opts) && (not . isAminoAcid $ genUnit) then convertToAminoAcidsCloneMap + (inputCodonTable opts) cloneMapNoEmptyClones else (cloneMapNoEmptyClones, Nothing) diff --git a/modify-fasta.cabal b/modify-fasta.cabal --- a/modify-fasta.cabal +++ b/modify-fasta.cabal @@ -2,7 +2,7 @@ -- documentation, see http://haskell.org/cabal/users-guide/ name: modify-fasta -version: 0.8.0.4 +version: 0.8.1.1 synopsis: Modify fasta (and CLIP) files in several optional ways -- description: homepage: https://github.com/GregorySchwartz/modify-fasta diff --git a/src/FilterCloneList.hs b/src/FilterCloneList.hs --- a/src/FilterCloneList.hs +++ b/src/FilterCloneList.hs @@ -24,8 +24,11 @@ -- | Remove highly mutated sequences (sequences with more than a third of -- their sequence being mutated). -filterHighlyMutatedEntry :: GeneticUnit -> CloneEntry -> CloneEntry -filterHighlyMutatedEntry !genUnit = newEntry +filterHighlyMutatedEntry :: GeneticUnit + -> CodonTable + -> CloneEntry + -> CloneEntry +filterHighlyMutatedEntry !genUnit !table = newEntry where newEntry (!germline, !fseqs) = ( germline , map snd @@ -68,4 +71,4 @@ | otherwise = False mutation x y = zip [1..] . T.zip x $ y readSeq Nucleotide x = Right x - readSeq AminoAcid x = translate 1 x + readSeq AminoAcid x = customTranslate table 1 x diff --git a/src/FilterCloneMap.hs b/src/FilterCloneMap.hs --- a/src/FilterCloneMap.hs +++ b/src/FilterCloneMap.hs @@ -37,8 +37,11 @@ -- | Remove highly mutated sequences (sequences with more than a third of -- their sequence being mutated). -filterHighlyMutated :: GeneticUnit -> CloneMap -> (CloneMap, Maybe String) -filterHighlyMutated !genUnit !cloneMap = (newCloneMap, errorString) +filterHighlyMutated :: GeneticUnit + -> CodonTable + -> CloneMap + -> (CloneMap, Maybe String) +filterHighlyMutated !genUnit !table !cloneMap = (newCloneMap, errorString) where newCloneMap = M.map (map snd . filter (not . fst) . rights) errorCloneMap @@ -84,7 +87,7 @@ | otherwise = False mutation x y = zip [1..] . T.zip x $ y readSeq Nucleotide x = Right x - readSeq AminoAcid x = translate 1 x + readSeq AminoAcid x = customTranslate table 1 x -- | Replace codons that have more than CodonMut mutations (make them "---" -- codons). @@ -115,11 +118,12 @@ -- | Remove clone sequences that have stop codons in the first stopRange -- codons removeStopsCloneMap :: GeneticUnit + -> CodonTable -> Int -> CloneMap -> (CloneMap, Maybe String) -removeStopsCloneMap !genUnit !stopRange !cloneMap = ( newCloneMap - , errorString ) +removeStopsCloneMap !genUnit !table !stopRange !cloneMap = ( newCloneMap + , errorString ) where errorString = listToMaybe' . unlines @@ -129,17 +133,17 @@ . M.map ( intercalate "\n" . map T.unpack . lefts - . map (translate 1) + . map (customTranslate table 1) ) $ cloneMap newCloneMap = M.map (filter (filterStops genUnit)) cloneMap - filterStops Nucleotide x = (isRight' . translate 1 $ x) + filterStops Nucleotide x = (isRight' . customTranslate table 1 $ x) && ( not . T.isInfixOf "*" . T.take stopRange . fastaSeq . fromEither - . translate 1 ) x + . customTranslate table 1 ) x filterStops AminoAcid x = not . T.isInfixOf "*" . T.take stopRange @@ -215,8 +219,10 @@ removeEmptyClone = M.filter (not . null) -- | Convert sequences to amino acids -convertToAminoAcidsCloneMap :: CloneMap -> (CloneMap, Maybe String) -convertToAminoAcidsCloneMap cloneMap = (newCloneMap, errorString) +convertToAminoAcidsCloneMap :: CodonTable + -> CloneMap + -> (CloneMap, Maybe String) +convertToAminoAcidsCloneMap table cloneMap = (newCloneMap, errorString) where newCloneMap = M.mapKeysWith (++) (\(!x, !y) -> (x, fromEither y)) . M.filterWithKey (\(_, !y) _ -> isRight' y) @@ -231,9 +237,9 @@ $ v ) $ errorCloneMap errorCloneMap = M.mapKeys keyMap - . M.map (map (translate 1)) + . M.map (map (customTranslate table 1)) $ cloneMap - keyMap (!x, !y) = (x, translate 1 y) + keyMap (!x, !y) = (x, customTranslate table 1 y) eitherToString (Right _) = "" eitherToString (Left x) = T.unpack x fromEither (Right x) = x diff --git a/src/FilterFastaList.hs b/src/FilterFastaList.hs --- a/src/FilterFastaList.hs +++ b/src/FilterFastaList.hs @@ -27,18 +27,19 @@ -- | Remove clone sequences that have stop codons in the first stopRange -- codons hasNoStops :: GeneticUnit + -> CodonTable -> Int -> FastaSequence -> Bool -hasNoStops genUnit stopRange = result . stop genUnit +hasNoStops genUnit table stopRange = result . stop genUnit where result (Right x) = x result (Left x) = error . T.unpack $ x - stop Nucleotide = fmap ( not - . T.isInfixOf "*" - . T.take stopRange - . fastaSeq ) - . translate 1 + stop Nucleotide = fmap ( not + . T.isInfixOf "*" + . T.take stopRange + . fastaSeq ) + . customTranslate table 1 stop AminoAcid = Right . not . T.isInfixOf "*" . T.take stopRange . fastaSeq -- | Remove out of frame sequences diff --git a/src/TransformFastaList.hs b/src/TransformFastaList.hs --- a/src/TransformFastaList.hs +++ b/src/TransformFastaList.hs @@ -32,8 +32,10 @@ import Utility -- | Convert sequences to amino acids -convertToAminoAcidsFastaSequence :: FastaSequence -> FastaSequence -convertToAminoAcidsFastaSequence = fromEither . translate 1 +convertToAminoAcidsFastaSequence :: CodonTable + -> FastaSequence + -> FastaSequence +convertToAminoAcidsFastaSequence table = fromEither . customTranslate table 1 where fromEither (Right x) = x fromEither (Left x) = error . T.unpack $ x diff --git a/src/Types.hs b/src/Types.hs --- a/src/Types.hs +++ b/src/Types.hs @@ -29,6 +29,7 @@ -- Advanced type CloneEntry = (Germline, [FastaSequence]) type FillInValue = (Field, Start, Char) +type CodonTable = [(T.Text, Char)] type Mutation = (Char, Char) type CountMap = Map.Map (Position, Mutation) Int type CodonMutations = [[(Position, Mutation)]]