sequenceTools (empty) → 1.4.0.1
raw patch · 13 files changed
+1213/−0 lines, 13 filesdep +ansi-wl-pprintdep +basedep +bytestringsetup-changed
Dependencies added: ansi-wl-pprint, base, bytestring, foldl, hspec, lens-family, optparse-applicative, pipes, pipes-group, pipes-ordered-zip, pipes-safe, random, rio, sequence-formats, sequenceTools, split, vector
Files
- Changelog.md +13/−0
- LICENSE +30/−0
- README.md +63/−0
- Setup.hs +2/−0
- sequenceTools.cabal +101/−0
- src-executables/genoStats.hs +157/−0
- src-executables/pileupCaller.hs +321/−0
- src-executables/vcf2eigenstrat.hs +130/−0
- src/SequenceTools/PileupCaller.hs +150/−0
- src/SequenceTools/Utils.hs +26/−0
- test/SequenceTools/PileupCallerSpec.hs +187/−0
- test/SequenceTools/UtilsSpec.hs +32/−0
- test/Spec.hs +1/−0
+ Changelog.md view
@@ -0,0 +1,13 @@+V 1.2.3 : Adapted to newest sequence-formats. Had to change all the chromosome-related code to the newType Chrom datatype. Also started implementing normaliseBimWithVCF.++V 1.2.4: normaliseBimWithVCF is ready.++V 1.3.0: Lots of refactoring. Lots of testing. Removed some features in vcf2eigenstrat and in pileupCaller, including+ the option in pileupCaller to call without a SNP file.++V 1.3.1: Bumped dependency on sequence-formats to new sequence-formats-1.4.0, which includes strand-information in pileup data, as well as + rsIds in freqSum to output the correct rsId, and an option to parse chromosomes X, Y and MT.++V 1.4.0: Added single strand mode, and new triallelic treatment.++V 1.4.0.1: Improved README, fixed output bug in genoStats.hs
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright Author name here (c) 2015++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++ * Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++ * Redistributions in binary form must reproduce the above+ copyright notice, this list of conditions and the following+ disclaimer in the documentation and/or other materials provided+ with the distribution.++ * Neither the name of Author name here nor the names of other+ contributors may be used to endorse or promote products derived+ from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ README.md view
@@ -0,0 +1,63 @@+# SequenceTools++[](https://anaconda.org/bioconda/sequencetools)++(bioconda package available thanks to [apeltzer](https://github.com/apeltzer)!)++This repository contains some programs that I use for processing sequencing data.++# Installation++1. Download stack (https://docs.haskellstack.org/en/stable/README/#how-to-install<Paste>).+2. Clone the repository using `git clone https://github.com/stschiff/sequenceTools.git; cd sequenceTools`.+3. Run `stack install` inside the sequenceTools directory. If everything goes as planned, you should end up having the executables defined in this package under `~/.local/bin`.+4. Add `~/.local/bin` to your PATH, for example by adding to your `~/.profile` or `~/.bash_profile` the line `PATH=$PATH:$HOME/.local/bin`. Run `source ~/.profile` or `source ~/.bash_profile`, respectively, to update your path.+5. Run `pileupCaller --version`. It should output `1.4.0`. You're all set.++## pileupCaller++The main tool in this repository is the program `pileupCaller` to sample alleles from low coverage sequence data. The first step is to generate a “pileup” file at all positions you wish to genotype. To do that, here is a typical command line, which restricts to mapping and base quality of 30:++ samtools mpileup -R -B -q30 -Q30 -l <list_of_positions.txt> \+ -f <reference_genome.fasta> \+ Sample1.bam Sample2.bam Sample3.bam > pileup.txt++Important Note: You should definitely use the `-B` flag, which disables base alignment quality recalibration. This mechanism is turned on by default and causes huge reference bias with low coverage ancient DNA data. This flag disables the mechanism.++In the above command line, the file "list_of_positions.txt" should either contain positions (0-based) or a bed file (see samtools manual for details). The output is a simple text file with all positions that could be genotyped in the three samples.++Next, you need to run my tool pileupCaller, which you run like this:++ pileupCaller --randomHaploid --sampleNames Sample1,Sample2,Sample3 \+ --samplePopName MyPop -f <Eigenstrat.snp> \+ -e <My_output_prefix> < pileup.txt++Here, options `--sampleNames` gives the names of the samples that is output in the Eigenstrat `*.ind` file, and and `-–samplePopName` is optional to also give the population names in that file (defaults to `Unknown`, you can also change it later in the output). Then, option `-f` gives an Eigenstrat positions file. This is important because the pileup file only contains sites which could be called in at least one of your samples. In order to later merge your dataset with another Eigenstrat file, pileupCaller will check every position in the other Eigenstrat file to make sure every position is output with the correct alleles and missing genotypes if appropriate. Finally, the `-e` option specifies Eigenstrat as output format and gives the prefix for the `*.ind`, `*.pos` and `*.geno` files. Without the `-e` option, pileupCaller will output in FreqSum format, described [here](https://rarecoal-docs.readthedocs.io/en/latest/rarecoal-tools.html#vcf2freqsum), which is useful for debugging your pipeline, since it's just a single file that is output into the terminal and can therefore easily be inspected.++You can also get some help by typing `pileupCaller -h`, which shows a lot more option, for example the sampling method, minimal coverage and other important options.++Note that you can also fuse the two steps above into one unix pipe:++ samtools mpileup -R -B -q30 -Q30 -l <list_of_positions.txt> \+ -f <reference_genome.fasta> \+ Sample1.bam Sample2.bam Sample3.bam | \+ pileupCaller --randomHaploid --sampleNames Sample1,Sample2,Sample3 \+ --samplePopName MyPop -f <Eigenstrat.snp> \+ -e <My_output_prefix>++There is however an issue here: If you have aligned your read data to a version of the reference genome that uses `chr1`, `chr2` and so on as chromosome names, the resulting Eigenstrat file will be valid, but won't merge with other Eigenstrat datasets that use chromosome names `1`, `2` and so on. I would therefore recommend to strip the `chr` from your chromosome names if necessary. You can do that easily using a little UNIX filter using the `sed` tool. In the full pipeline, it looks like this:++ samtools mpileup -R -B -q30 -Q30 -l <list_of_positions.txt> \+ -f <reference_genome.fasta> \+ Sample1.bam Sample2.bam Sample3.bam | sed 's/chr//' | \+ pileupCaller --sampleNames Sample1,Sample2,Sample3 \+ --samplePopName MyPop -f <Eigenstrat.snp> \+ -o EigenStrat -e <My_output_prefix>++## vcf2eigenstrat++Simple tool to convert a VCF file to an Eigenstrat file. Pretty self-explanatory. Please run `vcf2eigenstrat --help` to output some documentation.++## genoStats++A simple tool to get some per-individual statistics from an Eigenstrat or Freqsum-file. Run `genoStats --help` for documentation.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ sequenceTools.cabal view
@@ -0,0 +1,101 @@+cabal-version: >=1.10+name: sequenceTools+version: 1.4.0.1+license: GPL-3+license-file: LICENSE+maintainer: stephan.schiffels@mac.com+author: Stephan Schiffels+synopsis: A package with tools for processing next generation sequencing data, in particular for processing data from ancient DNA sequencing libraries.+description:+ Key tool in this package is pileupCaller, a tool to randomly sample genotypes from sequencing data.+category: Bioinformatics+build-type: Simple+extra-source-files:+ README.md+ Changelog.md++library+ exposed-modules:+ SequenceTools.Utils+ SequenceTools.PileupCaller+ hs-source-dirs: src+ other-modules:+ Paths_sequenceTools+ default-language: Haskell2010+ build-depends:+ base >=4.7 && <5,+ optparse-applicative >=0.14.3.0,+ random >=1.1,+ sequence-formats >=1.4.0,+ bytestring >=0.10.8.2,+ vector >=0.12.0.3,+ pipes >=4.3.11++executable pileupCaller+ main-is: pileupCaller.hs+ hs-source-dirs: src-executables+ default-language: Haskell2010+ ghc-options: -threaded -rtsopts -with-rtsopts=-N2+ build-depends:+ base >=4.12.0.0,+ sequenceTools -any,+ sequence-formats >=1.4.0,+ optparse-applicative >=0.14.3.0,+ pipes >=4.3.11,+ rio >=0.1.11.0,+ vector >=0.12.0.3,+ random >=1.1,+ bytestring >=0.10.8.2,+ pipes-safe >=2.3.1,+ pipes-ordered-zip >=1.1.0,+ split >=0.2.3.3,+ ansi-wl-pprint >=0.6.9++executable vcf2eigenstrat+ main-is: vcf2eigenstrat.hs+ hs-source-dirs: src-executables+ default-language: Haskell2010+ build-depends:+ base >=4.12.0.0,+ sequenceTools -any,+ pipes-ordered-zip >=1.1.0,+ sequence-formats >=1.4.0,+ bytestring >=0.10.8.2,+ vector >=0.12.0.3,+ optparse-applicative >=0.14.3.0,+ pipes >=4.3.11,+ pipes-safe >=2.3.1++executable genoStats+ main-is: genoStats.hs+ hs-source-dirs: src-executables+ default-language: Haskell2010+ build-depends:+ base >=4.12.0.0,+ sequence-formats >=1.4.0,+ sequenceTools -any,+ foldl >=1.4.5,+ bytestring >=0.10.8.2,+ vector >=0.12.0.3,+ lens-family >=1.2.3,+ optparse-applicative >=0.14.3.0,+ pipes >=4.3.11,+ pipes-group >=1.0.12,+ pipes-safe >=2.3.1++test-suite sequenceToolsTests+ type: exitcode-stdio-1.0+ main-is: Spec.hs+ hs-source-dirs: test+ other-modules:+ SequenceTools.UtilsSpec+ SequenceTools.PileupCallerSpec+ default-language: Haskell2010+ build-depends:+ base >=4.12.0.0,+ hspec >=2.7.1,+ sequenceTools -any,+ sequence-formats >=1.4.0,+ vector >=0.12.0.3,+ bytestring >=0.10.8.2,+ pipes >=4.3.11
+ src-executables/genoStats.hs view
@@ -0,0 +1,157 @@+{-# LANGUAGE OverloadedStrings #-}++import SequenceFormats.FreqSum (readFreqSumFile, readFreqSumStdIn, FreqSumHeader(..), + FreqSumEntry(..))+import SequenceFormats.Eigenstrat (readEigenstrat, GenoEntry(..), GenoLine, EigenstratSnpEntry(..), EigenstratIndEntry(..))+import SequenceFormats.Utils (Chrom(..))+import SequenceTools.Utils (versionInfoOpt, versionInfoText)++import Control.Applicative ((<|>))+import Control.Foldl (purely, Fold(..))+import Control.Monad (forM_)+import Control.Monad.IO.Class (MonadIO, liftIO)+import qualified Data.ByteString.Char8 as B+import qualified Data.Vector as V+import Lens.Family2 (view)+import qualified Options.Applicative as OP+import Pipes (for, Producer, (>->), yield, Consumer, cat)+import Pipes.Group (groupsBy, folds)+import Pipes.Safe (MonadSafe, runSafeT)+import qualified Pipes.Prelude as P+import System.IO (hPutStrLn, stderr)+import Text.Printf (printf)++data ProgOpt = ProgOpt InputOption++data InputOption = FreqsumInput (Maybe FilePath) | EigenstratInput FilePath FilePath FilePath++data InputEntry = InputEntry Chrom GenoLine deriving (Show)++type StatsReportAllSamples = [StatsReport]++data StatsReport = StatsReport {+ srNrSitesMissing :: Int,+ srNrSitesHomRef :: Int,+ srNrSitesHomAlt :: Int,+ srNrSitesHet :: Int+} deriving (Show)++main :: IO ()+main = OP.execParser optionSpec >>= runWithOpts++optionSpec :: OP.ParserInfo ProgOpt+optionSpec = OP.info (pure (.) <*> versionInfoOpt <*> OP.helper <*> argParser) (+ OP.progDesc ("A program \+ \to evaluate per-chromosome and total statistics \+ \of genotyping data, read either as Eigenstrat (by specifying options -g, -s and -i) or FreqSum (default, or by specifying option -f). " ++ versionInfoText))++argParser :: OP.Parser ProgOpt+argParser = ProgOpt <$> (parseFreqsumInput <|> parseEigenstratInput)++parseFreqsumInput :: OP.Parser InputOption+parseFreqsumInput =+ process <$> OP.strOption (OP.long "freqsum" <> OP.short 'f' <> OP.help "a freqsum file to read as input. Use - to read from stdin (the default)" <>+ OP.value "-" <> OP.showDefault <> OP.metavar "FILEPATH")+ where+ process p =+ if p == "-"+ then FreqsumInput Nothing+ else FreqsumInput (Just p)++parseEigenstratInput :: OP.Parser InputOption+parseEigenstratInput = EigenstratInput <$> parseGenoFile <*> parseSnpFile <*> parseIndFile+ where+ parseGenoFile = OP.strOption (OP.long "eigenstratGeno" <> OP.short 'g' <> OP.help "Eigenstrat Geno File" <> OP.metavar "FILEPATH")+ parseSnpFile = OP.strOption (OP.long "eigenstratSnp" <> OP.short 's' <> OP.help "Eigenstrat Snp File" <> OP.metavar "FILEPATH")+ parseIndFile = OP.strOption (OP.long "eigenstratInd" <> OP.short 'i' <> OP.help "Eigenstrat Ind File" <> OP.metavar "FILEPATH")++runWithOpts :: ProgOpt -> IO ()+runWithOpts (ProgOpt inputOpt) = runSafeT $ do+ (names, entryProducer) <- case inputOpt of + FreqsumInput fsFile -> runWithFreqSum fsFile+ EigenstratInput genoFile snpFile indFile -> runWithEigenstrat genoFile snpFile indFile+ liftIO $ hPutStrLn stderr ("processing samples: " <> show names)+ let p = runStats names entryProducer >-> P.tee (reportStats names)+ totalReport <- purely P.fold (accumulateAllChromStats names) p+ printReports names totalReport++runWithFreqSum :: (MonadSafe m) => Maybe FilePath -> m ([String], Producer InputEntry m ())+runWithFreqSum fsFile = do+ (FreqSumHeader names nrHaps, fsProd) <- case fsFile of+ Nothing -> readFreqSumStdIn+ Just fn -> readFreqSumFile fn+ let prod = for fsProd $ \(FreqSumEntry chrom _ _ _ _ counts) -> do+ let genotypes = V.fromList $ do+ (count', nrHap) <- zip counts nrHaps+ case count' of+ Just 0 -> return HomRef+ Just x | x == nrHap -> return HomAlt+ | x > 0 && x < nrHap -> return Het+ Nothing -> return Missing+ _ -> error "should not happen"+ yield $ InputEntry chrom genotypes+ return (map B.unpack names, prod)++runWithEigenstrat :: (MonadSafe m) =>+ FilePath -> FilePath -> FilePath -> m ([String], Producer InputEntry m ())+runWithEigenstrat genoFile snpFile indFile = do+ (indEntries, genoStream) <- readEigenstrat genoFile snpFile indFile+ let names = [name | EigenstratIndEntry name _ _ <- indEntries]+ let prod = genoStream >-> P.map (\(EigenstratSnpEntry chrom _ _ _ _ _, genoLine) -> InputEntry chrom genoLine)+ return (names, prod)++runStats :: (MonadIO m) => [String] -> Producer InputEntry m () ->+ Producer (Chrom, StatsReportAllSamples) m ()+runStats names entryProducer =+ let groupedProd = view (groupsBy (\(InputEntry c1 _) (InputEntry c2 _) -> c1 == c2)) + entryProducer+ in purely folds (runStatsPerChrom (length names)) groupedProd++runStatsPerChrom :: Int -> Fold InputEntry (Chrom, StatsReportAllSamples)+runStatsPerChrom nrSamples = (,) <$> getChrom <*>+ traverse runStatsPerChromPerSample [0..(nrSamples - 1)] + +getChrom :: Fold InputEntry Chrom+getChrom = Fold (\_ (InputEntry c _) -> c) (Chrom "") id++runStatsPerChromPerSample :: Int -> Fold InputEntry StatsReport+runStatsPerChromPerSample i = Fold step initial extract+ where+ step :: StatsReport -> InputEntry -> StatsReport+ step accum@(StatsReport miss homr homa het) (InputEntry _ line) = case (line V.! i) of+ Missing -> accum {srNrSitesMissing = miss + 1}+ HomRef -> accum {srNrSitesHomRef = homr + 1}+ HomAlt -> accum {srNrSitesHomAlt = homa + 1}+ Het -> accum {srNrSitesHet = het + 1}+ initial :: StatsReport+ initial = StatsReport 0 0 0 0+ extract :: StatsReport -> StatsReport+ extract = id+ +reportStats :: (MonadIO m) => [String] -> Consumer (Chrom, StatsReportAllSamples) m ()+reportStats names = do+ liftIO . putStrLn $ "Chrom\tSample\tMissing\tHomRef\tHomAlt\tHet"+ for cat $ \(chrom, reports) -> printReports names (chrom, reports)++printReports :: (MonadIO m) => [String] -> (Chrom, StatsReportAllSamples) -> m ()+printReports names (chrom, reports) =+ forM_ (zip names reports) $ \(n, StatsReport mis ref alt het) -> do+ let total = mis + ref + alt + het+ misPerc = round $ (fromIntegral mis / fromIntegral total) * (100.0 :: Double) :: Int+ liftIO . putStrLn $ printf "%s\t%s\t%d(%d%%)\t%d\t%d\t%d" (unChrom chrom) n mis misPerc ref alt het++accumulateAllChromStats :: [String] ->+ Fold (Chrom, StatsReportAllSamples) (Chrom, StatsReportAllSamples)+accumulateAllChromStats names = Fold step initial extract+ where+ step :: StatsReportAllSamples -> (Chrom, StatsReportAllSamples) -> StatsReportAllSamples+ step sumReports (_, newReports) = do+ (StatsReport smiss shomr shoma shet, StatsReport miss homr homa het) <-+ zip sumReports newReports+ return $ StatsReport (smiss + miss) (shomr + homr) (shoma + homa) (shet + het)+ initial :: StatsReportAllSamples+ initial = [StatsReport 0 0 0 0 | _ <- names]+ extract :: StatsReportAllSamples -> (Chrom, StatsReportAllSamples)+ extract r = (Chrom "Total", r)++
+ src-executables/pileupCaller.hs view
@@ -0,0 +1,321 @@+{-# LANGUAGE OverloadedStrings, NoImplicitPrelude #-}++import SequenceFormats.Eigenstrat (readEigenstratSnpFile, EigenstratSnpEntry(..), + EigenstratIndEntry(..), Sex(..), writeEigenstrat)+import SequenceFormats.FreqSum(FreqSumEntry(..), printFreqSumStdOut, FreqSumHeader(..))+import SequenceFormats.Pileup (PileupRow(..), readPileupFromStdIn)+import SequenceTools.Utils (versionInfoOpt, versionInfoText)+import SequenceTools.PileupCaller (CallingMode(..), callGenotypeFromPileup, callToDosage,+ freqSumToEigenstrat, filterTransitions, TransitionsMode(..), cleanSSdamageAllSamples)++import qualified Data.ByteString.Char8 as B+import Data.List.Split (splitOn)+import qualified Data.Vector.Unboxed.Mutable as V+import qualified Options.Applicative as OP+import Pipes (yield, (>->), runEffect, Producer, for)+import Pipes.OrderedZip (orderedZip, orderCheckPipe)+import qualified Pipes.Prelude as P+import Pipes.Safe (runSafeT, SafeT)+import RIO+import System.IO (readFile, hPutStrLn, stderr)+import System.Random (mkStdGen, setStdGen)+import Text.Printf (printf)+import qualified Text.PrettyPrint.ANSI.Leijen as PP++data OutFormat = EigenstratFormat String String | FreqSumFormat++data ProgOpt = ProgOpt CallingMode Bool (Maybe Int) Int TransitionsMode FilePath OutFormat (Either [String] FilePath)+ --optCallingMode :: CallingMode,+ --optKeepInCongruentReads :: Bool,+ --optSeed :: Maybe Int,+ --optMinDepth :: Int,+ --optTransitionsMode :: TransitionsMode,+ --optSnpFile :: FilePath,+ --optOutFormat :: OutFormat,+ --optSampleNames :: Either [String] FilePath++data ReadStats = ReadStats {+ rsTotalSites :: IORef Int,+ rsNonMissingSites :: V.IOVector Int,+ rsRawReads :: V.IOVector Int,+ rsReadsCleanedSS :: V.IOVector Int,+ rsReadsCongruent :: V.IOVector Int+}++data Env = Env {+ envCallingMode :: CallingMode,+ envKeepInCongruentReads :: Bool,+ envMinDepth :: Int,+ envTransitionsMode :: TransitionsMode,+ envOutFormat :: OutFormat,+ envSnpFile :: FilePath,+ envSampleNames :: [String],+ envStats :: ReadStats+}++type App = ReaderT Env (SafeT IO)++main :: IO ()+main = do+ args <- OP.execParser parserInfo+ env <- initialiseEnvironment args+ runSafeT $ runReaderT runMain env++parserInfo :: OP.ParserInfo ProgOpt+parserInfo = OP.info (pure (.) <*> versionInfoOpt <*> OP.helper <*> argParser)+ (OP.progDescDoc (Just programHelpDoc))++argParser :: OP.Parser ProgOpt+argParser = ProgOpt <$> parseCallingMode <*> parseKeepIncongruentReads <*> parseSeed <*> parseMinDepth <*>+ parseTransitionsMode <*> parseSnpFile <*> parseFormat <*> parseSampleNames+ where+ parseCallingMode = parseRandomCalling <|> parseMajorityCalling <|> parseRandomDiploidCalling+ parseRandomCalling = OP.flag' RandomCalling (OP.long "randomHaploid" <>+ OP.help "This method samples one read at random at each site, and uses the allele \+ \on that read as the one for the actual genotype. This results in a haploid \+ \call")+ parseRandomDiploidCalling = OP.flag' RandomDiploidCalling (OP.long "randomDiploid" <>+ OP.help "Sample two reads at random (without replacement) at each site and represent the \+ \individual by a diploid genotype constructed from those two random \+ \picks. This will always assign missing data to positions where only \+ \one read is present, even if minDepth=1. The main use case for this \+ \option is for estimating mean heterozygosity across sites.")+ parseMajorityCalling = MajorityCalling <$> (parseMajorityCallingFlag *> parseDownsamplingFlag)+ parseMajorityCallingFlag = OP.flag' True (OP.long "majorityCall" <> OP.help+ "Pick the allele supported by the \+ \most reads at a site. If an equal numbers of alleles fulfil this, pick one at \+ \random. This results in a haploid call. See --downSampling for best practices \+ \for calling rare variants")+ parseDownsamplingFlag = OP.switch (OP.long "downSampling" <> OP.help "When this switch is given, \+ \the MajorityCalling mode will downsample \+ \from the total number of reads a number of reads \+ \(without replacement) equal to the --minDepth given. This mitigates \+ \reference bias in the MajorityCalling model, which increases with higher coverage. \+ \The recommendation for rare-allele calling is --majorityCall --downsampling --minDepth 3")+ parseKeepIncongruentReads = OP.switch (OP.long "keepIncongruentReads" <> OP.help "By default, \+ \pileupCaller now removes reads with tri-allelic alleles that are neither of the two alleles specified in the SNP file. \+ \To keep those reads for sampling, set this flag. With this option given, if \+ \the sampled read has a tri-allelic allele that is neither of the two given alleles in the SNP file, a missing genotype is generated. \+ \IMPORTANT NOTE: The default behaviour has changed in pileupCaller version 1.4.0. If you want to emulate the previous \+ \behaviour, use this flag. I recommend now to NOT set this flag and use the new behaviour.")+ parseSeed = OP.option (Just <$> OP.auto) (OP.long "seed" <>+ OP.value Nothing <> OP.metavar "<RANDOM_SEED>" <>+ OP.help "random seed used for the random number generator. If not given, use \+ \system clock to seed the random number generator.")+ parseMinDepth = OP.option OP.auto (OP.long "minDepth" <> OP.short 'd' <>+ OP.value 1 <> OP.showDefault <> OP.metavar "<DEPTH>" <>+ OP.help "specify the minimum depth for a call. For sites with fewer \+ \reads than this number, declare Missing")+ parseTransitionsMode = parseSkipTransitions <|> parseTransitionsMissing <|> parseSingleStrandMode <|> pure AllSites+ parseSkipTransitions = OP.flag' SkipTransitions (OP.long "skipTransitions" <>+ OP.help "skip transition SNPs entirely in the output, resulting in a dataset with fewer sites.")+ parseTransitionsMissing = OP.flag' TransitionsMissing (OP.long "transitionsMissing" <>+ OP.help "mark transitions as missing in the output, but do output the sites.")+ parseSingleStrandMode = OP.flag' SingleStrandMode (OP.long "singleStrandMode" <>+ OP.help "At C/T polymorphisms, ignore reads aligning to the forward strand. \+ \At G/A polymorphisms, ignore reads aligning to the reverse strand. This should \+ \remove post-mortem damage in ancient DNA libraries prepared with the non-UDG single-stranded protocol.")+ parseSnpFile = OP.strOption (OP.long "snpFile" <> OP.short 'f' <>+ OP.metavar "<FILE>" <> OP.help "an Eigenstrat-formatted SNP list file for \+ \the positions and alleles to call. All \+ \positions in the SNP file will be output, adding missing data where \+ \there is no data. Note that pileupCaller automatically checks whether \+ \alleles in the SNP file are flipped with respect to the human \+ \reference, and in those cases flips the genotypes accordingly. \+ \But it assumes that the strand-orientation of the SNPs given in the SNP list is the one \+ \in the reference genome used in the BAM file underlying the pileup input. \+ \Note that both the SNP file and the incoming pileup data have to be \+ \ordered by chromosome and position, and this is checked. The chromosome order in humans is 1-22,X,Y,MT. \+ \Chromosome can generally begin with \"chr\". In case of non-human data with different chromosome \+ \names, you should convert all names to numbers. They will always considered to \+ \be numerically ordered, even beyond 22. Finally, I note that for internally, \+ \X is converted to 23, Y to 24 and MT to 90. This is the most widely used encoding in Eigenstrat \+ \databases for human data, so using a SNP file with that encoding will automatically be correctly aligned \+ \to pileup data with actual chromosome names X, Y and MT (or chrX, chrY and chrMT, respectively).")+ parseFormat = (EigenstratFormat <$> parseEigenstratPrefix <*> parseSamplePopName) <|> pure FreqSumFormat+ parseEigenstratPrefix = OP.strOption (OP.long "eigenstratOut" <> OP.short 'e' <>+ OP.metavar "<FILE_PREFIX>" <>+ OP.help "Set Eigenstrat as output format. Specify the filenames for the EigenStrat \+ \SNP and IND file outputs: <FILE_PREFIX>.snp.txt and <FILE_PREFIX>.ind.txt \+ \If not set, output will be FreqSum (Default). Note that freqSum format, described at \+ \https://rarecoal-docs.readthedocs.io/en/latest/rarecoal-tools.html#vcf2freqsum, \+ \is useful for testing your pipeline, since it's output to standard out")+ parseSampleNames = parseSampleNameList <|> parseSampleNameFile+ parseSampleNameList = OP.option (Left . splitOn "," <$> OP.str)+ (OP.long "sampleNames" <> OP.metavar "NAME1,NAME2,..." <>+ OP.help "give the names of the samples as comma-separated list (no spaces)")+ parseSampleNameFile = OP.option (Right <$> OP.str) (OP.long "sampleNameFile" <> OP.metavar "<FILE>" <>+ OP.help "give the names of the samples in a file with one name per \+ \line")+ parseSamplePopName = OP.strOption (OP.long "samplePopName" <> OP.value "Unknown" <> OP.showDefault <>+ OP.metavar "POP" <>+ OP.help "specify the population name of the samples, which is included\+ \ in the output *.ind.txt file in Eigenstrat output. This will be ignored if the output \+ \format is not Eigenstrat")++programHelpDoc :: PP.Doc+programHelpDoc =+ part1 PP.<$>+ PP.enclose PP.line PP.line (PP.indent 4 samtoolsExample) PP.<$>+ part2 PP.<$> PP.line PP.<$>+ PP.string versionInfoText+ where+ part1 = PP.fillSep . map PP.text . words $+ "PileupCaller is a tool to create genotype calls from bam files using read-sampling methods. \+ \To use this tool, you need to convert bam files into the mpileup-format, specified at \+ \http://www.htslib.org/doc/samtools.html (under \"mpileup\"). The recommended command line \+ \to create a multi-sample mpileup file to be processed with pileupCaller is"+ samtoolsExample = PP.hang 4 . PP.fillSep . map PP.text . words $+ "samtools mpileup -B -q30 -Q30 -l <BED_FILE> -R -f <FASTA_REFERENCE_FILE> \+ \Sample1.bam Sample2.bam Sample3.bam | pileupCaller ..."+ part2 = PP.fillSep . map PP.text . words $+ "You can lookup what these options do in the samtools documentation. \+ \Note that flag -B in samtools is very important to reduce reference \+ \bias in low coverage data."++initialiseEnvironment :: ProgOpt -> IO Env+initialiseEnvironment args = do+ let (ProgOpt callingMode keepInCongruentReads seed minDepth+ transitionsMode snpFile outFormat sampleNames) = args+ case seed of+ Nothing -> return ()+ Just seed_ -> liftIO . setStdGen $ mkStdGen seed_+ sampleNamesList <- case sampleNames of+ Left list -> return list+ Right fn -> lines <$> readFile fn+ let n = length sampleNamesList+ readStats <- ReadStats <$> newIORef 0 <*> makeVec n <*> makeVec n <*> makeVec n <*> makeVec n+ return $ Env callingMode keepInCongruentReads minDepth transitionsMode+ outFormat snpFile sampleNamesList readStats+ where+ makeVec n = do+ v <- V.new n + V.set v 0+ return v++runMain :: App ()+runMain = do+ let pileupProducer = readPileupFromStdIn+ snpFile <- asks envSnpFile+ freqSumProducer <- pileupToFreqSum snpFile pileupProducer+ outFormat <- asks envOutFormat+ case outFormat of+ FreqSumFormat -> outputFreqSum freqSumProducer+ EigenstratFormat outPrefix popName -> outputEigenStrat outPrefix popName freqSumProducer+ outputStats++pileupToFreqSum :: FilePath -> Producer PileupRow (SafeT IO) () ->+ App (Producer FreqSumEntry (SafeT IO) ())+pileupToFreqSum snpFileName pileupProducer = do+ nrSamples <- length <$> asks envSampleNames+ let snpProdOrderChecked =+ readEigenstratSnpFile snpFileName >-> orderCheckPipe cmpSnpPos+ pileupProdOrderChecked =+ pileupProducer >-> orderCheckPipe cmpPileupPos+ jointProd =+ orderedZip cmpSnpToPileupPos snpProdOrderChecked pileupProdOrderChecked+ mode <- asks envCallingMode+ keepInCongruentReads <- asks envKeepInCongruentReads+ transisitionsMode <- asks envTransitionsMode+ let singleStrandMode = (transisitionsMode == SingleStrandMode)+ minDepth <- asks envMinDepth+ readStats <- asks envStats+ let ret = Pipes.for jointProd $ \jointEntry ->+ case jointEntry of+ (Just esEntry, Nothing) -> do+ let (EigenstratSnpEntry chr pos _ id_ ref alt) = esEntry+ dosages = (replicate nrSamples Nothing) + liftIO $ addOneSite readStats+ yield $ FreqSumEntry chr pos (Just . B.unpack $ id_) ref alt dosages+ (Just esEntry, Just pRow) -> do+ let (EigenstratSnpEntry chr pos _ id_ ref alt) = esEntry+ (PileupRow _ _ _ rawPileupBasesPerSample rawStrandInfoPerSample) = pRow+ let cleanBasesPerSample =+ if singleStrandMode+ then cleanSSdamageAllSamples ref alt rawPileupBasesPerSample rawStrandInfoPerSample+ else rawPileupBasesPerSample+ let congruentBasesPerSample = + if keepInCongruentReads+ then cleanBasesPerSample+ else map (filter (\c -> c == ref || c == alt)) cleanBasesPerSample+ liftIO $ addOneSite readStats+ liftIO $ updateStatsAllSamples readStats (map length rawPileupBasesPerSample)+ (map length cleanBasesPerSample) (map length congruentBasesPerSample)+ calls <- liftIO $ mapM (callGenotypeFromPileup mode minDepth) congruentBasesPerSample+ let genotypes = map (callToDosage ref alt) calls+ yield (FreqSumEntry chr pos (Just . B.unpack $ id_) ref alt genotypes)+ _ -> return ()+ return (fst <$> ret)+ where+ cmpSnpPos :: EigenstratSnpEntry -> EigenstratSnpEntry -> Ordering+ cmpSnpPos es1 es2 = (snpChrom es1, snpPos es1) `compare` (snpChrom es2, snpPos es2)+ cmpPileupPos :: PileupRow -> PileupRow -> Ordering+ cmpPileupPos pr1 pr2 = (pileupChrom pr1, pileupPos pr1) `compare` (pileupChrom pr2, pileupPos pr2)+ cmpSnpToPileupPos :: EigenstratSnpEntry -> PileupRow -> Ordering+ cmpSnpToPileupPos es pr = (snpChrom es, snpPos es) `compare` (pileupChrom pr, pileupPos pr)++addOneSite :: ReadStats -> IO ()+addOneSite readStats = modifyIORef' (rsTotalSites readStats) (+1)++updateStatsAllSamples :: ReadStats -> [Int] -> [Int] -> [Int] -> IO ()+updateStatsAllSamples readStats rawBaseCounts damageCleanedBaseCounts congruencyCleanedBaseCounts = do+ sequence_ [V.modify (rsRawReads readStats) (+n) i | (i, n) <- zip [0..] rawBaseCounts]+ sequence_ [V.modify (rsReadsCleanedSS readStats) (+n) i | (i, n) <- zip [0..] damageCleanedBaseCounts]+ sequence_ [V.modify (rsReadsCongruent readStats) (+n) i | (i, n) <- zip [0..] congruencyCleanedBaseCounts]+ let nonMissingSites = [if n > 0 then 1 else 0 | n <- congruencyCleanedBaseCounts]+ sequence_ [V.modify (rsNonMissingSites readStats) (+n) i | (i, n) <- zip [0..] nonMissingSites]++outputFreqSum :: Producer FreqSumEntry (SafeT IO) () -> App ()+outputFreqSum freqSumProducer = do+ transitionsOnly <- asks envTransitionsMode+ callingMode <- asks envCallingMode+ sampleNames <- asks envSampleNames+ let nrHaplotypes = case callingMode of+ MajorityCalling _ -> 1 :: Int+ RandomCalling -> 1+ RandomDiploidCalling -> 2+ let header' = FreqSumHeader (map B.pack sampleNames) [nrHaplotypes | _ <- sampleNames]+ outProd = freqSumProducer >-> filterTransitions transitionsOnly+ lift . runEffect $ outProd >-> printFreqSumStdOut header'++outputEigenStrat :: FilePath -> String -> Producer FreqSumEntry (SafeT IO) () -> App ()+outputEigenStrat outPrefix popName freqSumProducer = do+ transitionsMode <- asks envTransitionsMode+ sampleNames <- asks envSampleNames+ callingMode <- asks envCallingMode+ let diploidizeCall = case callingMode of+ RandomCalling -> True+ MajorityCalling _ -> True+ RandomDiploidCalling -> False+ let snpOut = outPrefix <> "snp.txt"+ indOut = outPrefix <> "ind.txt"+ genoOut = outPrefix <> "geno.txt"+ let indEntries = [EigenstratIndEntry n Unknown popName | n <- sampleNames]+ lift . runEffect $ freqSumProducer >-> filterTransitions transitionsMode >->+ P.map (freqSumToEigenstrat diploidizeCall) >->+ writeEigenstrat genoOut snpOut indOut indEntries++outputStats :: App ()+outputStats = do+ ReadStats totalSites nonMissingSitesVec rawReadsVec damageCleanedReadsVec congruentReadsVec <- asks envStats+ sampleNames <- asks envSampleNames+ liftIO $ hPutStrLn stderr+ "# Summary Statistics per sample \n\+ \# SampleName: Name of the sample as given by the user \n\+ \# TotalSites: Total number of sites in the given Snp file (before transition filtering) \n\+ \# NonMissingCalls: Total number of sites output with a non-Missing call (before transition filtering) \n\+ \# avgRawReads: mean coverage of raw pileup input data across total sites (incl. missing sites) \n\+ \# avgDamageCleanedReads: mean coverage of pileup after single-stranded damage removal \n\+ \# avgSampledFrom: mean coverage of pileup after removing reads with tri-allelic alleles \n\+ \SampleName\tTotalSites\tNonMissingCalls\tavgRawReads\tavgDamageCleanedReads\tavgSampledFrom"+ forM_ (zip [0..] sampleNames) $ \(i, name) -> do+ totalS <- readIORef totalSites+ nonMissingSites <- V.read nonMissingSitesVec i+ rawReads <- V.read rawReadsVec i+ damageCleanedReads <- V.read damageCleanedReadsVec i+ congruentReads <- V.read congruentReadsVec i+ let avgRawReads = (fromIntegral rawReads / fromIntegral totalS) :: Double+ avgDamageCleanedReads = (fromIntegral damageCleanedReads / fromIntegral totalS) :: Double+ avgCongruentReads = (fromIntegral congruentReads / fromIntegral totalS) :: Double+ liftIO . hPutStrLn stderr $ printf "%s\t%d\t%d\t%g\t%g\t%g" name totalS nonMissingSites+ avgRawReads avgDamageCleanedReads avgCongruentReads
+ src-executables/vcf2eigenstrat.hs view
@@ -0,0 +1,130 @@+{-# LANGUAGE OverloadedStrings #-}++import Pipes.OrderedZip (orderedZip)+import SequenceFormats.VCF (readVCFfromStdIn, VCFheader(..), VCFentry(..),+ isBiallelicSnp, getDosages, vcfToFreqSumEntry)+import SequenceFormats.Eigenstrat (EigenstratSnpEntry(..), readEigenstratSnpFile, writeEigenstrat, + GenoLine, GenoEntry(..), Sex(..), EigenstratIndEntry(..))+import SequenceFormats.FreqSum (FreqSumEntry(..), freqSumEntryToText)+import SequenceFormats.Utils (Chrom(..))++import SequenceTools.Utils (versionInfoText, versionInfoOpt)++import Control.Exception.Base (throwIO, AssertionFailed(..))+import Control.Monad (when)+import Control.Monad.IO.Class (liftIO, MonadIO)+import qualified Data.ByteString.Char8 as B+import Data.Monoid ((<>))+-- import Debug.Trace (trace)+import Data.Vector (fromList)+import qualified Options.Applicative as OP+import Pipes (Pipe, yield, (>->), runEffect, Producer, Pipe, for, cat)+import qualified Pipes.Prelude as P+import Pipes.Safe (runSafeT, MonadSafe)++-- snpPosFile outPrefix+data ProgOpt = ProgOpt (Maybe FilePath) FilePath++main :: IO ()+main = readOptions >>= runMain++readOptions :: IO ProgOpt+readOptions = OP.execParser parserInfo+ where+ parserInfo = OP.info (pure (.) <*> versionInfoOpt <*> OP.helper <*> argParser)+ (OP.progDesc ("A program to convert a VCF file (stdin) to Eigenstrat. " ++ versionInfoText))++argParser :: OP.Parser ProgOpt+argParser = ProgOpt <$> parseSnpPosFile <*> parseOutPrefix+ where+ parseSnpPosFile = OP.option (Just <$> OP.str)+ (OP.long "snpFile" <> OP.short 'f' <> OP.value Nothing <> OP.metavar "<FILE>" <>+ OP.help "specify an Eigenstrat SNP file with the positions and alleles of a \+ \reference set. If this option is given, only positions that are both in the SNP file \+ \and in the VCF will be output. Without this option, all sites in the VCF will be output. \+ \WARNING: Sites that are not in the VCF will not be output, and this is new behaviour. \+ \Previously one could specify that they will be output as missing or hom-ref, but that \+ \feature was recently removed. I plan to implement this behaviour in the future in a new eigenstrat-merging tool.")+ parseOutPrefix = OP.strOption (OP.long "outPrefix" <> OP.short 'e' <>+ OP.metavar "<FILE_PREFIX>" <>+ OP.help "specify the filenames for the EigenStrat SNP and IND \+ \file outputs: <FILE_PREFIX>.snp.txt and <FILE_PREFIX>.ind.txt")++runMain :: ProgOpt -> IO ()+runMain (ProgOpt maybeSnpPosFile outPrefix) =+ runSafeT $ do+ (vcfHeader, vcfBody) <- readVCFfromStdIn+ let snpOut = outPrefix ++ ".snp.txt"+ indOut = outPrefix ++ ".ind.txt"+ genoOut = outPrefix ++ ".geno.txt"+ VCFheader _ sampleNames = vcfHeader+ nrInds = length sampleNames+ indEntries = [EigenstratIndEntry n Unknown "Unknown" | n <- sampleNames]+ let vcfBodyBiAllelic = vcfBody >-> P.filter (\e -> isBiallelicSnp (vcfRef e) (vcfAlt e))+ vcfProducer <- case maybeSnpPosFile of+ Just snpPosFile ->+ return $ runJointly vcfBodyBiAllelic nrInds snpPosFile+ Nothing -> return $ runSimple vcfBodyBiAllelic+ runEffect $ vcfProducer >-> eigenStratPipe >-> writeEigenstrat genoOut snpOut indOut indEntries++runJointly :: (MonadIO m, MonadSafe m) => Producer VCFentry m r -> Int -> FilePath -> Producer FreqSumEntry m r+runJointly vcfBody nrInds snpPosFile =+ let snpProd = readEigenstratSnpFile snpPosFile+ jointProd = snd <$> orderedZip cmp snpProd vcfBody+ in jointProd >-> processVcfWithSnpFile nrInds+ where+ cmp (EigenstratSnpEntry snpChrom' snpPos' _ _ _ _) vcfEntry = (snpChrom', snpPos') `compare` (vcfChrom vcfEntry, vcfPos vcfEntry)++processVcfWithSnpFile :: (MonadIO m) => Int -> Pipe (Maybe EigenstratSnpEntry, Maybe VCFentry) FreqSumEntry m r+processVcfWithSnpFile nrInds = for cat $ \jointEntry -> do+ case jointEntry of+ (Just (EigenstratSnpEntry snpChrom' snpPos' _ snpId' snpRef' snpAlt'), Nothing) -> do+ let dosages = replicate nrInds Nothing+ yield $ FreqSumEntry snpChrom' snpPos' (Just . B.unpack $ snpId') snpRef' snpAlt' dosages+ (Just (EigenstratSnpEntry snpChrom' snpPos' _ snpId' snpRef' snpAlt'), Just vcfEntry) -> do+ dosages <- case getDosages vcfEntry of+ Right dos -> return dos+ Left err -> liftIO . throwIO $ AssertionFailed err+ when (length dosages /= nrInds) $ (liftIO . throwIO) (AssertionFailed "inconsistent number of genotypes.")+ let normalizedDosages =+ case vcfAlt vcfEntry of+ [alt] -> if (vcfRef vcfEntry, alt) == (B.singleton snpRef', B.singleton snpAlt')+ then dosages+ else+ if (vcfRef vcfEntry, alt) == (B.singleton snpAlt', B.singleton snpRef')+ then map flipDosages dosages+ else replicate nrInds Nothing+ _ -> replicate nrInds Nothing+ yield $ FreqSumEntry snpChrom' snpPos' (Just . B.unpack $ snpId') snpRef' snpAlt' normalizedDosages+ _ -> return ()+ where+ flipDosages dos = case dos of+ Just 0 -> Just 2+ Just 1 -> Just 1+ Just 2 -> Just 0+ _ -> Nothing++runSimple :: (MonadIO m) => Producer VCFentry m r -> Producer FreqSumEntry m r+runSimple vcfBody = for vcfBody $ \e -> do+ case vcfToFreqSumEntry e of+ Right e' -> do+ liftIO . B.putStr . freqSumEntryToText $ e'+ yield e'+ Left err -> (liftIO . throwIO) (AssertionFailed err)++eigenStratPipe :: (MonadIO m) => Pipe FreqSumEntry (EigenstratSnpEntry, GenoLine) m r+eigenStratPipe = P.map vcfToEigenstrat+ where+ vcfToEigenstrat (FreqSumEntry chrom pos maybeSnpId ref alt dosages) =+ let snpId' = case maybeSnpId of+ Just i -> B.pack i+ Nothing -> B.pack (unChrom chrom <> "_" <> show pos)+ snpEntry = EigenstratSnpEntry chrom pos 0.0 snpId' ref alt+ genoLine = fromList [dosageToCall d | d <- dosages]+ in (snpEntry, genoLine)+ dosageToCall d = case d of+ Just 0 -> HomRef+ Just 1 -> Het+ Just 2 -> HomAlt+ Nothing -> Missing+ _ -> error ("unknown dosage " ++ show d)
+ src/SequenceTools/PileupCaller.hs view
@@ -0,0 +1,150 @@+module SequenceTools.PileupCaller (callToDosage, Call(..), callGenotypeFromPileup,+ callMajorityAllele, findMajorityAlleles, callRandomAllele,+ callRandomDiploid, dosageToEigenstratGeno, freqSumToEigenstrat, CallingMode(..),+ TransitionsMode(..), filterTransitions, cleanSSdamageAllSamples) where++import SequenceFormats.Utils (Chrom(..))+import SequenceFormats.FreqSum (FreqSumEntry(..))+import SequenceFormats.Eigenstrat (EigenstratSnpEntry(..), GenoEntry(..), GenoLine)+import SequenceFormats.Pileup (Strand(..))+import SequenceTools.Utils (sampleWithoutReplacement)++import qualified Data.ByteString.Char8 as B+import Data.List (sortOn, group, sort)+import Data.Vector (fromList)+import Pipes (Pipe, cat)+import qualified Pipes.Prelude as P++-- |A datatype to represent a single genotype call+data Call = HaploidCall Char | DiploidCall Char Char | MissingCall deriving (Show, Eq)++-- |A datatype to specify the calling mode+data CallingMode = MajorityCalling Bool | RandomCalling | RandomDiploidCalling++data TransitionsMode = TransitionsMissing | SkipTransitions | SingleStrandMode | AllSites deriving (Eq)++-- |a function to turn a call into the dosage of non-reference alleles+callToDosage :: Char -> Char -> Call -> Maybe Int+callToDosage refA altA call = case call of+ HaploidCall a | a == refA -> Just 0+ | a == altA -> Just 1+ | otherwise -> Nothing+ DiploidCall a1 a2 | (a1, a2) == (refA, refA) -> Just 0+ | (a1, a2) == (refA, altA) -> Just 1+ | (a1, a2) == (altA, refA) -> Just 1+ | (a1, a2) == (altA, altA) -> Just 2+ | otherwise -> Nothing+ MissingCall -> Nothing++-- |Make a call from alleles +callGenotypeFromPileup :: CallingMode -> Int -> String -> IO Call+callGenotypeFromPileup mode minDepth alleles =+ if length alleles < minDepth then return MissingCall else+ case mode of+ MajorityCalling withDownsampling -> callMajorityAllele withDownsampling minDepth alleles+ RandomCalling -> callRandomAllele alleles+ RandomDiploidCalling -> callRandomDiploid alleles++-- |Sample the majority allele, or one of the majority alleles+callMajorityAllele :: Bool -> Int-> String -> IO Call+callMajorityAllele withDownsampling minDepth alleles = do+ maybeAlleles <- if withDownsampling+ then sampleWithoutReplacement alleles minDepth+ else return (Just alleles)+ case maybeAlleles of+ Nothing -> return MissingCall+ Just alleles' -> do+ a <- case findMajorityAlleles alleles' of+ [] -> error "should not happen"+ [a'] -> return a'+ listA -> do+ r <- sampleWithoutReplacement listA 1+ case r of+ Just [r'] -> return r'+ _ -> error "should not happen"+ return $ HaploidCall a++-- |Find the majority allele(s)+findMajorityAlleles :: String -> [Char]+findMajorityAlleles alleles =+ let groupedAlleles = sortOn fst [(length g, head g) | g <- group . sort $ alleles]+ majorityCount = fst . last $ groupedAlleles+ in [a | (n, a) <- groupedAlleles, n == majorityCount]++-- |call a random allele+callRandomAllele :: String -> IO Call+callRandomAllele alleles = do+ res <- sampleWithoutReplacement alleles 1+ case res of+ Nothing -> return MissingCall+ Just [a] -> return $ HaploidCall a+ _ -> error "should not happen"++-- |call two random alleles+callRandomDiploid :: String -> IO Call+callRandomDiploid alleles = do+ res <- sampleWithoutReplacement alleles 2+ case res of+ Nothing -> return MissingCall+ Just [a1, a2] -> return $ DiploidCall a1 a2+ _ -> error "should not happen"++-- |convert a freqSum entry to an eigenstrat SNP entry+freqSumToEigenstrat :: Bool -> FreqSumEntry -> (EigenstratSnpEntry, GenoLine)+freqSumToEigenstrat diploidizeCall (FreqSumEntry chrom@(Chrom c) pos maybeSnpId ref alt calls) =+ let snpId_ = case maybeSnpId of + Just id_ -> B.pack id_+ Nothing -> B.pack $ c <> "_" <> show pos+ snpEntry = EigenstratSnpEntry chrom pos 0.0 snpId_ ref alt+ geno = fromList . map (dosageToEigenstratGeno diploidizeCall) $ calls+ in (snpEntry, geno)++-- |convert a Dosage to an eigenstrat-encoded genotype+dosageToEigenstratGeno :: Bool -> Maybe Int -> GenoEntry+dosageToEigenstratGeno diploidizeCall c =+ if diploidizeCall then+ case c of+ Just 0 -> HomRef+ Just 1 -> HomAlt+ Nothing -> Missing+ _ -> error "illegal call for pseudo-haploid Calling method"+ else+ case c of+ Just 0 -> HomRef+ Just 1 -> Het+ Just 2 -> HomAlt+ Nothing -> Missing+ _ -> error ("unknown genotype " ++ show c)++filterTransitions :: (Monad m) => TransitionsMode -> Pipe FreqSumEntry FreqSumEntry m ()+filterTransitions transversionsMode =+ case transversionsMode of+ SkipTransitions ->+ P.filter (\(FreqSumEntry _ _ _ ref alt _) -> isTransversion ref alt)+ TransitionsMissing ->+ P.map (\(FreqSumEntry chrom pos id_ ref alt calls) ->+ let calls' = if isTransversion ref alt then calls else [Nothing | _ <- calls]+ in FreqSumEntry chrom pos id_ ref alt calls')+ _ -> cat+ where+ isTransversion ref alt = not $ isTransition ref alt+ isTransition ref alt =+ ((ref == 'A') && (alt == 'G')) ||+ ((ref == 'G') && (alt == 'A')) ||+ ((ref == 'C') && (alt == 'T')) ||+ ((ref == 'T') && (alt == 'C'))++cleanSSdamageAllSamples :: Char -> Char -> [String] -> [[Strand]] -> [String]+cleanSSdamageAllSamples ref alt basesPerSample strandPerSample =+ if (ref, alt) == ('C', 'T') || (ref, alt) == ('T', 'C')+ then [removeForwardBases bases strands | (bases, strands) <- zip basesPerSample strandPerSample]+ else+ if (ref, alt) == ('G', 'A') || (ref, alt) == ('A', 'G')+ then [removeReverseBases bases strands | (bases, strands) <- zip basesPerSample strandPerSample]+ else basesPerSample+ where+ removeForwardBases = removeReads ForwardStrand+ removeReverseBases = removeReads ReverseStrand++removeReads :: Strand -> String -> [Strand] -> String+removeReads strand bases strands = [b | (b, s) <- zip bases strands, s /= strand]
+ src/SequenceTools/Utils.hs view
@@ -0,0 +1,26 @@+module SequenceTools.Utils (versionInfoOpt, versionInfoText, sampleWithoutReplacement) where ++import Paths_sequenceTools (version)+import Data.Version (showVersion)+import qualified Options.Applicative as OP+import System.Random (randomRIO)++versionInfoOpt :: OP.Parser (a -> a)+versionInfoOpt = OP.infoOption (showVersion version) (OP.long "version" <> OP.help "Print version and exit")++versionInfoText :: String+versionInfoText = "This tool is part of sequenceTools version " ++ showVersion version++sampleWithoutReplacement :: [a] -> Int -> IO (Maybe [a])+sampleWithoutReplacement = go []+ where+ go res _ 0 = return $ Just res+ go res xs n+ | n > length xs = return Nothing+ | n == length xs = return $ Just (xs ++ res)+ | otherwise = do+ rn <- randomRIO (0, length xs - 1)+ let a = xs !! rn+ xs' = remove rn xs+ go (a:res) xs' (n - 1)+ remove i xs = let (ys, zs) = splitAt i xs in ys ++ tail zs
+ test/SequenceTools/PileupCallerSpec.hs view
@@ -0,0 +1,187 @@+{-# LANGUAGE OverloadedStrings #-}+module SequenceTools.PileupCallerSpec (spec) where++import SequenceFormats.FreqSum (FreqSumEntry(..))+import SequenceFormats.Eigenstrat (GenoEntry(..), EigenstratSnpEntry(..))+import SequenceFormats.Utils (Chrom(..))+import SequenceFormats.Pileup (Strand(..))+import SequenceTools.PileupCaller (callToDosage, Call(..), callGenotypeFromPileup,+ callMajorityAllele, findMajorityAlleles, callRandomAllele,+ callRandomDiploid, dosageToEigenstratGeno, freqSumToEigenstrat, CallingMode(..),+ TransitionsMode(..), filterTransitions, cleanSSdamageAllSamples)++import Control.Monad (replicateM, forM_)+import qualified Data.ByteString.Char8 as B+import Data.List (sort)+import Data.Vector (fromList)+import Pipes (each, (>->))+import qualified Pipes.Prelude as P+import Test.Hspec++spec :: Spec+spec = do+ testCallToDosage+ testCallRandomDiploid+ testCallRandomAllele+ testCallMajorityAllele+ testFindMajorityAlleles+ testCallGenotypeFromPileup+ testDosageToEigenstratGeno+ testFreqSumToEigenstrat+ testFilterTransitions+ testCleanSSdamageAllSamples++testCallToDosage :: Spec+testCallToDosage = describe "callToDosage" $ do+ it "should return Nothing for missing call" $+ callToDosage 'A' 'C' MissingCall `shouldBe` Nothing+ it "should return Nothing for haploid non-congruent call" $+ callToDosage 'A' 'C' (HaploidCall 'G') `shouldBe` Nothing+ it "should return 0 for haploid ref call" $+ callToDosage 'A' 'C' (HaploidCall 'A') `shouldBe` Just 0+ it "should return 1 for haploid alt call" $+ callToDosage 'A' 'C' (HaploidCall 'C') `shouldBe` Just 1 + it "should return Nothing for diploid non-congruent call" $+ callToDosage 'A' 'C' (DiploidCall 'A' 'G') `shouldBe` Nothing+ it "should return 0 for diploid hom-ref call" $+ callToDosage 'A' 'C' (DiploidCall 'A' 'A') `shouldBe` Just 0+ it "should return 1 for diploid het call" $ do+ callToDosage 'A' 'C' (DiploidCall 'A' 'C') `shouldBe` Just 1+ callToDosage 'A' 'C' (DiploidCall 'C' 'A') `shouldBe` Just 1+ it "should return 2 for diploid hom-alt call" $+ callToDosage 'A' 'C' (DiploidCall 'C' 'C') `shouldBe` Just 2+++testCallGenotypeFromPileup :: Spec+testCallGenotypeFromPileup = describe "callGenotypeFromPileup" $ do+ it "should return missing if pileup below minDepth" $+ callGenotypeFromPileup RandomCalling 3 "A" `shouldReturn` MissingCall+ it "should not return missing if pileup above minDepth" $+ callGenotypeFromPileup RandomCalling 3 "AACCC" `shouldNotReturn` MissingCall+ ++testCallMajorityAllele :: Spec+testCallMajorityAllele = describe "callMajorityAllele" $ do+ it "should call A from AAA" $+ callMajorityAllele False 1 "AAA" `shouldReturn` HaploidCall 'A'+ it "should call A from AAAAA with ds" $+ callMajorityAllele True 3 "AAAAA" `shouldReturn` HaploidCall 'A'+ it "should call Missing from AA with ds 3" $+ callMajorityAllele True 3 "AA" `shouldReturn` MissingCall+ it "should call A from AAC" $+ callMajorityAllele False 1 "AAC" `shouldReturn` HaploidCall 'A'+ it "should call C from ACC" $+ callMajorityAllele False 1 "ACC" `shouldReturn` HaploidCall 'C'+ it "should call 50/50 from AACC" $ do+ r <- replicateM 1000 (callMajorityAllele False 1 "AACC")+ let c = [rr | rr <- r, rr == HaploidCall 'A']+ length c `shouldSatisfy` (\c' -> c' >= 418 && c' <= 582) --p < 1e-7++testFindMajorityAlleles :: Spec+testFindMajorityAlleles = describe "findMajorityAllele" $ do+ it "should return A for AAC" $+ findMajorityAlleles "AAC" `shouldBe` "A"+ it "should return C for ACC" $+ findMajorityAlleles "ACC" `shouldBe` "C"+ it "should return AC for AACC" $+ findMajorityAlleles "AACC" `shouldBe` "AC"++testCallRandomAllele :: Spec+testCallRandomAllele = describe "callRandomAllele" $ do+ it "should return A for AAA" $+ callRandomAllele "AAA" `shouldReturn` HaploidCall 'A'+ it "should return C for C" $+ callRandomAllele "C" `shouldReturn` HaploidCall 'C'+ it "should return A,C or G for ACG roughly with 30% each" $ do+ r <- replicateM 1000 (callRandomAllele "ACG")+ forM_ ['A', 'C', 'G'] $ \nuc -> do+ let n = length . filter (==HaploidCall nuc) $ r+ n `shouldSatisfy` (\nn -> nn >= 257 && nn <= 412) --p < 1e-7++testCallRandomDiploid :: Spec+testCallRandomDiploid = describe "callRandomDiploid" $ do+ it "should return Missing for A" $+ callRandomDiploid "A" `shouldReturn` MissingCall+ it "should return AC for AC" $ do+ DiploidCall a1 a2 <- callRandomDiploid "AC"+ sort [a1, a2] `shouldBe` "AC"+ it "should return 50% het for AACC" $ do+ r <- replicateM 1000 (callRandomDiploid "AACC")+ let n = length ['A' | DiploidCall a1 a2 <- r, a1 /= a2]+ n `shouldSatisfy` (\nn -> nn >= 588 && nn < 743) --p < 1e-7++testDosageToEigenstratGeno :: Spec+testDosageToEigenstratGeno = describe "dosageToEigenstratGeno" $ do+ it "should give Hom-Ref for 0 pseudo-haploid" $+ dosageToEigenstratGeno True (Just 0) `shouldBe` HomRef+ it "should give Hom-Alt for 1 pseudo-haploid" $+ dosageToEigenstratGeno True (Just 1) `shouldBe` HomAlt+ it "should give Missing for Nothing pseudo-haploid" $+ dosageToEigenstratGeno True Nothing `shouldBe` Missing+ it "should give Hom-Ref for 0 diploid" $+ dosageToEigenstratGeno False (Just 0) `shouldBe` HomRef+ it "should give Het for 1 diploid" $+ dosageToEigenstratGeno False (Just 1) `shouldBe` Het+ it "should give Hom-Alt for 2 diploid" $+ dosageToEigenstratGeno False (Just 2) `shouldBe` HomAlt+ it "should give Missing for Nothing diploid" $+ dosageToEigenstratGeno False Nothing `shouldBe` Missing++testFreqSumToEigenstrat :: Spec+testFreqSumToEigenstrat = describe "freqSumtoEigenstrat" $ do+ let fs = FreqSumEntry (Chrom "1") 1000 Nothing 'A' 'C' [Just 0, Just 1, Just 1, Nothing, Just 0]+ let es = EigenstratSnpEntry (Chrom "1") 1000 0.0 (B.pack "1_1000") 'A' 'C'+ genoLine = fromList [HomRef, HomAlt, HomAlt, Missing, HomRef]+ it "should convert a freqSum example correctly to eigenstrat" $+ freqSumToEigenstrat True fs `shouldBe` (es, genoLine)+ it "should convert a freqSum example with rsId correctly to eigenstrat" $+ freqSumToEigenstrat True (fs {fsSnpId = Just "rs123"}) `shouldBe` (es {snpId = "rs123"}, genoLine)+++mockFreqSumData :: [FreqSumEntry]+mockFreqSumData = [+ FreqSumEntry (Chrom "1") 1000 (Just "rs1") 'A' 'C' [Just 1, Just 2, Nothing, Just 0, Just 0],+ FreqSumEntry (Chrom "1") 2000 (Just "rs2") 'C' 'T' [Just 1, Just 2, Nothing, Just 0, Just 0],+ FreqSumEntry (Chrom "1") 3000 (Just "rs3") 'A' 'G' [Just 1, Just 2, Nothing, Just 0, Just 0],+ FreqSumEntry (Chrom "2") 1000 (Just "rs4") 'A' 'G' [Just 1, Just 2, Nothing, Just 0, Just 0],+ FreqSumEntry (Chrom "2") 2000 (Just "rs5") 'T' 'A' [Just 1, Just 2, Nothing, Just 0, Just 0],+ FreqSumEntry (Chrom "2") 3000 (Just "rs6") 'T' 'C' [Just 1, Just 2, Nothing, Just 0, Just 0]]++testFilterTransitions :: Spec+testFilterTransitions = describe "filterTransitions" $ do+ it "should remove transitions with SkipTransitions" $ do+ let r = P.toList $ each mockFreqSumData >-> filterTransitions SkipTransitions+ r `shouldBe` [+ FreqSumEntry (Chrom "1") 1000 (Just "rs1") 'A' 'C' [Just 1, Just 2, Nothing, Just 0, Just 0],+ FreqSumEntry (Chrom "2") 2000 (Just "rs5") 'T' 'A' [Just 1, Just 2, Nothing, Just 0, Just 0]]+ it "should mark transitions as missing with TransitionsMissing" $ do+ let r = P.toList $ each mockFreqSumData >-> filterTransitions TransitionsMissing+ r `shouldBe` [+ FreqSumEntry (Chrom "1") 1000 (Just "rs1") 'A' 'C' [Just 1, Just 2, Nothing, Just 0, Just 0],+ FreqSumEntry (Chrom "1") 2000 (Just "rs2") 'C' 'T' [Nothing, Nothing, Nothing, Nothing, Nothing],+ FreqSumEntry (Chrom "1") 3000 (Just "rs3") 'A' 'G' [Nothing, Nothing, Nothing, Nothing, Nothing],+ FreqSumEntry (Chrom "2") 1000 (Just "rs4") 'A' 'G' [Nothing, Nothing, Nothing, Nothing, Nothing],+ FreqSumEntry (Chrom "2") 2000 (Just "rs5") 'T' 'A' [Just 1, Just 2, Nothing, Just 0, Just 0],+ FreqSumEntry (Chrom "2") 3000 (Just "rs6") 'T' 'C' [Nothing, Nothing, Nothing, Nothing, Nothing]]+ it "should output all sites with AllSites" $ do+ let r = P.toList $ each mockFreqSumData >-> filterTransitions AllSites+ r `shouldBe` mockFreqSumData+ it "should output all sites with SingleStrandMode" $ do+ let r = P.toList $ each mockFreqSumData >-> filterTransitions SingleStrandMode+ r `shouldBe` mockFreqSumData++testCleanSSdamageAllSamples :: Spec+testCleanSSdamageAllSamples = describe "cleanSSdamageAllSamples" $ do+ let bases = ["AACATG", "AACATT", "AACTTG"]+ strands = [[f, r, r, f, r, r], [r, f, r, f, f, r], [f, f, r, f, f, r]]+ it "should not remove anything if not C/T or G/A SNP" $+ cleanSSdamageAllSamples 'C' 'A' bases strands `shouldBe` bases+ it "should remove forward reads from C/T SNPs" $ do+ cleanSSdamageAllSamples 'C' 'T' bases strands `shouldBe` ["ACTG", "ACT", "CG"]+ cleanSSdamageAllSamples 'T' 'C' bases strands `shouldBe` ["ACTG", "ACT", "CG"]+ it "should remove reverse reads from G/A SNPs" $ do+ cleanSSdamageAllSamples 'A' 'G' bases strands `shouldBe` ["AA", "AAT", "AATT"]+ cleanSSdamageAllSamples 'G' 'A' bases strands `shouldBe` ["AA", "AAT", "AATT"]+ where+ f = ForwardStrand+ r = ReverseStrand
+ test/SequenceTools/UtilsSpec.hs view
@@ -0,0 +1,32 @@+module SequenceTools.UtilsSpec (spec) where++import SequenceTools.Utils (sampleWithoutReplacement)++import Control.Monad (replicateM_)+import Data.List (sort, union, nub)+import Test.Hspec++spec :: Spec+spec = testSampleWithoutReplacement++testSampleWithoutReplacement :: Spec+testSampleWithoutReplacement = describe "sampleWithoutReplacement" $ do+ it "should return Nothing if sample 1 from empty list" $+ sampleWithoutReplacement ([] :: [Char]) 1 `shouldReturn` Nothing+ it "should return one item if sample 1 from 1" $+ sampleWithoutReplacement ['A'] 1 `shouldReturn` Just ['A']+ it "should return an empty list of 0 sampled from 0" $+ sampleWithoutReplacement ([] :: [Char]) 0 `shouldReturn` Just []+ it "should return an empty list if 0 sampled from 1" $+ sampleWithoutReplacement ['A'] 0 `shouldReturn` Just []+ it "should return Nothing if 2 sampled from 1" $+ sampleWithoutReplacement ['A'] 2 `shouldReturn` Nothing+ it "should return 2 items if 2 sampled from 2" $ do+ r <- sampleWithoutReplacement ['A', 'C'] 2+ fmap sort r `shouldBe` Just (sort ['A', 'C'])+ it "should return a non-duplicate subset of ABCDEFGHIJ if 4 are sampled" $ do+ replicateM_ 10 $ do+ Just r <- sampleWithoutReplacement "ABCDEFGHIJ" 4+ length (nub r) `shouldBe` 4+ length (union r "ABCDEFGHIJ") `shouldBe` 10+
+ test/Spec.hs view
@@ -0,0 +1,1 @@+{-# OPTIONS_GHC -F -pgmF hspec-discover #-}