diff --git a/app/Main.hs b/app/Main.hs --- a/app/Main.hs +++ b/app/Main.hs @@ -10,6 +10,7 @@ import qualified Data.Map as M import qualified System.IO as IO import Control.Monad +import Control.Monad.Trans.Maybe -- Cabal import qualified Data.Text as T @@ -67,8 +68,10 @@ , inputMutType :: String , inputChangeField :: String , inputCustomFilter :: String - , customGermlineFlag :: Bool , customRemoveFlag :: Bool + , customGermlineFlag :: Bool + , inputInsertionPosition :: Int + , inputReference :: Maybe String , inputGeneAlleleField :: Int , countFlag :: Bool , output :: String @@ -93,7 +96,7 @@ ( long "legacy" <> short 'L' <> help "Whether to use the legacy version with no pipes. Note: The\ - \ legacy version supports more features but is greedy\ + \ legacy version will be unsupported and is greedy\ \ in terms of speed and memory. Use only if really needed.\ \ Features that are legacy only are noted in this\ \ documentation" ) @@ -304,17 +307,33 @@ \ 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" ) + <*> switch + ( long "legacy-custom-germline" + <> short 'G' + <> help "Whether to apply the custom filter to germlines (>>)\ + \ instead of sequences (>). LEGACY ONLY" ) <*> option auto + ( long "input-insertion-position" + <> metavar "[1]|INT" + <> value 1 + <> help "The field (1 indexed) of the where to insert the sequence\ + \ into the reference from --input-insertion-reference.\ + \ For instance, with a reference of \"ATT\", a position of 2,\ + \ and a sequence of \"GC\", the result will be \"AGCTT\".\ + \ Keeps the header of the insertion sequence. Requires\ + \ --input-insertion-reference. No legacy." ) + <*> optional ( strOption + ( long "input-reference" + <> metavar "FILE" + <> help "The file containing a reference sequence.\ + \ Only the first sequence is used." ) + ) + <*> option auto ( long "legacy-gene-allele-field" <> short 'V' <> metavar "[1]|INT" @@ -356,6 +375,16 @@ hOut <- if null . output $ opts then return IO.stdout else IO.openFile (output opts) IO.WriteMode + + -- | Get a possible reference sequence. + let getRefSeq :: MaybeT IO FastaSequence + getRefSeq = do + file <- MaybeT . return $ inputReference opts + hInRef <- lift $ IO.openFile file IO.ReadMode + MaybeT . runEffect . P.head $ pipesFasta (PT.fromHandle hInRef) + + refSeq <- runMaybeT getRefSeq + let genUnit = aminoAcidsFlag opts stopRange = inputStopRange opts customFilters = fieldIntParser . inputCustomFilter $ opts @@ -363,6 +392,7 @@ codonMut = inputCodonMut opts codonMutType = T.pack . inputCodonMutType $ opts mutType = T.pack . inputMutType $ opts + insertPos = inputInsertionPosition opts -- Remove out of frame sequences seqInFrame x = not ( removeOutOfFrameFlag opts @@ -416,7 +446,7 @@ includeMutations fs = case trackMutations opts of Nothing -> fs - (Just x) -> + (Just x) -> addMutationsHeader (translateTrackMutations opts) x fs -- Fill in bad characters at the requested section with possible @@ -425,6 +455,15 @@ (-1, -1, 'X') -> id (f, s, c) -> fillInSequence f s c + -- Insert a sequence into a reference sequence. + insert fs = case (refSeq, insertPos) of + (Just r, p) -> + insertSequence + (read . T.unpack . getField p '|' $ fs) + r + fs + _ -> fs + -- Find the complement complement = if complementFlag opts then compl @@ -480,6 +519,7 @@ . trim . noNs . fillIn + . insert . 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 @@ -684,9 +724,9 @@ <> progDesc "Modify fasta (and CLIP) files in several optional ways.\ \ Order of transformation goes: seqInFrame -> customFilter\ \ -> noStops -> removeHighMutations -> getMutations ->\ - \ getFrequentMutations -> cutSequence -> fillIn -> noNs\ - \ -> changeHeader -> complement -> reverseComplement ->\ - \ ntToaa -> includeLength,\ + \ getFrequentMutations -> cutSequence -> insert -> fillIn\ + \ -> noNs -> changeHeader -> complement -> reverseComplement\ + \ -> 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." diff --git a/modify-fasta.cabal b/modify-fasta.cabal --- a/modify-fasta.cabal +++ b/modify-fasta.cabal @@ -2,14 +2,14 @@ -- documentation, see http://haskell.org/cabal/users-guide/ name: modify-fasta -version: 0.8.2.3 +version: 0.8.3.0 synopsis: Modify fasta (and CLIP) files in several optional ways -- description: homepage: https://github.com/GregorySchwartz/modify-fasta license: GPL-3 license-file: LICENSE author: GregorySchwartz -maintainer: gregory.schwartz@drexel.edu +maintainer: gsch@mail.med.upenn.edu -- copyright: category: Bioinformatics build-type: Simple @@ -29,12 +29,12 @@ , Utility build-depends: base >=4.6 && <5 , containers >=0.5 - , text - , text-show - , split >=0.2 + , fasta , regex-tdfa >=1.2 , regex-tdfa-text - , fasta + , split >=0.2 + , text + , text-show executable modify-fasta hs-source-dirs: app @@ -42,13 +42,14 @@ build-depends: modify-fasta , base >=4.6 && <5 , containers >=0.5 - , mtl >=2.1 - , text - , split >=0.2 , fasta + , mtl >=2.1 + , optparse-applicative >=0.13 , pipes >= 4.1 , pipes-text - , optparse-applicative >=0.11 , semigroups + , split >=0.2 + , text + , transformers -- Directories containing source files. ghc-options: -O2 diff --git a/src/TransformFastaList.hs b/src/TransformFastaList.hs --- a/src/TransformFastaList.hs +++ b/src/TransformFastaList.hs @@ -8,6 +8,7 @@ module TransformFastaList ( convertToAminoAcidsFastaSequence , replaceChars , fillInSequence + , insertSequence , changeField , changeAllFields , getRegionSequence @@ -46,6 +47,14 @@ new = (T.splitOn "|" . fastaHeader $ fs) !! (f - 1) (first, old) = T.splitAt (s - 1) . fastaSeq $ fs +-- | Insert a sequence into a reference. +insertSequence :: Start -> FastaSequence -> FastaSequence -> FastaSequence +insertSequence s ref fs = fs { fastaSeq = newFastaSeq } + where + newFastaSeq = before `mappend` ins `mappend` after + ins = fastaSeq fs + (before, after) = T.splitAt (s - 1) . fastaSeq $ ref + -- | 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 @@ -115,4 +124,3 @@ changeNuc x | toUpper x `elem` ("ATCGN.-" :: String) = x | otherwise = '-' -