sequenceTools 1.5.2 → 1.5.3.1
raw patch · 4 files changed
+140/−152 lines, 4 filesdep ~ansi-wl-pprintdep ~basedep ~bytestringPVP ok
version bump matches the API change (PVP)
Dependency ranges changed: ansi-wl-pprint, base, bytestring, foldl, hspec, lens-family, optparse-applicative, pipes, pipes-group, pipes-ordered-zip, pipes-safe, random, sequence-formats, split, transformers, vector
API changes (from Hackage documentation)
Files
- Changelog.md +17/−24
- README.md +37/−10
- sequenceTools.cabal +50/−94
- src-executables/pileupCaller.hs +36/−24
Changelog.md view
@@ -1,26 +1,19 @@-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--V 1.4.0.3: Updated to new sequence-formats version, now including reading of genetic position from eigenstrat files.--V 1.4.0.4:-* Fixed eigenstrat-output in pileupCaller to add a dot after the outputprefix before the file extensions.-* Updated haskell-stack wrapper scripts for EIGENSOFT and ADMIXTOOLS.-* Moved unmaintained scripts into unmaintained folder.--V 1.5.0: Added support for Plink output+# Changelog -V 1.5.1: Added automatic building+- 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.+- V 1.5.1: Added automatic building+- V 1.5.0: Added support for Plink output+- V 1.4.0.4:+ * Fixed eigenstrat-output in pileupCaller to add a dot after the outputprefix before the file extensions.+ * Updated haskell-stack wrapper scripts for EIGENSOFT and ADMIXTOOLS.+ * Moved unmaintained scripts into unmaintained folder.+- V 1.4.0.3: Updated to new sequence-formats version, now including reading of genetic position from eigenstrat files.+- V 1.4.0.1: Improved README, fixed output bug in genoStats.hs+- V 1.4.0: Added single strand mode, and new triallelic treatment.+- 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.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.2.4: normaliseBimWithVCF is ready.+- 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.5.2: Fixed a bug with --samplePopName having to be entered after -p or -e. Fixed a bug in the sequence-formats dependency.
README.md view
@@ -1,24 +1,34 @@ # SequenceTools -[](https://anaconda.org/bioconda/sequencetools)+[Install with Bioconda](https://anaconda.org/bioconda/sequencetools) This repository contains some programs that I use for processing sequencing data. # Installation-## Simple Installation via stack and hackage -This installation installs the latest version that is up on [hackage](https://hackage.haskell.org/package/sequenceTools):+## Installation via precompiled executables -1. Download stack (https://docs.haskellstack.org/en/stable/README/#how-to-install<Paste>).-2. Run `stack install sequenceTools --resolver nightly`. You should now have the executables from this package under `~/.local/bin`.-3. 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.+* [Download the latest Executable](https://github.com/stschiff/sequenceTools/releases/latest) that best matches your platform. -## Installation from source via stack+For example, to install `pileupCaller` in Linux, you can run the following commands to get started: -1. Download stack (https://docs.haskellstack.org/en/stable/README/#how-to-install<Paste>).-2. Clone this repository via `git clone https://github.com/stschiff/sequenceTools.git`-3. Install via `cd sequenceTools; stack install` +```bash+# download the current stable release binary+wget https://github.com/stschiff/sequenceTools/releases/latest/download/pileupCaller-linux+# make it executable+chmod +x pileupCaller-linux+# run it+./trident-linux -h+``` +## Installation from source++1. Clone the repository: `git clone https://github.com/stschiff/sequenceTools.git`+2. Go into the repository: `cd sequenceTools`+3. Compile the executables in the repository using `stack`: `stack install`++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@@ -62,6 +72,9 @@ --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.+ + 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> \@@ -70,6 +83,20 @@ 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.++### 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.++There is an important catch: If you have data from paired-end sequencing, and you are using _unmerged_ reads, then this approach will fail, as it will then _not_ discard potentially damaged reads.++So there are two options if you have Paired-end sequencing data:+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.+ ## vcf2eigenstrat
sequenceTools.cabal view
@@ -1,103 +1,59 @@-cabal-version: >=1.10-name: sequenceTools-version: 1.5.2-license: GPL-3-license-file: LICENSE-maintainer: stephan.schiffels@mac.com-author: Stephan Schiffels-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+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 library- exposed-modules:- SequenceTools.Utils- SequenceTools.PileupCaller+ 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 - hs-source-dirs: src- other-modules: Paths_sequenceTools- default-language: Haskell2010- build-depends:- base >=4.7 && <5,- optparse-applicative >=0.15.1.0,- random >=1.1,- sequence-formats >=1.6.3,- bytestring >=0.10.12.0,- vector >=0.12.1.2,- pipes >=4.3.14+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 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.14.1.0,- sequenceTools -any,- sequence-formats >=1.6.3,- optparse-applicative >=0.15.1.0,- pipes >=4.3.14,- vector >=0.12.1.2,- random >=1.1,- bytestring >=0.10.12.0,- pipes-safe >=2.3.2,- pipes-ordered-zip >=1.1.0,- split >=0.2.3.4,- ansi-wl-pprint >=0.6.9,- transformers >=0.5.6.2+ 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 executable vcf2eigenstrat- main-is: vcf2eigenstrat.hs- hs-source-dirs: src-executables- default-language: Haskell2010- build-depends:- base >=4.14.1.0,- sequenceTools -any,- pipes-ordered-zip >=1.1.0,- sequence-formats >=1.6.3,- bytestring >=0.10.12.0,- vector >=0.12.1.2,- optparse-applicative >=0.15.1.0,- pipes >=4.3.14,- pipes-safe >=2.3.2+ 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 executable genoStats- main-is: genoStats.hs- hs-source-dirs: src-executables- default-language: Haskell2010- build-depends:- base >=4.14.1.0,- sequence-formats >=1.6.3,- sequenceTools -any,- foldl >=1.4.10,- bytestring >=0.10.12.0,- vector >=0.12.1.2,- lens-family >=2.0.0,- optparse-applicative >=0.15.1.0,- pipes >=4.3.14,- pipes-group >=1.0.12,- pipes-safe >=2.3.2--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.14.1.0,- hspec >=2.7.8,- sequenceTools -any,- sequence-formats >=1.6.3,- vector >=0.12.1.2,- bytestring >=0.10.12.0,- pipes >=4.3.14+ 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
src-executables/pileupCaller.hs view
@@ -7,7 +7,7 @@ import SequenceTools.Utils (versionInfoOpt, versionInfoText, freqSumToEigenstrat) import SequenceTools.PileupCaller (CallingMode(..), callGenotypeFromPileup, callToDosage, filterTransitions, TransitionsMode(..), cleanSSdamageAllSamples)-import SequenceFormats.Plink (writePlink)+import SequenceFormats.Plink (writePlink, PlinkPopNameMode(..), eigenstratInd2PlinkFam) import Control.Applicative ((<|>)) import Control.Monad.Trans.Reader (ReaderT, asks, runReaderT)@@ -27,20 +27,20 @@ import Text.Printf (printf) import qualified Text.PrettyPrint.ANSI.Leijen as PP -data OutFormat = EigenstratFormat FilePath | PlinkFormat FilePath | FreqSumFormat deriving (Show)+data OutFormat = EigenstratFormat FilePath | PlinkFormat FilePath PlinkPopNameMode | FreqSumFormat deriving (Show) -data ProgOpt = ProgOpt {- optCallingMode :: CallingMode,- optKeepInCongruentReads :: Bool,- optSeed :: Maybe Int,- optMinDepth :: Int,- optTransitionsMode :: TransitionsMode,- optSnpFile :: FilePath,- optOutFormat :: OutFormat,- optSampleNames :: Either [String] FilePath,- optPopName :: String-}+data ProgOpt = ProgOpt+ CallingMode -- optCallingMode+ Bool -- optKeepInCongruentReads+ (Maybe Int) -- optSeed+ Int -- optMinDepth+ TransitionsMode -- optTransitionsMode+ FilePath -- optSnpFile+ OutFormat -- optOutFormat+ (Either [String] FilePath) -- optSampleNames+ String -- optPopName + data ReadStats = ReadStats { rsTotalSites :: IORef Int, rsNonMissingSites :: V.IOVector Int,@@ -150,7 +150,9 @@ \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) <|> pure FreqSumFormat+ parseFormat = (EigenstratFormat <$> parseEigenstratPrefix) <|>+ (PlinkFormat <$> parsePlinkPrefix <*> parsePlinkPopMode) <|>+ 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 \@@ -165,6 +167,13 @@ \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")+ parsePlinkPopMode = parsePlinkPopPhenotype <|> parsePlinkPopBoth <|> pure PlinkPopNameAsFamily+ parsePlinkPopPhenotype = OP.flag' PlinkPopNameAsPhenotype (OP.long "popNameAsPhenotype" <> OP.help "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)")+ 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)") parseSampleNames = parseSampleNameList <|> parseSampleNameFile parseSampleNameList = OP.option (Left . splitOn "," <$> OP.str) (OP.long "sampleNames" <> OP.metavar "NAME1,NAME2,..." <>@@ -227,8 +236,8 @@ popName <- asks envPopName case outFormat of FreqSumFormat -> outputFreqSum freqSumProducer- EigenstratFormat outPrefix -> outputEigenStratOrPlink outPrefix popName False freqSumProducer- PlinkFormat outPrefix -> outputEigenStratOrPlink outPrefix popName True freqSumProducer+ EigenstratFormat outPrefix -> outputEigenStratOrPlink outPrefix popName Nothing freqSumProducer+ PlinkFormat outPrefix popNameMode -> outputEigenStratOrPlink outPrefix popName (Just popNameMode) freqSumProducer outputStats pileupToFreqSum :: FilePath -> Producer PileupRow (SafeT IO) () ->@@ -305,8 +314,8 @@ outProd = freqSumProducer >-> filterTransitions transitionsOnly lift . runEffect $ outProd >-> printFreqSumStdOut header' -outputEigenStratOrPlink :: FilePath -> String -> Bool -> Producer FreqSumEntry (SafeT IO) () -> App ()-outputEigenStratOrPlink outPrefix popName formatIsPlink freqSumProducer = do+outputEigenStratOrPlink :: FilePath -> String -> Maybe PlinkPopNameMode -> Producer FreqSumEntry (SafeT IO) () -> App ()+outputEigenStratOrPlink outPrefix popName maybePlinkPopMode freqSumProducer = do transitionsMode <- asks envTransitionsMode sampleNames <- asks envSampleNames callingMode <- asks envCallingMode@@ -314,15 +323,18 @@ RandomCalling -> True MajorityCalling _ -> True RandomDiploidCalling -> False- let [snpOut, indOut, genoOut] =- if formatIsPlink- then map (outPrefix <>) [".bim", ".fam", ".bed"]- else map (outPrefix <>) [".snp", ".ind", ".geno"]- let writeFunc = if formatIsPlink then writePlink else writeEigenstrat+ 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 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 indEntries+ writeFunc genoOut snpOut indOut outputStats :: App () outputStats = do