diff --git a/src/MPileup.hs b/src/MPileup.hs --- a/src/MPileup.hs +++ b/src/MPileup.hs @@ -1,8 +1,8 @@ {-# Language BangPatterns #-} -module MPileup (Counts(..), readPile1, toList, major_allele, by_major_allele, showC, showV, sumList, MPileRecord(..)) where +module MPileup (Counts(..), readPile1, toList, major_allele, by_major_allele, sumList, MPileRecord(..)) where import Data.Char (toUpper) -import Data.List (intercalate,nub,elemIndex) +import Data.List (elemIndex) import qualified Data.ByteString.Lazy.Char8 as B import Variants hiding (parse) import Count @@ -69,16 +69,4 @@ var = (if c=='+' then Ins else Del) (B.unpack $ B.take (fromIntegral cnt) rest) in parse ref (addV cts var) (B.drop (fromIntegral cnt) rest) | otherwise = error ("Not a nucleotide: "++show c) - addRef = case ref of { 'A' -> addA_; 'C' -> addC_; 'G' -> addG_; 'T' -> addT_; _ -> addN_ } - --- | Show SNP counts and coverage -showC :: Counts -> (String,Int) -showC x = (" "++(intercalate ":" $ map show [getA_ x,getC_ x,getG_ x,getT_ x,getN_ x,getDel_ x]),covC x) - --- | Show structural variant count -showV :: [Counts] -> String -showV cs = let - vs = [[v | Ins v <- getV c] | c <- cs] - vuniq = nub $ concat vs - countV v = map (length . filter (==v)) vs - in intercalate "," $ [unwords (v:(map show $ countV v)) | v <- vuniq] + addRef = case toUpper ref of { 'A' -> addA_; 'C' -> addC_; 'G' -> addG_; 'T' -> addT_; _ -> addN_ } diff --git a/src/Options.hs b/src/Options.hs --- a/src/Options.hs +++ b/src/Options.hs @@ -13,7 +13,7 @@ , chi2, f_st, pi_k , conf, ds, dsw, esi, pconf, nd_all, maf :: Bool , input, output :: FilePath - , global :: Bool + , global, sync :: Bool , threads :: Int , min_cov, max_cov :: Int } deriving (Typeable,Data) @@ -31,6 +31,7 @@ , input = [] &= args &= typFile , global = False &= help "calculate global statistics" + , sync = False &= help "use 'sync'-compatible format (A:T:C:G:N:-)" -- Overall statistics , chi2 = False &= help "calculate chi² probability" &= ignore @@ -48,7 +49,7 @@ -- Statistics for all sample pairs , esi = False &= help "output conservative expected site information for SNPs using Agresti-Coull intervals" } - &= program "varan v0.5" + &= program "varan v0.5.1" &= summary "Identify genetic variants from pooled sequences." &= details ["Examples:", "" ,"Read input from a pipe, calculate site-wise Fst and confidence intervals, ignoring non-variant sites:" diff --git a/src/Process.hs b/src/Process.hs --- a/src/Process.hs +++ b/src/Process.hs @@ -2,16 +2,18 @@ module Process (proc_fused, run_procs, showPile') where -import Options +import qualified Options +import Options (Options, output, threads) import ParMap (parMap) -import MPileup (readPile1, counts, showC, showV, MPileRecord(..)) +import MPileup (readPile1, counts, MPileRecord(..)) import Metrics (pi_k, f_st, nd , conf_all, ds_all, dsw_all, maf , fst_params, ppi_params, dsconf_pairs) import Count (getV, covC, Counts(..), ptSum) +import Variants (Variant(..)) import ESIV (esiv) -import Data.List (tails) +import Data.List (tails, intercalate, nub) import qualified Data.ByteString.Char8 as B import qualified Data.ByteString.Lazy.Char8 as BL import Text.Printf @@ -89,7 +91,7 @@ where f (MPR sup _ _ _ cts) cur = let new = Metrics.fst_params cts cov = sum $ map covC cts - in if sup || (max_cov o > 0 && cov > max_cov o) || cov < min_cov o + in if sup || (Options.max_cov o > 0 && cov > Options.max_cov o) || cov < Options.min_cov o then cur else deepSeq $ zipWith (zipWith plus) cur new zero = repeat (repeat (0,0)) plus (a,c) (b,d) = (a+b,c+d) @@ -122,7 +124,7 @@ where f (MPR sup _ _ _ cts) (uv,cur) = let new = Metrics.ppi_params cts cov = sum $ map covC cts - ign = sup || (max_cov o > 0 && cov > max_cov o) || cov < min_cov o + ign = sup || (Options.max_cov o > 0 && cov > Options.max_cov o) || cov < Options.min_cov o nu = if ign then uv else add_uv uv (fromIntegral cov) nc = if ign then cur else (deepSeq (zipWith (zipWith plus) cur new)) in nu `seq` nc `seq` (nu,nc) @@ -194,8 +196,8 @@ showPile :: Options -> MPileRecord -> B.ByteString showPile _ (MPR _ _ _ _ []) = error "Pileup with no data?" -showPile o mpr = if suppress o && ignore mpr && (all null (map getV $ counts mpr) || not (variants o)) then B.empty else (B.concat - [ default_out mpr +showPile o mpr = if Options.suppress o && ignore mpr && (all null (map getV $ counts mpr) || not (Options.variants o)) then B.empty else (B.concat + [ if (Options.sync o) then sync_out mpr else default_out mpr , when (Options.f_st o) (printf "\t%.3f" (Metrics.f_st $ counts mpr)) , when (Options.pi_k o) (printf "\t%.3f" (Metrics.pi_k $ counts mpr)) -- , when (Options.chi2 o) (printf "\t%.3f" (Metrics.pearsons_chi² $ by_major_allele $ counts mpr)) @@ -220,10 +222,38 @@ default_out :: MPileRecord -> B.ByteString default_out (MPR _ chr pos ref stats) = B.concat ([chr',tab,pos',tab,B.singleton ref]++samples++fmtcounts) - where cnts = map showC stats + where cnts = map showC1 stats tab = B.pack "\t" samples = [B.append tab (B.pack s) | s <- map fst cnts] fmtcounts = [tab,B.pack $ show $ sum $ map snd cnts] -- todo: add indels? chr' = B.concat (BL.toChunks chr) pos' = B.concat (BL.toChunks pos) + +-- | Show SNP counts and coverage +showC1 :: Counts -> (String,Int) +showC1 x = (" "++(intercalate "/" $ map show [getA_ x,getC_ x,getG_ x,getT_ x]),covC x) + +-- | The default output, with only coverage statistics +sync_out :: MPileRecord -> B.ByteString +sync_out (MPR _ chr pos ref stats) = + B.concat ([chr',tab,pos',tab,B.singleton ref]++samples++fmtcounts) + where cnts = map showC2 stats + tab = B.pack "\t" + samples = [B.append tab (B.pack s) | s <- map fst cnts] + fmtcounts = [tab,B.pack $ show $ sum $ map snd cnts] -- todo: add indels? + chr' = B.concat (BL.toChunks chr) + pos' = B.concat (BL.toChunks pos) + +-- | Show SNP counts and coverage +showC2 :: Counts -> (String,Int) +showC2 x = (" "++(intercalate ":" $ map show [getA_ x,getC_ x,getG_ x,getT_ x,getN_ x,getDel_ x]),covC x) + + +-- | Show structural variant count +showV :: [Counts] -> String +showV cs = let + vs = [[v | Ins v <- getV c] | c <- cs] + vuniq = nub $ concat vs + countV v = map (length . filter (==v)) vs + in intercalate "," $ [unwords (v:(map show $ countV v)) | v <- vuniq] diff --git a/src/Sparks.hs b/src/Sparks.hs --- a/src/Sparks.hs +++ b/src/Sparks.hs @@ -33,7 +33,7 @@ main = do opts <- cmdArgsRun $ cmdArgsMode $ modes [test,disp,info,cite] &= summary "Visualize information from 'samtools mpileup' as sparklines" - &= program "sparks" + &= program "sparks v0.5.1" inp <- BL.getContents let ms = map readPile1 $ BL.lines inp mapM_ putStrLn $ case opts of diff --git a/src/VExtr.hs b/src/VExtr.hs --- a/src/VExtr.hs +++ b/src/VExtr.hs @@ -22,7 +22,7 @@ , fasta = False &= help "output FASTA header" , mincount = 1 &= help "ignore counts less than this" , minfreq = 5 &= help "ignore allele frequencies less than this" - } &= program "vextr v0.5" + } &= program "vextr v0.5.1" &= summary "Extract consensus sequence from pooled sequences" &= details ["Examples:", "" ,"Read input from a pipe, output IUPAC codes:" diff --git a/varan.cabal b/varan.cabal --- a/varan.cabal +++ b/varan.cabal @@ -1,5 +1,5 @@ Name: varan -Version: 0.5 +Version: 0.5.1 License: GPL Cabal-Version: >= 1.6 Build-Type: Simple @@ -7,10 +7,9 @@ Author: Ketil Malde Maintainer: Ketil Malde Synopsis: Process mpileup output to identify significant differences -Description: Using Agresti-Coull estimation of confidence interval, report - variant positions found in pooled samples along with significance of - the variant having different underlying allele frequency ('+' for 95%, - '*' for 99%). +Description: Post-processing output from `samtools mpileup` to extract various information, + including statistics (per-position or global), consensus sequence (in various + formats), and textual visualizations. Executable varan Hs-Source-Dirs: src @@ -19,7 +18,6 @@ Build-Depends: base >= 4 && < 5, random, mtl, parallel, statistics, cmdargs, bytestring Ghc-Options: -rtsopts -Wall -threaded --- this is just a quick hack, should be merged into the program proper Executable vextr Hs-Source-Dirs: src Main-Is: VExtr.hs @@ -33,8 +31,3 @@ Other-Modules: MPileup, Count Build-Depends: base >= 4 && < 5, bytestring, cmdargs, ansi-terminal Ghc-Options: -rtsopts -Wall -threaded - --- For parallel execution, we might want to add these, but they --- seem to degrade performance on older GHC --- Ghc-Options: -threaded -with-rtsopts=-N -