varan 0.5 → 0.5.1
raw patch · 6 files changed
+50/−38 lines, 6 files
Files
- src/MPileup.hs +3/−15
- src/Options.hs +3/−2
- src/Process.hs +38/−8
- src/Sparks.hs +1/−1
- src/VExtr.hs +1/−1
- varan.cabal +4/−11
src/MPileup.hs view
@@ -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_ }
src/Options.hs view
@@ -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:"
src/Process.hs view
@@ -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]
src/Sparks.hs view
@@ -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
src/VExtr.hs view
@@ -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:"
varan.cabal view
@@ -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 <ketil@malde.org> 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-