diff --git a/Changelog.md b/Changelog.md --- a/Changelog.md +++ b/Changelog.md @@ -18,3 +18,5 @@ * 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 diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -1,8 +1,6 @@ # SequenceTools -[![Anaconda package](https://anaconda.org/bioconda/sequencetools/badges/installer/conda.svg)](https://anaconda.org/bioconda/sequencetools) - -(bioconda package available thanks to [apeltzer](https://github.com/apeltzer)!) +[![Install with Bioconda](https://anaconda.org/bioconda/sequencetools/badges/installer/conda.svg)](https://anaconda.org/bioconda/sequencetools) This repository contains some programs that I use for processing sequencing data. @@ -11,7 +9,6 @@ 1. Download stack (https://docs.haskellstack.org/en/stable/README/#how-to-install). 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. -4. Run `pileupCaller --version`. It should output `1.4.0`. You're all set. ## pileupCaller @@ -31,7 +28,17 @@ --samplePopName MyPop -f \ -e < 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. +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, (required) option `-f` needs an Eigenstrat positions file. This is required for pileupCaller to know what is the reference and which the alternative allele in your reference dataset that you want to call. An Eigenstrat positions file is a line-based file format, where each line denotes a SNP position, and there are exactly six required columns, denoting in order i) SNP ID, ii) chromosome, iii) genetic position (can be set to zero), iv) physical position, v) reference allele, vi) alternate allele. Here is an example: + + rs0000 11 0.000000 0 A C + rs1111 11 0.001000 100000 A G + rs2222 11 0.002000 200000 A T + rs3333 11 0.003000 300000 C A + rs4444 11 0.004000 400000 G A + rs5555 11 0.005000 500000 T A + rs6666 11 0.006000 600000 G T + +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. diff --git a/sequenceTools.cabal b/sequenceTools.cabal --- a/sequenceTools.cabal +++ b/sequenceTools.cabal @@ -1,15 +1,16 @@ -cabal-version: >=1.10 -name: sequenceTools -version: 1.4.0.5 -license: GPL-3 -license-file: LICENSE -maintainer: stephan.schiffels@mac.com -author: Stephan Schiffels -synopsis: A package with tools for processing DNA sequencing data +cabal-version: >=1.10 +name: sequenceTools +version: 1.5.0 +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 + +category: Bioinformatics +build-type: Simple extra-source-files: README.md Changelog.md @@ -18,84 +19,85 @@ exposed-modules: SequenceTools.Utils SequenceTools.PileupCaller - hs-source-dirs: src - other-modules: - Paths_sequenceTools + + hs-source-dirs: src + other-modules: Paths_sequenceTools default-language: Haskell2010 build-depends: base >=4.7 && <5, - optparse-applicative >=0.14.3.0, + optparse-applicative >=0.15.1.0, random >=1.1, - sequence-formats >=1.4.1, - bytestring >=0.10.8.2, - vector >=0.12.0.3, - pipes >=4.3.11 + sequence-formats >=1.6.1, + bytestring >=0.10.12.0, + vector >=0.12.1.2, + pipes >=4.3.14 executable pileupCaller - main-is: pileupCaller.hs - hs-source-dirs: src-executables + main-is: pileupCaller.hs + hs-source-dirs: src-executables default-language: Haskell2010 - ghc-options: -threaded -rtsopts -with-rtsopts=-N2 + ghc-options: -threaded -rtsopts -with-rtsopts=-N2 build-depends: - base >=4.12.0.0, + base >=4.14.1.0, sequenceTools -any, - sequence-formats >=1.4.1, - optparse-applicative >=0.14.3.0, - pipes >=4.3.11, - rio >=0.1.11.0, - vector >=0.12.0.3, + sequence-formats >=1.6.1, + optparse-applicative >=0.15.1.0, + pipes >=4.3.14, + rio >=0.1.19.0, + vector >=0.12.1.2, random >=1.1, - bytestring >=0.10.8.2, - pipes-safe >=2.3.1, + bytestring >=0.10.12.0, + pipes-safe >=2.3.2, pipes-ordered-zip >=1.1.0, - split >=0.2.3.3, + split >=0.2.3.4, ansi-wl-pprint >=0.6.9 executable vcf2eigenstrat - main-is: vcf2eigenstrat.hs - hs-source-dirs: src-executables + main-is: vcf2eigenstrat.hs + hs-source-dirs: src-executables default-language: Haskell2010 build-depends: - base >=4.12.0.0, + base >=4.14.1.0, sequenceTools -any, pipes-ordered-zip >=1.1.0, - sequence-formats >=1.4.1, - 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 + sequence-formats >=1.6.1, + 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 executable genoStats - main-is: genoStats.hs - hs-source-dirs: src-executables + main-is: genoStats.hs + hs-source-dirs: src-executables default-language: Haskell2010 build-depends: - base >=4.12.0.0, - sequence-formats >=1.4.1, + base >=4.14.1.0, + sequence-formats >=1.6.1, 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, + 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.1 + pipes-safe >=2.3.2 test-suite sequenceToolsTests - type: exitcode-stdio-1.0 - main-is: Spec.hs - hs-source-dirs: test + 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, + base >=4.14.1.0, + hspec >=2.7.8, sequenceTools -any, - sequence-formats >=1.4.1, - vector >=0.12.0.3, - bytestring >=0.10.8.2, - pipes >=4.3.11 + sequence-formats >=1.6.1, + vector >=0.12.1.2, + bytestring >=0.10.12.0, + pipes >=4.3.14 diff --git a/src-executables/pileupCaller.hs b/src-executables/pileupCaller.hs --- a/src-executables/pileupCaller.hs +++ b/src-executables/pileupCaller.hs @@ -7,6 +7,7 @@ import SequenceTools.Utils (versionInfoOpt, versionInfoText, freqSumToEigenstrat) import SequenceTools.PileupCaller (CallingMode(..), callGenotypeFromPileup, callToDosage, filterTransitions, TransitionsMode(..), cleanSSdamageAllSamples) +import SequenceFormats.Plink (writePlink) import Data.List.Split (splitOn) import qualified Data.Vector.Unboxed.Mutable as V @@ -16,12 +17,12 @@ import qualified Pipes.Prelude as P import Pipes.Safe (runSafeT, SafeT) import RIO -import System.IO (readFile, hPutStrLn, stderr) +import System.IO (readFile, hPutStrLn, stderr, print) import System.Random (mkStdGen, setStdGen) import Text.Printf (printf) import qualified Text.PrettyPrint.ANSI.Leijen as PP -data OutFormat = EigenstratFormat String String | FreqSumFormat +data OutFormat = EigenstratFormat String String | PlinkFormat String String | FreqSumFormat deriving (Show) data ProgOpt = ProgOpt CallingMode Bool (Maybe Int) Int TransitionsMode FilePath OutFormat (Either [String] FilePath) --optCallingMode :: CallingMode, @@ -52,6 +53,9 @@ envStats :: ReadStats } +instance Show Env where + show (Env m r d t o sf sn _) = show (m, r, d, t, o, sf, sn) + type App = ReaderT Env (SafeT IO) main :: IO () @@ -131,14 +135,21 @@ \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 + parseFormat = (EigenstratFormat <$> parseEigenstratPrefix <*> parseSamplePopName) <|> (PlinkFormat <$> parsePlinkPrefix <*> parseSamplePopName) <|> pure FreqSumFormat parseEigenstratPrefix = OP.strOption (OP.long "eigenstratOut" <> OP.short 'e' <> OP.metavar "" <> OP.help "Set Eigenstrat as output format. Specify the filenames for the EigenStrat \ - \SNP and IND file outputs: .snp.txt and .ind.txt \ + \SNP, IND and GENO file outputs: .snp, .ind and .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") + parsePlinkPrefix = OP.strOption (OP.long "plinkOut" <> OP.short 'p' <> + OP.metavar "" <> + OP.help "Set Plink as output format. Specify the filenames for the Plink \ + \BIM, FAM and BED file outputs: .bim, .fam and .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") parseSampleNames = parseSampleNameList <|> parseSampleNameFile parseSampleNameList = OP.option (Left . splitOn "," <$> OP.str) (OP.long "sampleNames" <> OP.metavar "NAME1,NAME2,..." <> @@ -194,13 +205,16 @@ runMain :: App () runMain = do + env_ <- ask + -- liftIO $ print env_ 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 + EigenstratFormat outPrefix popName -> outputEigenStratOrPlink outPrefix popName False freqSumProducer + PlinkFormat outPrefix popName -> outputEigenStratOrPlink outPrefix popName True freqSumProducer outputStats pileupToFreqSum :: FilePath -> Producer PileupRow (SafeT IO) () -> @@ -277,8 +291,8 @@ outProd = freqSumProducer >-> filterTransitions transitionsOnly lift . runEffect $ outProd >-> printFreqSumStdOut header' -outputEigenStrat :: FilePath -> String -> Producer FreqSumEntry (SafeT IO) () -> App () -outputEigenStrat outPrefix popName freqSumProducer = do +outputEigenStratOrPlink :: FilePath -> String -> Bool -> Producer FreqSumEntry (SafeT IO) () -> App () +outputEigenStratOrPlink outPrefix popName formatIsPlink freqSumProducer = do transitionsMode <- asks envTransitionsMode sampleNames <- asks envSampleNames callingMode <- asks envCallingMode @@ -286,13 +300,15 @@ RandomCalling -> True MajorityCalling _ -> True RandomDiploidCalling -> False - let snpOut = outPrefix <> ".snp.txt" - indOut = outPrefix <> ".ind.txt" - genoOut = outPrefix <> ".geno.txt" + 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 indEntries = [EigenstratIndEntry n Unknown popName | n <- sampleNames] lift . runEffect $ freqSumProducer >-> filterTransitions transitionsMode >-> P.map (freqSumToEigenstrat diploidizeCall) >-> - writeEigenstrat genoOut snpOut indOut indEntries + writeFunc genoOut snpOut indOut indEntries outputStats :: App () outputStats = do diff --git a/src/SequenceTools/PileupCaller.hs b/src/SequenceTools/PileupCaller.hs --- a/src/SequenceTools/PileupCaller.hs +++ b/src/SequenceTools/PileupCaller.hs @@ -16,9 +16,9 @@ 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 CallingMode = MajorityCalling Bool | RandomCalling | RandomDiploidCalling deriving (Show) -data TransitionsMode = TransitionsMissing | SkipTransitions | SingleStrandMode | AllSites deriving (Eq) +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