diff --git a/src/Count.hs b/src/Count.hs --- a/src/Count.hs +++ b/src/Count.hs @@ -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 diff --git a/src/ESIV.hs b/src/ESIV.hs --- a/src/ESIV.hs +++ b/src/ESIV.hs @@ -20,15 +20,17 @@ esiv z epsilon c1 c2 = let t1 = covC c1 t2 = covC c2 - esiv1 (a1,b1) (a2,b2) - | j1+epsilonj2+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+epsj2+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) diff --git a/src/MPileup.hs b/src/MPileup.hs --- a/src/MPileup.hs +++ b/src/MPileup.hs @@ -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] diff --git a/src/Metrics.hs b/src/Metrics.hs --- a/src/Metrics.hs +++ b/src/Metrics.hs @@ -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? diff --git a/src/Options.hs b/src/Options.hs --- a/src/Options.hs +++ b/src/Options.hs @@ -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) diff --git a/src/Process.hs b/src/Process.hs --- a/src/Process.hs +++ b/src/Process.hs @@ -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) = diff --git a/src/Sparks.hs b/src/Sparks.hs new file mode 100644 --- /dev/null +++ b/src/Sparks.hs @@ -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 [] diff --git a/src/VExtr.hs b/src/VExtr.hs new file mode 100644 --- /dev/null +++ b/src/VExtr.hs @@ -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' diff --git a/varan.cabal b/varan.cabal --- a/varan.cabal +++ b/varan.cabal @@ -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