diff --git a/Changelog.md b/Changelog.md --- a/Changelog.md +++ b/Changelog.md @@ -20,3 +20,7 @@ * Moved unmaintained scripts into unmaintained folder. V 1.5.0: Added support for Plink output + +V 1.5.1: Added automatic building + +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. diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -4,11 +4,22 @@ This repository contains some programs that I use for processing sequencing data. -# Simple Installation of the main tools +# Installation +## Simple Installation via stack and hackage +This installation installs the latest version that is up on [hackage](https://hackage.haskell.org/package/sequenceTools): + 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. + +## Installation from source via stack + +1. Download stack (https://docs.haskellstack.org/en/stable/README/#how-to-install). +2. Clone this repository via `git clone https://github.com/stschiff/sequenceTools.git` +3. Install via `cd sequenceTools; stack install` + +# Commands ## pileupCaller diff --git a/sequenceTools.cabal b/sequenceTools.cabal --- a/sequenceTools.cabal +++ b/sequenceTools.cabal @@ -1,6 +1,6 @@ cabal-version: >=1.10 name: sequenceTools -version: 1.5.0 +version: 1.5.2 license: GPL-3 license-file: LICENSE maintainer: stephan.schiffels@mac.com @@ -27,7 +27,7 @@ base >=4.7 && <5, optparse-applicative >=0.15.1.0, random >=1.1, - sequence-formats >=1.6.1, + sequence-formats >=1.6.3, bytestring >=0.10.12.0, vector >=0.12.1.2, pipes >=4.3.14 @@ -40,17 +40,17 @@ build-depends: base >=4.14.1.0, sequenceTools -any, - sequence-formats >=1.6.1, + sequence-formats >=1.6.3, 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.12.0, pipes-safe >=2.3.2, pipes-ordered-zip >=1.1.0, split >=0.2.3.4, - ansi-wl-pprint >=0.6.9 + ansi-wl-pprint >=0.6.9, + transformers >=0.5.6.2 executable vcf2eigenstrat main-is: vcf2eigenstrat.hs @@ -60,7 +60,7 @@ base >=4.14.1.0, sequenceTools -any, pipes-ordered-zip >=1.1.0, - sequence-formats >=1.6.1, + sequence-formats >=1.6.3, bytestring >=0.10.12.0, vector >=0.12.1.2, optparse-applicative >=0.15.1.0, @@ -73,7 +73,7 @@ default-language: Haskell2010 build-depends: base >=4.14.1.0, - sequence-formats >=1.6.1, + sequence-formats >=1.6.3, sequenceTools -any, foldl >=1.4.10, bytestring >=0.10.12.0, @@ -97,7 +97,7 @@ base >=4.14.1.0, hspec >=2.7.8, sequenceTools -any, - sequence-formats >=1.6.1, + sequence-formats >=1.6.3, 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 @@ -1,4 +1,4 @@ -{-# LANGUAGE OverloadedStrings, NoImplicitPrelude #-} +{-# LANGUAGE OverloadedStrings #-} import SequenceFormats.Eigenstrat (readEigenstratSnpFile, EigenstratSnpEntry(..), EigenstratIndEntry(..), Sex(..), writeEigenstrat) @@ -9,30 +9,37 @@ filterTransitions, TransitionsMode(..), cleanSSdamageAllSamples) import SequenceFormats.Plink (writePlink) +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 Pipes.OrderedZip (orderedZip, orderCheckPipe) import qualified Pipes.Prelude as P +import Pipes.OrderedZip (orderedZip, orderCheckPipe) import Pipes.Safe (runSafeT, SafeT) -import RIO -import System.IO (readFile, hPutStrLn, stderr, print) +import System.IO (hPutStrLn, stderr) import System.Random (mkStdGen, setStdGen) import Text.Printf (printf) import qualified Text.PrettyPrint.ANSI.Leijen as PP -data OutFormat = EigenstratFormat String String | PlinkFormat String String | FreqSumFormat deriving (Show) +data OutFormat = EigenstratFormat FilePath | PlinkFormat FilePath | FreqSumFormat deriving (Show) -data ProgOpt = ProgOpt CallingMode Bool (Maybe Int) Int TransitionsMode FilePath OutFormat (Either [String] FilePath) - --optCallingMode :: CallingMode, - --optKeepInCongruentReads :: Bool, - --optSeed :: Maybe Int, - --optMinDepth :: Int, - --optTransitionsMode :: TransitionsMode, - --optSnpFile :: FilePath, - --optOutFormat :: OutFormat, - --optSampleNames :: Either [String] FilePath +data ProgOpt = ProgOpt { + optCallingMode :: CallingMode, + optKeepInCongruentReads :: Bool, + optSeed :: Maybe Int, + optMinDepth :: Int, + optTransitionsMode :: TransitionsMode, + optSnpFile :: FilePath, + optOutFormat :: OutFormat, + optSampleNames :: Either [String] FilePath, + optPopName :: String +} data ReadStats = ReadStats { rsTotalSites :: IORef Int, @@ -50,11 +57,12 @@ envOutFormat :: OutFormat, envSnpFile :: FilePath, envSampleNames :: [String], + envPopName :: String, envStats :: ReadStats } instance Show Env where - show (Env m r d t o sf sn _) = show (m, r, d, t, o, sf, sn) + show (Env m r d t o sf sn pn _) = show (m, r, d, t, o, sf, sn, pn) type App = ReaderT Env (SafeT IO) @@ -69,8 +77,15 @@ (OP.progDescDoc (Just programHelpDoc)) argParser :: OP.Parser ProgOpt -argParser = ProgOpt <$> parseCallingMode <*> parseKeepIncongruentReads <*> parseSeed <*> parseMinDepth <*> - parseTransitionsMode <*> parseSnpFile <*> parseFormat <*> parseSampleNames +argParser = ProgOpt <$> parseCallingMode + <*> parseKeepIncongruentReads + <*> parseSeed + <*> parseMinDepth + <*> parseTransitionsMode + <*> parseSnpFile + <*> parseFormat + <*> parseSampleNames + <*> parsePopName where parseCallingMode = parseRandomCalling <|> parseMajorityCalling <|> parseRandomDiploidCalling parseRandomCalling = OP.flag' RandomCalling (OP.long "randomHaploid" <> @@ -135,7 +150,7 @@ \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) <|> (PlinkFormat <$> parsePlinkPrefix <*> parseSamplePopName) <|> pure FreqSumFormat + parseFormat = (EigenstratFormat <$> parseEigenstratPrefix) <|> (PlinkFormat <$> parsePlinkPrefix) <|> 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 \ @@ -157,7 +172,7 @@ parseSampleNameFile = OP.option (Right <$> OP.str) (OP.long "sampleNameFile" <> OP.metavar "" <> OP.help "give the names of the samples in a file with one name per \ \line") - parseSamplePopName = OP.strOption (OP.long "samplePopName" <> OP.value "Unknown" <> OP.showDefault <> + 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\ \ in the output *.ind.txt file in Eigenstrat output. This will be ignored if the output \ @@ -186,7 +201,7 @@ initialiseEnvironment :: ProgOpt -> IO Env initialiseEnvironment args = do let (ProgOpt callingMode keepInCongruentReads seed minDepth - transitionsMode snpFile outFormat sampleNames) = args + transitionsMode snpFile outFormat sampleNames popName) = args case seed of Nothing -> return () Just seed_ -> liftIO . setStdGen $ mkStdGen seed_ @@ -196,7 +211,7 @@ let n = length sampleNamesList readStats <- ReadStats <$> newIORef 0 <*> makeVec n <*> makeVec n <*> makeVec n <*> makeVec n return $ Env callingMode keepInCongruentReads minDepth transitionsMode - outFormat snpFile sampleNamesList readStats + outFormat snpFile sampleNamesList popName readStats where makeVec n = do v <- V.new n @@ -205,16 +220,15 @@ runMain :: App () runMain = do - env_ <- ask - -- liftIO $ print env_ let pileupProducer = readPileupFromStdIn snpFile <- asks envSnpFile freqSumProducer <- pileupToFreqSum snpFile pileupProducer outFormat <- asks envOutFormat + popName <- asks envPopName case outFormat of FreqSumFormat -> outputFreqSum freqSumProducer - EigenstratFormat outPrefix popName -> outputEigenStratOrPlink outPrefix popName False freqSumProducer - PlinkFormat outPrefix popName -> outputEigenStratOrPlink outPrefix popName True freqSumProducer + EigenstratFormat outPrefix -> outputEigenStratOrPlink outPrefix popName False freqSumProducer + PlinkFormat outPrefix -> outputEigenStratOrPlink outPrefix popName True freqSumProducer outputStats pileupToFreqSum :: FilePath -> Producer PileupRow (SafeT IO) () -> @@ -324,7 +338,7 @@ \# avgSampledFrom: mean coverage of pileup after removing reads with tri-allelic alleles \n\ \SampleName\tTotalSites\tNonMissingCalls\tavgRawReads\tavgDamageCleanedReads\tavgSampledFrom" forM_ (zip [0..] sampleNames) $ \(i, name) -> do - totalS <- readIORef totalSites + totalS <- liftIO $ readIORef totalSites nonMissingSites <- V.read nonMissingSitesVec i rawReads <- V.read rawReadsVec i damageCleanedReads <- V.read damageCleanedReadsVec i diff --git a/src-executables/vcf2eigenstrat.hs b/src-executables/vcf2eigenstrat.hs --- a/src-executables/vcf2eigenstrat.hs +++ b/src-executables/vcf2eigenstrat.hs @@ -13,7 +13,6 @@ import Control.Monad (when) import Control.Monad.IO.Class (liftIO, MonadIO) import qualified Data.ByteString.Char8 as B -import Data.Monoid ((<>)) -- import Debug.Trace (trace) import qualified Options.Applicative as OP import Pipes (Pipe, yield, (>->), runEffect, Producer, Pipe, for, cat) diff --git a/src/SequenceTools/PileupCaller.hs b/src/SequenceTools/PileupCaller.hs --- a/src/SequenceTools/PileupCaller.hs +++ b/src/SequenceTools/PileupCaller.hs @@ -105,13 +105,13 @@ ((ref == 'T') && (alt == 'C')) cleanSSdamageAllSamples :: Char -> Char -> [String] -> [[Strand]] -> [String] -cleanSSdamageAllSamples ref alt basesPerSample strandPerSample = - if (ref, alt) == ('C', 'T') || (ref, alt) == ('T', 'C') - then [removeForwardBases bases strands | (bases, strands) <- zip basesPerSample strandPerSample] - else - if (ref, alt) == ('G', 'A') || (ref, alt) == ('A', 'G') - then [removeReverseBases bases strands | (bases, strands) <- zip basesPerSample strandPerSample] - else basesPerSample +cleanSSdamageAllSamples ref alt basesPerSample strandPerSample + | (ref, alt) == ('C', 'T') || (ref, alt) == ('T', 'C') = + [removeForwardBases bases strands | (bases, strands) <- zip basesPerSample strandPerSample] + | (ref, alt) == ('G', 'A') || (ref, alt) == ('A', 'G') = + [removeReverseBases bases strands | (bases, strands) <- zip basesPerSample strandPerSample] + | otherwise = + basesPerSample where removeForwardBases = removeReads ForwardStrand removeReverseBases = removeReads ReverseStrand