varan 0.3 → 0.5
raw patch · 9 files changed
+441/−77 lines, 9 filesdep +ansi-terminalnew-component:exe:sparksnew-component:exe:vextr
Dependencies added: ansi-terminal
Files
- src/Count.hs +23/−10
- src/ESIV.hs +7/−5
- src/MPileup.hs +13/−11
- src/Metrics.hs +98/−24
- src/Options.hs +18/−10
- src/Process.hs +31/−16
- src/Sparks.hs +124/−0
- src/VExtr.hs +111/−0
- varan.cabal +16/−1
src/Count.hs view
@@ -1,9 +1,10 @@ -- High performance counting data structure module Count (Counts(..)- , addA, addC, addG, addT, addV- , addA_, addC_, addG_, addT_ - , getA, getC, getG, getT+ , addA, addC, addG, addT, addN, addDel, addV+ , addA_, addC_, addG_, addT_, addN_, addDel_+ , getA, getC, getG, getT, getN, getDel+ , czero , covC, ptAdd, ptSum , toList, sumList ) where@@ -12,46 +13,58 @@ import Data.Int import Data.List (foldl1') -data Counts = C { getA_, getC_, getG_, getT_ :: {-# UNPACK #-} !Int32+data Counts = C { getA_, getC_, getG_, getT_, getN_, getDel_ :: {-# UNPACK #-} !Int32 , getV :: ![Variant] }+czero :: Counts+czero = C 0 0 0 0 0 0 [] -getA, getC, getG, getT :: Integral i => Counts -> i+getA, getC, getG, getT, getN, getDel :: Integral i => Counts -> i getA = fromIntegral . getA_ getC = fromIntegral . getC_ getG = fromIntegral . getG_ getT = fromIntegral . getT_+getN = fromIntegral . getN_+getDel = fromIntegral . getDel_ {-# INLINE getA #-} {-# INLINE getC #-} {-# INLINE getG #-} {-# INLINE getT #-}+{-# INLINE getN #-}+{-# INLINE getDel #-} -addA, addC, addG, addT :: Integral i => Counts -> i -> Counts+addA, addC, addG, addT, addN, addDel :: Integral i => Counts -> i -> Counts addA c i = c { getA_ = getA_ c + fromIntegral i } addC c i = c { getC_ = getC_ c + fromIntegral i } addG c i = c { getG_ = getG_ c + fromIntegral i } addT c i = c { getT_ = getT_ c + fromIntegral i }+addN c i = c { getN_ = getN_ c + fromIntegral i }+addDel c i = c { getDel_ = getDel_ c + fromIntegral i } {-# INLINE addA #-} {-# INLINE addC #-} {-# INLINE addG #-} {-# INLINE addT #-}+{-# INLINE addN #-}+{-# INLINE addDel #-} -addA_, addC_, addG_, addT_ :: Counts -> Int32 -> Counts+addA_, addC_, addG_, addT_, addN_, addDel_ :: Counts -> Int32 -> Counts addA_ c i = c { getA_ = getA_ c + i } addC_ c i = c { getC_ = getC_ c + i } addG_ c i = c { getG_ = getG_ c + i } addT_ c i = c { getT_ = getT_ c + i }-+addN_ c i = c { getN_ = getN_ c + i }+addDel_ c i = c { getDel_ = getDel_ c + i } addV :: Counts -> Variant -> Counts addV c v = c { getV = v : getV c } {-# INLINE addV #-} covC :: Counts -> Int-covC c = fromIntegral (getA_ c + getC_ c + getG_ c + getT_ c)+covC c = fromIntegral (getA_ c + getC_ c + getG_ c + getT_ c) -- + getN_ c + getDel_ c)+ -- mostly, we want to compare the identified bases. ptAdd :: Counts -> Counts -> Counts-ptAdd a b = a `addA` (getA_ b) `addC` (getC_ b) `addG` (getG_ b) `addT` (getT_ b)+ptAdd a b = a `addA` (getA_ b) `addC` (getC_ b) `addG` (getG_ b) `addT` (getT_ b) `addN` (getN_ b) `addDel` (getDel_ b) -- todo: concat variants? ptSum :: [Counts] -> Counts ptSum = foldl1' ptAdd
src/ESIV.hs view
@@ -20,15 +20,17 @@ esiv z epsilon c1 c2 = let t1 = covC c1 t2 = covC c2- esiv1 (a1,b1) (a2,b2) - | j1+epsilon<i2 = esiv_score (j1+epsilon/2) (i2-epsilon/2)- | i1>j2+epsilon = esiv_score (j2+epsilon/2) (i1-epsilon/2)+ in sum [abs $ esiv1 z epsilon (x,t1-x) (y,t2-y) | (x,y) <- zip (toList c1) (toList c2)]++esiv1 :: Double -> Double -> (Int, Int) -> (Int, Int) -> Double+esiv1 z eps (a1,b1) (a2,b2) + | j1+eps<i2 = esiv_score (j1+eps/2) (i2-eps/2)+ | i1>j2+eps = esiv_score (i1-eps/2) (j2+eps/2) | otherwise = 0 where (i1,j1) = confidenceInterval z a1 b1 (i2,j2) = confidenceInterval z a2 b2- in sum [esiv1 (x,t1-x) (y,t2-y) | (x,y) <- zip (toList c1) (toList c2)] esiv_score :: Double -> Double -> Double-esiv_score p1 p2 = abs ((p1+p2)/2*logBase 2 (p1/p2))+esiv_score p1 p2 = (p1+p2)/2*logBase 2 (p1/p2)
src/MPileup.hs view
@@ -38,13 +38,14 @@ parse1 (chr:pos:r:rest) = let trs = triples ref rest ref = B.head r in MPR (ign trs) chr pos ref trs- parse1 xs = error ("parse1: insufficiently long line:"++show xs)+ parse1 xs = error ("parse1: line too short: "++show xs) -- set the ignore flag if we only see one or zero alleles ign xs = (length $ filter (/=(0::Int)) $ toList $ ptSum xs) <= 1+ -- variants are checked for in the showPile' function triples _ [] = []- triples ref (_cnt:bases:_quals:rest) = let this = parse ref (C 0 0 0 0 []) (B.map toUpper bases) + triples ref (_cnt:bases:_quals:rest) = let this = parse ref czero (B.map toUpper bases) in this `seq` this : triples ref rest triples _ _ = error "triples: incorrect number of columns" @@ -60,23 +61,24 @@ | c == 'C' = parse ref (addC_ cts 1) str | c == 'G' = parse ref (addG_ cts 1) str | c == 'T' = parse ref (addT_ cts 1) str- | c == 'N' = parse ref cts str - | c == '^' = parse ref (addRef cts 1) $ B.drop 1 str- | c == '*' || c == '$' = parse ref cts str -- * is a deletion, also reported as variant...+ | c == 'N' = parse ref (addN_ cts 1) str+ | c == '^' = parse ref cts $ B.drop 1 str+ | c == '*' = parse ref (addDel_ cts 1) str+ | c == '$' = parse ref cts str | c == '-' || c == '+' = let Just (cnt,rest) = B.readInt str 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_; _ -> const }+ 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 (toList x :: [Int])),covC x)+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 = nub $ concatMap getV cs- countV :: Variant -> Counts -> Int- countV v c = length . filter (==v) $ getV c- in intercalate "\t" (show vs:[unwords $ map (\v -> show $ countV v c) vs | c <- cs])+ 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/Metrics.hs view
@@ -1,12 +1,18 @@+{-# Options_Ghc -fno-warn-unused-binds #-}+ -- | Calculate various metrics/statistics-module Metrics where+module Metrics+ ( pi_k, f_st, nd, maf+ , conf_all, ds_all, dsw_all+ , fst_params, ppi_params, dsconf_pairs)+ where import AgrestiCoull import MPileup (by_major_allele) import Count import Statistics.Distribution import Statistics.Distribution.ChiSquared-import Data.List (foldl1')+import Data.List (foldl1', tails, sort) -- | Calculate vector angle between allele frequencies. This is -- similar to `dist`, but from 1 (equal) to 0 (orthogonal)@@ -18,18 +24,26 @@ -- Calculate pairwise nucleotide diversities ppi_params :: [Counts] -> [[Double]]-ppi_params (c:cs) = map (\x -> pi_k [c,x]) (c:cs) : ppi_params cs+ppi_params (c:cs) = map (\x -> nd2 c x) (c:cs) : ppi_params cs ppi_params [] = [] -- calculate diversity within and between sample pairs fst_params :: [Counts] -> [[(Double,Double)]] fst_params (x:xs) = go (x:xs)- where go (y:ys) = map (heteroz $ y) ys : go ys+ where go (y:ys) = map (heteroz y) ys : go ys go [] = [] fst_params [] = [] -heteroz :: Counts -> Counts -> (Double,Double)-heteroz c1 c2 = let+-- | Calculate heterozyogisity total, and within groups+-- Not weighting by coverage.+heteroz :: Counts -> Counts -> (Double, Double)+heteroz c1 c2 = let nd_tot = nd c1 + nd c2+ in if covC c1 == 0 || covC c2 == 0 then (0,0) + else (nd (c1 `ptAdd` c2), nd_tot/2)++-- Weighted heterozygosity+heteroz_w :: Counts -> Counts -> (Double,Double)+heteroz_w c1 c2 = let c1s = fromIntegral $ covC c1 c2s = fromIntegral $ covC c2 total = c1s + c2s@@ -39,8 +53,8 @@ in if c1s == 0 || c2s == 0 || h_tot == 0.0 then (0,0) else (h_tot,h_subs) -heteroz_ :: [Double] -> [Double] -> (Double, Double)-heteroz_ c1 c2 = let+heteroz_w2 :: [Double] -> [Double] -> (Double, Double)+heteroz_w2 c1 c2 = let c1s = sum c1 c2s = sum c2 total = c1s + c2s@@ -53,17 +67,27 @@ in if c1s == 0 || c2s == 0 || h_tot == 0 then (0,0) else (h_tot,h_subs) +-- | Simple calculation, samples represent equal populations (weights = 1/n) f_st :: [Counts] -> Double f_st xs = let- hz x = 1 - fromIntegral (sq (getA x::Int) + sq (getC x) + sq (getG x) + sq (getT x))/fromIntegral (sq $ covC x) where sq z = z*z+ l = fromIntegral (length xs)+ nd_sub = map ((/l) . nd) xs+ -- total heterozygosity based on average allele frequencies over populations+ nd_tot = 1 - sum (map (\x->x*x) $ map (/l) $ sumList (map pi_freqs xs))+ in if nd_tot == 0 then 0.0 else (nd_tot - sum nd_sub) / nd_tot++-- | Calculate F_ST. Note that this is weighted by the number of sequences (coverage)+-- this is not what we want for sequencing data!+f_st__ :: [Counts] -> Double+f_st__ xs = let h_subs, weights :: [Double]- h_tot = hz (ptSum xs)- h_subs = map hz xs+ h_tot = nd (ptSum xs)+ h_subs = map nd xs weights = let t = fromIntegral $ covC (ptSum xs) in [fromIntegral (covC x)/t | x <- xs] in if h_tot == 0 then 0.0 else (h_tot - sum (zipWith (*) h_subs weights)) / h_tot --- | Calculate F_ST+-- | Calculate F_ST - equivalent to the above f_st_ :: [Counts] -> Double f_st_ cs = let cs' = map toList cs@@ -80,17 +104,42 @@ in if h_tot == 0 then 0.0 -- no heterozygosity in the population! else (h_tot - sum (zipWith (*) h_subs weights)) / h_tot --- | Calculate Pi (my version), the expected number of differences--- between two random samples from the populations. I.e. the--- probability that sampling once from each population will not be all--- the same. One weakness is that if one population has fifty-fifty--- allele frequencies, the result is always exactly 0.5. I.e. it--- can't identify divergent allele frequencies in that case. Like Fst, this--- also is indifferent to the actual counts, so reliability depends on coverage.+-- | Caluclate MAF (minor allele frequency)+maf :: Counts -> Double+maf c = if cv == 0 then 0.0+ else (head . tail . reverse . sort . toList $ c)/fromIntegral cv+ where cv = covC c +-- | Calculate nucleotide diversity, the probability that sampling+-- twice will give you two different results. Should we correct by+-- a factor of c/(c-1) here? Note this gets weird with e.g. allele count of +-- 1 and 1 (nd=1, rather than 0.5)+nd :: Counts -> Double+nd cs = let fs = pi_freqs cs+ in 1-sum (zipWith (*) fs fs)++-- | Nucleotide diversity between - the probability of selecting+-- one from each will differ.+nd2 :: Counts -> Counts -> Double+nd2 cs1 cs2 = let fs1 = pi_freqs cs1+ fs2 = pi_freqs cs2+ in 1 - sum (zipWith (*) fs1 fs2)++-- | Calculate Pi, the probability of getting different nucleotides by+-- sampling from two different populations. pi_k :: [Counts] -> Double-pi_k cs = let fs = map pi_freqs cs- c = fromIntegral $ covC $ ptSum $ cs+pi_k cs = sum [nd2 x y | (x:x2:xs) <- tails cs, y <- (x2:xs)] * 2/(lcs*(lcs-1))+ where lcs = fromIntegral (length cs)++-- pi_k_WRONG is the probability that sampling once from each+-- population will not be all the same. One weakness is that if one+-- population has fifty-fifty allele frequencies, the result is always+-- exactly 0.5. I.e. it can't identify divergent allele frequencies+-- in that case. Like Fst, this also is indifferent to the actual+-- counts, so reliability depends on coverage.+pi_k_WRONG :: [Counts] -> Double+pi_k_WRONG cs = let fs = map pi_freqs cs+ c = fromIntegral $ covC $ ptSum $ cs in if c>1 then c/(c-1)*(1 - (sum $ foldl1' (zipWith (*)) fs)) else 0 pi_freqs :: Counts -> [Double]@@ -117,7 +166,7 @@ in if any (==0) (cols t) || any (==0) (rows t) then 1.0 else complCumulative (chiSquared df) chi -- | Use AgrestiCoull to calculate significant differences between--- allele frequency spectra+-- allele frequency spectra. Output * or +, depending on significance for each allele. conf :: Counts -> Counts -> String conf x y = let s1 = covC x -- don't count structural variants@@ -140,7 +189,7 @@ else '.' -- | Use AgrestiCoull to calculate significant difference from--- a combined distribution, with error.+-- the combined distribution, with error. -- FIXME: use an error threshold min dist between major allele frequency conf intervals conf_all :: [Counts] -> String conf_all cs' = let@@ -155,7 +204,8 @@ go [] = "" ds x y = if delta_sigma 2.326 x y > e then '*' else if delta_sigma 1.65 x y > e then '+' else '.' -- | Calculate distance (in absolute numbers) between confidence intervals --- with the given z-score+-- with the given z-score. This is overly conservative, and it is better to+-- use wald or wald_p below (or some more complex method, like Newcomb). delta_sigma :: Double -> (Int,Int) -> (Int,Int) -> Double delta_sigma z (s1,f1) (s2,f2) = let (i1,j1) = confidenceInterval z s1 f1@@ -173,6 +223,30 @@ (bs,bf) = (sum (map fst xs), sum (map snd xs)) pairs = [((s,f),(bs-s,bf-f)) | (s,f) <- xs ] in map (uncurry (delta_sigma sig)) pairs++-- as above, only using wald_p instead of agresti-coull+dsw_all :: Double -> [Counts] -> [Double]+dsw_all sig counts = let+ xs = by_major_allele counts+ (bs,bf) = (sum (map fst xs), sum (map snd xs))+ pairs = [((s,f),(bs-s,bf-f)) | (s,f) <- xs ]+ in map (uncurry (wald_p sig)) pairs++-- | Wald intervals with pseudocounts+-- See Agresti and Caffo, 2000.+wald_p :: Double -> (Int, Int) -> (Int, Int) -> Double+wald_p z (s1,f1) (s2,f2) = wald z (s1+1,f1+1) (s2+1,f2+1)++-- | Calculate lower bound of Wald confidence interval for difference between frequencies+wald :: Double -> (Int, Int) -> (Int, Int) -> Double+wald z (s1,f1) (s2,f2) = let+ -- estimated success frequencies+ (//) x y = fromIntegral x / fromIntegral y+ n1 = s1+f1+ n2 = s2+f2+ p1 = s1//n1+ p2 = s2//n2+ in abs (p1-p2) - z*sqrt ((s1*f1)//(n1*n1*n1)+(s2*f2)//(n2*n2*n2)) -- | Calculate distance between approximate distributions -- in terms of their standard deviation. Perhaps use binomial distribution directly?
src/Options.hs view
@@ -10,11 +10,12 @@ data Options = Opts { suppress, variants- , chi2, f_st, pi_k, conf, ds, esi, pconf :: Bool+ , chi2, f_st, pi_k+ , conf, ds, dsw, esi, pconf, nd_all, maf :: Bool , input, output :: FilePath , global :: Bool , threads :: Int- , min_cov, max_cov :: Int+ , min_cov, max_cov :: Int } deriving (Typeable,Data) defopts :: Options@@ -39,17 +40,24 @@ -- Per sample statistics (vs. pool of all populations) , conf = False &= help "check if major allele frequency confidence intervals overlap" -- uses conf_all (for each allele) , pconf = False &= help "pairwise major allele confidence" -- uses dsconf_pairs (by major allele)- , ds = False &= help "output distance between major allele frequency confidence intervals" -- uses ds_all (by major allele)+ , ds = False &= help "distance between major allele frequency confidence intervals, using Agresti-Coull" -- uses ds_all (by major allele)+ , dsw = False &= help "lower bound for distance between major allele frequencies, using Wald"+ , nd_all = False &= help "nucleotide diversity (unadjusted), per sample and overall"+ , maf = False &= help "minor allele frequency per position" -- Statistics for all sample pairs , esi = False &= help "output conservative expected site information for SNPs using Agresti-Coull intervals"- } &= program "varan"- &= 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:",""- ," samtools mpileup -f ref.fasta reads.bam | varan --fst --conf -s -o snps.txt",""- ,"Read input from a file, send the site-wise output to /dev/null, and only output global statistics to standard output:",""- ," varan --global -o /dev/null input.mpile",""+ }+ &= program "varan v0.5"+ &= 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:"+ ,"", " samtools mpileup -f ref.fasta reads.bam | varan --fst --conf -s -o snps.txt", ""+ ,"Read input from a file, send the site-wise output to /dev/null, and only output global statistics to standard output:"+ ,"", " varan --global -o /dev/null input.mpile", ""+ ,"If you use this program, please cite:"+ ," BMC Genomics 2014, 15(Suppl 6):S20"+ ," http://www.biomedcentral.com/1471-2164/15/S6/S20" ] getArgs :: IO (IO BL.ByteString,Options)
src/Process.hs view
@@ -3,11 +3,13 @@ module Process (proc_fused, run_procs, showPile') where import Options-import ParMap-import MPileup-import Metrics-import Count-import ESIV+import ParMap (parMap)+import MPileup (readPile1, counts, showC, showV, 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 ESIV (esiv) import Data.List (tails) import qualified Data.ByteString.Char8 as B@@ -36,7 +38,7 @@ outh <- if use_stdout then return stdout else openFile (output o) WriteMode B.hPutStr outh $ gen_header o r1 (gi,gfin) <- start_proc proc_gpi- (\x -> putStrLn ("Global pi_k (nucleotide diversity): "++show x++"\n"))+ (\x -> printf "Global pi_k (nucleotide diversity): %d\n" (round x :: Integer)) (ppi,pfin) <- start_proc (proc_gppi o) out_gppi (fi,ffin) <- start_proc (proc_gfst o) out_gfst let run (M r s:rs) = do@@ -92,6 +94,7 @@ zero = repeat (repeat (0,0)) plus (a,c) (b,d) = (a+b,c+d) deepSeq x | x == x = x+ | True = error (show x) -- | Calculate and print global pairwise Fst out_gfst :: [[(Double,Double)]] -> IO ()@@ -104,11 +107,12 @@ putStrLn $ unwords $ map (\(t,w) -> printf "%.3f" ((t-w)/t)) l go (i+1) ls go _ [] = return ()+ -- outputs one line too many? -- -------------------------------------------------- data UniVar = UV { _count :: {-# UNPACK #-} !Int, _sum1, _sum2 :: {-# UNPACK #-} !Double }-+ deriving Show add_uv :: UniVar -> Double -> UniVar add_uv (UV c s s2) d = UV (c+1) (s+d) (s2+d*d) @@ -118,13 +122,14 @@ where f (MPR sup _ _ _ cts) (uv,cur) = let new = Metrics.ppi_params cts cov = sum $ map covC cts- nu = add_uv uv $ fromIntegral cov- nc = if sup || (max_cov o > 0 && cov > max_cov o) || cov < min_cov o - then cur else deepSeq $ zipWith (zipWith plus) cur new + ign = sup || (max_cov o > 0 && cov > max_cov o) || cov < 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) zero = (UV 0 0 0, repeat (repeat 0)) plus a b = if isNaN b then a else a+b deepSeq x | x == x = x+ | True = error (show x) -- | Calculate and print global pairwise ND -- Todo: divide by genome size @@ -137,7 +142,7 @@ go (i+1) ls go _ [] = return () putStrLn "Coverage statistics:"- _ <- printf " covered sites: %d\n avg. cover: %.2f\n std. dev.: %.2f\n\n" n (s1/n') (sqrt (s2/n'-(s1*s1)/(n'*n')))+ _ <- printf " observed variant sites: %d\n avg. cover: %.2f\n std. dev.: %.2f\n\n" n (s1/n') (sqrt (s2/n'-(s1*s1)/(n'*n'))) putStrLn "Pairwise Nucleotide Diversities:" putStrLn (" "++ concat [ " s"++show i | i <- [1..length (head xs)]]) go 1 xs@@ -171,7 +176,10 @@ ,if Options.chi2 o then "\tChi²" else "" ,if Options.conf o then concat ["\tCI "++show n | n <- [1..(length cs)]] else "" ,if Options.pconf o then "\tpconf" else ""- ,if Options.ds o then "\tdelta-sigma" else ""+ ,if Options.ds o then "\tds-agresti" else ""+ ,if Options.dsw o then "\tds-wald" else ""+ ,if Options.nd_all o then "\tNuc divs\tNd tot" else ""+ ,if Options.maf o then "\tMAF" else "" ,if Options.esi o then "\tESI" else "" ,if Options.variants o then "\tVariants" else "" ,"\n" @@ -186,7 +194,7 @@ showPile :: Options -> MPileRecord -> B.ByteString showPile _ (MPR _ _ _ _ []) = error "Pileup with no data?"-showPile o mpr = if suppress o && ignore mpr then B.empty else (B.concat+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 , 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))@@ -194,13 +202,20 @@ , when (Options.conf o) (conf_all $ counts mpr) , when (Options.pconf o) ("\t"++dsconf_pairs 0.01 (counts mpr)) , when (Options.ds o) ("\t"++(unwords $ map (\x -> if x>=0 then printf "%.2f" x else " - ") $ ds_all 2.326 $ counts mpr))- , when (Options.esi o) ("\t"++concat [ concat [printf " %2.2f" (ESIV.esiv 1.64 0.01 c1 c2) | c2 <- rest] | (c1:rest) <- Data.List.tails (counts mpr)])--- , when (Options.mafci o) ("\t"++concat counts mpr)...something+ , when (Options.dsw o) ("\t"++(unwords $ map (\x -> if x>=0 then printf "%.2f" x else " - ") $ dsw_all 2.326 $ counts mpr))+ , when (Options.nd_all o) ("\t"++(unwords $ map (printf "%.2f" . nd) $ counts mpr)++"\t"++printf "%.3f" (nd $ ptSum $ counts mpr))+ , when (Options.maf o) ("\t"++unwords (map (printf "%.2f" . Metrics.maf) (counts mpr)))++ -- Between pairs of samples+ , when (Options.esi o) (pairwise (ESIV.esiv 1.64 0.01) (counts mpr))+ , when (Options.variants o) ("\t"++showV (counts mpr)) , B.pack "\n" ]) where when p s = if p then B.pack s else B.empty- + pairwise :: (Counts -> Counts -> Double) -> [Counts] -> String+ pairwise f cs = "\t"++concat [ concat [printf " %2.2f" (f c1 c2) | c2 <- rest] | (c1:rest) <- Data.List.tails cs]+ -- | The default output, with only coverage statistics default_out :: MPileRecord -> B.ByteString default_out (MPR _ chr pos ref stats) =
+ src/Sparks.hs view
@@ -0,0 +1,124 @@+{-# Language DeriveDataTypeable #-}++-- Draw Unicode sparklines, using code points 0x2581 to 0x2588++module Main where++import System.Console.ANSI+import Data.Char (chr)+import Data.List (sort)+import System.Console.CmdArgs++import ESIV+import Count+import VExtr (makeConsensus, Format(IUPAC))+import MPileup+-- import Count+import qualified Data.ByteString.Lazy.Char8 as BL++data Options = Test+ | Disp+ | Info -- maybe select samples to contrast? color?+ | Cite+ deriving (Typeable,Data)++test, disp, info, cite :: Options+test = Test &= details ["Prints a test string.",+ "This is useful for checking that your terminal supports the graphical characters needed to generate sparklines. It should look like gradual transitions from one color to the next."]+disp = Disp &= details ["Show per sample consensus sequences.","Reads mpile-formatted input, and shows the consensus for each sample, one line each. This can potentially result in very long lines, use 'samtools mpileup' with the '-r' option to restrict output."]+info = Info &= details ["Show expected information values.","This shows the information value for observing each allele, i.e. how diagnostic each site is between the two samples."]+cite = Cite &= details ["Output citation information."]++main :: IO ()+main = do+ opts <- cmdArgsRun $ cmdArgsMode $ modes [test,disp,info,cite]+ &= summary "Visualize information from 'samtools mpileup' as sparklines"+ &= program "sparks"+ inp <- BL.getContents+ let ms = map readPile1 $ BL.lines inp+ mapM_ putStrLn $ case opts of+ Test -> teststr+ Disp -> sparklines ms+ Info -> infoline ms+ Cite -> citestr++citestr :: [String]+citestr = ["If you use this program, please cite:"+ ," BMC Genomics 2014, 15(Suppl 6):S20"+ ," http://www.biomedcentral.com/1471-2164/15/S6/S20"+ ]++teststr :: [String]+teststr = [ concat [count2char [5,b,0,0] | b <- [0..10]]+ ++concat [count2char [0,5,b,0] | b <- [0..10]]+ ++ concat [count2char [0,0,5,b] | b <- [0..10]]+ ++ count2char [0,0,0,10] ++ setSGRCode []]++infoline :: [MPileRecord] -> [String]+infoline ms = (cons:esivs)+ where+ cons = makeConsensus (IUPAC,1,5) ms+ esivs = map (concat . (++ rst)) $ transpose $ map esivstr ms+ rst = [setSGRCode []]++esivstr :: MPileRecord -> [String]+esivstr m = let+ (s1:s2:_) = counts m+ [t1,t2] = map covC [s1,s2]+ [l1,l2] = map toList [s1,s2]+ [maxpos1,maxpos2] = map (snd . last . sort) [zip xs [0..3] | xs <- [l1,l2]]+ es = [esiv1 1.64 0.01 (x,t1-x) (y,t2-y) | (x,y) <- zip l1 l2]+ sorted = sort (zip es [0..3::Int])+ echar :: Color -> Color -> Double -> String+ echar bg fg f = setSGRCode [SetColor Foreground Dull fg,SetColor Background Dull bg]+ ++ if f<0 then [chr (0x2589-max 1 (min 8 (round (20*abs f))))]+ else [chr (0x2580+max 1 (min 8 (round (10*f))))]+ toCol = ([Red,Blue,Green,Yellow]!!)+ in if sum (map abs es) > 0.1+ then let+ (plus,ppos) = last sorted+ (minus,mpos) = head sorted+ in [echar Black (toCol ppos) plus, echar (toCol mpos) Black minus]+ else [echar Black (toCol maxpos1) 0, echar (toCol maxpos2) Black 99]++sparklines :: [MPileRecord] -> [String]+sparklines = map sparkline . transpose . map counts++transpose :: [[a]] -> [[a]]+transpose xs+ | any null xs = []+ | otherwise = map head xs : transpose (map tail xs)++sparkline :: [Counts] -> String+sparkline = (++setSGRCode []) . concatMap (count2char . toList)++-- display a frequency spectrum as a sparkline char++count2char :: [Int] -> String+count2char [0,0,0,0] = setSGRCode [SetColor Background Dull Black] ++ " "+count2char cts = let+ ((a,b):(c,d):_) = reverse $ sort $ zip cts [0..3]+ [(majpos,major),(minpos,minor)] = sort [(b,a),(d,c)]+ freq = fromIntegral major / fromIntegral (major + minor)+ in f2char ("acgt"!!majpos) ("acgt"!!minpos) freq+ -- show (majpos,minpos,majfreq,counts) -- ++-- From colors and frequency, calculate the appropriate sparkline char+-- f should range be in the range [0..9]+f2char :: Char -> Char -> Double -> String+f2char c1 c2 f+ | f<=0 = setSGRCode [SetColor Foreground Dull Black,SetColor Background Dull bg] ++ [c2]+ | f>=1 = setSGRCode [SetColor Foreground Dull Black,SetColor Background Dull fg] ++ [c1]+ | otherwise = cols ++ [chr (0x2580+round (9*f))]+ where+ cols = setSGRCode [SetColor Foreground Dull fg,SetColor Background Dull bg]+ fg = toCol c1+ bg = toCol c2+ toCol 'a' = Red -- A+ toCol 'c' = Blue -- C+ toCol 'g' = Green -- G + toCol 't' = Yellow -- T+ toCol x = error ("illegal character "++show x)++reset :: IO ()+reset = setSGR []
+ src/VExtr.hs view
@@ -0,0 +1,111 @@+{-# Language DeriveDataTypeable #-}++module VExtr where++import MPileup+import Count+import qualified Data.ByteString.Lazy.Char8 as BL++import System.Console.CmdArgs++data Options = Opts { infile, outfile :: Maybe FilePath+ , format :: Format+ , fasta :: Bool+ , mincount :: Int+ , minfreq :: Int } deriving (Data,Typeable)++opts :: Options+opts = Opts + { infile = Nothing &= args &= typFile+ , outfile = Nothing &= help "output file"+ , format = IUPAC &= help "output X, N, or [a/b] instead of IUPAC codes for variable sites"+ , 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"+ &= summary "Extract consensus sequence from pooled sequences"+ &= details ["Examples:", ""+ ,"Read input from a pipe, output IUPAC codes:"+ ,"", " samtools mpileup -f ref.fasta reads.bam | vextr --format=iupac", ""+ ,"Read input from a file, create consensus FASTA sequence:"+ ,"", " vextr input.mpile --fasta -o output.fasta", ""+ ,"If you use this program, please cite:"+ ," BMC Genomics 2014, 15(Suppl 6):S20"+ ," http://www.biomedcentral.com/1471-2164/15/S6/S20"+ ]++data Format = Xs | IUPAC | Regex deriving (Data,Typeable,Show)++main :: IO ()+main = do+ o <- cmdArgs opts+ inp <- case infile o of Nothing -> BL.getContents+ Just f -> BL.readFile f+ let ms = map readPile1 $ BL.lines inp+ outf = case outfile o of Nothing -> putStr+ Just f -> writeFile f+ gen = if fasta o then makeFasta else makeConsensus+ outf $ gen (format o,mincount o,minfreq o) ms+++makeFasta :: (Format,Int,Int) -> [MPileRecord] -> String+makeFasta fi ms = let+ header = case ms of+ (m1:_) -> '>':BL.unpack (chrom m1)++":"++BL.unpack (cpos m1)+ [] -> ""+ breaks str = case splitAt 60 str of+ (rest,"") -> [rest]+ (this,more) -> this : breaks more+ in unlines (header:breaks (makeConsensus fi ms))++makeConsensus :: (Format,Int,Int) -> [MPileRecord] -> String+makeConsensus (iup,mct,mfq) = concatMap (fixiup iup . selectChar mct mfq . ptSum . counts)+ -- todo: include variants++-- this doesn't work so well with high coverage/many libraries++-- | Optionally change from IUPAC code to X or regex+fixiup :: Format -> Char -> String+fixiup iup c | c `elem` "ACGTacgtNn" = [c]+ | otherwise = case iup of+ Xs -> "X"+ IUPAC -> [c]+ Regex -> case c of+ 'R' -> "[A/G]"+ 'Y' -> "[C/T]"+ 'S' -> "[C/G]"+ 'W' -> "[A/T]"+ 'K' -> "[G/T]"+ 'M' -> "[A/C]"+ 'B' -> "[C/G/T]"+ 'D' -> "[A/G/T]"+ 'H' -> "[A/C/T]"+ 'V' -> "[A/C/G]"+ x -> [x]++-- | Convert allele counts into IUPAC character+selectChar :: Int -> Int -> Counts -> Char+selectChar mct mfq ss = case toList ss of+ [0,0,0,0] -> 'n'+ [_,0,0,0] -> 'A'+ [0,_,0,0] -> 'C'+ [0,0,_,0] -> 'G'+ [0,0,0,_] -> 'T'++ [x,0,y,0] -> maybeWild x y 'a' 'g' 'R'+ [0,x,0,y] -> maybeWild x y 'c' 't' 'Y'+ [0,x,y,0] -> maybeWild x y 'c' 'g' 'S'+ [x,0,0,y] -> maybeWild x y 'a' 't' 'W'+ [0,0,x,y] -> maybeWild x y 'g' 't' 'K'+ [x,y,0,0] -> maybeWild x y 'a' 'c' 'M'++ [0,_,_,_] -> 'B'+ [_,0,_,_] -> 'D'+ [_,_,0,_] -> 'H'+ [_,_,_,0] -> 'V'++ _ -> 'N'+ where maybeWild x y c1 c2 c3 =+ let xok = x > mct && x > (x+y)*mfq`div`100+ yok = y > mct && y > (x+y)*mfq`div`100+ in if xok && yok then c3 else if xok then c1 else if yok then c2 else 'n'
varan.cabal view
@@ -1,5 +1,5 @@ Name: varan-Version: 0.3+Version: 0.5 License: GPL Cabal-Version: >= 1.6 Build-Type: Simple@@ -17,6 +17,21 @@ Main-Is: Varan.hs Other-Modules: AgrestiCoull, MPileup, RandomSelect, Variants, Metrics, Options, ParMap, Process, Count, ESIV 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+ Other-Modules: MPileup, Count+ Build-Depends: base >= 4 && < 5, bytestring, cmdargs+ Ghc-Options: -rtsopts -Wall -threaded -main-is VExtr++Executable sparks+ Hs-Source-Dirs: src+ Main-Is: Sparks.hs+ 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