modify-fasta 0.8.0.2 → 0.8.0.3
raw patch · 12 files changed
+1564/−637 lines, 12 filesdep +modify-fasta
Dependencies added: modify-fasta
Files
- app/Main.hs +627/−0
- modify-fasta.cabal +26/−10
- src/Diversity.hs +57/−0
- src/FilterCloneList.hs +71/−0
- src/FilterCloneMap.hs +240/−0
- src/FilterFastaList.hs +83/−0
- src/Main.hs +0/−627
- src/Print.hs +97/−0
- src/TransformCloneList.hs +145/−0
- src/TransformFastaList.hs +119/−0
- src/Types.hs +34/−0
- src/Utility.hs +65/−0
+ app/Main.hs view
@@ -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" )
modify-fasta.cabal view
@@ -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
+ src/Diversity.hs view
@@ -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+
+ src/FilterCloneList.hs view
@@ -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
+ src/FilterCloneMap.hs view
@@ -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)
+ src/FilterFastaList.hs view
@@ -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
− src/Main.hs
@@ -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" )
+ src/Print.hs view
@@ -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
+ src/TransformCloneList.hs view
@@ -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+
+ src/TransformFastaList.hs view
@@ -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 = '-'+
+ src/Types.hs view
@@ -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)]]
+ src/Utility.hs view
@@ -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)