diff --git a/app/Main.hs b/app/Main.hs new file mode 100644 --- /dev/null +++ b/app/Main.hs @@ -0,0 +1,627 @@ +-- modify-fasta +-- By Gregory W. Schwartz + +-- Takes a fasta file filters the fasta file in several ways. + +{-# LANGUAGE BangPatterns #-} + +-- Built-in +import Data.Maybe +import qualified Data.Map as M +import qualified System.IO as IO +import qualified Data.Text.IO as IO +import Control.Monad + +-- Cabal +import qualified Data.Text as T +import qualified Data.Text.IO as T +import Options.Applicative +import Data.Fasta.Text +import Pipes +import qualified Pipes.Prelude as P +import qualified Pipes.Text as PT +import qualified Pipes.Text.IO as PT +import qualified Data.List.Split as Split + +-- Local +import Types +import Utility +import FilterCloneMap +import FilterFastaList +import FilterCloneList +import TransformFastaList +import TransformCloneList +import Print + +-- Command line arguments +data Options = Options { input :: String + , aminoAcidsFlag :: GeneticUnit + , legacyFlag :: Bool + , clipFastaFlag :: Bool + , convertToAminoAcidsFlag :: Bool + , inputFillIn :: FillInValue + , inputStart :: Maybe Int + , inputStop :: Maybe Int + , inputMinSequenceMutation :: Maybe Int + , inputMutationCount :: Maybe Int + , inputMutationPercent :: Maybe Double + , addLengthFlag :: Bool + , removeTheNsFlag :: Bool + , removeGermlinesPreFlag :: Bool + , removeHighlyMutatedFlag :: Bool + , removeStopsFlag :: Bool + , removeDuplicatesFlag :: Bool + , removeOutOfFrameFlag :: Bool + , removeUnknownNuc :: Bool + , inputInFrame :: Maybe Field + , inputOutFrame :: Maybe Field + , trimFrame :: Bool + , inputStopRange :: Int + , inputCodonMut :: CodonMut + , inputCodonMutType :: String + , inputMutType :: String + , inputChangeField :: String + , inputCustomFilter :: String + , customGermlineFlag :: Bool + , customRemoveFlag :: Bool + , inputGeneAlleleField :: Int + , countFlag :: Bool + , output :: String + } + +-- Command line options +options :: Parser Options +options = Options + <$> strOption + ( long "input" + <> short 'i' + <> metavar "FILE" + <> value "" + <> help "The input fasta file or CLIP fasta file" ) + <*> option auto + ( long "unit" + <> short 'u' + <> metavar "AminoAcid | Nucleotide" + <> help "Whether these sequences are composed of\ + \ amino acids (AminoAcid) or nucleotides (Nucleotide)" ) + <*> switch + ( long "legacy" + <> short 'L' + <> help "Whether to use the legacy version with no pipes. Note: The\ + \ legacy version supports more features but is greedy\ + \ in terms of speed and memory. Use only if really needed.\ + \ Features that are legacy only are noted in this\ + \ documentation" ) + <*> switch + ( long "clip-fasta" + <> short 'A' + <> help "Whether the input is a clip fasta file (has germline >>\ + \ sequences)." ) + <*> switch + ( long "convert-to-amino-acids" + <> short 'C' + <> help "Whether to convert the filtered sequences to amino acids\ + \ in the output. Applied last, even after add length." ) + <*> option auto + ( long "fill-in" + <> short 'F' + <> metavar "(FIELD, START, 'CHARACTER')" + <> value (-1, -1, 'X') + <> help "Use the FIELD index (1 indexed, split by '|') in the\ + \ header to fill in the unknown character CHARACTER with the\ + \ respective character in that field. Use the START value\ + \ to let the program know where the string in the field starts\ + \ from in the sequence. For instance, a header of >H2O|HELLO\ + \ for the sequence HIMANHELXODUDE and a value of fill-in of\ + \ (2, 6, 'X') would change the sequence to HIMANHELLODUDE. If\ + \ the string in the field has the unknown character as well,\ + \ i.e. >H2O|HELXO, then the sequence is considered bad and is\ + \ removed" ) + <*> optional ( option auto + ( long "start" + <> short 't' + <> metavar "[ ] | INT" + <> help "Remove everything before this position (1 indexed).\ + \ Done first just after filtering." ) ) + <*> optional ( option auto + ( long "stop" + <> short 'p' + <> metavar "[ ] | INT" + <> help "Remove everything after this position (1 indexed).\ + \ Done first just after filtering." ) ) + <*> optional ( option auto + ( long "frequent-mutation-min" + <> metavar "[ ] | INT" + <> help "Minimum number of sequences required for a clone to be valid\ + \ in the calculation of frequent mutations, replaces\ + \ with gaps otherwise" ) ) + <*> optional ( option auto + ( long "frequent-mutation-count" + <> metavar "[ ] | INT" + <> help "Only include codons containing a mutation present in this\ + \ many sequences in the clone or more. 0 is all sequences.\ + \ Converts the unincluded codons to gaps." ) ) + <*> optional ( option auto + ( long "frequent-mutation-percent" + <> metavar "[ ] | PERCENT" + <> help "Only include codons containing a mutation present in this\ + \ percentage of sequences in the clone or more.\ + \ 0 is all sequences. Converts the unincluded codons\ + \ to gaps." ) ) + <*> switch + ( long "add-length" + <> short 'l' + <> help "Whether to append the length of the sequence to the end of\ + \ the header, calculated after --convert-to-amino-acids\ + \ if enabled" ) + <*> switch + ( long "remove-N" + <> short 'N' + <> help "Whether to replace N or n in the sequence with a gap, '-'" ) + <*> switch + ( long "remove-germlines" + <> short 'g' + <> help "Whether to remove germlines." ) + <*> switch + ( long "remove-highly-mutated" + <> short 'h' + <> help "Whether to remove highly mutated clone sequences (a third\ + \ of their sequence are different amino acids)." ) + <*> switch + ( long "remove-stops" + <> short 's' + <> help "Whether to remove sequences with stop codons" ) + <*> switch + ( long "legacy-remove-duplicates" + <> short 'd' + <> help "Whether to remove duplicate sequences. LEGACY ONLY" ) + <*> switch + ( long "remove-out-of-frame" + <> short 'O' + <> help "Whether to remove sequences that are out of frame--if the\ + \ sequences or number of gaps is not divisible by 3" ) + <*> switch + ( long "remove-unknown-nucleotides" + <> short 'Y' + <> help "Convert unknown nucleotides (not ACGTN-.) to gaps (-)" ) + <*> optional ( option auto + ( long "input-inframe" + <> metavar "[ ] | FIELD" + <> help "Represents the 1 indexed field split by '|'\ + \ containing the inframe\ + \ value (frames are 0, 1, or 2 like\ + \ USCS definitions). For use with trim-frame." ) + ) + <*> optional ( option auto + ( long "input-outframe" + <> metavar "[ ] | FIELD" + <> help "Represents the 1 indexed field split by '|'\ + \ containing the outframe\ + \ value (frames are 0, 1, or 2 like\ + \ USCS definitions). For use with trim-frame." ) + ) + <*> switch + ( long "trim-frame" + <> short 'y' + <> help "Trim each sequence to be in frame by remove extra nucleotides\ + \ at the end. If input-inframe or input-outframe is\ + \ specified, follow those rules instead." ) + <*> option auto + ( long "input-stop-range" + <> short 'r' + <> metavar "[106]|INT" + <> value 106 + <> help "Only search for stops with remove-stops up to this\ + \ amino acid position" ) + <*> option auto + ( long "input-codon-mut" + <> short 'c' + <> metavar "[-1]|0|1|2|3" + <> value (-1) + <> help "Only include codons with this many mutations or less or more,\ + \ depending on input-codon-mut-type (-1 is the same as include\ + \ all codons). Converts the unincluded codon to gaps." ) + <*> strOption + ( long "input-codon-mut-type" + <> short 'T' + <> metavar "[=]|>|<" + <> value "=" + <> help "Only include codons with this many mutations (=)\ + \ (or lesser (<) or greater (>), depending on\ + \ input-codon-mut). Converts the unincluded codon to gaps.\ + \ For use with input-codon-mut." ) + <*> strOption + ( long "input-mut-type" + <> short 'M' + <> metavar "[All]|Silent|Replacement" + <> value "All" + <> help "Only include codons with all mutations (All),\ + \ (or silent (Silent) or replacement (Replacement)).\ + \ For use with input-codon-mut." ) + <*> strOption + ( long "input-change-field" + <> short 'e' + <> metavar "((FIELD (Int), VALUE (String))" + <> value "" + <> help "Change a field to a match, so a regex \"ch.*_\" to field 2\ + \ of \">abc|brie_cheese_dude\" would result in\ + \ \">abc|cheese_\". Useful for getting specific properties\ + \ from a field. Can take a list of format\ + \ \"(Int, String)&&(Int, String)&& ...\" and so on. The String\ + \ is in regex format (POSIX extended).\ + \ The first in the tuple is the location of the field\ + \ (1 indexed, split by '|')." ) + <*> strOption + ( long "input-custom-filter" + <> short 'f' + <> metavar "((FIELD (Int), VALUE (String))" + <> value "" + <> help "A custom filter. Can take a list of format\ + \ \"(Int, String)&&(Int, String)&& ...\" and so on. The String\ + \ is in regex format (POSIX extended), so if the entire\ + \ string 'a' is in the whole field, then you need to input\ + \ '^a$' for the beginning to the end!\ + \ The first in the tuple is the location of the field\ + \ (1 indexed, split by '|'). If you want to apply to\ + \ the entire header, either have the location as 0 or\ + \ exclude the location altogether (, Day 3|IGHV3) for instance\ + \ will match if the entire header is '>Day 3|IGHV3'.\ + \ This list will be filtered one at a time, so you cannot\ + \ get multiple filters, but you can remove multiple filters." ) + <*> switch + ( long "legacy-custom-germline" + <> short 'G' + <> help "Whether to apply the custom filter to germlines (>>)\ + \ instead of sequences (>). LEGACY ONLY" ) + <*> switch + ( long "custom-remove" + <> short 'm' + <> help "Whether to remove the sequences containing the custom filter\ + \ as opposed to remove the sequences that don't contain the\ + \ filter" ) + <*> option auto + ( long "legacy-gene-allele-field" + <> short 'V' + <> metavar "[1]|INT" + <> value 1 + <> help "The field (1 indexed) of the gene allele name. LEGACY ONLY" ) + <*> switch + ( long "legacy-count" + <> short 'v' + <> help "Do not save output, just count genes and alleles from\ + \ the results. Requires gene-allele-field. LEGACY ONLY" ) + <*> strOption + ( long "output" + <> short 'o' + <> metavar "FILE" + <> value "" + <> help "The output fasta file" ) + +-- | Parse the argument fieldInt +fieldIntParser :: String -> [(Maybe Int, T.Text)] +fieldIntParser "" = [] +fieldIntParser s = map (\x -> (first x, second x)) . Split.splitOn "&&" $ s + where + first x + | (head . Split.splitOn "," $ x) == "(" = Nothing + | otherwise = + Just (read (tail . head . Split.splitOn "," $ x) :: Int) + second = T.pack . init . dropWhile (== ' ') . last . Split.splitOn "," + +-- | Check for amino acid +isAminoAcid :: GeneticUnit -> Bool +isAminoAcid AminoAcid = True +isAminoAcid _ = False + +modifyFastaList :: Options -> IO () +modifyFastaList opts = do + hIn <- if null . input $ opts + then return IO.stdin + else IO.openFile (input opts) IO.ReadMode + hOut <- if null . output $ opts + then return IO.stdout + else IO.openFile (output opts) IO.WriteMode + let genUnit = aminoAcidsFlag opts + stopRange = inputStopRange opts + customFilters = fieldIntParser . inputCustomFilter $ opts + changeFields = fieldIntParser . inputChangeField $ opts + codonMut = inputCodonMut opts + codonMutType = T.pack . inputCodonMutType $ opts + mutType = T.pack . inputMutType $ opts + + -- Remove out of frame sequences + seqInFrame x = not ( removeOutOfFrameFlag opts + && (not . isAminoAcid $ genUnit) + ) + || isInFrame x + + -- Start filtering out sequences + -- Include only custom filter sequences + customFilter x = null customFilters + || hasAllCustomFilters + (customRemoveFlag opts) + customFilters + x + + -- Remove clones with stops in the range + noStops x = not (removeStopsFlag opts) || hasNoStops genUnit stopRange x + + -- Remove Ns from fasta list + noNs = if removeTheNsFlag opts && (not . isAminoAcid $ genUnit) + then removeN + else id + + -- Change fasta headers with match + changeHeader x = if not . null $ changeFields + then changeAllFields x changeFields + else x + + -- Get a specific region of the sequence + cutSequence = case (inputStart opts, inputStop opts) of + (Nothing, Nothing) -> id + (start, stop) -> getRegionSequence start stop + + -- Remove non standard nucleotides + removeUnknown = if removeUnknownNuc opts + then removeUnknownNucs + else id + + -- Trim sequence + trim fs = + if trimFrame opts + then trimFasta + genUnit + ((read . T.unpack . flip getField fs) <$> inputInFrame opts) + ((read . T.unpack . flip getField fs) <$> inputOutFrame opts) + fs + else fs + + -- Fill in bad characters at the requested section with possible + -- replacements + fillIn = case inputFillIn opts of + (-1, -1, 'X') -> id + (f, s, c) -> fillInSequence f s c + + -- Convert to amino acids + ntToaa = if convertToAminoAcidsFlag opts + then convertToAminoAcidsFastaSequence + else id + + -- Include sequence length in header at the end + includeLength = if addLengthFlag opts + then addLengthHeader + else id + + -- CLIP fasta specific filters and transformations + + -- Remove highly mutated sequences + removeHighMutations = if removeHighlyMutatedFlag opts + then filterHighlyMutatedEntry genUnit + else id + + -- Extract mutations to a certain degree + getMutations = if codonMut > -1 + then onlyMutations codonMut codonMutType mutType + else id + + -- Extract mutations to a certain degree + getFrequentMutations = if isJust (inputMutationCount opts) + || isJust (inputMutationPercent opts) + then frequentMutations + (inputMinSequenceMutation opts) + (inputMutationCount opts) + (inputMutationPercent opts) + else id + + -- Final order + filterOrder x = seqInFrame x && customFilter x && noStops x + transformOrder = includeLength + . ntToaa + . changeHeader + . removeUnknown + . trim + . noNs + . fillIn + . cutSequence + -- Specifically for germlines, as we don't want to change header or + -- fill in the germline because that would make no sense in this + -- case + transformGermline = includeLength + . ntToaa + . removeUnknown + . trim + . noNs + . cutSequence + -- Specifically for CLIP fasta files + transformOrderCLIP = getFrequentMutations + . getMutations + . removeHighMutations + + executePrintFasta x = mappend (showFasta x) (T.pack "\n") + executePrintCLIPFasta x = + mappend + (printCloneEntry (removeGermlinesPreFlag opts) x) + (T.pack "\n") + + -- Execute pipes + if clipFastaFlag opts + then + runEffect $ ( ( pipesCLIPFasta (PT.fromHandle hIn) + >-> P.map ( \(!germline, !fseqs) -> + (germline, filter filterOrder fseqs) + ) -- Filter sequences + >-> P.map transformOrderCLIP -- Transform specifically for CLIP fasta + >-> P.map ( \(!germline, !fseqs) -> + ( transformGermline germline + , map transformOrder fseqs + ) + ) -- Transform sequences, only do some for germline + >-> P.map ( \(!germline, !fs) + -> ( germline + , filter (not . T.null . fastaSeq) fs + ) + ) -- Remove empty sequences + >-> P.filter (not . null . snd) -- Remove empty clones + >-> P.map executePrintCLIPFasta ) -- Print the results + >> yield (T.pack "\n") ) -- want that newline at the end + >-> PT.toHandle hOut + else + runEffect $ ( ( pipesFasta (PT.fromHandle hIn) + >-> P.filter filterOrder -- Filter + >-> P.map transformOrder -- Transform + >-> P.filter (not . T.null . fastaSeq) -- Remove empty sequences + >-> P.map executePrintFasta ) -- Print the results + >> 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) + +-- Legacy function +modifyFastaCloneMap :: Options -> IO () +modifyFastaCloneMap opts = do + contents <- if null . input $ opts + then T.getContents + else T.readFile . input $ opts + -- No redundant newlines in sequence + let genUnit = aminoAcidsFlag opts + stopRange = inputStopRange opts + codonMut = inputCodonMut opts + codonMutType = T.pack . inputCodonMutType $ opts + mutType = T.pack . inputMutType $ opts + customFilters = fieldIntParser $ inputCustomFilter opts + removeGermlinesFlag = if not . clipFastaFlag $ opts + then True + else removeGermlinesPreFlag opts + + -- Initiate CloneMap + cloneMapFrames = if not . clipFastaFlag $ opts + then addFillerGermlines + . parsecFasta + $ contents + else parsecCLIPFasta contents + + -- Remove out of frame sequences + cloneMapInFrame = if (removeOutOfFrameFlag opts && ( not + . isAminoAcid + $ genUnit ) ) + then removeOutOfFrameSeqs cloneMapFrames + else cloneMapFrames + + -- Remove Ns from CloneMap + cloneMap = if (removeTheNsFlag opts && ( not + . isAminoAcid + $ genUnit ) ) + then removeCLIPNs cloneMapInFrame + else cloneMapInFrame + + -- Start filtering out sequences + -- Include only custom filter sequences + cloneMapCustom = if not . null $ customFilters + then removeAllCustomFilters + (customGermlineFlag opts) + (customRemoveFlag opts) + cloneMap + customFilters + else cloneMap + -- Remove clones with stops in the range + (cloneMapNoStops, errorString) = if removeStopsFlag opts + then removeStopsCloneMap + genUnit + stopRange + cloneMapCustom + else (cloneMapCustom, Nothing) + -- Output Error if necessary + case errorString of + Nothing -> return () + Just x -> error x + + -- Remove duplicate sequences + let cloneMapNoDuplicates = if removeDuplicatesFlag opts + then removeDuplicatesCloneMap + cloneMapNoStops + else cloneMapNoStops + + -- Remove clones that are highly mutated + (cloneMapLowMutation, errorString2) = if (removeHighlyMutatedFlag opts) + && ( not + . isAminoAcid + $ genUnit ) + then filterHighlyMutated + genUnit + cloneMapNoDuplicates + else ( cloneMapNoDuplicates + , Nothing ) + + -- Output Error if necessary + case errorString2 of + Nothing -> return () + Just x -> error x + + -- Remove codons with codons with a certain number of mutations + let cloneMapNoCodonMut = if codonMut > -1 + then removeCodonMutCount codonMut + codonMutType + mutType + cloneMapLowMutation + else cloneMapLowMutation + -- Remove empty clones + cloneMapNoEmptyClones = removeEmptyClone cloneMapNoCodonMut + + -- Convert sequences to amino acids + (cloneMapAA, errorString3) = if (convertToAminoAcidsFlag opts) + && (not . isAminoAcid $ genUnit) + then convertToAminoAcidsCloneMap + cloneMapNoEmptyClones + else (cloneMapNoEmptyClones, Nothing) + + -- Output Error if necessary + case errorString3 of + Nothing -> return () + Just x -> error x + + -- Break if there are no sequences to output + case M.null cloneMapAA of + True -> error "No sequences left! Nothing written." + False -> return () + + -- What to do with results + case countFlag opts of + True -> do + -- Print results + let outputText = printSequenceCount + (clipFastaFlag opts) + (inputGeneAlleleField opts) + cloneMapAA + -- Print results to stdout + T.putStrLn outputText + False -> do + -- Print results + let outputText = if removeGermlinesFlag + then printFastaNoGermline cloneMapAA + else printFasta cloneMapAA + + -- Save results + if null . output $ opts + then T.putStrLn outputText + else T.writeFile (output opts) outputText + +modifyFasta :: Options -> IO () +modifyFasta opts = if legacyFlag opts + then modifyFastaCloneMap opts + else modifyFastaList opts + +main :: IO () +main = execParser opts >>= modifyFasta + where + opts = info (helper <*> options) + ( fullDesc + <> progDesc "Modify fasta (and CLIP) files in several optional ways.\ + \ Order of transformation goes: seqInFrame -> customFilter\ + \ -> noStops -> removeHighMutations -> getMutations ->\ + \ getFrequentMutations -> cutSequence -> fillIn -> noNs\ + \ -> changeHeader -> ntToaa -> includeLength,\ + \ so if you require a different\ + \ order (which can change results dramatically), then do\ + \ so one at a time through the wonderful world of piping." + <> header "modify-fasta, Gregory W. Schwartz" ) 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.2 +version: 0.8.0.3 synopsis: Modify fasta (and CLIP) files in several optional ways -- description: homepage: https://github.com/GregorySchwartz/modify-fasta @@ -15,23 +15,39 @@ build-type: Simple cabal-version: >=1.8 +library + ghc-options: -O2 + hs-source-dirs: src + exposed-modules: Diversity + , FilterCloneList + , FilterCloneMap + , FilterFastaList + , Print + , TransformCloneList + , TransformFastaList + , Types + , Utility + build-depends: base >=4.6 && <5 + , containers >=0.5 + , text + , text-show + , split >=0.2 + , regex-tdfa >=1.2 + , regex-tdfa-text + , fasta + executable modify-fasta + hs-source-dirs: app main-is: Main.hs - -- other-modules: - build-depends: base >=4.6 && <5 + build-depends: modify-fasta + , base >=4.6 && <5 , containers >=0.5 , mtl >=2.1 , text - , text-show , split >=0.2 - , optparse-applicative >=0.11 , fasta - , regex-tdfa >=1.2 - , regex-tdfa-text , pipes >= 4.1 , pipes-text - + , optparse-applicative >=0.11 -- Directories containing source files. - hs-source-dirs: src - ghc-options: -O2 diff --git a/src/Diversity.hs b/src/Diversity.hs new file mode 100644 --- /dev/null +++ b/src/Diversity.hs @@ -0,0 +1,57 @@ +-- Diversity module. +-- By G.W. Schwartz +-- +-- Collection of functions pertaining to finding the diversity of samples. + +module Diversity where + +-- Built-in +import Data.List +import qualified Data.Text as T + +-- Takes two strings, returns Hamming distance +hamming :: T.Text -> T.Text -> Int +hamming xs ys = length $ filter (not . uncurry (==)) $ T.zip xs ys + +-- Returns the diversity of a list of things +diversity :: (Ord b) => Double -> [b] -> Double +diversity order sample + | length sample == 0 = 0 + | order == 1 = exp . h $ speciesList + | otherwise = (sum . map ((** order) . p_i) $ speciesList) ** pow + where + pow = 1 / (1 - order) + h = negate . sum . map (\x -> (p_i x) * (log (p_i x))) + p_i x = ((fromIntegral . length $ x) :: Double) / + ((fromIntegral . length $ sample) :: Double) + speciesList = group . sort $ sample + +-- Calculates the binary coefficient +choose :: (Integral a) => a -> a -> a +choose _ 0 = 1 +choose 0 _ = 0 +choose n k = choose (n - 1) (k - 1) * n `div` k + +-- Returns the rarefaction curve for each position in a list +rarefactionCurve :: (Eq a, Ord a) => [a] -> [Double] +rarefactionCurve xs = map rarefact [1..n_total] + where + rarefact n + | n == 0 = 0 + | n == 1 = 1 + | n == n_total = k + | otherwise = k - ((1 / (fromIntegral (choose n_total n))) * inner n) + inner n = fromIntegral . + sum . + map (\g -> choose (n_total - length g) n) $ + grouped + n_total = length xs + k = genericLength grouped + grouped = group . sort $ xs + +-- Calculates the percent of the curve that is above 95% of height of the curve +rarefactionViable :: [Double] -> Double +rarefactionViable xs = (genericLength valid / genericLength xs) * 100 + where + valid = dropWhile (< (0.95 * last xs)) xs + diff --git a/src/FilterCloneList.hs b/src/FilterCloneList.hs new file mode 100644 --- /dev/null +++ b/src/FilterCloneList.hs @@ -0,0 +1,71 @@ +-- FilterFastaList module. +-- By Gregory W. Schwartz +-- +-- Collection of functions for the filtering of a pipesFasta + +{-# LANGUAGE BangPatterns, OverloadedStrings, FlexibleContexts #-} + +module FilterCloneList ( filterHighlyMutatedEntry + ) where + +-- Built in +import Data.List +import Data.Maybe +import Data.Either +import Text.Regex.TDFA +import Text.Regex.TDFA.Text +import qualified Data.Text as T + +-- Cabal +import Data.Fasta.Text + +-- Local +import Types + +-- Remove highly mutated sequences (sequences with more than a third of +-- their sequence being mutated). +filterHighlyMutatedEntry :: GeneticUnit -> CloneEntry -> CloneEntry +filterHighlyMutatedEntry !genUnit = newEntry + where + newEntry (!germline, !fseqs) = ( germline + , map snd + . filter (not . fst) + . rights + . assignMutated germline + $ fseqs + ) + assignMutated k = map (isHighlyMutated k) + isHighlyMutated !k !x = + case (readSeq genUnit k, readSeq genUnit x) of + ((Right a), (Right b)) -> (\n -> Right (n, b)) + $ ( (fromIntegral (T.length (fastaSeq a)) :: Double) + / 3 ) + <= ( ( genericLength + . realMutations (fastaSeq a) + $ fastaSeq b ) ) + ((Left a), (Right _)) -> error ("Error in germline: " ++ T.unpack a) + ((Right _), (Left b)) -> error ("Error in sequence: " ++ T.unpack b) + ((Left a), (Left b)) -> error (unwords [ "Error in sequence:" + , T.unpack b + , "with germline:" + , T.unpack a ] ) + realMutations k x = filterCodonMutStab (\(!y, !z) -> y /= z) + . map snd + . mutation k + $ x + filterCodonMutStab isWhat = filter (filterRules genUnit isWhat) + filterRules AminoAcid isWhat x = isWhat x + && not (inTuple '-' x) + && not (inTuple '.' x) + && not (inTuple '~' x) + filterRules Nucleotide isWhat x = isWhat x + && not (inTuple '-' x) + && not (inTuple '.' x) + && not (inTuple '~' x) + && not (inTuple 'N' x) + inTuple c (x, y) + | c == x || c == y = True + | otherwise = False + mutation x y = zip [1..] . T.zip x $ y + readSeq Nucleotide x = Right x + readSeq AminoAcid x = translate 1 x diff --git a/src/FilterCloneMap.hs b/src/FilterCloneMap.hs new file mode 100644 --- /dev/null +++ b/src/FilterCloneMap.hs @@ -0,0 +1,240 @@ +-- FilterCloneMap module. +-- By G.W. Schwartz +-- +-- Collection of functions for the filtering of a CloneMap + +{-# LANGUAGE BangPatterns, OverloadedStrings, FlexibleContexts #-} + +module FilterCloneMap where + +-- Built in +import Data.List +import Data.Char +import Data.Maybe +import Data.Either +import qualified Data.Set as S +import qualified Data.Map as M +import Text.Regex.TDFA +import Text.Regex.TDFA.Text +import qualified Data.Text as T + +-- Cabal +import Data.Fasta.Text + +-- Local +import Types +import Diversity + +-- Check if the data structure is Right +isRight' :: Either a b -> Bool +isRight' (Right _) = True +isRight' _ = False + +-- Altered version of listToMaybe +listToMaybe' :: [a] -> Maybe [a] +listToMaybe' [] = Nothing +listToMaybe' x = Just x + +-- 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) + where + newCloneMap = M.map (map snd . filter (not . fst) . rights) + errorCloneMap + errorString = listToMaybe' + . unlines + . filter (not . null) + . map snd + . M.toAscList + . M.map (intercalate "\n" . lefts) + $ errorCloneMap + errorCloneMap = M.mapWithKey assignMutated cloneMap + assignMutated k = map (isHighlyMutated (snd k)) + isHighlyMutated !k !x = + case (readSeq genUnit k, readSeq genUnit x) of + ((Right a), (Right b)) -> (\n -> Right (n, b)) + $ ( (fromIntegral (T.length (fastaSeq a)) :: Double) + / 3 ) + <= ( ( genericLength + . realMutations (fastaSeq a) + $ fastaSeq b ) ) + ((Left a), (Right _)) -> Left (unwords ["Germline: ", T.unpack a]) + ((Right _), (Left b)) -> Left (unwords ["Sequence: ", T.unpack b]) + ((Left a), (Left b)) -> Left (unwords [ "Sequence:" + , T.unpack b + , "with Germline:" + , T.unpack a ] ) + realMutations k x = filterCodonMutStab (\(!y, !z) -> y /= z) + . map snd + . mutation k + $ x + filterCodonMutStab isWhat = filter (filterRules genUnit isWhat) + filterRules AminoAcid isWhat x = isWhat x + && not (inTuple '-' x) + && not (inTuple '.' x) + && not (inTuple '~' x) + filterRules Nucleotide isWhat x = isWhat x + && not (inTuple '-' x) + && not (inTuple '.' x) + && not (inTuple '~' x) + && not (inTuple 'N' x) + inTuple c (x, y) + | c == x || c == y = True + | otherwise = False + mutation x y = zip [1..] . T.zip x $ y + readSeq Nucleotide x = Right x + readSeq AminoAcid x = translate 1 x + +-- Replace codons that have more than CodonMut mutations (make them "---" +-- codons). +removeCodonMutCount :: CodonMut -> T.Text -> T.Text -> CloneMap -> CloneMap +removeCodonMutCount codonMut codonMutType mutType = M.mapWithKey mapRemove + where + mapRemove (_, germ) = map (removeCodon germ) + removeCodon germ clone = clone { fastaSeq + = remove (fastaSeq germ) + . fastaSeq $ clone } + remove germSeq = mconcat + . map (snd . replaceCodon) + . zip (codonSplit germSeq) + . codonSplit + replaceCodon (x, y) + | (codonMutOp codonMutType) (hamming x y) codonMut + && isMutType (T.toUpper mutType) x y = (x, y) + | otherwise = ("---", "---") + codonSplit = fullCodon . T.chunksOf 3 + fullCodon = filter ((== 3) . T.length) + codonMutOp ">" = (>) + codonMutOp "<" = (<) + codonMutOp "=" = (==) + isMutType "REPLACEMENT" x y = codon2aa x /= codon2aa y + isMutType "SILENT" x y = codon2aa x == codon2aa y + isMutType _ _ _ = True + +-- Remove clone sequences that have stop codons in the first stopRange +-- codons +removeStopsCloneMap :: GeneticUnit + -> Int + -> CloneMap + -> (CloneMap, Maybe String) +removeStopsCloneMap !genUnit !stopRange !cloneMap = ( newCloneMap + , errorString ) + where + errorString = listToMaybe' + . unlines + . filter (not . null) + . map snd + . M.toAscList + . M.map ( intercalate "\n" + . map T.unpack + . lefts + . map (translate 1) + ) + $ cloneMap + newCloneMap = M.map (filter (filterStops genUnit)) cloneMap + filterStops Nucleotide x = (isRight' . translate 1 $ x) + && ( not + . T.isInfixOf "*" + . T.take stopRange + . fastaSeq + . fromEither + . translate 1 ) x + filterStops AminoAcid x = not + . T.isInfixOf "*" + . T.take stopRange + . fastaSeq + $ x + fromEither (Right x) = x + fromEither (Left x) = error (T.unpack x) + +-- Remove duplicate sequences +removeDuplicatesCloneMap :: CloneMap -> CloneMap +removeDuplicatesCloneMap cloneMap = M.map + (filter (`S.member` duplicateSet)) + cloneMap + where + duplicateSet = S.fromList + . nubBy (\x y -> fastaSeq x == fastaSeq y) + . concatMap snd + . M.toAscList + $ cloneMap + +-- Remove out of frame sequences +removeOutOfFrameSeqs :: CloneMap -> CloneMap +removeOutOfFrameSeqs = M.map (filter isInFrame) + where + isInFrame = (== 0) + . mod 3 + . T.length + . T.filter (\x -> not $ T.isInfixOf (T.singleton x) ".-") + . fastaSeq + +-- Remove sequences that do not contain the string customFilter in the +-- customField location, split by "|". Note that this is 1 indexed and +-- 0 means to search the entire header for the customFilter. If the +-- customRemove option is enabled, this function will instead remove +-- sequences that have headers which match the custom filter, as opposed to +-- the other way around (this is defined in the "equal" function). Also +-- takes into account whether to filter on the germline versus the actual +-- sequences. +removeCustomFilter :: Bool + -> Bool + -> Maybe Int + -> T.Text + -> CloneMap + -> CloneMap +removeCustomFilter germ rm customField customFilter cloneMap + | germ && ((customField == Just 0) || (isNothing customField)) + = M.filterWithKey (\(_, k) _ -> inField k) cloneMap + | germ && customField > Just 0 + = M.filterWithKey (\(_, k) _ -> inCustomField k) cloneMap + | (customField == Just 0) || (isNothing customField) + = M.map (filter inField) cloneMap + | customField > Just 0 = + M.map (filter inCustomField) cloneMap + where + inField = equal rm customFilter . fastaHeader + inCustomField x = equal rm customFilter + . (!!) (T.splitOn "|" . fastaHeader $ x) + $ (fromJust customField - 1) + equal False x y = y =~ x :: Bool + equal True x y = not . equal False x $ y + +removeAllCustomFilters :: Bool + -> Bool + -> CloneMap + -> [(Maybe Int, T.Text)] + -> CloneMap +removeAllCustomFilters germ rm = foldl' filterMap + where + filterMap acc (x, y) = removeCustomFilter germ rm x y acc + +-- Remove clones that do not have any sequences after the filtrations +removeEmptyClone :: CloneMap -> CloneMap +removeEmptyClone = M.filter (not . null) + +-- Convert sequences to amino acids +convertToAminoAcidsCloneMap :: CloneMap -> (CloneMap, Maybe String) +convertToAminoAcidsCloneMap cloneMap = (newCloneMap, errorString) + where + newCloneMap = M.mapKeysWith (++) (\(!x, !y) -> (x, fromEither y)) + . M.filterWithKey (\(_, !y) _ -> isRight' y) + . M.map rights + $ errorCloneMap + errorString = listToMaybe' + . concatMap snd + . M.toAscList + . M.mapWithKey (\(_, !y) v -> (++) (eitherToString y) + . concatMap T.unpack + . lefts + $ v ) + $ errorCloneMap + errorCloneMap = M.mapKeys keyMap + . M.map (map (translate 1)) + $ cloneMap + keyMap (!x, !y) = (x, translate 1 y) + eitherToString (Right _) = "" + eitherToString (Left x) = T.unpack x + fromEither (Right x) = x + fromEither (Left x) = error (T.unpack x) diff --git a/src/FilterFastaList.hs b/src/FilterFastaList.hs new file mode 100644 --- /dev/null +++ b/src/FilterFastaList.hs @@ -0,0 +1,83 @@ +-- FilterFastaList module. +-- By Gregory W. Schwartz +-- +-- Collection of functions for the filtering of a pipesFasta + +{-# LANGUAGE OverloadedStrings, FlexibleContexts #-} + +module FilterFastaList ( hasNoStops + , isInFrame + , hasCustomFilter + , hasAllCustomFilters + ) where + +-- Built in +import Data.List +import Data.Maybe +import Text.Regex.TDFA +import Text.Regex.TDFA.Text +import qualified Data.Text as T + +-- Cabal +import Data.Fasta.Text + +-- Local +import Types + +-- Remove clone sequences that have stop codons in the first stopRange +-- codons +hasNoStops :: GeneticUnit + -> Int + -> FastaSequence + -> Bool +hasNoStops genUnit 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 AminoAcid = Right . not . T.isInfixOf "*" . T.take stopRange . fastaSeq + +-- Remove out of frame sequences +isInFrame :: FastaSequence -> Bool +isInFrame = (== 0) + . mod 3 + . T.length + . T.filter (\x -> not . T.isInfixOf (T.singleton x) $ ".-") + . fastaSeq + +-- Remove sequences that do not contain the string customFilter in the +-- customField location, split by "|". Note that this is 1 indexed and +-- 0 means to search the entire header for the customFilter. If the +-- customRemove option is enabled, this function will instead remove +-- sequences that have headers which match the custom filter, as opposed to +-- the other way around (this is defined in the "equal" function). Also +-- takes into account whether to filter on the germline versus the actual +-- sequences. +hasCustomFilter :: Bool + -> Maybe Int + -> T.Text + -> FastaSequence + -> Bool +hasCustomFilter rm customField customFilter fasta + | customField == Just 0 || isNothing customField = inField fasta + | customField > Just 0 = inCustomField fasta + where + inField = equal rm customFilter . fastaHeader + inCustomField x = equal rm customFilter + . (!!) (T.splitOn "|" . fastaHeader $ x) + $ (fromJust customField - 1) + equal :: Bool -> T.Text -> T.Text -> Bool + equal False x y = y =~ x :: Bool + equal True x y = not . equal False x $ y + +hasAllCustomFilters :: Bool + -> [(Maybe Int, T.Text)] + -> FastaSequence + -> Bool +hasAllCustomFilters rm filters f = all filterMap filters + where + filterMap (x, y) = hasCustomFilter rm x y f diff --git a/src/Main.hs b/src/Main.hs deleted file mode 100644 --- a/src/Main.hs +++ /dev/null @@ -1,627 +0,0 @@ --- modify-fasta --- By Gregory W. Schwartz - --- Takes a fasta file filters the fasta file in several ways. - -{-# LANGUAGE BangPatterns #-} - --- Built-in -import Data.Maybe -import qualified Data.Map as M -import qualified System.IO as IO -import qualified Data.Text.IO as IO -import Control.Monad - --- Cabal -import qualified Data.Text as T -import qualified Data.Text.IO as T -import Options.Applicative -import Data.Fasta.Text -import Pipes -import qualified Pipes.Prelude as P -import qualified Pipes.Text as PT -import qualified Pipes.Text.IO as PT -import qualified Data.List.Split as Split - --- Local -import Types -import Utility -import FilterCloneMap -import FilterFastaList -import FilterCloneList -import TransformFastaList -import TransformCloneList -import Print - --- Command line arguments -data Options = Options { input :: String - , aminoAcidsFlag :: GeneticUnit - , legacyFlag :: Bool - , clipFastaFlag :: Bool - , convertToAminoAcidsFlag :: Bool - , inputFillIn :: FillInValue - , inputStart :: Maybe Int - , inputStop :: Maybe Int - , inputMinSequenceMutation :: Maybe Int - , inputMutationCount :: Maybe Int - , inputMutationPercent :: Maybe Double - , addLengthFlag :: Bool - , removeTheNsFlag :: Bool - , removeGermlinesPreFlag :: Bool - , removeHighlyMutatedFlag :: Bool - , removeStopsFlag :: Bool - , removeDuplicatesFlag :: Bool - , removeOutOfFrameFlag :: Bool - , removeUnknownNuc :: Bool - , inputInFrame :: Maybe Field - , inputOutFrame :: Maybe Field - , trimFrame :: Bool - , inputStopRange :: Int - , inputCodonMut :: CodonMut - , inputCodonMutType :: String - , inputMutType :: String - , inputChangeField :: String - , inputCustomFilter :: String - , customGermlineFlag :: Bool - , customRemoveFlag :: Bool - , inputGeneAlleleField :: Int - , countFlag :: Bool - , output :: String - } - --- Command line options -options :: Parser Options -options = Options - <$> strOption - ( long "input" - <> short 'i' - <> metavar "FILE" - <> value "" - <> help "The input fasta file or CLIP fasta file" ) - <*> option auto - ( long "unit" - <> short 'u' - <> metavar "AminoAcid | Nucleotide" - <> help "Whether these sequences are composed of\ - \ amino acids (AminoAcid) or nucleotides (Nucleotide)" ) - <*> switch - ( long "legacy" - <> short 'L' - <> help "Whether to use the legacy version with no pipes. Note: The\ - \ legacy version supports more features but is greedy\ - \ in terms of speed and memory. Use only if really needed.\ - \ Features that are legacy only are noted in this\ - \ documentation" ) - <*> switch - ( long "clip-fasta" - <> short 'A' - <> help "Whether the input is a clip fasta file (has germline >>\ - \ sequences)." ) - <*> switch - ( long "convert-to-amino-acids" - <> short 'C' - <> help "Whether to convert the filtered sequences to amino acids\ - \ in the output. Applied last, even after add length." ) - <*> option auto - ( long "fill-in" - <> short 'F' - <> metavar "(FIELD, START, 'CHARACTER')" - <> value (-1, -1, 'X') - <> help "Use the FIELD index (1 indexed, split by '|') in the\ - \ header to fill in the unknown character CHARACTER with the\ - \ respective character in that field. Use the START value\ - \ to let the program know where the string in the field starts\ - \ from in the sequence. For instance, a header of >H2O|HELLO\ - \ for the sequence HIMANHELXODUDE and a value of fill-in of\ - \ (2, 6, 'X') would change the sequence to HIMANHELLODUDE. If\ - \ the string in the field has the unknown character as well,\ - \ i.e. >H2O|HELXO, then the sequence is considered bad and is\ - \ removed" ) - <*> optional ( option auto - ( long "start" - <> short 't' - <> metavar "[ ] | INT" - <> help "Remove everything before this position (1 indexed).\ - \ Done first just after filtering." ) ) - <*> optional ( option auto - ( long "stop" - <> short 'p' - <> metavar "[ ] | INT" - <> help "Remove everything after this position (1 indexed).\ - \ Done first just after filtering." ) ) - <*> optional ( option auto - ( long "frequent-mutation-min" - <> metavar "[ ] | INT" - <> help "Minimum number of sequences required for a clone to be valid\ - \ in the calculation of frequent mutations, replaces\ - \ with gaps otherwise" ) ) - <*> optional ( option auto - ( long "frequent-mutation-count" - <> metavar "[ ] | INT" - <> help "Only include codons containing a mutation present in this\ - \ many sequences in the clone or more. 0 is all sequences.\ - \ Converts the unincluded codons to gaps." ) ) - <*> optional ( option auto - ( long "frequent-mutation-percent" - <> metavar "[ ] | PERCENT" - <> help "Only include codons containing a mutation present in this\ - \ percentage of sequences in the clone or more.\ - \ 0 is all sequences. Converts the unincluded codons\ - \ to gaps." ) ) - <*> switch - ( long "add-length" - <> short 'l' - <> help "Whether to append the length of the sequence to the end of\ - \ the header, calculated after --convert-to-amino-acids\ - \ if enabled" ) - <*> switch - ( long "remove-N" - <> short 'N' - <> help "Whether to replace N or n in the sequence with a gap, '-'" ) - <*> switch - ( long "remove-germlines" - <> short 'g' - <> help "Whether to remove germlines." ) - <*> switch - ( long "remove-highly-mutated" - <> short 'h' - <> help "Whether to remove highly mutated clone sequences (a third\ - \ of their sequence are different amino acids)." ) - <*> switch - ( long "remove-stops" - <> short 's' - <> help "Whether to remove sequences with stop codons" ) - <*> switch - ( long "legacy-remove-duplicates" - <> short 'd' - <> help "Whether to remove duplicate sequences. LEGACY ONLY" ) - <*> switch - ( long "remove-out-of-frame" - <> short 'O' - <> help "Whether to remove sequences that are out of frame--if the\ - \ sequences or number of gaps is not divisible by 3" ) - <*> switch - ( long "remove-unknown-nucleotides" - <> short 'Y' - <> help "Convert unknown nucleotides (not ACGTN-.) to gaps (-)" ) - <*> optional ( option auto - ( long "input-inframe" - <> metavar "[ ] | FIELD" - <> help "Represents the 1 indexed field split by '|'\ - \ containing the inframe\ - \ value (frames are 0, 1, or 2 like\ - \ USCS definitions). For use with trim-frame." ) - ) - <*> optional ( option auto - ( long "input-outframe" - <> metavar "[ ] | FIELD" - <> help "Represents the 1 indexed field split by '|'\ - \ containing the outframe\ - \ value (frames are 0, 1, or 2 like\ - \ USCS definitions). For use with trim-frame." ) - ) - <*> switch - ( long "trim-frame" - <> short 'y' - <> help "Trim each sequence to be in frame by remove extra nucleotides\ - \ at the end. If input-inframe or input-outframe is\ - \ specified, follow those rules instead." ) - <*> option auto - ( long "input-stop-range" - <> short 'r' - <> metavar "[106]|INT" - <> value 106 - <> help "Only search for stops with remove-stops up to this\ - \ amino acid position" ) - <*> option auto - ( long "input-codon-mut" - <> short 'c' - <> metavar "[-1]|0|1|2|3" - <> value (-1) - <> help "Only include codons with this many mutations or less or more,\ - \ depending on input-codon-mut-type (-1 is the same as include\ - \ all codons). Converts the unincluded codon to gaps." ) - <*> strOption - ( long "input-codon-mut-type" - <> short 'T' - <> metavar "[=]|>|<" - <> value "=" - <> help "Only include codons with this many mutations (=)\ - \ (or lesser (<) or greater (>), depending on\ - \ input-codon-mut). Converts the unincluded codon to gaps.\ - \ For use with input-codon-mut." ) - <*> strOption - ( long "input-mut-type" - <> short 'M' - <> metavar "[All]|Silent|Replacement" - <> value "All" - <> help "Only include codons with all mutations (All),\ - \ (or silent (Silent) or replacement (Replacement)).\ - \ For use with input-codon-mut." ) - <*> strOption - ( long "input-change-field" - <> short 'e' - <> metavar "((FIELD (Int), VALUE (String))" - <> value "" - <> help "Change a field to a match, so a regex \"ch.*_\" to field 2\ - \ of \">abc|brie_cheese_dude\" would result in\ - \ \">abc|cheese_\". Useful for getting specific properties\ - \ from a field. Can take a list of format\ - \ \"(Int, String)&&(Int, String)&& ...\" and so on. The String\ - \ is in regex format (POSIX extended).\ - \ The first in the tuple is the location of the field\ - \ (1 indexed, split by '|')." ) - <*> strOption - ( long "input-custom-filter" - <> short 'f' - <> metavar "((FIELD (Int), VALUE (String))" - <> value "" - <> help "A custom filter. Can take a list of format\ - \ \"(Int, String)&&(Int, String)&& ...\" and so on. The String\ - \ is in regex format (POSIX extended), so if the entire\ - \ string 'a' is in the whole field, then you need to input\ - \ '^a$' for the beginning to the end!\ - \ The first in the tuple is the location of the field\ - \ (1 indexed, split by '|'). If you want to apply to\ - \ the entire header, either have the location as 0 or\ - \ exclude the location altogether (, Day 3|IGHV3) for instance\ - \ will match if the entire header is '>Day 3|IGHV3'.\ - \ This list will be filtered one at a time, so you cannot\ - \ get multiple filters, but you can remove multiple filters." ) - <*> switch - ( long "legacy-custom-germline" - <> short 'G' - <> help "Whether to apply the custom filter to germlines (>>)\ - \ instead of sequences (>). LEGACY ONLY" ) - <*> switch - ( long "custom-remove" - <> short 'm' - <> help "Whether to remove the sequences containing the custom filter\ - \ as opposed to remove the sequences that don't contain the\ - \ filter" ) - <*> option auto - ( long "legacy-gene-allele-field" - <> short 'V' - <> metavar "[1]|INT" - <> value 1 - <> help "The field (1 indexed) of the gene allele name. LEGACY ONLY" ) - <*> switch - ( long "legacy-count" - <> short 'v' - <> help "Do not save output, just count genes and alleles from\ - \ the results. Requires gene-allele-field. LEGACY ONLY" ) - <*> strOption - ( long "output" - <> short 'o' - <> metavar "FILE" - <> value "" - <> help "The output fasta file" ) - --- | Parse the argument fieldInt -fieldIntParser :: String -> [(Maybe Int, T.Text)] -fieldIntParser "" = [] -fieldIntParser s = map (\x -> (first x, second x)) . Split.splitOn "&&" $ s - where - first x - | (head . Split.splitOn "," $ x) == "(" = Nothing - | otherwise = - Just (read (tail . head . Split.splitOn "," $ x) :: Int) - second = T.pack . init . dropWhile (== ' ') . last . Split.splitOn "," - --- | Check for amino acid -isAminoAcid :: GeneticUnit -> Bool -isAminoAcid AminoAcid = True -isAminoAcid _ = False - -modifyFastaList :: Options -> IO () -modifyFastaList opts = do - hIn <- if null . input $ opts - then return IO.stdin - else IO.openFile (input opts) IO.ReadMode - hOut <- if null . output $ opts - then return IO.stdout - else IO.openFile (output opts) IO.WriteMode - let genUnit = aminoAcidsFlag opts - stopRange = inputStopRange opts - customFilters = fieldIntParser . inputCustomFilter $ opts - changeFields = fieldIntParser . inputChangeField $ opts - codonMut = inputCodonMut opts - codonMutType = T.pack . inputCodonMutType $ opts - mutType = T.pack . inputMutType $ opts - - -- Remove out of frame sequences - seqInFrame x = not ( removeOutOfFrameFlag opts - && (not . isAminoAcid $ genUnit) - ) - || isInFrame x - - -- Start filtering out sequences - -- Include only custom filter sequences - customFilter x = null customFilters - || hasAllCustomFilters - (customRemoveFlag opts) - customFilters - x - - -- Remove clones with stops in the range - noStops x = not (removeStopsFlag opts) || hasNoStops genUnit stopRange x - - -- Remove Ns from fasta list - noNs = if removeTheNsFlag opts && (not . isAminoAcid $ genUnit) - then removeN - else id - - -- Change fasta headers with match - changeHeader x = if not . null $ changeFields - then changeAllFields x changeFields - else x - - -- Get a specific region of the sequence - cutSequence = case (inputStart opts, inputStop opts) of - (Nothing, Nothing) -> id - (start, stop) -> getRegionSequence start stop - - -- Remove non standard nucleotides - removeUnknown = if removeUnknownNuc opts - then removeUnknownNucs - else id - - -- Trim sequence - trim fs = - if trimFrame opts - then trimFasta - genUnit - ((read . T.unpack . flip getField fs) <$> inputInFrame opts) - ((read . T.unpack . flip getField fs) <$> inputOutFrame opts) - fs - else fs - - -- Fill in bad characters at the requested section with possible - -- replacements - fillIn = case inputFillIn opts of - (-1, -1, 'X') -> id - (f, s, c) -> fillInSequence f s c - - -- Convert to amino acids - ntToaa = if convertToAminoAcidsFlag opts - then convertToAminoAcidsFastaSequence - else id - - -- Include sequence length in header at the end - includeLength = if addLengthFlag opts - then addLengthHeader - else id - - -- CLIP fasta specific filters and transformations - - -- Remove highly mutated sequences - removeHighMutations = if removeHighlyMutatedFlag opts - then filterHighlyMutatedEntry genUnit - else id - - -- Extract mutations to a certain degree - getMutations = if codonMut > -1 - then onlyMutations codonMut codonMutType mutType - else id - - -- Extract mutations to a certain degree - getFrequentMutations = if isJust (inputMutationCount opts) - || isJust (inputMutationPercent opts) - then frequentMutations - (inputMinSequenceMutation opts) - (inputMutationCount opts) - (inputMutationPercent opts) - else id - - -- Final order - filterOrder x = seqInFrame x && customFilter x && noStops x - transformOrder = includeLength - . ntToaa - . changeHeader - . removeUnknown - . trim - . noNs - . fillIn - . cutSequence - -- Specifically for germlines, as we don't want to change header or - -- fill in the germline because that would make no sense in this - -- case - transformGermline = includeLength - . ntToaa - . removeUnknown - . trim - . noNs - . cutSequence - -- Specifically for CLIP fasta files - transformOrderCLIP = getFrequentMutations - . getMutations - . removeHighMutations - - executePrintFasta x = mappend (showFasta x) (T.pack "\n") - executePrintCLIPFasta x = - mappend - (printCloneEntry (removeGermlinesPreFlag opts) x) - (T.pack "\n") - - -- Execute pipes - if clipFastaFlag opts - then - runEffect $ ( ( pipesCLIPFasta (PT.fromHandle hIn) - >-> P.map ( \(!germline, !fseqs) -> - (germline, filter filterOrder fseqs) - ) -- Filter sequences - >-> P.map transformOrderCLIP -- Transform specifically for CLIP fasta - >-> P.map ( \(!germline, !fseqs) -> - ( transformGermline germline - , map transformOrder fseqs - ) - ) -- Transform sequences, only do some for germline - >-> P.map ( \(!germline, !fs) - -> ( germline - , filter (not . T.null . fastaSeq) fs - ) - ) -- Remove empty sequences - >-> P.filter (not . null . snd) -- Remove empty clones - >-> P.map executePrintCLIPFasta ) -- Print the results - >> yield (T.pack "\n") ) -- want that newline at the end - >-> PT.toHandle hOut - else - runEffect $ ( ( pipesFasta (PT.fromHandle hIn) - >-> P.filter filterOrder -- Filter - >-> P.map transformOrder -- Transform - >-> P.filter (not . T.null . fastaSeq) -- Remove empty sequences - >-> P.map executePrintFasta ) -- Print the results - >> 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) - --- Legacy function -modifyFastaCloneMap :: Options -> IO () -modifyFastaCloneMap opts = do - contents <- if null . input $ opts - then T.getContents - else T.readFile . input $ opts - -- No redundant newlines in sequence - let genUnit = aminoAcidsFlag opts - stopRange = inputStopRange opts - codonMut = inputCodonMut opts - codonMutType = T.pack . inputCodonMutType $ opts - mutType = T.pack . inputMutType $ opts - customFilters = fieldIntParser $ inputCustomFilter opts - removeGermlinesFlag = if not . clipFastaFlag $ opts - then True - else removeGermlinesPreFlag opts - - -- Initiate CloneMap - cloneMapFrames = if not . clipFastaFlag $ opts - then addFillerGermlines - . parsecFasta - $ contents - else parsecCLIPFasta contents - - -- Remove out of frame sequences - cloneMapInFrame = if (removeOutOfFrameFlag opts && ( not - . isAminoAcid - $ genUnit ) ) - then removeOutOfFrameSeqs cloneMapFrames - else cloneMapFrames - - -- Remove Ns from CloneMap - cloneMap = if (removeTheNsFlag opts && ( not - . isAminoAcid - $ genUnit ) ) - then removeCLIPNs cloneMapInFrame - else cloneMapInFrame - - -- Start filtering out sequences - -- Include only custom filter sequences - cloneMapCustom = if not . null $ customFilters - then removeAllCustomFilters - (customGermlineFlag opts) - (customRemoveFlag opts) - cloneMap - customFilters - else cloneMap - -- Remove clones with stops in the range - (cloneMapNoStops, errorString) = if removeStopsFlag opts - then removeStopsCloneMap - genUnit - stopRange - cloneMapCustom - else (cloneMapCustom, Nothing) - -- Output Error if necessary - case errorString of - Nothing -> return () - Just x -> error x - - -- Remove duplicate sequences - let cloneMapNoDuplicates = if removeDuplicatesFlag opts - then removeDuplicatesCloneMap - cloneMapNoStops - else cloneMapNoStops - - -- Remove clones that are highly mutated - (cloneMapLowMutation, errorString2) = if (removeHighlyMutatedFlag opts) - && ( not - . isAminoAcid - $ genUnit ) - then filterHighlyMutated - genUnit - cloneMapNoDuplicates - else ( cloneMapNoDuplicates - , Nothing ) - - -- Output Error if necessary - case errorString2 of - Nothing -> return () - Just x -> error x - - -- Remove codons with codons with a certain number of mutations - let cloneMapNoCodonMut = if codonMut > -1 - then removeCodonMutCount codonMut - codonMutType - mutType - cloneMapLowMutation - else cloneMapLowMutation - -- Remove empty clones - cloneMapNoEmptyClones = removeEmptyClone cloneMapNoCodonMut - - -- Convert sequences to amino acids - (cloneMapAA, errorString3) = if (convertToAminoAcidsFlag opts) - && (not . isAminoAcid $ genUnit) - then convertToAminoAcidsCloneMap - cloneMapNoEmptyClones - else (cloneMapNoEmptyClones, Nothing) - - -- Output Error if necessary - case errorString3 of - Nothing -> return () - Just x -> error x - - -- Break if there are no sequences to output - case M.null cloneMapAA of - True -> error "No sequences left! Nothing written." - False -> return () - - -- What to do with results - case countFlag opts of - True -> do - -- Print results - let outputText = printSequenceCount - (clipFastaFlag opts) - (inputGeneAlleleField opts) - cloneMapAA - -- Print results to stdout - T.putStrLn outputText - False -> do - -- Print results - let outputText = if removeGermlinesFlag - then printFastaNoGermline cloneMapAA - else printFasta cloneMapAA - - -- Save results - if null . output $ opts - then T.putStrLn outputText - else T.writeFile (output opts) outputText - -modifyFasta :: Options -> IO () -modifyFasta opts = if legacyFlag opts - then modifyFastaCloneMap opts - else modifyFastaList opts - -main :: IO () -main = execParser opts >>= modifyFasta - where - opts = info (helper <*> options) - ( fullDesc - <> progDesc "Modify fasta (and CLIP) files in several optional ways.\ - \ Order of transformation goes: seqInFrame -> customFilter\ - \ -> noStops -> removeHighMutations -> getMutations ->\ - \ getFrequentMutations -> cutSequence -> fillIn -> noNs\ - \ -> changeHeader -> ntToaa -> includeLength,\ - \ so if you require a different\ - \ order (which can change results dramatically), then do\ - \ so one at a time through the wonderful world of piping." - <> header "modify-fasta, Gregory W. Schwartz" ) diff --git a/src/Print.hs b/src/Print.hs new file mode 100644 --- /dev/null +++ b/src/Print.hs @@ -0,0 +1,97 @@ +-- Print module +-- By Gregory W. Schwartz +-- +-- Collection of functions for the printing of data (converting data +-- structures into strings for use with writing to output files). + +{-# LANGUAGE BangPatterns, OverloadedStrings #-} + +module Print where + +-- Built-in +import Data.List +import qualified Data.Map as M +import qualified Data.Text as T + +-- Cabal +import Data.Fasta.Text +import TextShow + +-- Local +import Types + +-- Return the results of the filtration in text form for saving +-- to a file +printFasta :: CloneMap -> T.Text +printFasta = body + where + body = T.unlines + . map mapGerm + . M.toAscList + . M.map (T.intercalate "\n" . map showFasta) + mapGerm ((_, y), z) = mconcat [ ">>" + , fastaHeader y + , "\n" + , fastaSeq y + , "\n" + , z + ] + +-- Return the results of the filtration in text form for saving +-- to a file and excluding germline +printFastaNoGermline :: CloneMap -> T.Text +printFastaNoGermline = body + where + body = T.unlines + . map mapGerm + . M.toAscList + . M.map (T.intercalate "\n" . map showFasta) + mapGerm ((_, _), z) = z + +printSequenceCount :: Bool -> Int -> CloneMap -> T.Text +printSequenceCount clip idx s = body + where + body = T.unlines [ "Allele List: " + , "" + , "----------------------------------------------------" + , "" + , alleleCounts + , "" + , "----------------------------------------------------" + , "" + , "Gene List: " + , "" + , "----------------------------------------------------" + , "" + , geneCounts + , "" + , "----------------------------------------------------" + , "" + , mappend "Number of genes: " (showt . M.size $ germlineMap) + , mappend "Number of alleles: " (showt . M.size $ alleleMap) + ] + geneCounts = T.intercalate "\n" . map toLine . M.toAscList $ germlineMap + alleleCounts = T.intercalate "\n" . map fst . M.toAscList $ alleleMap + toLine (x, y) = x + `mappend` (T.replicate (30 - T.length x) " ") + `mappend` showt y + germlineMap = M.fromListWith (+) + . map (\(x, y) -> (head . T.splitOn "*" $ x, y)) + $ geneAlleleList + alleleMap = M.fromListWith (+) geneAlleleList + geneAlleleList = map (countProp clip) . M.toAscList $ s + countProp True ((_, x), y) = (getField idx x, length y) + countProp False ((_, _), y) = (getField idx . head $ y, 1) + getField f h = splitHeader h !! (f - 1) + splitHeader = T.splitOn "|" . fastaHeader + +-- | Takes a clone entry and returns a formatted text with or without +-- germline +printCloneEntry :: Bool -> CloneEntry -> T.Text +printCloneEntry False (!germline, !fseqs) = + T.intercalate "\n" [ T.cons '>' $ showFasta germline + , T.intercalate "\n" . map showFasta $ fseqs + ] +printCloneEntry True (!germline, !fseqs) = T.intercalate "\n" + . map showFasta + $ fseqs diff --git a/src/TransformCloneList.hs b/src/TransformCloneList.hs new file mode 100644 --- /dev/null +++ b/src/TransformCloneList.hs @@ -0,0 +1,145 @@ +-- TransformCloneList module. +-- By Gregory W. Schwartz +-- +-- Collection of functions that transform a clone entry in some way + +{-# LANGUAGE BangPatterns, OverloadedStrings #-} + +module TransformCloneList ( onlyMutations + , frequentMutations + ) where + +-- Built-in +import Data.Maybe +import Data.List +import qualified Data.Text as T +import qualified Data.Map.Strict as Map +import Data.Tuple +import Data.Monoid +import Control.Arrow (second) + +-- Cabal +import Data.Fasta.Text +import qualified Data.List.Split as Split + +-- Local +import Types +import Utility +import Diversity + +-- | Return True if there are no gap characters in the text +noGaps :: T.Text -> Bool +noGaps = not . any (\x -> x == '.' || x == '-') . T.unpack + +-- Replace codons that have more than CodonMut mutations (make them "---" +-- codons) and don't have gaps in them. +onlyMutations :: CodonMut -> T.Text -> T.Text -> CloneEntry -> CloneEntry +onlyMutations codonMut codonMutType mutType = newEntry + where + newEntry (!germ, !fseqs) = (germ, map (removeCodon germ) fseqs) + removeCodon germ fseq = fseq { fastaSeq = remove (fastaSeq germ) + . fastaSeq + $ fseq + } + remove germSeq = mconcat + . map (snd . replaceCodon) + . zip (codonSplit germSeq) + . codonSplit + replaceCodon (x, y) + | (codonMutOp codonMutType) (hamming x y) codonMut + && isMutType (T.toUpper mutType) x y + && noGaps x + && noGaps y = (x, y) + | otherwise = ("---", "---") + codonSplit = fullCodon . T.chunksOf 3 + fullCodon = filter ((== 3) . T.length) + codonMutOp ">" = (>) + codonMutOp "<" = (<) + codonMutOp "=" = (==) + isMutType "REPLACEMENT" x y = codon2aa x /= codon2aa y + isMutType "SILENT" x y = codon2aa x == codon2aa y + isMutType "ALL" _ _ = True + isMutType _ _ _ = error "Unknown mutation type" + +-- Only include codons containing mutations found in a certain number of +-- mutants +frequentMutations :: Maybe Int + -> Maybe Int + -> Maybe Double + -> CloneEntry + -> CloneEntry +frequentMutations minSeqs mutCount mutPercent entry@(!germline, !fseqs) = + newEntry + where + newEntry = (germline, map removeCodon fseqs) + removeCodon fseq = fseq { fastaSeq = replace fseq } + replace = rejoinCodons + . replaceCodon minSeqs mutCount mutPercent numSeqs countMap + . positionalCodons germline + rejoinCodons = T.pack . concatMap (map (snd . snd)) + countMap = getCountMap germline fseqs + numSeqs = length fseqs + +-- | Replace codons that are not valid with a gap "---". Important to note +-- that the predicates return False if we have what we want because the +-- monoid instance of Any has True trumping False, but we want False to +-- trump True, so invert it all and invert it back at the end. +replaceCodon :: Maybe Int + -> Maybe Int + -> Maybe Double + -> Int + -> CountMap + -> CodonMutations + -> CodonMutations +replaceCodon minSeqs mutCount mutPercent numSeqs countMap = map replace + where + replace x = if (any isValid . filter isMutation $ x) + && (not . any (isGapMut . snd) $ x) + then x + else map ((second . second) (const '-')) x + isGapMut (x, y) = x == '.' || x == '-' || y == '.' || y == '-' + isValid x = case Map.lookup x countMap of + (Just count) -> not + . fromMaybe False + . fmap getAny + $ minSeqsValid + <> mutCountValid count + <> mutPercentValid count + Nothing -> error ("Mutation not found: " ++ show x) + minSeqsValid = fmap (Any . not . (>=) numSeqs) minSeqs + mutCountValid count = case mutCount of + (Just 0) -> fmap (Any . not . (==) count) + $ Just numSeqs + (Just x) -> Just . Any . not . (>=) count $ x + Nothing -> Nothing + mutPercentValid count = + fmap + (Any . not . (>=) ((fromIntegral count / fromIntegral numSeqs) * 100)) + mutPercent + isMutation (_, (x, y)) = x /= y + +-- | Get the complete mutation codons with a position provided for each +-- nucleotide +positionalCodons :: Germline -> FastaSequence -> CodonMutations +positionalCodons germline = fullCodon . codonSplit . T.unpack . fastaSeq + where + codonSplit = fullCodon + . Split.chunksOf 3 + . zip [1..] + . zip (T.unpack . fastaSeq $ germline) + fullCodon = filter ((== 3) . length) + +-- | Get the number of times each mutation appears +getCountMap :: Germline + -> [FastaSequence] + -> CountMap +getCountMap (FastaSequence { fastaSeq = germlineSeq }) = + Map.unionsWith (+) . map initializeMaps + where + initializeMaps = Map.fromList + . map swap + . zip [1,1..] + . zip [1..] + . T.zip germlineSeq + . fastaSeq + diff --git a/src/TransformFastaList.hs b/src/TransformFastaList.hs new file mode 100644 --- /dev/null +++ b/src/TransformFastaList.hs @@ -0,0 +1,119 @@ +-- TransformFastaList module. +-- By Gregory W. Schwartz +-- +-- Collection of functions that transform a fasta sequence in some way + +{-# LANGUAGE BangPatterns, OverloadedStrings #-} + +module TransformFastaList ( convertToAminoAcidsFastaSequence + , replaceChars + , fillInSequence + , changeField + , changeAllFields + , getRegionSequence + , trimFasta + , removeUnknownNucs + ) where + +-- Built-in +import Data.List +import Data.Char +import qualified Data.Sequence as Seq +import qualified Data.Foldable as F +import qualified Data.Text as T + +-- Cabal +import Data.Fasta.Text +import Text.Regex.TDFA +import Text.Regex.TDFA.Text + +-- Local +import Types +import Utility + +-- | Convert sequences to amino acids +convertToAminoAcidsFastaSequence :: FastaSequence -> FastaSequence +convertToAminoAcidsFastaSequence = fromEither . translate 1 + where + fromEither (Right x) = x + fromEither (Left x) = error . T.unpack $ x + +-- | Fill in the sequence with corrected nucleotides or amino acids +fillInSequence :: Field -> Start -> Char -> FastaSequence -> FastaSequence +fillInSequence f s c fs = fs { fastaSeq = newFastaSeq } + where + newFastaSeq = (\xs -> if T.singleton c `T.isInfixOf` xs then "" else xs) + $ first `mappend` replaceChars c old new + new = (T.splitOn "|" . fastaHeader $ fs) !! (f - 1) + (first, old) = T.splitAt (s - 1) . fastaSeq $ fs + +-- | Change a field to a match, so a regex "ch.*_" to field 2 of +-- ">abc|brie_cheese_dude" would result in ">abc|cheese_". Useful for +-- getting specific properties from a field +changeField :: Maybe Field -> T.Text -> FastaSequence -> FastaSequence +changeField Nothing _ fs = fs +changeField (Just field) regex fs = fs { fastaHeader = newFastaHeader } + where + newFastaHeader = T.intercalate "|" + . F.toList + . Seq.update (field - 1) newField + $ splitField + newField = (Seq.index splitField (field - 1)) =~ regex :: T.Text + splitField = Seq.fromList . T.splitOn "|" . fastaHeader $ fs + +-- | Change all fields to their matches based on changeField +changeAllFields :: FastaSequence -> [(Maybe Int, T.Text)] -> FastaSequence +changeAllFields = foldl' (\fs (!x, !y) -> changeField x y fs) + +-- | Get a region of a text, 0 indexed +getRegion :: Maybe Start -> Maybe Stop -> T.Text -> T.Text +getRegion Nothing Nothing = id +getRegion Nothing (Just stop) = T.take stop +getRegion (Just start) Nothing = T.drop start +getRegion (Just start) (Just stop) = T.take (stop - start) . T.drop start + +-- | Get a region of a sequence, 1 indexed +getRegionSequence :: Maybe Start -> Maybe Stop -> FastaSequence -> FastaSequence +getRegionSequence start0 stop fs = fs { fastaSeq = newFastaSeq } + where + newFastaSeq = getRegion start stop . fastaSeq $ fs + start = fmap (flip (-) 1) start0 + +-- | Trim the sequence. For default, trim the extra nucleotides off the +-- end. Otherwise, use UCSC frames (0, 1, 2) to cut off the beginning (0 +-- means in frame) or the end (0 means the nucleotide AFTER the end of the +-- sequence is in frame). For amino acids, if it's not in frame, just cut +-- that amino acid. +trim :: GeneticUnit -> Maybe FrameType -> Maybe Frame -> T.Text -> T.Text +trim Nucleotide Nothing Nothing x = T.dropEnd + (T.length x `mod` 3) + x +trim _ _ Nothing x = x +trim _ _ (Just 0) x = x +trim Nucleotide (Just InFrame) (Just frame) x = T.drop (3 - frame) x +trim Nucleotide (Just OutFrame) (Just frame) x = T.dropEnd frame x +trim AminoAcid (Just InFrame) (Just _) x = T.drop 1 x +trim AminoAcid (Just OutFrame) (Just _) x = T.dropEnd 1 x + +-- | Trim off extra nucleotides (or amino acids) from a fasta sequence. If +-- inframe and outframe are specified, instead cut off based on those frames. +trimFasta :: GeneticUnit + -> Maybe Frame + -> Maybe Frame + -> FastaSequence + -> FastaSequence +trimFasta Nucleotide Nothing Nothing fs = + fs { fastaSeq = trim Nucleotide Nothing Nothing . fastaSeq $ fs } +trimFasta genU inF outF fs = fs { fastaSeq = trim genU (Just InFrame) inF + . trim genU (Just OutFrame) outF + . fastaSeq $ fs + } + +-- | Convert non standard nucleotides to gaps +removeUnknownNucs :: FastaSequence -> FastaSequence +removeUnknownNucs fs = fs { fastaSeq = T.map changeNuc . fastaSeq $ fs } + where + changeNuc x + | toUpper x `elem` ("ATCGN.-" :: String) = x + | otherwise = '-' + diff --git a/src/Types.hs b/src/Types.hs new file mode 100644 --- /dev/null +++ b/src/Types.hs @@ -0,0 +1,34 @@ +-- Types module. +-- By Gregory W. Schwartz +-- +-- Collects all application specific types. + +module Types where + +-- Built-in +import qualified Data.Text as T +import qualified Data.Map.Strict as Map + +-- Cabal +import Data.Fasta.Text.Types + +-- Algebraic +data GeneticUnit = AminoAcid | Nucleotide deriving (Read, Show) +data FrameType = InFrame | OutFrame deriving (Read, Show) + +-- Basic +type ID = Int +type Codon = T.Text +type CodonMut = Int +type Field = Int +type Start = Int +type Stop = Int +type Position = Int +type Frame = Int + +-- Advanced +type CloneEntry = (Germline, [FastaSequence]) +type FillInValue = (Field, Start, Char) +type Mutation = (Char, Char) +type CountMap = Map.Map (Position, Mutation) Int +type CodonMutations = [[(Position, Mutation)]] diff --git a/src/Utility.hs b/src/Utility.hs new file mode 100644 --- /dev/null +++ b/src/Utility.hs @@ -0,0 +1,65 @@ +-- Utility module +-- By Gregory W. Schwartz +-- +-- Collects utility functions for the main files + +{-# LANGUAGE OverloadedStrings, ViewPatterns #-} + +module Utility ( addLengthHeader + , addFillerGermlines + , replaceChars + , getField + ) where + +-- Built-in +import qualified Data.Map as M +import qualified Data.Text as T + +-- Cabal +import Data.Fasta.Text +import TextShow + +-- | Adds the length of a sequence to the header of that sequence +addLengthHeader :: FastaSequence -> FastaSequence +addLengthHeader fSeq = fSeq { fastaHeader = fastaHeader fSeq + `mappend` "|" + `mappend` (showt . T.length . fastaSeq $ fSeq) + } + +-- | Adds filler germlines to normal fasta files +addFillerGermlines :: [FastaSequence] -> CloneMap +addFillerGermlines = M.fromList . labelGermlines . map insertDummy + where + labelGermlines = map (\(x, (y, z)) -> ((x, y), z)) . zip [0..] + insertDummy x = (dummy, [x]) + dummy = FastaSequence {fastaHeader = "filler", fastaSeq = "---"} + +-- Like zipWith, but if one if one list is longer than the other than use +-- the remaining, needs to be the same type +zipWithRetain :: (a -> a -> a) -> [a] -> [a] -> [a] +zipWithRetain _ [] [] = [] +zipWithRetain _ xs [] = xs +zipWithRetain _ [] ys = ys +zipWithRetain f (x:xs) (y:ys) = f x y : zipWithRetain f xs ys + +-- Like zipWithRetain, but for text +zipWithRetainText :: (Char -> Char -> Char) -> T.Text -> T.Text -> T.Text +zipWithRetainText _ (T.uncons -> Nothing) (T.uncons -> Nothing) = T.empty +zipWithRetainText _ xs (T.uncons -> Nothing) = xs +zipWithRetainText _ (T.uncons -> Nothing) ys = ys +zipWithRetainText f (T.uncons -> Just (x, xs)) (T.uncons -> Just (y, ys)) + = f x y `T.cons` zipWithRetainText f xs ys + +-- Replace characters in the first string with another in the second string +-- if they are equal to a certain character and they aren't replaced with +-- a gap. +replaceChars :: Char -> T.Text -> T.Text -> T.Text +replaceChars c = zipWithRetainText changeChar + where + changeChar a b = if a == c && (not . T.isInfixOf (T.singleton b)) ".-" + then b + else a + +-- | Get the field of a fasta sequence, 1 indexed split by "|" +getField :: Int -> FastaSequence -> T.Text +getField f fs = (T.splitOn "|" . fastaHeader $ fs) !! (f - 1)