{-# Options_Ghc -fno-warn-unused-binds #-} -- | Calculate various metrics/statistics 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', tails, sort) -- | 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 -> 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 go [] = [] fst_params [] = [] -- | 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 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_w2 :: [Double] -> [Double] -> (Double, Double) heteroz_w2 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) -- | Simple calculation, samples represent equal populations (weights = 1/n) f_st :: [Counts] -> Double f_st xs = let 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 = 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 - equivalent to the above 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 -- | 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 = 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] 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. Output * or +, depending on significance for each allele. 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 -- 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 -- 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. 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 (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 -- 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? -- 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)