sequenceTools 1.5.3.1 → 1.6.0.0
raw patch · 12 files changed
+752/−360 lines, 12 filesdep +prettyprinterdep +textdep −ansi-wl-pprintdep ~basedep ~bytestringdep ~foldlsetup-changed
Dependencies added: prettyprinter, text
Dependencies removed: ansi-wl-pprint
Dependency ranges changed: base, bytestring, foldl, hspec, lens-family, optparse-applicative, pipes, pipes-group, pipes-ordered-zip, pipes-safe, random, sequence-formats, split, transformers, vector
Files
- Changelog.md +8/−0
- LICENSE +1/−1
- README.md +200/−17
- Setup.hs +1/−1
- sequenceTools.cabal +98/−50
- src-executables/genoStats.hs +40/−42
- src-executables/pileupCaller.hs +236/−105
- src-executables/vcf2eigenstrat.hs +37/−37
- src/SequenceTools/PileupCaller.hs +35/−24
- src/SequenceTools/Utils.hs +29/−31
- test/SequenceTools/PileupCallerSpec.hs +63/−48
- test/SequenceTools/UtilsSpec.hs +4/−4
Changelog.md view
@@ -1,5 +1,13 @@ # Changelog +- V 1.6.0.0: + - pileupCaller now supports writing of VCF (Variant Call Format), via the standard output. See option --vcf and the Readme.+- V 1.5.4.0:+ - updated sequence-formats dependency allows more lenient parsing of pileup-data, now also allowing for arbitrary reference alleles (not just ACTGN). This won't affect calling (reads that support an allele that is not in the SNP-file input are treated as before), but will be less disruptive when parsing pileup-input, for example without a bed-file in samtools.+ - improved error output for parsing problems with pileup-format data. Now only a small part of the problematic chunk is output, hopefully easing error interpretation in such cases+ - output a useful error message if the number of samples passed in --sampleNames is inconsistent with the pileup-input+ - `--samplePopName` now accepts multiple pop-names, separated by comma. The number of pop-names must then match the number of samples.+- V 1.5.3.2: fixed a bug in vcf2eigenstrat that would fail on VCFs with missing Quality values. - V 1.5.3.1: updated to latest GHC pedantic compilation - V 1.5.3: Upgraded to sequence-formats 1.7.0 introducing an option for plink popName encoding, and improved pileup-Parsing to allow for skip-reference characters - V 1.5.2: Fixed a bug with --samplePopName having to be entered after -p or -e. Fixed a bug in the sequence-formats dependency.
LICENSE view
@@ -1,4 +1,4 @@-Copyright Author name here (c) 2015+Copyright Stephan Schiffels (c) 2024 All rights reserved.
README.md view
@@ -29,11 +29,11 @@ This last step will take a while, as it not only compiles the source, but also first downloads the correct version of the Haskell compiler. -# Commands+# PileupCaller -## pileupCaller+## Overview -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:+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 and uses a predefined set of positions to generate the pileup for (optional, see below): samtools mpileup -R -B -q30 -Q30 -l <list_of_positions.txt> \ -f <reference_genome.fasta> \@@ -41,9 +41,9 @@ 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.+In the above command line, if you use a positions-file, it 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:+Next, you need to run pileupCaller, which you run like this: pileupCaller --randomHaploid --sampleNames Sample1,Sample2,Sample3 \ --samplePopName MyPop -f <Eigenstrat.snp> \@@ -61,32 +61,192 @@ 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> \+ samtools mpileup -R -B -q30 -Q30 \ -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> -Note that `--randomHaploid` is only one way to call genotypes. If you need stricter calling, you may want to try `--majorityCall --downSampling --minDepth 3`, which calls genotypes only on sites with at least three reads, downsamples to three if there are more, and then calls whatever of the two alleles has the majority. This will reduce errors, but also yield less data in case of lower coverage.- +Here, I omitted the positions-file in the samtools command, because pileupCaller itself will ensure filtering for the positions listed in the Eigenstrat-Positions file given via option `-f`. Note that `--randomHaploid` is only one way to call genotypes. If you need stricter calling, you may want to try `--majorityCall --downSampling --minDepth 3`, which calls genotypes only on sites with at least three reads, downsamples to three if there are more, and then calls whatever of the two alleles has the majority. This will reduce errors, but also yield less data in case of lower coverage. -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:+You will possibly encounter an issue: 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> \+ samtools mpileup -R -B -q30 -Q30 \ -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> - -Note: You do not have to use a positions file in your `samtools` step. You can also just generate pileup-data for every covered position (default without passing a positions file via `-l`) and have pileupCaller filter the sites for you. This makes sense for dense genotyping, but a positions file might speed up the process for sparser genotyping.+## Options +You can see all options via `pileupCaller --help`, which outputs:++```+Usage: pileupCaller [--version] + (--randomHaploid | --majorityCall [--downSampling] | + --randomDiploid) [--keepIncongruentReads] + [--seed <RANDOM_SEED>] [-d|--minDepth <DEPTH>] + [--skipTransitions | --transitionsMissing | + --singleStrandMode] (-f|--snpFile <FILE>) + [(-e|--eigenstratOut <FILE_PREFIX>) [-z|----zip] | + (-p|--plinkOut <FILE_PREFIX>) + [--popNameAsPhenotype | --popNameAsBoth] [-z|----zip] |+ --vcf] + (--sampleNames NAME1,NAME2,... | --sampleNameFile <FILE>) + [--samplePopName POP(s)]++ 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++ samtools mpileup -B -q30 -Q30 -l <BED_FILE> -R -f <FASTA_REFERENCE_FILE>+ Sample1.bam Sample2.bam Sample3.bam | pileupCaller ...++ 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.+++ This tool is part of sequenceTools version 1.6.0.0++Available options:+ --version Print version and exit+ -h,--help Show this help text+ --randomHaploid 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+ --majorityCall 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+ --downSampling 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+ --randomDiploid 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.+ --keepIncongruentReads 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.+ --seed <RANDOM_SEED> random seed used for the random number generator. If+ not given, use system clock to seed the random number+ generator.+ -d,--minDepth <DEPTH> specify the minimum depth for a call. For sites with+ fewer reads than this number, declare Missing+ (default: 1)+ --skipTransitions skip transition SNPs entirely in the output,+ resulting in a dataset with fewer sites.+ --transitionsMissing mark transitions as missing in the output, but do+ output the sites.+ --singleStrandMode [THIS IS CURRENTLY AN EXPERIMENTAL FEATURE]. 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.+ -f,--snpFile <FILE> 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).+ -e,--eigenstratOut <FILE_PREFIX>+ Set Eigenstrat as output format. Specify the+ filenames for the EigenStrat SNP, IND and GENO file+ outputs: <FILE_PREFIX>.snp, <FILE_PREFIX>.ind and+ <FILE_PREFIX>.geno. 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+ -z,----zip GZip the output Eigenstrat or Plink genotype and SNP+ files. Filenames will be appended with '.gz'. To zip+ FreqSum or VCF output, just zip the standard output+ of this program, for example `pileupCaller ... --vcf+ | gzip -c > out.vcf.gz+ -p,--plinkOut <FILE_PREFIX>+ Set Plink as output format. Specify the filenames for+ the Plink BIM, FAM and BED file outputs:+ <FILE_PREFIX>.bim, <FILE_PREFIX>.fam and+ <FILE_PREFIX>.bed. 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+ --popNameAsPhenotype Only valid for Plink Output: Write the population+ name into the last column of the fam file, as a+ Phenotype according to the Plink Spec. By default,+ the population name is specified as the first column+ only (family name in the Plink spec)+ --popNameAsBoth Only valid for Plink Output: Write the population+ name into both the first and last column of the fam+ file, so both as Family-ID and as a Phenotype+ according to the Plink Spec. By default, the+ population name is specified only as the first column+ (family name in the Plink spec)+ -z,----zip GZip the output Eigenstrat or Plink genotype and SNP+ files. Filenames will be appended with '.gz'. To zip+ FreqSum or VCF output, just zip the standard output+ of this program, for example `pileupCaller ... --vcf+ | gzip -c > out.vcf.gz+ --vcf output VCF format to stdout+ --sampleNames NAME1,NAME2,...+ give the names of the samples as comma-separated list+ (no spaces)+ --sampleNameFile <FILE> give the names of the samples in a file with one name+ per line+ --samplePopName POP(s) specify the population name(s) of the samples, which+ are included in the output *.ind.txt file in+ Eigenstrat output. This will be ignored if the output+ format is not Eigenstrat. If a single name is given,+ it is applied to all samples, if multiple are given,+ their number must match the the number of samples+ (default: Left "Unknown")+```+ ### SingleStrandMode pileupCaller supports a special calling mode (`--singleStrandMode`) for sequencing data generated from single-stranded libraries (Gansauge, Marie-Theres, and Matthias Meyer. 2013. “Single-Stranded DNA Library Preparation for the Sequencing of Ancient or Damaged DNA.” Nature Protocols 8 (4): 737–48.). The idea is that at C/T SNPs, forward mapping reads are discarded, and at G/A SNPs, reverse mapping reads are discarded. This will get rid of post-mortem ancient DNA damage in a conservative way, i.e. it will remove more than necessary and make sure that the remainder of the data is clean of DNA damage, improving the overall calling quality.@@ -97,16 +257,39 @@ 1) Use only merged reads and `--singleStrandMode` 2) Use all reads but do _not_ use `--singleStrandMode`. Instead, in such cases I recommend to trim reads from both ends to remove ancient DNA damage. Depending on the details of the library construction, you may have UDG-treated data, in which case fewer basepairs would have to be trimmed. +### VCF output+VCF output was added in version 1.6.0.0. The VCF format is specified in detail at https://samtools.github.io/hts-specs/VCFv4.5.pdf. I just mention two specifics. First, with calling modes `--randomHaploid` and `--majorityCall`, the output genotypes will be haploid. This means that instead of genotypes like `0/0`, `0/1`, `1/1` or `./.`, you will instead just see `0`, `1` or `.`. Second, I added some possibly useful filters and statistics to the output, which are described in the header of the VCF: -## vcf2eigenstrat+```+##fileformat=VCFv4.2+##source=pileupCaller_v1.6.0.0+##command_line=pileupCaller --randomHaploid -f 1240k_eigenstrat_snp_short.snp.txt --sampleNames 1,2,3,4 --vcf+##group_names=Unknown,Unknown,Unknown,Unknown+##INFO=<ID=NS,Number=1,Type=Integer,Description="Number of Samples With Data">+##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">+##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency">+##FILTER=<ID=s50,Description="Less than 50% of samples have data">+##FILTER=<ID=s10,Description="Less than 10% of samples have data">+##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">+##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read Depth">+##FORMAT=<ID=DP8,Number=8,Type=Integer,Description="Nr of Reads supporting A,C,G,T in forward strand, followed by the same quartet in reverse strand">+``` +As you can see, Info fields NS, DP and AF are added and defined, as well as two filters which might come in handy. Also, beyond the required genotype `GT` tag, I added two per-sample tags `DP` and `DP8` as defined in the header. ++### Summary Statistics++PileupCaller automatically outputs a few lines of summary statistics, including the number of sites called for each sample, and the average read depth. These are output to the stderr, so do not affect stdout or file output.++# vcf2eigenstrat+ Simple tool to convert a VCF file to an Eigenstrat file. Pretty self-explanatory. Please run `vcf2eigenstrat --help` to output some documentation. -## genoStats+# genoStats A simple tool to get some per-individual statistics from an Eigenstrat or Freqsum-file. Run `genoStats --help` for documentation. -## Scripts+# Scripts This package also contains several haskell wrapper scripts for the following [ADMIXTOOLS and EIGENSOFT](https://reich.hms.harvard.edu/software) commands: convertf, mergeit, qp3Pop, qpDstat and smartPCA. The original tools require parameter files as input, which I find tedious to use in bioinformatics pipelines. I wrote those wrapper scripts to be able to start the tools with a simple command line option interface. If you have `stack` installed your system (see above), you should be able to run those scripts on your machine without any difficult setup. Simply clone this repository, navigate to the `scripts` subfolder and invoke any script using standard bash execution, for example
Setup.hs view
@@ -1,2 +1,2 @@-import Distribution.Simple+import Distribution.Simple main = defaultMain
sequenceTools.cabal view
@@ -1,59 +1,107 @@-name: sequenceTools-version: 1.5.3.1-synopsis: A package with tools for processing DNA sequencing data-description: The tools in this package process sequencing Data, in particular from ancient DNA sequencing libraries. Key tool in this package is pileupCaller, a tool to randomly sample genotypes from sequencing data.-license: GPL-3-license-file: LICENSE-author: Stephan Schiffels-maintainer: stephan.schiffels@mac.com-category: Bioinformatics-build-type: Simple-cabal-version: >=1.10-Homepage: https://github.com/stschiff/sequenceTools-Bug-Reports: https://github.com/stschiff/sequenceTools/issues-extra-source-files: README.md,- Changelog.md+cabal-version: >=1.10+name: sequenceTools+version: 1.6.0.0+license: GPL-3+license-file: LICENSE+maintainer: stephan.schiffels@mac.com+author: Stephan Schiffels+homepage: https://github.com/stschiff/sequenceTools+bug-reports: https://github.com/stschiff/sequenceTools/issues+synopsis: A package with tools for processing DNA sequencing data+description:+ The tools in this package process sequencing Data, in particular from ancient DNA sequencing libraries. 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- build-depends: base >= 4.7 && < 5, optparse-applicative, random,- sequence-formats, bytestring, vector, pipes- other-modules: Paths_sequenceTools- default-language: Haskell2010+ exposed-modules:+ SequenceTools.Utils+ SequenceTools.PileupCaller -Test-Suite sequenceToolsTests- type: exitcode-stdio-1.0- main-is: Spec.hs- hs-source-dirs: test- build-depends: base, hspec, sequenceTools, sequence-formats,- vector, bytestring, pipes- other-modules: SequenceTools.UtilsSpec,- SequenceTools.PileupCallerSpec- default-language: Haskell2010+ hs-source-dirs: src+ other-modules: Paths_sequenceTools+ default-language: Haskell2010+ build-depends:+ base >=4.7 && <5,+ optparse-applicative >=0.18.1.0,+ random >=1.2.1.2,+ sequence-formats >=1.11.0.0,+ bytestring >=0.11.5.3,+ vector >=0.13.1.0,+ pipes >=4.3.16 executable pileupCaller- main-is: pileupCaller.hs- hs-source-dirs: src-executables- build-depends: base, sequenceTools, sequence-formats,- optparse-applicative, pipes, vector, random, bytestring,- pipes-safe, pipes-ordered-zip, split, ansi-wl-pprint, transformers- default-language: Haskell2010- ghc-options: -threaded -rtsopts -with-rtsopts=-N2+ main-is: pileupCaller.hs+ hs-source-dirs: src-executables+ other-modules: Paths_sequenceTools+ default-language: Haskell2010+ ghc-options: -threaded -rtsopts -with-rtsopts=-N2+ build-depends:+ base >=4.18.2.1,+ sequenceTools,+ sequence-formats >=1.11.0.0,+ optparse-applicative >=0.18.1.0,+ pipes >=4.3.16,+ vector >=0.13.1.0,+ random >=1.2.1.2,+ bytestring >=0.11.5.3,+ pipes-safe >=2.3.5,+ pipes-ordered-zip >=1.2.1,+ split >=0.2.5,+ prettyprinter >=1.7.1,+ transformers >=0.6.1.0,+ text >=2.0.2 executable vcf2eigenstrat- main-is: vcf2eigenstrat.hs- hs-source-dirs: src-executables- build-depends: base, sequenceTools, pipes-ordered-zip,- sequence-formats, bytestring, vector,- optparse-applicative, pipes, pipes-safe- default-language: Haskell2010+ main-is: vcf2eigenstrat.hs+ hs-source-dirs: src-executables+ default-language: Haskell2010+ build-depends:+ base >=4.18.2.1,+ sequenceTools,+ pipes-ordered-zip >=1.2.1,+ sequence-formats >=1.11.0.0,+ bytestring >=0.11.5.3,+ vector >=0.13.1.0,+ optparse-applicative >=0.18.1.0,+ pipes >=4.3.16,+ pipes-safe >=2.3.5 executable genoStats- main-is: genoStats.hs- hs-source-dirs: src-executables- build-depends: base, sequence-formats, sequenceTools,- foldl, bytestring, vector, lens-family, optparse-applicative,- pipes, pipes-group, pipes-safe- default-language: Haskell2010+ main-is: genoStats.hs+ hs-source-dirs: src-executables+ default-language: Haskell2010+ build-depends:+ base >=4.18.2.1,+ sequence-formats >=1.11.0.0,+ sequenceTools,+ foldl >=1.4.17,+ bytestring >=0.11.5.3,+ vector >=0.13.1.0,+ lens-family >=2.1.3,+ optparse-applicative >=0.18.1.0,+ pipes >=4.3.16,+ pipes-group >=1.0.12,+ pipes-safe >=2.3.5++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.18.2.1,+ hspec >=2.11.10,+ sequenceTools,+ sequence-formats >=1.11.0.0,+ vector >=0.13.1.0,+ bytestring >=0.11.5.3,+ pipes >=4.3.16
src-executables/genoStats.hs view
@@ -1,24 +1,31 @@ {-# 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 SequenceFormats.Eigenstrat (EigenstratIndEntry (..),+ EigenstratSnpEntry (..),+ GenoEntry (..), GenoLine,+ readEigenstrat)+import SequenceFormats.FreqSum (FreqSumEntry (..),+ FreqSumHeader (..),+ readFreqSumFile, readFreqSumStdIn)+import SequenceFormats.Utils (Chrom (..))+import SequenceTools.Utils (dosageToEigenstratGeno,+ versionInfoOpt, versionInfoText) -import Control.Applicative ((<|>))-import Control.Foldl (purely, Fold(..))-import Control.Monad (forM_)-import Control.Monad.IO.Class (MonadIO, liftIO)-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)+import Control.Applicative ((<|>))+import Control.Foldl (Fold (..), purely)+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 (Consumer, Producer, cat, for,+ (>->))+import Pipes.Group (folds, groupsBy)+import qualified Pipes.Prelude as P+import Pipes.Safe (MonadSafe, runSafeT)+import System.IO (hPutStrLn, stderr)+import Text.Printf (printf) data ProgOpt = ProgOpt InputOption @@ -30,9 +37,9 @@ data StatsReport = StatsReport { srNrSitesMissing :: Int,- srNrSitesHomRef :: Int,- srNrSitesHomAlt :: Int,- srNrSitesHet :: Int+ srNrSitesHomRef :: Int,+ srNrSitesHomAlt :: Int,+ srNrSitesHet :: Int } deriving (Show) main :: IO ()@@ -66,7 +73,7 @@ runWithOpts :: ProgOpt -> IO () runWithOpts (ProgOpt inputOpt) = runSafeT $ do- (names, entryProducer) <- case inputOpt of + (names, entryProducer) <- case inputOpt of FreqsumInput fsFile -> runWithFreqSum fsFile EigenstratInput genoFile snpFile indFile -> runWithEigenstrat genoFile snpFile indFile liftIO $ hPutStrLn stderr ("processing samples: " <> show names)@@ -76,19 +83,10 @@ runWithFreqSum :: (MonadSafe m) => Maybe FilePath -> m ([String], Producer InputEntry m ()) runWithFreqSum fsFile = do- (FreqSumHeader names nrHaps, fsProd) <- case fsFile of+ (FreqSumHeader names _, 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+ let prod = fsProd >-> P.map (\fs -> InputEntry (fsChrom fs) (V.fromList . map dosageToEigenstratGeno . fsCounts $ fs)) return (names, prod) runWithEigenstrat :: (MonadSafe m) =>@@ -97,19 +95,19 @@ (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)+ return (map B.unpack 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)) + 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)] - + traverse runStatsPerChromPerSample [0..(nrSamples - 1)]+ getChrom :: Fold InputEntry Chrom getChrom = Fold (\_ (InputEntry c _) -> c) (Chrom "") id @@ -119,14 +117,14 @@ 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}+ 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"@@ -153,4 +151,4 @@ extract :: StatsReportAllSamples -> (Chrom, StatsReportAllSamples) extract r = (Chrom "Total", r) - +
src-executables/pileupCaller.hs view
@@ -1,68 +1,111 @@ {-# LANGUAGE OverloadedStrings #-} -import SequenceFormats.Eigenstrat (readEigenstratSnpFile, EigenstratSnpEntry(..), - EigenstratIndEntry(..), Sex(..), writeEigenstrat)-import SequenceFormats.FreqSum(FreqSumEntry(..), printFreqSumStdOut, FreqSumHeader(..))-import SequenceFormats.Pileup (PileupRow(..), readPileupFromStdIn)-import SequenceTools.Utils (versionInfoOpt, versionInfoText, freqSumToEigenstrat)-import SequenceTools.PileupCaller (CallingMode(..), callGenotypeFromPileup, callToDosage,- filterTransitions, TransitionsMode(..), cleanSSdamageAllSamples)-import SequenceFormats.Plink (writePlink, PlinkPopNameMode(..), eigenstratInd2PlinkFam)+import SequenceFormats.Eigenstrat (EigenstratIndEntry (..),+ EigenstratSnpEntry (..),+ Sex (..),+ readEigenstratSnpFile,+ writeEigenstrat)+import SequenceFormats.FreqSum (FreqSumEntry (..),+ FreqSumHeader (..),+ printFreqSumStdOut)+import SequenceFormats.Pileup (PileupRow (..),+ readPileupFromStdIn,+ Strand(..))+import SequenceFormats.Plink (PlinkPopNameMode (..),+ eigenstratInd2PlinkFam,+ writePlink)+import SequenceFormats.Utils (SeqFormatException (..))+import SequenceFormats.VCF (VCFentry (..), VCFheader (..),+ printVCFtoStdOut)+import SequenceTools.PileupCaller (CallingMode (..),+ TransitionsMode (..),+ callGenotypeFromPileup,+ callToDosage,+ cleanSSdamageAllSamples,+ computeAlleleFreq,+ filterTransitions)+import SequenceTools.Utils (UserInputException (..),+ freqSumToEigenstrat,+ versionInfoOpt,+ versionInfoText) -import Control.Applicative ((<|>))-import Control.Monad.Trans.Reader (ReaderT, asks, runReaderT)-import Control.Monad.Trans.Class (lift)-import Control.Monad.IO.Class (liftIO)-import Control.Monad (forM_)-import Data.IORef (IORef, readIORef, modifyIORef', newIORef)-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 qualified Pipes.Prelude as P-import Pipes.OrderedZip (orderedZip, orderCheckPipe)-import Pipes.Safe (runSafeT, SafeT)-import System.IO (hPutStrLn, stderr)-import System.Random (mkStdGen, setStdGen)-import Text.Printf (printf)-import qualified Text.PrettyPrint.ANSI.Leijen as PP+import Control.Applicative ((<|>))+import Control.Exception (catch, throwIO)+import Control.Monad (forM_, when)+import Control.Monad.IO.Class (liftIO)+import Control.Monad.Trans.Class (lift)+import Control.Monad.Trans.Reader (ReaderT, asks,+ runReaderT)+import qualified Data.ByteString.Char8 as B+import Data.IORef (IORef, modifyIORef', newIORef,+ readIORef)+import Data.List (intercalate)+import Data.List.Split (splitOn)+import qualified Data.Text as T+import qualified Data.Vector.Unboxed.Mutable as V+import Data.Version (Version, showVersion)+import qualified Options.Applicative as OP+import qualified Options.Applicative.Help.Pretty as PP+import Paths_sequenceTools (version)+import Pipes (Producer, for, runEffect,+ yield, (>->))+import Pipes.OrderedZip (orderCheckPipe, orderedZip)+import qualified Pipes.Prelude as P+import Pipes.Safe (SafeT, runSafeT)+import qualified Prettyprinter.Util+import System.Environment (getArgs, getProgName)+import System.IO (hPutStrLn, stderr)+import System.Random (mkStdGen, setStdGen)+import Text.Printf (printf) -data OutFormat = EigenstratFormat FilePath | PlinkFormat FilePath PlinkPopNameMode | FreqSumFormat deriving (Show)+data OutFormat = EigenstratFormat FilePath Bool+ | PlinkFormat FilePath PlinkPopNameMode Bool+ | VCFformat+ | FreqSumFormat deriving (Show) data ProgOpt = ProgOpt- CallingMode -- optCallingMode- Bool -- optKeepInCongruentReads- (Maybe Int) -- optSeed- Int -- optMinDepth- TransitionsMode -- optTransitionsMode- FilePath -- optSnpFile- OutFormat -- optOutFormat- (Either [String] FilePath) -- optSampleNames- String -- optPopName-+ CallingMode --optCallingMode+ Bool --optKeepIncongruentReads+ (Maybe Int) --optSeed+ Int --optMinDepth+ TransitionsMode --optTransitionsMode+ FilePath --optSnpFile+ OutFormat --optOutFormat+ (Either [String] FilePath) --optSampleNames+ (Either String [String]) --optPopName data ReadStats = ReadStats {- rsTotalSites :: IORef Int,+ rsTotalSites :: IORef Int, rsNonMissingSites :: V.IOVector Int,- rsRawReads :: V.IOVector Int,- rsReadsCleanedSS :: V.IOVector Int,- rsReadsCongruent :: V.IOVector Int+ rsRawReads :: V.IOVector Int,+ rsReadsCleanedSS :: V.IOVector Int,+ rsReadsCongruent :: V.IOVector Int } data Env = Env {- envCallingMode :: CallingMode,+ envCallingMode :: CallingMode, envKeepInCongruentReads :: Bool,- envMinDepth :: Int,- envTransitionsMode :: TransitionsMode,- envOutFormat :: OutFormat,- envSnpFile :: FilePath,- envSampleNames :: [String],- envPopName :: String,- envStats :: ReadStats+ envMinDepth :: Int,+ envTransitionsMode :: TransitionsMode,+ envOutFormat :: OutFormat,+ envSnpFile :: FilePath,+ envSampleNames :: [String],+ envPopName :: Either String [String],+ envVersion :: Version,+ envStats :: ReadStats } instance Show Env where- show (Env m r d t o sf sn pn _) = show (m, r, d, t, o, sf, sn, pn)+ show e = show (+ envCallingMode e,+ envKeepInCongruentReads e,+ envMinDepth e,+ envTransitionsMode e,+ envOutFormat e,+ envSnpFile e,+ envSampleNames e,+ envPopName e,+ envVersion e) type App = ReaderT Env (SafeT IO) @@ -70,7 +113,9 @@ main = do args <- OP.execParser parserInfo env <- initialiseEnvironment args- runSafeT $ runReaderT runMain env+ let handler = \(SeqFormatException msg) -> do+ throwIO $ SeqFormatException (take 200 msg)+ catch (runSafeT $ runReaderT runMain env) handler parserInfo :: OP.ParserInfo ProgOpt parserInfo = OP.info (pure (.) <*> versionInfoOpt <*> OP.helper <*> argParser)@@ -150,9 +195,13 @@ \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) <|>- (PlinkFormat <$> parsePlinkPrefix <*> parsePlinkPopMode) <|>- pure FreqSumFormat+ parseFormat = (EigenstratFormat <$> parseEigenstratPrefix <*> parseZipOut) <|>+ (PlinkFormat <$> parsePlinkPrefix <*> parsePlinkPopMode <*> parseZipOut) <|>+ parseVCFformat <|> pure FreqSumFormat+ parseZipOut = OP.switch (OP.long "--zip" <> OP.short 'z' <>+ OP.help "GZip the output Eigenstrat or Plink genotype and SNP files. \+ \Filenames will be appended with '.gz'. To zip FreqSum or VCF output, just zip the standard output of this\+ \ program, for example `pileupCaller ... --vcf | gzip -c > out.vcf.gz") 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 \@@ -174,6 +223,7 @@ parsePlinkPopBoth = OP.flag' PlinkPopNameAsBoth (OP.long "popNameAsBoth" <> OP.help "Only valid for Plink Output: \ \Write the population name into both the first and last column of the fam file, so both as Family-ID and as a \ \Phenotype according to the Plink Spec. By default, the population name is specified only as the first column (family name in the Plink spec)")+ parseVCFformat = OP.flag' VCFformat (OP.long "vcf" <> OP.help "output VCF format to stdout") parseSampleNames = parseSampleNameList <|> parseSampleNameFile parseSampleNameList = OP.option (Left . splitOn "," <$> OP.str) (OP.long "sampleNames" <> OP.metavar "NAME1,NAME2,..." <>@@ -181,28 +231,27 @@ 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")- parsePopName = 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\+ parsePopName = OP.option (processPopNames . splitOn "," <$> OP.str) (OP.long "samplePopName" <> OP.value (Left "Unknown") <> OP.showDefault <>+ OP.metavar "POP(s)" <>+ OP.help "specify the population name(s) of the samples, which are included\ \ in the output *.ind.txt file in Eigenstrat output. This will be ignored if the output \- \format is not Eigenstrat")+ \format is not Eigenstrat. If a single name is given, it is applied to all samples, if multiple are given, their number must match the \+ \the number of samples")+ processPopNames names = if length names == 1 then Left (head names) else Right names 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+programHelpDoc = PP.vsep [part1, PP.enclose PP.line PP.line (PP.indent 4 samtoolsExample),+ part2, PP.line, PP.fillSep . Prettyprinter.Util.words . T.pack $ versionInfoText] where- part1 = PP.fillSep . map PP.text . words $+ part1 = PP.fillSep . Prettyprinter.Util.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 $+ samtoolsExample = PP.hang 4 . PP.fillSep . Prettyprinter.Util.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 $+ part2 = PP.fillSep . Prettyprinter.Util.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."@@ -212,37 +261,48 @@ let (ProgOpt callingMode keepInCongruentReads seed minDepth transitionsMode snpFile outFormat sampleNames popName) = args case seed of- Nothing -> return ()+ Nothing -> return () Just seed_ -> liftIO . setStdGen $ mkStdGen seed_ sampleNamesList <- case sampleNames of Left list -> return list- Right fn -> lines <$> readFile fn+ 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 popName readStats+ outFormat snpFile sampleNamesList popName version readStats where makeVec n = do- v <- V.new n + 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- popName <- asks envPopName+ popNameSpec <- asks envPopName+ n <- length <$> asks envSampleNames+ let popNames = case popNameSpec of+ Left singlePopName -> replicate n singlePopName+ Right p -> if length p /= n then error "number of specified populations must equal sample size" else p case outFormat of- FreqSumFormat -> outputFreqSum freqSumProducer- EigenstratFormat outPrefix -> outputEigenStratOrPlink outPrefix popName Nothing freqSumProducer- PlinkFormat outPrefix popNameMode -> outputEigenStratOrPlink outPrefix popName (Just popNameMode) freqSumProducer+ FreqSumFormat -> do+ freqSumProducer <- pileupToFreqSum pileupProducer+ outputFreqSum (freqSumProducer >-> P.map snd)+ EigenstratFormat outPrefix zipOut -> do+ freqSumProducer <- pileupToFreqSum pileupProducer+ outputEigenStratOrPlink outPrefix zipOut popNames Nothing (freqSumProducer >-> P.map snd)+ PlinkFormat outPrefix popNameMode zipOut -> do+ freqSumProducer <- pileupToFreqSum pileupProducer+ outputEigenStratOrPlink outPrefix zipOut popNames (Just popNameMode) (freqSumProducer >-> P.map snd)+ VCFformat -> do+ freqSumProducer <- pileupToFreqSum pileupProducer+ outputVCF popNames freqSumProducer outputStats -pileupToFreqSum :: FilePath -> Producer PileupRow (SafeT IO) () ->- App (Producer FreqSumEntry (SafeT IO) ())-pileupToFreqSum snpFileName pileupProducer = do+pileupToFreqSum :: Producer PileupRow (SafeT IO) () -> App (Producer (Maybe PileupRow, FreqSumEntry) (SafeT IO) ())+pileupToFreqSum pileupProducer = do+ snpFileName <- asks envSnpFile nrSamples <- length <$> asks envSampleNames let snpProdOrderChecked = readEigenstratSnpFile snpFileName >-> orderCheckPipe cmpSnpPos@@ -252,17 +312,17 @@ orderedZip cmpSnpToPileupPos snpProdOrderChecked pileupProdOrderChecked mode <- asks envCallingMode keepInCongruentReads <- asks envKeepInCongruentReads- transisitionsMode <- asks envTransitionsMode- let singleStrandMode = (transisitionsMode == SingleStrandMode)+ transitionsMode <- asks envTransitionsMode+ let singleStrandMode = (transitionsMode == SingleStrandMode) minDepth <- asks envMinDepth readStats <- asks envStats let ret = Pipes.for jointProd $ \jointEntry -> case jointEntry of (Just esEntry, Nothing) -> do let (EigenstratSnpEntry chr pos gpos id_ ref alt) = esEntry- dosages = (replicate nrSamples Nothing) + dosages = (replicate nrSamples Nothing) liftIO $ addOneSite readStats- yield $ FreqSumEntry chr pos (Just id_) (Just gpos) ref alt dosages+ yield $ (Nothing, FreqSumEntry chr pos (Just id_) (Just gpos) ref alt dosages) (Just esEntry, Just pRow) -> do let (EigenstratSnpEntry chr pos gpos id_ ref alt) = esEntry (PileupRow _ _ _ rawPileupBasesPerSample rawStrandInfoPerSample) = pRow@@ -270,7 +330,7 @@ if singleStrandMode then cleanSSdamageAllSamples ref alt rawPileupBasesPerSample rawStrandInfoPerSample else rawPileupBasesPerSample- let congruentBasesPerSample = + let congruentBasesPerSample = if keepInCongruentReads then cleanBasesPerSample else map (filter (\c -> c == ref || c == alt)) cleanBasesPerSample@@ -279,9 +339,9 @@ (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 id_) (Just gpos) ref alt genotypes)+ yield $ (Just pRow, FreqSumEntry chr pos (Just id_) (Just gpos) ref alt genotypes) _ -> return ()- return (fst <$> ret)+ return $ (fst <$> ret) >-> filterTransitions transitionsMode where cmpSnpPos :: EigenstratSnpEntry -> EigenstratSnpEntry -> Ordering cmpSnpPos es1 es2 = (snpChrom es1, snpPos es1) `compare` (snpChrom es2, snpPos es2)@@ -295,6 +355,11 @@ updateStatsAllSamples :: ReadStats -> [Int] -> [Int] -> [Int] -> IO () updateStatsAllSamples readStats rawBaseCounts damageCleanedBaseCounts congruencyCleanedBaseCounts = do+ let nSamples = V.length (rsRawReads readStats)+ when (length rawBaseCounts /= nSamples) . throwIO . UserInputException $+ "number of individuals specified (" ++ show nSamples +++ ") differs from number of individuals in the pileup input (" +++ show (length rawBaseCounts) ++ ")" 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]@@ -303,38 +368,104 @@ 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+ MajorityCalling _ -> 1 :: Int+ RandomCalling -> 1 RandomDiploidCalling -> 2 let header' = FreqSumHeader sampleNames [nrHaplotypes | _ <- sampleNames]- outProd = freqSumProducer >-> filterTransitions transitionsOnly- lift . runEffect $ outProd >-> printFreqSumStdOut header'+ lift . runEffect $ freqSumProducer >-> printFreqSumStdOut header' -outputEigenStratOrPlink :: FilePath -> String -> Maybe PlinkPopNameMode -> Producer FreqSumEntry (SafeT IO) () -> App ()-outputEigenStratOrPlink outPrefix popName maybePlinkPopMode freqSumProducer = do- transitionsMode <- asks envTransitionsMode+outputEigenStratOrPlink :: FilePath -> Bool -> [String] -> Maybe PlinkPopNameMode -> Producer FreqSumEntry (SafeT IO) () -> App ()+outputEigenStratOrPlink outPrefix zipOut popNames maybePlinkPopMode freqSumProducer = do sampleNames <- asks envSampleNames- callingMode <- asks envCallingMode- let diploidizeCall = case callingMode of- RandomCalling -> True- MajorityCalling _ -> True- RandomDiploidCalling -> False- let (snpOut, indOut, genoOut) = case maybePlinkPopMode of- Just _ -> (outPrefix <> ".bim", outPrefix <> ".fam", outPrefix <> ".bed")- Nothing -> (outPrefix <> ".snp", outPrefix <> ".ind", outPrefix <> ".geno")- let indEntries = [EigenstratIndEntry n Unknown popName | n <- sampleNames]+ let (snpOut, indOut, genoOut) = case (maybePlinkPopMode, zipOut) of+ (Just _, False) -> (outPrefix <> ".bim", outPrefix <> ".fam", outPrefix <> ".bed")+ (Just _, True) -> (outPrefix <> ".bim.gz", outPrefix <> ".fam", outPrefix <> ".bed.gz")+ (Nothing, False) -> (outPrefix <> ".snp", outPrefix <> ".ind", outPrefix <> ".geno")+ (Nothing, True) -> (outPrefix <> ".snp.gz", outPrefix <> ".ind", outPrefix <> ".geno.gz")+ let indEntries = [EigenstratIndEntry (B.pack n) Unknown (B.pack p) | (n, p) <- zip sampleNames popNames] let writeFunc = case maybePlinkPopMode of Nothing -> (\g s i -> writeEigenstrat g s i indEntries) Just popMode -> let famEntries = map (eigenstratInd2PlinkFam popMode) indEntries in (\g s i -> writePlink g s i famEntries)- lift . runEffect $ freqSumProducer >-> filterTransitions transitionsMode >->- P.map (freqSumToEigenstrat diploidizeCall) >->- writeFunc genoOut snpOut indOut+ lift . runEffect $ freqSumProducer >-> P.map freqSumToEigenstrat >-> writeFunc genoOut snpOut indOut++outputVCF :: [String] -> Producer (Maybe PileupRow, FreqSumEntry) (SafeT IO) () -> App ()+outputVCF popNames freqSumProd = do+ sampleNames <- map B.pack <$> asks envSampleNames+ ver <- asks envVersion+ prog_name <- liftIO getProgName+ prog_args <- liftIO getArgs+ let command_line = prog_name ++ " " ++ intercalate " " prog_args+ let metaInfoLines = map B.pack [+ "##fileformat=VCFv4.2",+ "##source=pileupCaller_v" ++ showVersion ver,+ "##command_line=" ++ command_line,+ "##group_names=" ++ intercalate "," popNames,+ "##INFO=<ID=NS,Number=1,Type=Integer,Description=\"Number of Samples With Data\">",+ "##INFO=<ID=DP,Number=1,Type=Integer,Description=\"Total Depth\">",+ "##INFO=<ID=AF,Number=A,Type=Float,Description=\"Allele Frequency\">",+ "##FILTER=<ID=s50,Description=\"Less than 50% of samples have data\">",+ "##FILTER=<ID=s10,Description=\"Less than 10% of samples have data\">",+ "##FORMAT=<ID=GT,Number=1,Type=String,Description=\"Genotype\">",+ "##FORMAT=<ID=DP,Number=1,Type=Integer,Description=\"Read Depth\">",+ "##FORMAT=<ID=DP8,Number=8,Type=Integer,Description=\"Nr of Reads supporting A,C,G,T in forward strand, followed by the same quartet in reverse strand\">"]+ vcfh = VCFheader metaInfoLines sampleNames+ lift . runEffect $ freqSumProd >-> P.map (\(mpr, fse) -> createVcfEntry mpr fse) >-> printVCFtoStdOut vcfh++createVcfEntry :: Maybe PileupRow -> FreqSumEntry -> VCFentry+createVcfEntry maybePileupRow (FreqSumEntry chrom pos maybeSnpId _ ref alt calls) =+ VCFentry chrom pos maybeSnpId (B.pack [ref]) [B.pack [alt]] Nothing (Just filterString) infoFields genotypeInfos+ where+ nrMissing = length . filter (==Nothing) $ calls+ nrSamples = length calls+ filterString =+ if nrMissing * 10 > 9 * nrSamples then "s10;s50"+ else if nrMissing * 2 > nrSamples then "s50"+ else "PASS"+ totalDepth = case maybePileupRow of+ Nothing -> 0+ Just pr -> sum . map length . pileupBases $ pr+ nrSamplesWithData = nrSamples - nrMissing+ alleleFreq = computeAlleleFreq calls+ infoFields = [+ B.pack $ "NS=" ++ show nrSamplesWithData,+ B.pack $ "DP=" ++ show totalDepth] ++ + case alleleFreq of+ Just f ->+ let roundedFreq = fromIntegral (round (f * 100.0) :: Int) / 100.0 :: Double+ in [B.pack $ "AF=" ++ show roundedFreq]+ Nothing -> []+ formatField = case maybePileupRow of+ Nothing -> ["GT"]+ Just _ -> ["GT", "DP", "DP8"]+ genotypeFields = do -- list monad over samples+ i <- [0 .. (nrSamples - 1)]+ let ca = calls !! i+ let gt = case ca of+ Nothing -> "."+ Just (0, 1) -> "0"+ Just (1, 1) -> "1"+ Just (0, 2) -> "0/0"+ Just (1, 2) -> "0/1"+ Just (2, 2) -> "1/1"+ _ -> error "should never happen"+ case maybePileupRow of+ Nothing -> return [gt]+ Just pr -> do+ let bases = pileupBases pr !! i+ let strands = pileupStrandInfo pr !! i+ let dp = length bases+ let dp8 = do -- list monad+ strand <- [ForwardStrand, ReverseStrand] -- outer loop+ allele <- ['A', 'C', 'G', 'T'] -- inner loop+ return . show . length . filter (\(a, s) -> a == allele && s == strand) $ zip bases strands+ return [gt, B.pack $ show dp, B.pack $ intercalate "," dp8]+ genotypeInfos = Just (formatField, genotypeFields)+ outputStats :: App () outputStats = do
src-executables/vcf2eigenstrat.hs view
@@ -1,23 +1,29 @@ {-# LANGUAGE OverloadedStrings #-} -import Pipes.OrderedZip (orderedZip)-import SequenceFormats.VCF (readVCFfromStdIn, VCFheader(..), VCFentry(..),- isBiallelicSnp, getDosages, vcfToFreqSumEntry)-import SequenceFormats.Eigenstrat (EigenstratSnpEntry(..), readEigenstratSnpFile, writeEigenstrat,- Sex(..), EigenstratIndEntry(..))-import SequenceFormats.FreqSum (FreqSumEntry(..), freqSumEntryToText)+import Pipes.OrderedZip (orderedZip)+import SequenceFormats.Eigenstrat (EigenstratIndEntry (..),+ EigenstratSnpEntry (..), Sex (..),+ readEigenstratSnpFile,+ writeEigenstrat)+import SequenceFormats.FreqSum (FreqSumEntry (..))+import SequenceFormats.VCF (VCFentry (..), VCFheader (..),+ getDosages, isBiallelicSnp,+ readVCFfromStdIn,+ vcfToFreqSumEntry) -import SequenceTools.Utils (versionInfoText, versionInfoOpt, freqSumToEigenstrat)+import SequenceTools.Utils (freqSumToEigenstrat,+ versionInfoOpt, versionInfoText) -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 Control.Exception.Base (AssertionFailed (..), throwIO)+import Control.Monad (when)+import Control.Monad.IO.Class (MonadIO, liftIO)+import qualified Data.ByteString.Char8 as B -- import Debug.Trace (trace)-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)+import qualified Options.Applicative as OP+import Pipes (Pipe, Producer, cat, for,+ runEffect, yield, (>->))+import qualified Pipes.Prelude as P+import Pipes.Safe (MonadSafe, runSafeT) -- snpPosFile outPrefix data ProgOpt = ProgOpt (Maybe FilePath) FilePath@@ -38,22 +44,19 @@ (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.")+ \and in the VCF will be output. Without this option, all sites in the VCF will be output") 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")+ OP.help "specify the prefix for the EigenStrat output files. Output files will then be \+ \<FILE_PREFIX>.geno, <FILE_PREFIX>.snp and <FILE_PREFIX>.ind") 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"+ let snpOut = outPrefix ++ ".snp"+ indOut = outPrefix ++ ".ind"+ genoOut = outPrefix ++ ".geno" VCFheader _ sampleNames = vcfHeader nrInds = length sampleNames indEntries = [EigenstratIndEntry n Unknown "Unknown" | n <- sampleNames]@@ -62,7 +65,7 @@ Just snpPosFile -> return $ runJointly vcfBodyBiAllelic nrInds snpPosFile Nothing -> return $ runSimple vcfBodyBiAllelic- runEffect $ vcfProducer >-> P.map (freqSumToEigenstrat False) >-> writeEigenstrat genoOut snpOut indOut indEntries+ runEffect $ vcfProducer >-> P.map freqSumToEigenstrat >-> writeEigenstrat genoOut snpOut indOut indEntries runJointly :: (MonadIO m, MonadSafe m) => Producer VCFentry m r -> Int -> FilePath -> Producer FreqSumEntry m r runJointly vcfBody nrInds snpPosFile =@@ -79,9 +82,7 @@ let dosages = replicate nrInds Nothing yield $ FreqSumEntry snpChrom' snpPos' (Just snpId') (Just gpos) snpRef' snpAlt' dosages (Just (EigenstratSnpEntry snpChrom' snpPos' gpos snpId' snpRef' snpAlt'), Just vcfEntry) -> do- dosages <- case getDosages vcfEntry of- Right dos -> return dos- Left err -> liftIO . throwIO $ AssertionFailed err+ dosages <- liftIO $ getDosages vcfEntry when (length dosages /= nrInds) $ (liftIO . throwIO) (AssertionFailed "inconsistent number of genotypes.") let normalizedDosages = case vcfAlt vcfEntry of@@ -96,15 +97,14 @@ _ -> return () where flipDosages dos = case dos of- Just 0 -> Just 2- Just 1 -> Just 1- Just 2 -> Just 0- _ -> Nothing+ Just (0, 2) -> Just (2, 2)+ Just (1, 2) -> Just (1, 2)+ Just (2, 2) -> Just (0, 2)+ Just (0, 1) -> Just (1, 1)+ Just (1, 1) -> Just (0, 1)+ _ -> 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)+ fs <- liftIO $ vcfToFreqSumEntry e+ yield fs
src/SequenceTools/PileupCaller.hs view
@@ -2,15 +2,16 @@ module SequenceTools.PileupCaller (callToDosage, Call(..), callGenotypeFromPileup, callMajorityAllele, findMajorityAlleles, callRandomAllele, callRandomDiploid, CallingMode(..),- TransitionsMode(..), filterTransitions, cleanSSdamageAllSamples) where+ TransitionsMode(..), filterTransitions, cleanSSdamageAllSamples,+ computeAlleleFreq) where -import SequenceFormats.FreqSum (FreqSumEntry(..))-import SequenceFormats.Pileup (Strand(..))-import SequenceTools.Utils (sampleWithoutReplacement)+import SequenceFormats.FreqSum (FreqSumEntry (..))+import SequenceFormats.Pileup (PileupRow (..), Strand (..))+import SequenceTools.Utils (sampleWithoutReplacement) -import Data.List (sortOn, group, sort)-import Pipes (Pipe, cat)-import qualified Pipes.Prelude as P+import Data.List (group, sort, sortOn)+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)@@ -21,19 +22,19 @@ data TransitionsMode = TransitionsMissing | SkipTransitions | SingleStrandMode | AllSites deriving (Eq, Show) -- |a function to turn a call into the dosage of non-reference alleles-callToDosage :: Char -> Char -> Call -> Maybe Int+callToDosage :: Char -> Char -> Call -> Maybe (Int, Int) callToDosage refA altA call = case call of- HaploidCall a | a == refA -> Just 0- | a == altA -> Just 1+ HaploidCall a | a == refA -> Just (0, 1)+ | a == altA -> Just (1, 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+ DiploidCall a1 a2 | (a1, a2) == (refA, refA) -> Just (0, 2)+ | (a1, a2) == (refA, altA) -> Just (1, 2)+ | (a1, a2) == (altA, refA) -> Just (1, 2)+ | (a1, a2) == (altA, altA) -> Just (2, 2) | otherwise -> Nothing MissingCall -> Nothing --- |Make a call from alleles +-- |Make a call from alleles callGenotypeFromPileup :: CallingMode -> Int -> String -> IO Call callGenotypeFromPileup mode minDepth alleles = if length alleles < minDepth then return MissingCall else@@ -58,7 +59,7 @@ r <- sampleWithoutReplacement listA 1 case r of Just [r'] -> return r'- _ -> error "should not happen"+ _ -> error "should not happen" return $ HaploidCall a -- |Find the majority allele(s)@@ -73,28 +74,31 @@ callRandomAllele alleles = do res <- sampleWithoutReplacement alleles 1 case res of- Nothing -> return MissingCall+ Nothing -> return MissingCall Just [a] -> return $ HaploidCall a- _ -> error "should not happen"+ _ -> 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+ Nothing -> return MissingCall Just [a1, a2] -> return $ DiploidCall a1 a2- _ -> error "should not happen"+ _ -> error "should not happen" -filterTransitions :: (Monad m) => TransitionsMode -> Pipe FreqSumEntry FreqSumEntry m ()+-- the basic information stream is a tuple of a PileupRow (if data is present at a SNP), and a FreqSumEntry that contains the calls.+-- For Eigenstrat and Plink we don't need the PileupRow, but for VCF, we can store additional information beyond the mere calls,+-- that's why we're streaming both, to have an output-agnostic stream.+filterTransitions :: (Monad m) => TransitionsMode -> Pipe (Maybe PileupRow, FreqSumEntry) (Maybe PileupRow, FreqSumEntry) m () filterTransitions transversionsMode = case transversionsMode of SkipTransitions ->- P.filter (\(FreqSumEntry _ _ _ _ ref alt _) -> isTransversion ref alt)+ P.filter (\(_, FreqSumEntry _ _ _ _ ref alt _) -> isTransversion ref alt) TransitionsMissing ->- P.map (\(FreqSumEntry chrom pos id_ gpos ref alt calls) ->+ P.map (\(mp, FreqSumEntry chrom pos id_ gpos ref alt calls) -> let calls' = if isTransversion ref alt then calls else [Nothing | _ <- calls]- in FreqSumEntry chrom pos id_ gpos ref alt calls')+ in (mp, FreqSumEntry chrom pos id_ gpos ref alt calls')) _ -> cat where isTransversion ref alt = not $ isTransition ref alt@@ -118,3 +122,10 @@ removeReads :: Strand -> String -> [Strand] -> String removeReads strand bases strands = [b | (b, s) <- zip bases strands, s /= strand]++computeAlleleFreq :: [Maybe (Int, Int)] -> Maybe Double+computeAlleleFreq dosages =+ let nrTotalAlleles = sum . map (maybe 0 snd) $ dosages+ nrNonRefAlleles = sum . map (maybe 0 fst) $ dosages+ in if nrTotalAlleles == 0 then Nothing else+ Just (fromIntegral nrNonRefAlleles / fromIntegral nrTotalAlleles)
src/SequenceTools/Utils.hs view
@@ -1,18 +1,23 @@ {-# LANGUAGE OverloadedStrings #-} module SequenceTools.Utils (versionInfoOpt, versionInfoText, sampleWithoutReplacement,- freqSumToEigenstrat, dosageToEigenstratGeno) where + freqSumToEigenstrat, dosageToEigenstratGeno, UserInputException(..)) where -import SequenceFormats.FreqSum (FreqSumEntry(..))-import SequenceFormats.Eigenstrat (EigenstratSnpEntry(..), GenoLine, GenoEntry(..))-import SequenceFormats.Utils (Chrom(..))+import SequenceFormats.Eigenstrat (EigenstratSnpEntry (..),+ GenoEntry (..), GenoLine)+import SequenceFormats.FreqSum (FreqSumEntry (..))+import SequenceFormats.Utils (Chrom (..)) -import qualified Data.ByteString.Char8 as B-import Data.Vector (fromList)-import Data.Version (showVersion)-import qualified Options.Applicative as OP-import Paths_sequenceTools (version)-import System.Random (randomRIO)+import Control.Exception (Exception)+import qualified Data.ByteString.Char8 as B+import Data.Vector (fromList)+import Data.Version (showVersion)+import qualified Options.Applicative as OP+import Paths_sequenceTools (version)+import System.Random (randomRIO) +data UserInputException = UserInputException String deriving (Show)+instance Exception UserInputException+ versionInfoOpt :: OP.Parser (a -> a) versionInfoOpt = OP.infoOption (showVersion version) (OP.long "version" <> OP.help "Print version and exit") @@ -34,31 +39,24 @@ remove i xs = let (ys, zs) = splitAt i xs in ys ++ tail zs -- |convert a freqSum entry to an eigenstrat SNP entry-freqSumToEigenstrat :: Bool -> FreqSumEntry -> (EigenstratSnpEntry, GenoLine)-freqSumToEigenstrat diploidizeCall (FreqSumEntry chrom@(Chrom c) pos maybeSnpId maybeGeneticPos ref alt calls) =- let snpId_ = case maybeSnpId of +freqSumToEigenstrat :: FreqSumEntry -> (EigenstratSnpEntry, GenoLine)+freqSumToEigenstrat (FreqSumEntry chrom@(Chrom c) pos maybeSnpId maybeGeneticPos ref alt calls) =+ let snpId_ = case maybeSnpId of Just id_ -> id_- Nothing -> c <> "_" <> B.pack (show pos)+ Nothing -> c <> "_" <> B.pack (show pos) geneticPos = case maybeGeneticPos of- Just p -> p+ Just p -> p Nothing -> 0.0 snpEntry = EigenstratSnpEntry chrom pos geneticPos snpId_ ref alt- geno = fromList . map (dosageToEigenstratGeno diploidizeCall) $ calls+ geno = fromList . map dosageToEigenstratGeno $ 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)+dosageToEigenstratGeno :: Maybe (Int, Int) -> GenoEntry+dosageToEigenstratGeno Nothing = Missing+dosageToEigenstratGeno (Just (0, 1)) = HomRef+dosageToEigenstratGeno (Just (1, 1)) = HomAlt+dosageToEigenstratGeno (Just (0, 2)) = HomRef+dosageToEigenstratGeno (Just (1, 2)) = Het+dosageToEigenstratGeno (Just (2, 2)) = HomAlt+dosageToEigenstratGeno c = error ("unknown genotype " ++ show c)
test/SequenceTools/PileupCallerSpec.hs view
@@ -1,23 +1,31 @@ {-# LANGUAGE OverloadedStrings #-} module SequenceTools.PileupCallerSpec (spec) where -import SequenceTools.Utils (dosageToEigenstratGeno, freqSumToEigenstrat)-import SequenceTools.PileupCaller (callToDosage, Call(..), callGenotypeFromPileup,- callMajorityAllele, findMajorityAlleles, callRandomAllele,- callRandomDiploid, CallingMode(..),- TransitionsMode(..), filterTransitions, cleanSSdamageAllSamples)+import SequenceFormats.Pileup (PileupRow)+import SequenceTools.PileupCaller (Call (..), CallingMode (..),+ TransitionsMode (..),+ callGenotypeFromPileup,+ callMajorityAllele,+ callRandomAllele,+ callRandomDiploid, callToDosage,+ cleanSSdamageAllSamples,+ filterTransitions,+ findMajorityAlleles)+import SequenceTools.Utils (dosageToEigenstratGeno,+ freqSumToEigenstrat) -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 SequenceFormats.FreqSum (FreqSumEntry(..))-import SequenceFormats.Eigenstrat (GenoEntry(..), EigenstratSnpEntry(..))-import SequenceFormats.Utils (Chrom(..))-import SequenceFormats.Pileup (Strand(..))-import Test.Hspec+import Control.Monad (forM_, replicateM)+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 SequenceFormats.Eigenstrat (EigenstratSnpEntry (..),+ GenoEntry (..))+import SequenceFormats.FreqSum (FreqSumEntry (..))+import SequenceFormats.Pileup (Strand (..))+import SequenceFormats.Utils (Chrom (..))+import Test.Hspec spec :: Spec spec = do@@ -39,18 +47,18 @@ 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+ callToDosage 'A' 'C' (HaploidCall 'A') `shouldBe` Just (0, 1) it "should return 1 for haploid alt call" $- callToDosage 'A' 'C' (HaploidCall 'C') `shouldBe` Just 1 + callToDosage 'A' 'C' (HaploidCall 'C') `shouldBe` Just (1, 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+ callToDosage 'A' 'C' (DiploidCall 'A' 'A') `shouldBe` Just (0, 2) 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+ callToDosage 'A' 'C' (DiploidCall 'A' 'C') `shouldBe` Just (1, 2)+ callToDosage 'A' 'C' (DiploidCall 'C' 'A') `shouldBe` Just (1, 2) it "should return 2 for diploid hom-alt call" $- callToDosage 'A' 'C' (DiploidCall 'C' 'C') `shouldBe` Just 2+ callToDosage 'A' 'C' (DiploidCall 'C' 'C') `shouldBe` Just (2, 2) testCallGenotypeFromPileup :: Spec@@ -59,8 +67,8 @@ 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" $@@ -114,55 +122,62 @@ testDosageToEigenstratGeno :: Spec testDosageToEigenstratGeno = describe "dosageToEigenstratGeno" $ do it "should give Hom-Ref for 0 pseudo-haploid" $- dosageToEigenstratGeno True (Just 0) `shouldBe` HomRef+ dosageToEigenstratGeno (Just (0, 1)) `shouldBe` HomRef it "should give Hom-Alt for 1 pseudo-haploid" $- dosageToEigenstratGeno True (Just 1) `shouldBe` HomAlt+ dosageToEigenstratGeno (Just (1, 1)) `shouldBe` HomAlt it "should give Missing for Nothing pseudo-haploid" $- dosageToEigenstratGeno True Nothing `shouldBe` Missing+ dosageToEigenstratGeno Nothing `shouldBe` Missing it "should give Hom-Ref for 0 diploid" $- dosageToEigenstratGeno False (Just 0) `shouldBe` HomRef+ dosageToEigenstratGeno (Just (0, 2)) `shouldBe` HomRef it "should give Het for 1 diploid" $- dosageToEigenstratGeno False (Just 1) `shouldBe` Het+ dosageToEigenstratGeno (Just (1, 2)) `shouldBe` Het it "should give Hom-Alt for 2 diploid" $- dosageToEigenstratGeno False (Just 2) `shouldBe` HomAlt+ dosageToEigenstratGeno (Just (2, 2)) `shouldBe` HomAlt it "should give Missing for Nothing diploid" $- dosageToEigenstratGeno False Nothing `shouldBe` Missing+ dosageToEigenstratGeno Nothing `shouldBe` Missing testFreqSumToEigenstrat :: Spec testFreqSumToEigenstrat = describe "freqSumtoEigenstrat" $ do- let fs = FreqSumEntry (Chrom "1") 1000 Nothing Nothing 'A' 'C' [Just 0, Just 1, Just 1, Nothing, Just 0]+ let fs = FreqSumEntry (Chrom "1") 1000 Nothing Nothing 'A' 'C' [Just (0, 1), Just (1, 1), Just (1, 1), Nothing, Just (0, 1)] 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)+ freqSumToEigenstrat 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)+ freqSumToEigenstrat (fs {fsSnpId = Just "rs123"}) `shouldBe` (es {snpId = "rs123"}, genoLine) -mockFreqSumData :: [FreqSumEntry]+mockFreqSumData :: [(Maybe PileupRow, FreqSumEntry)] mockFreqSumData = [- FreqSumEntry (Chrom "1") 1000 (Just "rs1") Nothing 'A' 'C' [Just 1, Just 2, Nothing, Just 0, Just 0],- FreqSumEntry (Chrom "1") 2000 (Just "rs2") Nothing 'C' 'T' [Just 1, Just 2, Nothing, Just 0, Just 0],- FreqSumEntry (Chrom "1") 3000 (Just "rs3") Nothing 'A' 'G' [Just 1, Just 2, Nothing, Just 0, Just 0],- FreqSumEntry (Chrom "2") 1000 (Just "rs4") Nothing 'A' 'G' [Just 1, Just 2, Nothing, Just 0, Just 0],- FreqSumEntry (Chrom "2") 2000 (Just "rs5") Nothing 'T' 'A' [Just 1, Just 2, Nothing, Just 0, Just 0],- FreqSumEntry (Chrom "2") 3000 (Just "rs6") Nothing 'T' 'C' [Just 1, Just 2, Nothing, Just 0, Just 0]]+ (Nothing, FreqSumEntry (Chrom "1") 1000 (Just "rs1") Nothing 'A' 'C' [Just (1, 2), Just (2, 2), Nothing, Just (0, 2), Just (0, 2)]),+ (Nothing, FreqSumEntry (Chrom "1") 2000 (Just "rs2") Nothing 'C' 'T' [Just (1, 2), Just (2, 2), Nothing, Just (0, 2), Just (0, 2)]),+ (Nothing, FreqSumEntry (Chrom "1") 3000 (Just "rs3") Nothing 'A' 'G' [Just (1, 2), Just (2, 2), Nothing, Just (0, 2), Just (0, 2)]),+ (Nothing, FreqSumEntry (Chrom "2") 1000 (Just "rs4") Nothing 'A' 'G' [Just (1, 2), Just (2, 2), Nothing, Just (0, 2), Just (0, 2)]),+ (Nothing, FreqSumEntry (Chrom "2") 2000 (Just "rs5") Nothing 'T' 'A' [Just (1, 2), Just (2, 2), Nothing, Just (0, 2), Just (0, 2)]),+ (Nothing, FreqSumEntry (Chrom "2") 3000 (Just "rs6") Nothing 'T' 'C' [Just (1, 2), Just (2, 2), Nothing, Just (0, 2), Just (0, 2)])] 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") Nothing 'A' 'C' [Just 1, Just 2, Nothing, Just 0, Just 0],- FreqSumEntry (Chrom "2") 2000 (Just "rs5") Nothing 'T' 'A' [Just 1, Just 2, Nothing, Just 0, Just 0]]+ let r = P.toList $ each mockFreqSumData >->+ filterTransitions SkipTransitions+ map snd r `shouldBe` [+ FreqSumEntry (Chrom "1") 1000 (Just "rs1") Nothing 'A' 'C'+ [Just (1, 2), Just (2, 2), Nothing, Just (0, 2), Just (0, 2)],+ FreqSumEntry (Chrom "2") 2000 (Just "rs5") Nothing 'T' 'A'+ [Just (1, 2), Just (2, 2), Nothing, Just (0, 2), Just (0, 2)]] 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") Nothing 'A' 'C' [Just 1, Just 2, Nothing, Just 0, Just 0],+ let r = P.toList $+ each mockFreqSumData >->+ filterTransitions TransitionsMissing+ map snd r `shouldBe` [+ FreqSumEntry (Chrom "1") 1000 (Just "rs1") Nothing 'A' 'C'+ [Just (1, 2), Just (2, 2), Nothing, Just (0, 2), Just (0, 2)], FreqSumEntry (Chrom "1") 2000 (Just "rs2") Nothing 'C' 'T' [Nothing, Nothing, Nothing, Nothing, Nothing], FreqSumEntry (Chrom "1") 3000 (Just "rs3") Nothing 'A' 'G' [Nothing, Nothing, Nothing, Nothing, Nothing], FreqSumEntry (Chrom "2") 1000 (Just "rs4") Nothing 'A' 'G' [Nothing, Nothing, Nothing, Nothing, Nothing],- FreqSumEntry (Chrom "2") 2000 (Just "rs5") Nothing 'T' 'A' [Just 1, Just 2, Nothing, Just 0, Just 0],+ FreqSumEntry (Chrom "2") 2000 (Just "rs5") Nothing 'T' 'A'+ [Just (1, 2), Just (2, 2), Nothing, Just (0, 2), Just (0, 2)], FreqSumEntry (Chrom "2") 3000 (Just "rs6") Nothing 'T' 'C' [Nothing, Nothing, Nothing, Nothing, Nothing]] it "should output all sites with AllSites" $ do let r = P.toList $ each mockFreqSumData >-> filterTransitions AllSites
test/SequenceTools/UtilsSpec.hs view
@@ -1,10 +1,10 @@ module SequenceTools.UtilsSpec (spec) where -import SequenceTools.Utils (sampleWithoutReplacement)+import SequenceTools.Utils (sampleWithoutReplacement) -import Control.Monad (replicateM_)-import Data.List (sort, union, nub)-import Test.Hspec+import Control.Monad (replicateM_)+import Data.List (nub, sort, union)+import Test.Hspec spec :: Spec spec = testSampleWithoutReplacement