packages feed

varan (empty) → 0.3

raw patch · 13 files changed

+855/−0 lines, 13 filesdep +basedep +bytestringdep +cmdargssetup-changed

Dependencies added: base, bytestring, cmdargs, mtl, parallel, random, statistics

Files

+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ src/AgrestiCoull.hs view
@@ -0,0 +1,15 @@+module AgrestiCoull where++-- zscores: 1.654 is one-sided 95%+--          1.96 is two-sided+--          2.326 is one-sided 99%++confidenceInterval :: Double -> Int -> Int -> (Double,Double)+confidenceInterval z succs fails = let+  s = fromIntegral succs+  f = fromIntegral fails+  z2 = z * z+  nest = s + f + z2+  pest = (s + z2/2)/nest+  delta = z*sqrt(pest*(1-pest)/nest)+  in (pest-delta,pest+delta)
+ src/Count.hs view
@@ -0,0 +1,67 @@+-- High performance counting data structure++module Count (Counts(..)+             , addA, addC, addG, addT, addV+             , addA_, addC_, addG_, addT_                                       +             , getA, getC, getG, getT+             , covC, ptAdd, ptSum+             , toList, sumList +             ) where++import Variants       +import Data.Int+import Data.List (foldl1')++data Counts = C { getA_, getC_, getG_, getT_ :: {-# UNPACK #-} !Int32+                , getV :: ![Variant]+                }++getA, getC, getG, getT :: Integral i => Counts -> i+getA = fromIntegral . getA_+getC = fromIntegral . getC_+getG = fromIntegral . getG_+getT = fromIntegral . getT_+{-# INLINE getA #-}+{-# INLINE getC #-}+{-# INLINE getG #-}+{-# INLINE getT #-}++addA, addC, addG, addT :: 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 }+{-# INLINE addA #-}+{-# INLINE addC #-}+{-# INLINE addG #-}+{-# INLINE addT #-}++addA_, addC_, addG_, addT_ :: 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 }+++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)++ptAdd :: Counts -> Counts -> Counts+ptAdd a b = a `addA` (getA_ b) `addC` (getC_ b) `addG` (getG_ b) `addT` (getT_ b)++ptSum :: [Counts] -> Counts+ptSum = foldl1' ptAdd++-- Convert 'Counts' to a list of allele counts+toList :: Num a => Counts -> [a]+toList x = map fromIntegral [getA_ x,getC_ x,getG_ x,getT_ x]++-- | Pointwise summation of the input lists+sumList :: Num a => [[a]] -> [a]+sumList = foldr (zipWith (+)) [0,0,0,0]+{-# DEPRECATED sumList "use ptSum" #-}+
+ src/ESIV.hs view
@@ -0,0 +1,34 @@+{- +  ESIV - Expected Site Information Value from SNPs++  Information value is the log odds for each allele, times the probablity of observing it+  not using any prior (or prior is 50/50)++  Biallelic:+  abs (avg(p1,p2)*log(p1/p2) - (1-avg(p1,p2))*log((1-p1)/(1-p2))+-}++module ESIV where++import Count+import AgrestiCoull++-- given a z-score `z` for confidence interval, and a+-- minimum error rate `epsilon`, calculate the ESIV conservatively+-- using the conf interval boundaries as frequencies.+esiv :: Double -> Double -> Counts -> Counts -> Double+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)+    | 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))+
+ src/MPileup.hs view
@@ -0,0 +1,82 @@+{-# Language BangPatterns #-}+module MPileup (Counts(..), readPile1, toList, major_allele, by_major_allele, showC, showV, sumList, MPileRecord(..)) where++import Data.Char (toUpper)+import Data.List (intercalate,nub,elemIndex)+import qualified Data.ByteString.Lazy.Char8 as B+import Variants hiding (parse)+import Count++data MPileRecord = MPR { ignore :: !Bool+                       , chrom, cpos :: !B.ByteString+                       , refnuc :: !Char+                       , counts :: ![Counts]+                       }++-- convert counts to major/non-major allele counts+by_major_allele :: [Counts] -> [(Int,Int)] -- always length 2+by_major_allele cs = let+  ls = map toList cs+  s  = toList $ ptSum cs :: [Int]+  Just i = elemIndex (maximum s) s+  in map (\l -> (l!!i,sum l-l!!i)) ls++-- pick out major allele in first count, and output number of same/different+major_allele :: Counts -> Counts -> ((Int,Int),(Int,Int))+major_allele x y =+  let s1 = toList x+      s2 = toList y+      m = maximum s1+      Just i = elemIndex m s1+  in ((m, sum s1-m),(s2!!i,sum [s2!!j | j <- [0..3], j /= i]))++-- count major allele in first sample+-- return flag whether informative, chrom, pos, ref, +readPile1 :: B.ByteString -> MPileRecord +readPile1 = parse1 . B.split '\t'  -- later samtools sometimes outputs empty strings in columns+  where+    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)+    +    -- set the ignore flag if we only see one or zero alleles    +    ign xs = (length $ filter (/=(0::Int)) $ toList $ ptSum xs) <= 1++    triples _ [] = []+    triples ref (_cnt:bases:_quals:rest) = let this = parse ref (C 0 0 0 0 []) (B.map toUpper bases) +                                           in this `seq` this : triples ref rest+    triples _ _ = error "triples: incorrect number of columns"+    +    -- note: does not (yet) deal with '<' and '>', which apparently can occur+    -- perhaps it would be faster to do this as a sequence of BS ops, which would fuse?+    parse :: Char -> Counts ->  B.ByteString -> Counts+    parse ref !cts bs = case B.uncons bs of+        Nothing -> cts+        Just (c_,str_) -> p c_ str_+        where p !c !str +                | c == '.' || c == ',' = parse ref (addRef cts 1) str+                | c == 'A'             = parse ref (addA_ cts 1) str+                | 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 == '-' || 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 }++-- | Show SNP counts and coverage+showC :: Counts -> (String,Int)+showC x = (" "++(intercalate ":" $ map show (toList x :: [Int])),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])
+ src/Metrics.hs view
@@ -0,0 +1,188 @@+-- | Calculate various metrics/statistics+module Metrics where++import AgrestiCoull+import MPileup (by_major_allele)+import Count+import Statistics.Distribution+import Statistics.Distribution.ChiSquared+import Data.List (foldl1')++-- | Calculate vector angle between allele frequencies.  This is +--   similar to `dist`, but from 1 (equal) to 0 (orthogonal)+angle :: Counts -> Counts -> Double+angle c1' c2' = let+  (c1,c2) = (toList c1', toList c2')+  vnorm = sqrt . sum . map ((**2))+  in sum $ zipWith (*) (map (/vnorm c1) c1) (map (/vnorm c2) c2)++-- 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 [] = []++-- 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+        go [] = []+fst_params [] = []++heteroz :: Counts -> Counts -> (Double,Double)+heteroz c1 c2 = let+  c1s = fromIntegral $ covC c1+  c2s = fromIntegral $ covC c2+  total = c1s + c2s+  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+  h_tot = hz (c1 `ptAdd` c2)+  h_subs = (hz c1*c1s + hz c2*c2s)/total+  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+  c1s = sum c1+  c2s = sum c2+  total = c1s + c2s+  hz :: [Double] -> Double+  hz xs' = let s = sum xs' in 1 - sum (map ((**2) . (/s)) xs')+  h_tot, h_subs :: Double++  h_tot = hz $ zipWith (+) c1 c2+  h_subs = (hz c1*c1s + hz c2*c2s)/total+  in if c1s == 0 || c2s == 0 || h_tot == 0 then (0,0) +     else (h_tot,h_subs)++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+  h_subs, weights :: [Double]+  h_tot = hz (ptSum xs)+  h_subs = map hz 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+f_st_ :: [Counts] -> Double+f_st_ cs = let+  cs' = map toList cs+  -- hm, er ikke dette bare nuc div?+  hz :: [Double] -> Double+  hz xs' = let s = sum xs'+           in 1 - sum (map ((**2) . (/s)) xs')+  h_tot :: Double+  h_tot = hz $ sumList $ cs'+  h_subs, weights :: [Double]+  h_subs = map hz cs'+  weights = let total = sum $ concat cs'+            in [sum c/total | c <- cs']+  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.++pi_k :: [Counts] -> Double+pi_k 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]+pi_freqs x = let s = fromIntegral $ covC x+             in [fromIntegral (getA_ x)/s,fromIntegral (getC_ x)/s,fromIntegral (getG_ x)/s,fromIntegral (getT_ x)/s]++-- Or, equivalent+pi_k_alt :: [Counts] -> Double+pi_k_alt cs' = let+  cs = map toList cs'+  no_diff = sum $ foldr1 (zipWith (*)) cs+  all_pairs = product $ map sum cs+  in (all_pairs - no_diff) / all_pairs++-- should probably include a warning if more than 20% of cells < 5 expected or some such+pearsons_chi² :: [[Int]] -> Double+pearsons_chi² t = let+  cols   = map sum+  rows x = if all null x then []+           else sum (map head x) : rows (map tail x)+  exps   = [[ fromIntegral (r*c) / fromIntegral (sum $ rows t) | r <- rows t] | c <- cols t ]+  chi    = sum [ (fromIntegral a-b)^(2::Int)/b | (a,b) <- zip (concat t) (concat exps) ]+  df     = (length t-1)*(length (head t) - 1)+  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+conf :: Counts -> Counts -> String+conf x y = let+  s1 = covC x  -- don't count structural variants+  s2 = covC y+  in [overlap (getA x,s1-getA x) (getA y,s2-getA y)+     ,overlap (getC x,s1-getC x) (getC y,s2-getC y)+     ,overlap (getG x,s1-getG x) (getG y,s2-getG y)+     ,overlap (getT x,s1-getT x) (getT y,s2-getT y)+     ]++-- | Helper function for conf+overlap :: (Int,Int) -> (Int,Int) -> Char+overlap (succ1,fail1) (succ2,fail2) =+  let (i1,j1) = confidenceInterval 1.65 succ1 fail1+      (i2,j2) = confidenceInterval 1.65 succ2 fail2+  in if i2>=j1 || i1>=j2 then+       let (k1,l1) = confidenceInterval 2.326 succ1 fail1+           (k2,l2) = confidenceInterval 2.326 succ2 fail2+       in if k2>=l1 || k1>=l2 then '*' else '+'+     else '.'++-- | Use AgrestiCoull to calculate significant difference from+--   a 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+  -- attempt to smooth errors by subtracting one, seems to work:+  -- m x = max (x-1) 0+  -- rm_err x = C (0 `addA` (m $ getA x) `addC` (m $ getC x) `addG` (m $ getG x) `addT` (m $ getT x)) []+  in concat ["\t"++x | x <- map (conf (ptSum cs')) cs']++dsconf_pairs :: Double -> [Counts] -> String+dsconf_pairs e cs = go $ by_major_allele cs+  where go (x:xs) = " " ++ map (ds x) xs ++ go xs+        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+delta_sigma :: Double -> (Int,Int) -> (Int,Int) -> Double+delta_sigma z (s1,f1) (s2,f2) =+  let (i1,j1) = confidenceInterval z s1 f1+      (i2,j2) = confidenceInterval z s2 f2+      mu1 = i1+j1 -- all values are times two (so it cancels out)+      mu2 = i2+j2+      sd1 = j1-i1+      sd2 = j2-i2+  in (abs (mu2-mu1) - (sd1+sd2))/2++-- use on output from by_major_allele+ds_all :: Double -> [Counts] -> [Double]+ds_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 (delta_sigma sig)) pairs++-- | Calculate distance between approximate distributions+-- in terms of their standard deviation.  Perhaps use binomial distribution directly?+-- This is a z-score, i.e. score of 2 means that the 95% CIs barely overlap.+ci_dist :: Double -> (Int,Int) -> (Int,Int) -> Double+ci_dist z (s1,f1) (s2,f2) =+  let (i1,j1) = confidenceInterval z s1 f1+      (i2,j2) = confidenceInterval z s2 f2+      mu1 = i1+j1 -- all values are times two (so it cancels out)+      mu2 = i2+j2+      sd1 = j1-i1+      sd2 = j2-i2+  in if sd1+sd2 == 0 then 0 else abs (mu2-mu1)/(sd1+sd2)
+ src/Options.hs view
@@ -0,0 +1,69 @@+{-# Language DeriveDataTypeable #-}++module Options (Options(..), getArgs) where++import System.Console.CmdArgs+import qualified Data.ByteString.Lazy.Char8 as BL+import System.IO++-- data SnpMode = Default | Pairwise | Groups deriving (Typeable,Data)++data Options = Opts +  { suppress, variants+  , chi2, f_st, pi_k, conf, ds, esi, pconf :: Bool+  , input, output :: FilePath+  , global :: Bool+  , threads :: Int+    , min_cov, max_cov :: Int+  } deriving (Typeable,Data)++defopts :: Options+defopts = Opts +  { +  -- General parameters+    output = "" &= help "output file name" &= typFile+  , suppress = False &= help "omit non-variant lines from output"+  , variants = False &= help "output list of non-SNP variants"+  , threads = 100 &= help "queue lenght for parallel execution"+  , min_cov = 0     &= help "minimum coverage to include"+  , max_cov = 32000 &= help "maximum coverage to include"+  , input  = [] &= args &= typFile+  +  , global = False &= help "calculate global statistics"++  -- Overall statistics+  , chi2   = False &= help "calculate chi² probability" &= ignore+  , f_st   = False &= help "estimate fixation index, F_st" &= name "fst"+  , pi_k   = False &= help "estimate nucleotide diversity, Pi_k"++  -- 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)++  -- 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",""+               ]    ++getArgs :: IO (IO BL.ByteString,Options)+getArgs = do+  opts <- cmdArgs defopts+  inp  <- get_input opts+  return (inp,opts)+  +-- | Determine where and how to read the input mpile data.+get_input :: Options -> IO (IO BL.ByteString)+get_input o = do+  interactive <- hIsTerminalDevice stdin+  return $ if null (input o) then +             if interactive then error "Cowardly refusing to read input from terminal.\nIf you really want this, specify '-' as input.\n(or use --help for help)."+             else BL.getContents+           else if input o == "-" then BL.getContents +                else BL.readFile (input o)
+ src/ParMap.hs view
@@ -0,0 +1,25 @@+module ParMap (parMap) where++import Control.Concurrent+import Control.Monad+import System.IO.Unsafe++parMap :: Int -> (b -> a) -> [b] -> IO [a]+parMap t f xs = do+  outs <- replicateM t newEmptyMVar+  ins  <- replicateM t newEmptyMVar+  sequence_ [forkIO (worker i o) | (i,o) <- zip ins outs]+  _ <- forkIO $ sequence_ [putMVar i (f x) | (i,x) <- zip (cycle ins) xs]+  sequence' [takeMVar o | (o,_) <- zip (cycle outs) xs]++worker :: MVar a -> MVar a -> IO b+worker imv omv = forever $ do +  ex <- takeMVar imv +  ex `seq` putMVar omv ex++sequence' :: [IO a] -> IO [a]+sequence' [] = return []+sequence' (m:ms) = do+  x <- m+  xs <- unsafeInterleaveIO $ sequence' ms+  return (x:xs)
+ src/Process.hs view
@@ -0,0 +1,214 @@+{-# Language BangPatterns #-}++module Process (proc_fused, run_procs, showPile') where++import Options+import ParMap+import MPileup+import Metrics+import Count+import ESIV++import Data.List (tails)+import qualified Data.ByteString.Char8 as B+import qualified Data.ByteString.Lazy.Char8 as BL+import Text.Printf+import System.IO+import Control.Concurrent++proc_fused :: Options -> [BL.ByteString] -> IO ()+proc_fused o (l:ls) = do+  outh <- if null (output o) || output o == "-" then return stdout else openFile (output o) WriteMode+  B.hPutStr outh $ gen_header o $ readPile1 l+  if threads o > 1 +    then mapM_ (B.hPutStr outh) =<< parMap (threads o) (showPile o . readPile1) (l:ls)+    else mapM_ (B.hPutStr outh) $ map (showPile o . readPile1) (l:ls)+  hClose outh+proc_fused _ [] = error "No input?"++-- | Runs a set of processes, distributes each MPileRecord to them+--   and runs the finalizer (collecting and outputting the results)  +run_procs :: Options -> [MPileRecord'] -> IO ()+run_procs _ [] = putStrLn "No input records found! (Or try: varan --help)"+run_procs o recs@(M r1 _:_) = do+  -- initialize default output+  let use_stdout = null (output o) || output o == "-"+  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"))+  (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+        push_procs (Just r) [gi,ppi,fi]+        B.hPutStr outh s+        run rs+      run [] = do+        push_procs Nothing [gi,ppi,fi]+        putStrLn ""+        sequence_ [gfin,pfin,ffin]+  run recs++-- | The main process (in the first parameter) reads from 'inv' and puts the result+--   in 'out' when it's done. It returns the input MVar, and the IO action that+--   runs the finalizer (second parameter) on the final result. +start_proc :: (MVar inv -> MVar out -> IO ()) -> (out -> IO ()) -> IO (MVar inv, IO ())+start_proc f g = do+  imv <- newEmptyMVar+  omv <- newEmptyMVar+  _ <- forkIO $ f imv omv+  return (imv, takeMVar omv >>= g)++-- | Shortcut for feeding a datum to a list of process inputs.+push_procs :: a -> [MVar a] -> IO ()+push_procs r mvs = mapM_ (\m -> putMVar m r) mvs++proc_fold :: a -> (MPileRecord -> a -> a) -> MVar (Maybe MPileRecord) -> MVar a -> IO ()+proc_fold z f mv e = go z +  where go !c = do+          v <- takeMVar mv+          case v of+            Nothing -> putMVar e c +            Just x  -> go (f x c)++-- --------------------------------------------------+-- Actual output calculations++-- | Calculate global nucleotide diversity+proc_gpi :: MVar (Maybe MPileRecord) -> MVar Double -> IO ()+proc_gpi = proc_fold 0.0 f+  where f mpr cur =+          let p = Metrics.pi_k (counts mpr)+          in if ignore mpr || isNaN p then cur else cur+p++-- | Collect variation within and between for global pairwise Fst+proc_gfst :: Options -> MVar (Maybe MPileRecord) -> MVar [[(Double, Double)]] -> IO ()+proc_gfst o = proc_fold zero f+  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 +             then cur else deepSeq $ zipWith (zipWith plus) cur new+        zero = repeat (repeat (0,0))+        plus (a,c) (b,d) = (a+b,c+d)+        deepSeq x | x == x = x++-- | Calculate and print global pairwise Fst+out_gfst :: [[(Double,Double)]] -> IO ()+out_gfst xs = do+  putStrLn "Pairwise FST values:"+  putStrLn (" "++ concat [ "    s"++show (i+1) | i <- [1..length $ head xs]])+  go 1 xs+  where go i (l:ls) = do +          putStr ("s"++show i++replicate (6*i-4) ' ')+          putStrLn $ unwords $ map (\(t,w) -> printf "%.3f" ((t-w)/t)) l+          go (i+1) ls+        go _ [] = return ()++-- --------------------------------------------------++data UniVar = UV { _count :: {-# UNPACK #-} !Int, _sum1, _sum2 :: {-# UNPACK #-} !Double }++add_uv :: UniVar -> Double -> UniVar+add_uv (UV c s s2) d = UV (c+1) (s+d) (s2+d*d)++-- | Collect nucleotide diversity within and between for global pairwise ND (pi)+proc_gppi :: Options -> MVar (Maybe MPileRecord) -> MVar (UniVar,[[Double]]) -> IO ()+proc_gppi o = proc_fold zero f+  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 +          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++-- | Calculate and print global pairwise ND+--   Todo: divide by genome size        +out_gppi :: (UniVar,[[Double]]) -> IO ()+out_gppi (UV n s1 s2,xs) = do+  let n' = fromIntegral n+      go i (l:ls) = do +          let s = show i in putStr ("s"++s++replicate (7*i-3-length s) ' ')+          putStrLn $ unwords $ map (\t -> printf "%.4f" (t/n')) l+          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')))+  putStrLn "Pairwise Nucleotide Diversities:"+  putStrLn (" "++ concat [ "     s"++show i | i <- [1..length (head xs)]])+  go 1 xs+  putStrLn ""++-- | Calculating per site statistics writing to a file or standard out.+proc_default :: Options -> MVar (Maybe MPileRecord) -> MVar () -> IO ()+proc_default o imv omv = do+  let use_stdout = null (output o) || output o == "-"+  outh <- if use_stdout then return stdout else openFile (output o) WriteMode+  Just l <- takeMVar imv -- or fail!+  B.hPutStr outh $ gen_header o l+  B.hPutStr outh $ showPile o l+  let run = do+        ml <- takeMVar imv+        case ml of+          Nothing -> do+            if (not use_stdout) then hClose outh else return ()+            putMVar omv ()+          Just x -> do+            B.hPutStr outh $ showPile o x+            run+  run++-- generate the appropriate header, based on number of input pools+gen_header :: Options -> MPileRecord -> B.ByteString+gen_header o (MPR _ _ _ _ cs) = B.pack $ concat [+  standard+  ,if Options.f_st o then "\tF_st" else ""+  ,if Options.pi_k o then "\tPi_k" else ""+  ,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.esi  o then "\tESI" else ""+  ,if Options.variants o then "\tVariants" else ""+  ,"\n"  +  ]+  where+    standard = "#Target seq.\tpos\tref"++concat ["\tsample "++show x | x <- [1..(length cs)]]++"\tcover"++data MPileRecord' = M !MPileRecord !B.ByteString++showPile' :: Options -> MPileRecord -> MPileRecord'+showPile' o m = M m (showPile o m)++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+          [ 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))+          , 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.variants o) ("\t"++showV (counts mpr))+          , B.pack "\n"+          ])+  where when p s = if p then B.pack s else B.empty+        +-- | The default output, with only coverage statistics+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+          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)+
+ src/RandomSelect.hs view
@@ -0,0 +1,88 @@+-- Calculates p-values using a permutation test (not currently used)++module RandomSelect where++import System.Random+import Control.Monad.State+import Variants+import AgrestiCoull (confidenceInterval)++pval_max_count :: Int+pval_max_count = 1000++pval_acc :: Double+pval_acc = 0.01++pval :: RandomGen g => g -> ([Counts] -> Double) -> [Counts] -> (Double,Double)+pval g f cs = let thresh = f cs+                  xs = take pval_max_count $ rselect g cs+                  go suc tot (y:ys) = let (a,b) = confidenceInterval 1.0 suc (tot-suc) +                                      in if b-a <= pval_acc then fromIntegral suc/fromIntegral tot +                                         else go (suc + if f y >= thresh then 1 else 0) (tot+1) ys+                  go suc tot [] = fromIntegral suc/fromIntegral tot+              in if thresh > 0 then (thresh, go 0 0 xs) else (0,1)++type AlleleSample = (Int,Int,Int,Int)++-- generate an infinite stream of sampled allele distributions+rselect :: RandomGen g => g -> [Counts] -> [[Counts]]+rselect g cs = map fromLists $ fst $ runState sel (g,sum_al,pool_sizes)+  where (_tot,sum_al,pool_sizes) = count cs+        +fromLists :: [AlleleSample] -> [Counts]+fromLists = map f+  where f (a,b,c,d) = C a b c d []+        +sel :: RandomGen g => State (g,AlleleSample,[Int]) [[AlleleSample]]+sel = do+  a <- select1+  as <- sel+  return (a:as)+        +select1 :: RandomGen g => State (g,AlleleSample,[Int]) [AlleleSample]+select1 = do+  -- select randomly for each p_sz+  (g,s_al,p_sz) <- get+  let (g',g'') = split g+  put (g'',s_al,p_sz) -- restore state+  return $ pickNs g' p_sz s_al++-- pick a random sampling of alleles given pool sizes+pickNs :: RandomGen g => g -> [Int] -> AlleleSample -> [AlleleSample]+pickNs _ [] _als = [] -- verify that als is empty+pickNs g (p:ps) als = let+  (g1,g2) = split g+  in pickN' g1 als p : pickNs g2 ps als++-- sample alleles given randomgen, count, and sum alleles (i.e. distribution)+pickN' :: RandomGen g => g -> AlleleSample -> Int -> AlleleSample+pickN' g als = go (0,0,0,0) g+  where+        go acc _ 0 = acc+        go acc g1 cnt = +          let (i,g2) = randomR (1,sum' als) g1+          in go (add i als acc) g2 (cnt-1)+        sum' (a,b,c,d) = a+b+c+d+        +-- count total alleles, number of each variant, and number in each count+count :: [Counts] -> (Int,AlleleSample,[Int])+count cs = let xs = sumList . map toList $ cs+               ys = map (sum . toList) cs+               tuple [a,b,c,d] = (a,b,c,d)+           in (sum xs,tuple xs,ys)++add :: Int -> AlleleSample -> AlleleSample -> AlleleSample+add i (a,b,c,d) (x,y,z,w)+  | i <= a = (x+1,y,z,w)+  | i <= a+b = (x,y+1,z,w)+  | i <= a+b+c = (x,y,z+1,w)               +  | i <= a+b+c+d = (x,y,z,w+1)               +  | otherwise = error "too high i"++{-+-- remove allele number x from input+pick :: Int -> AlleleSample -> AlleleSample+pick x (a:as) | x <= a = 1 : map (const 0) as+              | otherwise = 0 : pick (x-a) as+pick _ [] = error "pick ran out of input, x too large"+-}
+ src/Varan.hs view
@@ -0,0 +1,23 @@+module Main where++import MPileup (readPile1)+import Process (proc_fused, run_procs, showPile')+import Options+import Control.Monad (when)+import ParMap++import qualified Data.ByteString.Lazy.Char8 as BL++main :: IO ()+main = do+  (inp,o) <- Options.getArgs+  lns <- BL.lines `fmap` inp+  when (null lns) $ error "No lines in input!"+  if not (global o)+    then proc_fused o lns -- faster?+    else do+      -- seems faster for few CPUs?+      let comb = showPile' o . readPile1+      recs <- if threads o > 1 then parMap (threads o) comb lns +              else return $ map comb lns+      run_procs o recs
+ src/Variants.hs view
@@ -0,0 +1,23 @@+{-# Language GADTs #-}++module Variants where+import Data.Char (isDigit)++data Variant where+  Nuc :: Char -> Variant+  Ins, Del :: String -> Variant +  deriving (Show,Eq)++-- manually inlined in the `MPileup` module.  Don't use this one.+parse :: Char -> String -> [Variant]+parse _ [] = []+parse ref (c:str) +      | c == '^'             = parse ref $ drop 1 str+      | c == '*' || c == '$' = parse ref str                               +      | c == '.' || c == ',' = Nuc ref : parse ref str+      | c == '-' || c == '+' = let (x,rest) = span isDigit str+                                   cnt = read x+                               in (if c=='+' then Ins else Del) (take cnt rest) +                                  : parse ref (drop cnt rest)+      | otherwise            = Nuc c : parse ref str+
+ varan.cabal view
@@ -0,0 +1,25 @@+Name:          varan+Version:       0.3+License:       GPL+Cabal-Version: >= 1.6+Build-Type:    Simple+Category:      Bioinformatics+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%).++Executable varan+    Hs-Source-Dirs:  src+    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++-- For parallel execution, we might want to add these, but they+-- seem to degrade performance on older GHC+-- Ghc-Options: -threaded -with-rtsopts=-N+