diff --git a/biohazard.cabal b/biohazard.cabal --- a/biohazard.cabal +++ b/biohazard.cabal @@ -1,5 +1,5 @@ Name: biohazard -Version: 0.6.5 +Version: 0.6.6 Synopsis: bioinformatics support library Description: This is a collection of modules I separated from various bioinformatics tools. The hope is to make @@ -14,10 +14,8 @@ Author: Udo Stenzel Maintainer: udo.stenzel@eva.mpg.de -Copyright: (C) 2010-2015 Udo Stenzel +Copyright: (C) 2010-2016 Udo Stenzel -Tested-With: GHC == 7.4.2, GHC == 7.6.3, GHC == 7.8.4, - GHC == 7.10.3, GHC == 8.0.1 Extra-Source-Files: man/man7/biohazard.7 man/man1/bam-meld.1 man/man1/bam-rewrap.1 @@ -31,19 +29,21 @@ Cabal-version: >=1.9.2 Build-type: Custom +-- Tested-With: GHC == 7.4.2, GHC == 7.6.3, GHC == 7.8.4, GHC == 7.10.3 source-repository head type: git - location: git://github.com/udo-stenzel/biohazard.git + location: https://bitbucket.org/ustenzel/biohazard.git source-repository this type: git - location: git://github.com/udo-stenzel/biohazard.git - tag: 0.6.5 + location: https://bitbucket.org/ustenzel/biohazard.git + tag: 0.6.6 Library Exposed-modules: Bio.Base, + Bio.Adna, Bio.Align, Bio.Bam, Bio.Bam.Evan, @@ -59,7 +59,6 @@ Bio.Bam.Trim, Bio.Bam.Writer, Bio.Genocall, - Bio.Genocall.Adna, Bio.Genocall.AvroFile, Bio.Genocall.Metadata, Bio.Iteratee, @@ -107,7 +106,7 @@ vector-th-unbox == 0.2.*, zlib >= 0.5 && < 0.7 - Ghc-options: -Wall -auto-all + Ghc-options: -Wall Hs-source-dirs: src Install-Includes: src/cbits/myers_align.h C-sources: src/cbits/myers_align.c @@ -120,15 +119,14 @@ -- Build-tools: -- Test-Suite test-biohazard - -- Ghc-options: -Wall -auto-all -threaded -rtsopts -with-rtsopts=-N + -- Ghc-options: -Wall -threaded -rtsopts -with-rtsopts=-N -- Type: exitcode-stdio-1.0 -- Main-is: test-biohazard.hs Executable redeye-dar Main-is: redeye-dar.hs - Ghc-options: -Wall -auto-all Hs-Source-Dirs: tools - Ghc-options: -Wall -auto-all -threaded -rtsopts -with-rtsopts=-N + Ghc-options: -Wall -threaded -rtsopts -with-rtsopts=-N Build-depends: async, base, biohazard, @@ -139,7 +137,7 @@ Executable redeye-div Main-is: redeye-div.hs - Ghc-options: -Wall -auto-all -threaded -rtsopts -with-rtsopts=-N + Ghc-options: -Wall -threaded -rtsopts -with-rtsopts=-N Hs-Source-Dirs: tools Build-depends: async, base, @@ -152,7 +150,7 @@ Executable redeye-pileup Main-is: redeye-pileup.hs - Ghc-options: -Wall -auto-all -threaded -rtsopts -with-rtsopts=-N + Ghc-options: -Wall -threaded -rtsopts -with-rtsopts=-N Hs-Source-Dirs: tools Build-depends: aeson, base, @@ -169,7 +167,7 @@ Executable redeye-single Main-is: redeye-single.hs - Ghc-options: -Wall -auto-all -rtsopts + Ghc-options: -Wall -rtsopts Hs-Source-Dirs: tools Build-depends: aeson, base, @@ -186,7 +184,7 @@ Executable gt-scan Main-is: gt-scan.hs - Ghc-options: -Wall -auto-all + Ghc-options: -Wall Hs-Source-Dirs: tools Build-depends: aeson, base, @@ -205,10 +203,10 @@ Executable afroengineer Main-is: afroengineer.hs - Ghc-options: -Wall -auto-all + Ghc-options: -Wall Hs-source-dirs: tools - -- Ghc-options: -Wall -auto-all -threaded -rtsopts -with-rtsopts=-N - Ghc-options: -Wall -auto-all -rtsopts + -- Ghc-options: -Wall -threaded -rtsopts -with-rtsopts=-N + Ghc-options: -Wall -rtsopts Other-Modules: Align, SimpleSeed Build-Depends: base, biohazard, @@ -221,10 +219,10 @@ Executable bam-fixpair Main-is: bam-fixpair.hs - Ghc-options: -Wall -auto-all + Ghc-options: -Wall Hs-Source-Dirs: tools - -- Ghc-options: -Wall -auto-all -threaded -rtsopts -with-rtsopts=-N - Ghc-options: -Wall -auto-all -rtsopts + -- Ghc-options: -Wall -threaded -rtsopts -with-rtsopts=-N + Ghc-options: -Wall -rtsopts Build-depends: base, binary, biohazard, @@ -235,10 +233,10 @@ Executable bam-meld Main-is: bam-meld.hs - Ghc-options: -Wall -auto-all + Ghc-options: -Wall Hs-Source-Dirs: tools - -- Ghc-options: -Wall -auto-all -threaded -rtsopts -with-rtsopts=-N - Ghc-options: -Wall -auto-all -rtsopts + -- Ghc-options: -Wall -threaded -rtsopts -with-rtsopts=-N + Ghc-options: -Wall -rtsopts Build-depends: base, biohazard, bytestring, @@ -246,10 +244,10 @@ Executable bam-resample Main-is: bam-resample.hs - Ghc-options: -Wall -auto-all + Ghc-options: -Wall Hs-Source-Dirs: tools - -- Ghc-options: -Wall -auto-all -threaded -rtsopts -with-rtsopts=-N - Ghc-options: -Wall -auto-all -rtsopts + -- Ghc-options: -Wall -threaded -rtsopts -with-rtsopts=-N + Ghc-options: -Wall -rtsopts Build-depends: base, biohazard, bytestring, @@ -257,10 +255,10 @@ Executable bam-rewrap Main-is: bam-rewrap.hs - Ghc-options: -Wall -auto-all + Ghc-options: -Wall Hs-Source-Dirs: tools - -- Ghc-options: -Wall -auto-all -threaded -rtsopts -with-rtsopts=-N - Ghc-options: -Wall -auto-all -rtsopts + -- Ghc-options: -Wall -threaded -rtsopts -with-rtsopts=-N + Ghc-options: -Wall -rtsopts Build-depends: base, biohazard, bytestring, @@ -268,10 +266,10 @@ Executable bam-rmdup Main-is: bam-rmdup.hs - Ghc-options: -Wall -auto-all + Ghc-options: -Wall Hs-Source-Dirs: tools - -- Ghc-options: -Wall -auto-all -threaded -rtsopts -with-rtsopts=-N - Ghc-options: -Wall -auto-all -rtsopts + -- Ghc-options: -Wall -threaded -rtsopts -with-rtsopts=-N + Ghc-options: -Wall -rtsopts Build-depends: base, biohazard, bytestring, @@ -283,20 +281,20 @@ Executable bam-trim Main-is: bam-trim.hs - Ghc-options: -Wall -auto-all + Ghc-options: -Wall Hs-Source-Dirs: tools - -- Ghc-options: -Wall -auto-all -threaded -rtsopts -with-rtsopts=-N - Ghc-options: -Wall -auto-all -rtsopts + -- Ghc-options: -Wall -threaded -rtsopts -with-rtsopts=-N + Ghc-options: -Wall -rtsopts Build-depends: base, biohazard, bytestring Executable fastq2bam Main-is: fastq2bam.hs - Ghc-options: -Wall -auto-all + Ghc-options: -Wall Hs-Source-Dirs: tools - -- Ghc-options: -Wall -auto-all -threaded -rtsopts -with-rtsopts=-N - Ghc-options: -Wall -auto-all -rtsopts + -- Ghc-options: -Wall -threaded -rtsopts -with-rtsopts=-N + Ghc-options: -Wall -rtsopts Build-depends: base, biohazard, bytestring, @@ -309,8 +307,8 @@ Hs-Source-Dirs: tools C-sources: src/cbits/jive.c CC-options: -std=c99 -ffast-math - -- Ghc-options: -Wall -auto-all -threaded -rtsopts -with-rtsopts=-N - Ghc-options: -Wall -auto-all -rtsopts + -- Ghc-options: -Wall -threaded -rtsopts -with-rtsopts=-N + Ghc-options: -Wall -rtsopts Other-modules: Index Build-depends: aeson, base, @@ -329,10 +327,10 @@ Executable mt-anno Main-is: mt-anno.hs - Ghc-options: -Wall -auto-all + Ghc-options: -Wall Hs-Source-Dirs: tools - -- Ghc-options: -Wall -auto-all -threaded -rtsopts -with-rtsopts=-N - Ghc-options: -Wall -auto-all -rtsopts + -- Ghc-options: -Wall -threaded -rtsopts -with-rtsopts=-N + Ghc-options: -Wall -rtsopts Other-Modules: Anno, Seqs, Xlate Build-Depends: base, bytestring, @@ -341,10 +339,10 @@ Executable mt-ccheck Main-is: mt-ccheck.hs - Ghc-options: -Wall -auto-all + Ghc-options: -Wall Hs-Source-Dirs: tools - -- Ghc-options: -Wall -auto-all -threaded -rtsopts -with-rtsopts=-N - Ghc-options: -Wall -auto-all -rtsopts + -- Ghc-options: -Wall -threaded -rtsopts -with-rtsopts=-N + Ghc-options: -Wall -rtsopts Build-Depends: base, bytestring, biohazard, diff --git a/doc/genotyping.tex b/doc/genotyping.tex --- a/doc/genotyping.tex +++ b/doc/genotyping.tex @@ -9,6 +9,8 @@ \newcommandx{\oops}[2][1=]{\todo[nolist,noline,inline,linecolor=yellow,backgroundcolor=yellow!25,bordercolor=yellow,#1]{#2}} \newcommandx{\result}[2][1=]{\todo[nolist,noline,inline,linecolor=green,backgroundcolor=green!25,bordercolor=green,#1]{#2}} +\newcommand{\argmax}{\operatornamewithlimits{arg\,max}} + \title{Deamination Aware Genotype Calling} \author{Udo Stenzel} @@ -117,9 +119,9 @@ single stranded or double stranded model being correct. (For brevity, all parameters not mentioned here are assumed to be zero.) \begin{align*} -\hat{\sigma_s}, \hat{\delta_s}, \hat{\lambda_s}, \hat{\kappa_s} &:= \arg \max_{\sigma_s, \delta_s, \lambda_s, \kappa_s} +\hat{\sigma_s}, \hat{\delta_s}, \hat{\lambda_s}, \hat{\kappa_s} &:= \argmax_{\sigma_s, \delta_s, \lambda_s, \kappa_s} P(D | \sigma_s, \delta_s, \lambda_s, \kappa_s ) \\ -\hat{\sigma_d}, \hat{\delta_d}, \hat{\lambda_d} &:= \arg \max_{\sigma_d, \delta_d, \lambda_d} +\hat{\sigma_d}, \hat{\delta_d}, \hat{\lambda_d} &:= \argmax_{\sigma_d, \delta_d, \lambda_d} P(D | \sigma_d, \delta_d, \lambda_d ) \\ P_s(D) &:= \frac{ P(D | \hat{\sigma_s}, \hat{\delta_s}, \hat{\lambda_s}, \hat{\kappa_s}) }{ @@ -528,30 +530,49 @@ When co-calling individuals from multiple populations, the correct prior for the genotypes would be based on a covariance matrix\footnote{We restrict to site-by-site analysis and a couple more simplifying -assumptions.}. Estimating -that matrix allows Treemix, Patterson's~D and Pruefer's Divergence, -possibly more. % And the replacement for TreeMix is "Blandskog" (== - % mixed forest) +assumptions.}. Estimating that matrix allows Treemix, Patterson's~D, +Pr\:ufer's Divergence, and PCA; possibly more. We also get a clean +method for imputation for free. +% And the replacement for TreeMix is "Blandskog" (== mixed forest) Conceptually, it's easy: the covariance matrix serves as prior for the -allele frequencies in multiple populations, the allele frequency -(together with a small term for new mutations) serves as prior for the -genotypes\todo{Equation!}. Maximizing the covariance matrix is -straight forward, but it would require integrating over the space of -possible combinations of allele frequencies, which sounds impractical, -and becomes exponentially more impractical the more samples are considered. Even if -symbolic integration was possible (I don't think it is), the resulting -term would grow with the number of genotype \emph{combinations}, which -is still exponential in the number of samples. +allele frequencies in multiple populations, the allele frequency serves +as prior for the genotype. The goal is to maximize the likelihood with +respect to the covariance matrix. The likelihood function looks more or +less like this: -Instead, we take inspiration from SeqEm\cite{seqem}. We treat the -allele frequencies $\mathcal{f}$ as hidden parameters, marginalize over -the genotypes, and maximize the joint probability $P(\mathcal{D}, -\mathcal{f} | \Sigma)$ with respect to the variance-covariance matrix -$\Sigma$ using an EM algorithm\footnote{Effectively, we estimate the -allele frequency at every position for every sample. Literally, this is of course -impossible, but in aggregate makes sense for populations.}. +\begin{equation*} +L = \frac{1}{\sqrt{2^k \pi^k |\Sigma|}} \int_{0^k}^{1^k} f^T G f + \exp \left( -\frac{1}{2}(f-\mu)^T \Sigma^{-1} (f-\mu) \right) df +\end{equation*} +This appears to require integration over the space of possible +combinations of allele frequencies, which sounds impractical. If this +has to be done numerically (Monte-Carlo integration seems at least +feasible), it will be slow and horrible. It's not obvious how to do it +symbolically (it would only make sense if the genotypes can be +separated, otherwise marginalization over all \emph{combinations} of +genotypes remains), but even then, a few problems remain: + +\begin{itemize} +\item The formula lacks a normalization factor arising from the +integration bounds. +\item PopGen models based on allele frequencies fail catastrophically if +the frequencies approach 1 or 0. It seems disingenious to integrate out +to these boundaries. +\item How to deal with $\mu$ is not obvious. +\end{itemize} + +Instead, we sidestep these problems by taking inspiration from +SeqEm\cite{seqem}. We treat the allele frequencies $\bar{f}$ as +hidden parameters, marginalize over the genotypes(!), and maximize the +joint probability $P(\mathcal{D}, \bar{f} | \Sigma, \mu)$ with +respect to the variance-covariance matrix $\Sigma$ and the mean +distances $\mu$ using an EM algorithm\footnote{Effectively, we estimate +the allele frequency at every position for every sample. Literally, +this is of course impossible, but in aggregate makes sense for +populations.}. + The expectation step, which is finding estimates for the allele frequencies, is much easier, as it requires only summation over the possible genotypes of \emph{one} individual at a time and optimizes @@ -561,11 +582,42 @@ matrix from allele frequencies, and this is again easy. The obvious moment based estimator should work. +Maximization (actually done by the moment based estimator for sample +mean and sample covariances): + +\begin{equation*} +\hat{\mu}, \hat{\Sigma} := \argmax_{\mu,\Sigma} + P(\bar{f}_1,\ldots,\bar{f}_k| \mu,\Sigma) +\end{equation*} + +Expectation (for each site): +\begin{align*} +\hat{f_1} \lell \hat{f_k} &= \argmax_{f_1,\ldots,f_k} + P(\mathcal{D} | f_1,\ldots,f_k) P(f_1,\ldots,f_k | \mu, \Sigma) \\ + &= \argmax_{f_1,\ldots,f_k} + \prod_{i=1}^k \sum_{G_i} P(\mathcal{D} | G_i) P( G_i | f_i ) + e^{ -\frac{1}{2} ((f_1,\ldots,f_k) - \mu)^T \Sigma^{-1} ((f_1,\ldots,f_k) - \mu) } +\end{align*} + +Here, $P(\mathcal{D}|G_i)$ is a genotype likelihood, acting as a +constant, and $P( G_i | f_i)$ follows by assuming +Hardy-Weinberg-Equilibrium\footnote{Even if HWE is violated in practice, +it will be absorbed by the estimation of more drift along the affected +lineage.}: + +\begin{equation*} +P(G=RR) = (1-f)^2 \quad P(G=RA) = 2f(1-f) \quad P(AA) = f^2 +\end{equation*} + +This \emph{looks} as if it can be optimized analytically. The summation +over the genotypes requires no nested sums, so it \emph{looks} tractable. + \idea{The practical implementation should probably not store the aforementioned 600GB of likelihoods, but only 6GB or thereabout of allele frequency data per sample. The likelihoods can be generated from the smaller BAM files on the fly. That probably means 'pileup' needs to -get a lot faster.} +get a lot faster. On second thought, maybe the frequencies don't need +to be stored at all.} \idea{Dealing with contamination becomes obvious in this framework. For a contaminated sample, every position has two allele frequencies (or two @@ -584,11 +636,11 @@ its empirical standard deviation $\sqrt{f_j (1 - f_j)}$. For multiallelic sites, one frequency is used per variant. -Since $m < n$, matrix $X = \frac{1}{n} M M^T$ is small and \sc{Lapack} -can trivially perform PCA on it. It just so happens that $X$ is the -same covariance matrix we estimated above. (Before David and Nick got -involved, everybody PCA's $\frac{1}{m} M^T M$ instead, which is much -bigger. I don't know why.) +Since $m < n$, matrix $X = \frac{1}{n} M M^T$ is small and +\textsc{Lapack} can trivially perform PCA on it. It just so happens +that $X$ is the same covariance matrix we estimated above. (Before +David and Nick got involved, everybody PCA'd $\frac{1}{m} M^T M$ +instead, which is much bigger. I don't know why.) diff --git a/src/Bio/Adna.hs b/src/Bio/Adna.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Adna.hs @@ -0,0 +1,474 @@ +{-# LANGUAGE BangPatterns, RecordWildCards, FlexibleContexts #-} +module Bio.Adna ( + DmgStats(..), + damagePatternsIter, + damagePatternsIterMD, + damagePatternsIter2Bit, + + DamageParameters(..), + DamageModel, + bang, + Alignment(..), + FragType(..), + To(..), + NPair, + + noDamage, + univDamage, + memoDamageModel, + Mat44, + Mat44D + ) where + +import Bio.Bam +import Bio.Base +import Bio.TwoBit +import Control.Arrow ( (***) ) +import Control.Monad +import Data.Bits ( xor ) +import Data.Vec ( Mat44, Mat44D, identity, getElem, Vec4, (:.)((:.)) ) + +import qualified Data.Vector.Generic as G +import qualified Data.Vector as V +import qualified Data.Vector.Storable as VS +import qualified Data.Vector.Unboxed as U +import qualified Data.Vector.Unboxed.Mutable as UM + +-- ^ Things specific to ancient DNA, e.g. damage models. +-- +-- For aDNA, we need a substitution probability. We have three options: +-- use an empirically determined PSSM, use an arithmetically defined +-- PSSM based on the /Johnson/ model, use a context sensitive PSSM based +-- on the /Johnson/ model and an alignment. Using /Dindel/, actual +-- substitutions relative to a called haplotype would be taken into +-- account. Since we're not going to do that, taking alignments into +-- account is difficult, somewhat approximate, and therefore not worth +-- the hassle. +-- +-- We represent substitution matrices by the type 'Mat44D'. Internally, +-- this is a vector of packed vectors. Conveniently, each of the packed +-- vectors represents all transition /into/ the given nucleotide. + + +-- | A 'DamageModel' is a function that gives substitution matrices for +-- each position in a read. The 'DamageModel' can depend on the length +-- of the read and whether its alignment is reversed. In practice, we +-- should probably memoize precomputed damage models somehow. + +type DamageModel a = Bool -> Int -> V.Vector (Mat44 a) +data To = Nucleotide :-> Nucleotide + +infix 9 :-> +infix 8 `bang` + +-- | Convenience function to access a substitution matrix that has a +-- mnemonic reading. +{-# INLINE bang #-} +bang :: Mat44D -> To -> Double +bang m (N x :-> N y) = getElem (fromIntegral x) $ getElem (fromIntegral y) m + +-- | 'DamageModel' for undamaged DNA. The likelihoods follow directly +-- from the quality score. This needs elaboration to see what to do +-- with amibiguity codes (even though those haven't actually been +-- observed in the wild). + +{-# SPECIALIZE noDamage :: DamageModel Double #-} +noDamage :: Num a => DamageModel a +noDamage _ l = V.replicate l identity + + +-- | Parameters for the universal damage model. +-- +-- We assume the correct model is either no damage, or single strand +-- damage, or double strand damage. Each of them comes with a +-- probability. It turns out that blending them into one is simply +-- accomplished by multiplying these probabilities onto the deamination +-- probabilities. +-- +-- For single stranded library prep, only one kind of damage occurs (C +-- to T), it occurs at low frequency ('ssd_delta') everywhere, at high +-- frequency ('ssd_sigma') in single stranded parts, and the overhang +-- length is distributed exponentially with parameter 'ssd_lambda' at +-- the 5' end and 'ssd_kappa' at the 3' end. (Without UDG treatment, +-- those will be equal. With UDG, those are much smaller and in fact +-- don't literally represent overhangs.) +-- +-- For double stranded library prep, we get C->T damage at the 5' end +-- and G->A at the 3' end with rate 'dsd_sigma' and both in the interior +-- with rate 'dsd_delta'. Everything is symmetric, and therefore the +-- orientation of the aligned read doesn't matter either. Both +-- overhangs follow a distribution with parameter 'dsd_lambda'. + +data DamageParameters float = DP { ssd_sigma :: !float -- deamination rate in ss DNA, SS model + , ssd_delta :: !float -- deamination rate in ds DNA, SS model + , ssd_lambda :: !float -- param for geom. distribution, 5' end, SS model + , ssd_kappa :: !float -- param for geom. distribution, 3' end, SS model + , dsd_sigma :: !float -- deamination rate in ss DNA, DS model + , dsd_delta :: !float -- deamination rate in ds DNA, DS model + , dsd_lambda :: !float } -- param for geom. distribution, DS model + deriving (Read, Show) + +-- | Generic substitution matrix, has C->T and G->A deamination as +-- parameters. Setting 'p' or 'q' to 0 as appropriate makes this apply +-- to the single stranded or undamaged case. + +{-# INLINE genSubstMat #-} +genSubstMat :: Fractional a => a -> a -> Mat44 a +genSubstMat p q = vec4 ( vec4 1 0 q 0 ) + ( vec4 0 (1-p) 0 0 ) + ( vec4 0 0 (1-q) 0 ) + ( vec4 0 p 0 1 ) + where + vec4 :: a -> a -> a -> a -> Vec4 a + vec4 a b c d = a :. b :. c :. d :. () + +memoDamageModel :: DamageModel a -> DamageModel a +memoDamageModel f = \r l -> if l > 512 || l < 0 then f r l + else if r then V.unsafeIndex rev l + else V.unsafeIndex fwd l + where + rev = V.generate 512 $ f True + fwd = V.generate 512 $ f False + +{-# SPECIALIZE univDamage :: DamageParameters Double -> DamageModel Double #-} +univDamage :: Fractional a => DamageParameters a -> DamageModel a +univDamage DP{..} r l = V.generate l mat + where + mat i = genSubstMat (p1+p2) (q1+q2) + where + (p1, q1) = if r then let lam5 = ssd_lambda ^ (l-i) + lam3 = ssd_kappa ^ (1+i) + lam = lam3 + lam5 - lam3 * lam5 + p = ssd_sigma * lam + ssd_delta * (1-lam) + in (0,p) + else let lam5 = ssd_lambda ^ (1+i) + lam3 = ssd_kappa ^ (l-i) + lam = lam3 + lam5 - lam3 * lam5 + p = ssd_sigma * lam + ssd_delta * (1-lam) + in (p,0) + + p2 = dsd_sigma * lam5_ds + dsd_delta * (1-lam5_ds) + q2 = dsd_sigma * lam3_ds + dsd_delta * (1-lam3_ds) + lam5_ds = dsd_lambda ^ (1+i) + lam3_ds = dsd_lambda ^ (l-i) + + + +-- | Collected \"traditional\" statistics: +-- +-- * Base composition near 5' end and near 3' end. Each consists of +-- five vectors of counts of A,C,G,T, and everything else. 'basecompo5' +-- begins with 'context' bases to the left of the 5' end, 'basecompo3' +-- ends with 'context' bases to the right of the 3' end. +-- +-- * Substitutions. Counted from the reconstructed alignment, once +-- around the 5' end and once around the 3' end. For a total of 2*4*4 +-- different substitutions. +-- +-- * Substitutions at CpG motivs. Also counted from the reconstructed +-- alignment, and a CpG site is simply the sequence CG in the reference. +-- Gaps may interfere with that, but that might actually be desirable. +-- +-- * Conditional substitutions. Need an exact definition, and then a +-- couple additional tables. +-- +-- XXX This got kind of ugly. We'll see where this goes... + +data DmgStats a = DmgStats { + basecompo5 :: [( Maybe Nucleotide, U.Vector Int )], + basecompo3 :: [( Maybe Nucleotide, U.Vector Int )], + substs5 :: [( To, U.Vector Int )], + substs3 :: [( To, U.Vector Int )], + substs5d5 :: [( To, U.Vector Int )], + substs3d5 :: [( To, U.Vector Int )], + substs5d3 :: [( To, U.Vector Int )], + substs3d3 :: [( To, U.Vector Int )], + substs5dd :: [( To, U.Vector Int )], + substs3dd :: [( To, U.Vector Int )], + substs5cpg :: [( To, U.Vector Int )], + substs3cpg :: [( To, U.Vector Int )], + stats_more :: a } + +data FragType = Complete | Leading | Trailing deriving (Show, Eq) +type NPair = ( Nucleotides, Nucleotides ) + +-- Alignment record, might have been gotten from practically anywhere +-- with varying completeness. We record anything we can get, most is +-- optional. Reference sequence is filled with Ns if missing. +data Alignment = ALN + { a_sequence :: !(U.Vector NPair) -- the alignment proper + , a_fragment_type :: !FragType } -- was the adapter trimmed? + + +-- | Enumeratee (almost) that computes some statistics from plain BAM +-- (no MD field needed) and a 2bit file. The 'Alignment' is also +-- reconstructed and passed downstream. The result of any downstream +-- processing is available in the 'stats_more' field of the result. +-- +-- * Get the reference sequence including both contexts once. If this +-- includes invalid sequence (negative coordinate), pad suitably. +-- * Accumulate counts for the valid parts around 5' and 3' ends as +-- appropriate from flags and config. +-- * Combine the part that was aligned to (so no context) with the read +-- to reconstruct the alignment. +-- +-- Arguments are the table of reference names, the 2bit file with the +-- reference, the amount of context outside the alignment desired, the +-- amount of context inside desired, and whether to compensate for +-- leeHom breakage. +-- +-- For 'Complete' fragments, we cut the read in the middle, so the 5' +-- and 3' plots stay clean from each other's influence. 'Leading' and +-- 'Trailing' fragments count completely towards the appropriate end. + +damagePatternsIter2Bit :: MonadIO m + => Refs -> TwoBitFile -> Int -> Int -> Bool + -> Iteratee [Alignment] m b + -> Iteratee [BamRaw] m (DmgStats b) +damagePatternsIter2Bit refs tbf ctx rng = + damagePatternsIter ctx rng $ \br -> + let b@BamRec{..} = unpackBam br + ref_nm = sq_name $ getRef refs b_rname + ref = getFragment tbf ref_nm (b_pos - ctx) (alignedLength b_cigar + 2*ctx) + pairs = aln_from_ref (U.drop ctx ref) b_seq b_cigar + in if isUnmapped b + then Nothing + else Just (b, ref, pairs) + + +-- | Enumeratee (almost) that computes some statistics from plain BAM +-- with a valid MD field. The 'Alignment' is also reconstructed and +-- passed downstream. The result of any downstream processing is +-- available in the 'stats_more' field of the result. +-- +-- * Reconstruct the alignment from CIGAR, SEQ, and MD. +-- * Filter the alignment to get the reference sequence, accumulate it. +-- * Accumulate everything over the alignment. +-- +-- Arguments are just the amount of context inside desired. +-- Arguments are the amount of context inside desired, and whether to +-- compensate for leeHom breakage. +-- +-- For 'Complete' fragments, we cut the read in the middle, so the 5' +-- and 3' plots stay clean from each other's influence. 'Leading' and +-- 'Trailing' fragments count completely towards the appropriate end. + +damagePatternsIterMD :: MonadIO m + => Int -> Bool + -> Iteratee [Alignment] m b + -> Iteratee [BamRaw] m (DmgStats b) +damagePatternsIterMD rng = + damagePatternsIter 0 rng $ \br -> do + let b@BamRec{..} = unpackBam br + guard (not $ isUnmapped b) + md <- getMd b + let pairs = aln_from_md b_seq b_cigar md + ref = U.map fromN $ U.filter ((/=) gap . fst) pairs + return (b, ref, pairs) + where + fromN (ns,_) | ns == nucsA = 2 + | ns == nucsC = 1 + | ns == nucsG = 3 + | ns == nucsT = 0 + | otherwise = 4 + +-- | Common logic for statistics. The function 'get_ref_and_aln' +-- reconstructs reference sequence and alignment from a Bam record. It +-- is expected to construct the alignment with respect to the forwards +-- strand of the reference; we reverse-complement it if necessary. +damagePatternsIter :: MonadIO m + => Int -> Int + -> (BamRaw -> Maybe (BamRec, U.Vector Word8, U.Vector NPair)) + -> Bool -> Iteratee [Alignment] m b + -> Iteratee [BamRaw] m (DmgStats b) +damagePatternsIter ctx rng get_ref_and_aln leeHom it = do + let maxwidth = ctx + rng + acc_bc <- liftIO $ UM.replicate (2 * 5 * maxwidth) (0::Int) + acc_st <- liftIO $ UM.replicate (2 * 4 * 4 * 4 * rng) (0::Int) + acc_cg <- liftIO $ UM.replicate (2 * 2 * 4 * rng) (0::Int) + + let do_bc br = case fmap revcom_both $ get_ref_and_aln br of + Nothing -> return [] + Just (b@BamRec{..}, ref, a_sequence) -> liftIO $ do + + let good_pairs = U.indexed a_sequence + good_pairs_rev = U.indexed $ U.reverse a_sequence + a_fragment_type | isFirstMate b && isPaired b = Leading + | isSecondMate b && isPaired b = Trailing + | leeHom = Complete + | isFirstMate b || isSecondMate b = Complete -- old style flagging + | isTrimmed b || isMerged b = Complete -- new style flagging + | otherwise = Leading + + -- basecompositon near 5' end, near 3' end + let width = ctx + min rng (alignedLength b_cigar `div` 2) + mapM_ (\i -> bump (fromIntegral (ref U.! i ) * maxwidth + i) acc_bc) [0 .. width-1] + mapM_ (\i -> bump (fromIntegral (ref U.! (i + U.length ref) +6) * maxwidth + i) acc_bc) [-width .. -1] + + -- For substitutions, decide what damage class we're in: + -- 0 - no damage, 1 - damaged 5' end, 2 - damaged 3' end, 3 - both + let dmgbase = 2*4*4*rng * + ( (if U.null a_sequence || U.head a_sequence /= (nucsC,nucsT) then 1 else 0) + + (if U.null a_sequence || U.last a_sequence /= (nucsC,nucsT) then 2 else 0) ) + + -- substitutions near 5' end + let len_at_5 = case a_fragment_type of Leading -> min rng (G.length b_seq) + Complete -> min rng (G.length b_seq `div` 2) + Trailing -> 0 + U.forM_ (U.take len_at_5 good_pairs) $ + \(i,uv) -> withPair uv $ \j -> bump (j * rng + i + dmgbase) acc_st + + -- substitutions at CpG sites near 5' end + U.zipWithM_ + (\(i,(u,v)) (_,(w,z)) -> + when (u == nucsC && w == nucsG) $ do + withNs v $ \y -> bump ( y * rng + i ) acc_cg + withNs z $ \y -> bump ((y+4) * rng + i+1) acc_cg) + (U.take len_at_5 good_pairs) (U.drop 1 good_pairs) + + -- substitutions near 3' end + let len_at_3 = case a_fragment_type of Leading -> 0 + Complete -> min rng (G.length b_seq `div` 2) + Trailing -> min rng (G.length b_seq) + U.forM_ (U.take len_at_3 good_pairs_rev) $ + \(i,uv) -> withPair uv $ \j -> bump ((17+j) * rng -i -1 + dmgbase) acc_st + + -- substitutions at CpG sites near 3' end + U.zipWithM_ + (\(_,(u,v)) (i,(w,z)) -> + when (u == nucsC && w == nucsG) $ do + withNs v $ \y -> bump ((y+ 9) * rng - i-2) acc_cg + withNs z $ \y -> bump ((y+13) * rng - i-1) acc_cg) + (U.drop 1 good_pairs_rev) (U.take len_at_3 good_pairs_rev) + + + return [ ALN{..} ] + where + {-# INLINE withPair #-} + withPair (Ns u, Ns v) k = case pairTab `U.unsafeIndex` fromIntegral (16*u+v) of + j -> if j >= 0 then k j else return () + + !pairTab = U.replicate 256 (-1) U.// + [ (fromIntegral $ 16*u+v, x*4+y) | (Ns u,x) <- zip [nucsA, nucsC, nucsG, nucsT] [0,1,2,3] + , (Ns v,y) <- zip [nucsA, nucsC, nucsG, nucsT] [0,1,2,3] ] + + {-# INLINE bump #-} + bump i v = UM.unsafeRead v i >>= UM.unsafeWrite v i . succ + + {-# INLINE withNs #-} + withNs ns k | ns == nucsA = k 0 + | ns == nucsC = k 1 + | ns == nucsG = k 2 + | ns == nucsT = k 3 + | otherwise = return () + + + it' <- concatMapStreamM do_bc it + + let nsubsts = 2*4*4*rng + mk_substs off = sequence [ (,) (n1 :-> n2) <$> U.unsafeFreeze (UM.slice ((4*i+j)*rng + off*nsubsts) rng acc_st) + | (i,n1) <- zip [0..] [nucA..nucT] + , (j,n2) <- zip [0..] [nucA..nucT] ] + + accs <- liftIO $ DmgStats <$> sequence [ (,) nuc <$> U.unsafeFreeze (UM.slice (i*maxwidth) maxwidth acc_bc) + | (i,nuc) <- zip [2,1,3,0,4] [Just nucA, Just nucC, Just nucG, Just nucT, Nothing] ] + <*> sequence [ (,) nuc <$> U.unsafeFreeze (UM.slice (i*maxwidth) maxwidth acc_bc) + | (i,nuc) <- zip [7,6,8,5,9] [Just nucA, Just nucC, Just nucG, Just nucT, Nothing] ] + + <*> mk_substs 0 + <*> mk_substs 1 + <*> mk_substs 2 + <*> mk_substs 3 + <*> mk_substs 4 + <*> mk_substs 5 + <*> mk_substs 6 + <*> mk_substs 7 + + <*> sequence [ (,) (n1 :-> n2) <$> U.unsafeFreeze (UM.slice ((i+j)*rng) rng acc_cg) + | (i,n1) <- [(0,nucC), (4,nucG)] + , (j,n2) <- zip [0..] [nucA..nucT] ] + + <*> sequence [ (,) (n1 :-> n2) <$> U.unsafeFreeze (UM.slice ((i+j)*rng) rng acc_cg) + | (i,n2) <- [(8,nucC), (12,nucG)] + , (j,n1) <- zip [0..] [nucA..nucT] ] + + accs' <- accs `liftM` lift (run it') + let vsum = foldl1 (zipWith s1) . map ($ accs') + s1 (x :-> y, u) (z :-> w, v) | x == z && y == w = (x :-> y, U.zipWith (+) u v) + | otherwise = error $ "Mismatch in zip. This is a bug." + + return $ accs' { substs5 = vsum [ substs5, substs5d5, substs5d3, substs5dd ] + , substs3 = vsum [ substs3, substs3d5, substs3d3, substs3dd ] + , substs5d5 = vsum [ substs5d5, substs5dd ] + , substs3d5 = vsum [ substs3d5, substs3dd ] + , substs5d3 = vsum [ substs5d3, substs5dd ] + , substs3d3 = vsum [ substs3d3, substs3dd ] } + + +revcom_both :: ( BamRec, U.Vector Word8, U.Vector (Nucleotides, Nucleotides) ) + -> ( BamRec, U.Vector Word8, U.Vector (Nucleotides, Nucleotides) ) +revcom_both (b, ref, pairs) + | isReversed b = ( b, revcom_ref ref, revcom_pairs pairs ) + | otherwise = ( b, ref, pairs ) + where + revcom_ref = U.reverse . U.map (xor 2) + revcom_pairs = U.reverse . U.map (compls *** compls) + + +-- | Reconstructs the alignment from reference, query, and cigar. Only +-- positions where the query is not gapped are produced. +aln_from_ref :: U.Vector Word8 -> Vector_Nucs_half Nucleotides -> VS.Vector Cigar -> U.Vector NPair +aln_from_ref ref0 qry0 cig0 = U.fromList $ step ref0 qry0 cig0 + where + step ref qry cig1 + | U.null ref || G.null qry || G.null cig1 = [] + | otherwise = case G.unsafeHead cig1 of { op :* n -> + case G.unsafeTail cig1 of { cig -> + case op of { + + Mat -> zipWith (\r q -> (nn r,q)) (G.toList (G.take n ref)) + (G.toList (G.take n qry)) ++ step (G.drop n ref) (G.drop n qry) cig ; + Del -> step (G.drop n ref) qry cig ; + Ins -> map (\q -> ( gap, q )) (G.toList (G.take n qry)) ++ step ref (G.drop n qry) cig ; + SMa -> map (\q -> ( gap, q )) (G.toList (G.take n qry)) ++ step ref (G.drop n qry) cig ; + HMa -> replicate n (gap, nucsN) ++ step ref qry cig ; + Nop -> step ref qry cig ; + Pad -> step ref qry cig }}} + + nn 0 = nucsT + nn 1 = nucsC + nn 2 = nucsA + nn 3 = nucsG + nn _ = nucsN + + +-- | Reconstructs the alignment from query, cigar, and md. Only +-- positions where the query is not gapped are produced. +aln_from_md :: Vector_Nucs_half Nucleotides -> VS.Vector Cigar -> [MdOp] -> U.Vector NPair +aln_from_md qry0 cig0 md0 = U.fromList $ step qry0 cig0 md0 + where + step qry cig1 md + | G.null qry || G.null cig1 || null md = [] + | otherwise = case G.unsafeHead cig1 of op :* n -> step' qry op n (G.unsafeTail cig1) md + + step' qry _ 0 cig md = step qry cig md + step' qry op n cig (MdNum 0 : md) = step' qry op n cig md + step' qry op n cig (MdDel [] : md) = step' qry op n cig md + + step' qry Mat n cig (MdNum m : md) + | n < m = map (\q -> (q,q)) (G.toList (G.take n qry)) ++ step (G.drop n qry) cig (MdNum (m-n) : md) + | n > m = map (\q -> (q,q)) (G.toList (G.take m qry)) ++ step' (G.drop m qry) Mat (n-m) cig md + | n == m = map (\q -> (q,q)) (G.toList (G.take n qry)) ++ step (G.drop n qry) cig md + step' qry Mat n cig (MdRep c : md) = ( c, G.head qry ) : step' (G.tail qry) Mat (n-1) cig md + step' _ Mat _ _ _ = [] + + step' qry Del n cig (MdDel (_:ss) : md) = step' qry Del (n-1) cig (MdDel ss : md) + step' _ Del _ _ _ = [] + + step' qry Ins n cig md = map ((,) gap) (G.toList (G.take n qry)) ++ step (G.drop n qry) cig md + step' qry SMa n cig md = map ((,) gap) (G.toList (G.take n qry)) ++ step (G.drop n qry) cig md + step' qry HMa n cig md = replicate n (gap, nucsN) ++ step qry cig md + step' qry Nop _ cig md = step qry cig md + step' qry Pad _ cig md = step qry cig md + + diff --git a/src/Bio/Genocall.hs b/src/Bio/Genocall.hs --- a/src/Bio/Genocall.hs +++ b/src/Bio/Genocall.hs @@ -1,9 +1,9 @@ {-# LANGUAGE BangPatterns #-} module Bio.Genocall where +import Bio.Adna import Bio.Bam.Pileup import Bio.Base -import Bio.Genocall.Adna import Control.Applicative import Data.Foldable hiding ( sum, product ) import Data.List ( inits, tails, sortBy ) @@ -44,7 +44,7 @@ | (_q,(d,i)) <- vars , not (null d) || not (null i) ] - match = zipWith $ \(DB b q _ m) n -> let p = m ! n :-> b + match = zipWith $ \(DB b q _ m) n -> let p = m `bang` n :-> b p' = fromQual q in toProb $ p + p' - p * p' @@ -63,12 +63,12 @@ simple_snp_call from_qual ploidy vars = snp_gls (simple_call ploidy mkpls vars) ref where ref = case vars of (_, DB _ _ r _) : _ -> r ; _ -> nucsN - mkpls (q, DB b qq _ m) = [ toProb $ x + pe*(s-x) | n <- [0..3], let x = m ! N n :-> b ] + mkpls (q, DB b qq _ m) = [ toProb $ x + pe*(s-x) | n <- [0..3], let x = m `bang` N n :-> b ] where !p1 = from_qual q !p2 = from_qual qq !pe = p1 + p2 - p1*p2 - !s = sum [ m ! N n :-> b | n <- [0..3] ] / 4 + !s = sum [ m `bang` N n :-> b | n <- [0..3] ] / 4 -- | Compute @GL@ values for the simple case. The simple case is where -- we sample 'ploidy' alleles with equal probability and assume that @@ -162,7 +162,7 @@ let -- P(X|Q,H), a vector of four (x is fixed, h is not) -- this is the simple form where we set all w to 1/4 - p_x__q_h_ = Vec.map (\h -> 0.25 * fromQualRaised (theta ** (k ! h :-> db_call x)) (db_qual x)) everynuc + p_x__q_h_ = Vec.map (\h -> 0.25 * fromQualRaised (theta ** (k `bang` h :-> db_call x)) (db_qual x)) everynuc -- eh, this is cumbersome... what was I thinking?! p_x__q_h = Vec.zipWith (\p h -> if db_call x == h then 1 + p - Vec.sum p_x__q_h_ else p) p_x__q_h_ everynuc @@ -176,7 +176,7 @@ k' = Vec.setElem (fromIntegral . unN $ db_call x) kk k acc' = acc * toProb p_x__q - meh = Vec.map (\h -> k ! h :-> db_call x) everynuc -- XXX + meh = Vec.map (\h -> k `bang` h :-> db_call x) everynuc -- XXX in {- trace (unlines ["gt " ++ show gt ,"p(x|q,h) " ++ show p_x__q_h ,"dg " ++ show dg ++ ", call = " ++ show (db_call x) diff --git a/src/Bio/Genocall/Adna.hs b/src/Bio/Genocall/Adna.hs deleted file mode 100644 --- a/src/Bio/Genocall/Adna.hs +++ /dev/null @@ -1,126 +0,0 @@ -{-# LANGUAGE BangPatterns, RecordWildCards #-} -module Bio.Genocall.Adna where - -import Bio.Base -import Data.Vec -import qualified Data.Vector as V - --- ^ Things specific to ancient DNA, e.g. damage models. --- --- For aDNA, we need a substitution probability. We have three options: --- use an empirically determined PSSM, use an arithmetically defined --- PSSM based on the /Johnson/ model, use a context sensitive PSSM based --- on the /Johnson/ model and an alignment. Using /Dindel/, actual --- substitutions relative to a called haplotype would be taken into --- account. Since we're not going to do that, taking alignments into --- account is difficult, somewhat approximate, and therefore not worth --- the hassle. --- --- We represent substitution matrices by the type 'Mat44D'. Internally, --- this is a vector of packed vectors. Conveniently, each of the packed --- vectors represents all transition /into/ the given nucleotide. - - --- | A 'DamageModel' is a function that gives substitution matrices for --- each position in a read. The 'DamageModel' can depend on the length --- of the read and whether its alignment is reversed. In practice, we --- should probably memoize precomputed damage models somehow. - -type DamageModel a = Bool -> Int -> V.Vector (Mat44 a) - -data To = Nucleotide :-> Nucleotide - -infix 9 :-> -infix 8 ! - --- | Convenience function to access a substitution matrix that has a --- mnemonic reading. -{-# INLINE (!) #-} -(!) :: Mat44D -> To -> Double -(!) m (N x :-> N y) = getElem (fromIntegral x) $ getElem (fromIntegral y) m - --- | 'DamageModel' for undamaged DNA. The likelihoods follow directly --- from the quality score. This needs elaboration to see what to do --- with amibiguity codes (even though those haven't actually been --- observed in the wild). - -{-# SPECIALIZE noDamage :: DamageModel Double #-} -noDamage :: Num a => DamageModel a -noDamage _ l = V.replicate l identity - - --- | Parameters for the universal damage model. --- --- We assume the correct model is either no damage, or single strand --- damage, or double strand damage. Each of them comes with a --- probability. It turns out that blending them into one is simply --- accomplished by multiplying these probabilities onto the deamination --- probabilities. --- --- For single stranded library prep, only one kind of damage occurs (C --- to T), it occurs at low frequency ('ssd_delta') everywhere, at high --- frequency ('ssd_sigma') in single stranded parts, and the overhang --- length is distributed exponentially with parameter 'ssd_lambda' at --- the 5' end and 'ssd_kappa' at the 3' end. (Without UDG treatment, --- those will be equal. With UDG, those are much smaller and in fact --- don't literally represent overhangs.) --- --- For double stranded library prep, we get C->T damage at the 5' end --- and G->A at the 3' end with rate 'dsd_sigma' and both in the interior --- with rate 'dsd_delta'. Everything is symmetric, and therefore the --- orientation of the aligned read doesn't matter either. Both --- overhangs follow a distribution with parameter 'dsd_lambda'. - -data DamageParameters float = DP { ssd_sigma :: !float -- deamination rate in ss DNA, SS model - , ssd_delta :: !float -- deamination rate in ds DNA, SS model - , ssd_lambda :: !float -- param for geom. distribution, 5' end, SS model - , ssd_kappa :: !float -- param for geom. distribution, 3' end, SS model - , dsd_sigma :: !float -- deamination rate in ss DNA, DS model - , dsd_delta :: !float -- deamination rate in ds DNA, DS model - , dsd_lambda :: !float } -- param for geom. distribution, DS model - deriving (Read, Show) - --- | Generic substitution matrix, has C->T and G->A deamination as --- parameters. Setting 'p' or 'q' to 0 as appropriate makes this apply --- to the single stranded or undamaged case. - -{-# INLINE genSubstMat #-} -genSubstMat :: Fractional a => a -> a -> Mat44 a -genSubstMat p q = vec4 ( vec4 1 0 q 0 ) - ( vec4 0 (1-p) 0 0 ) - ( vec4 0 0 (1-q) 0 ) - ( vec4 0 p 0 1 ) - where - vec4 :: a -> a -> a -> a -> Vec4 a - vec4 a b c d = a :. b :. c :. d :. () - -memoDamageModel :: DamageModel a -> DamageModel a -memoDamageModel f = \r l -> if l > 512 || l < 0 then f r l - else if r then V.unsafeIndex rev l - else V.unsafeIndex fwd l - where - rev = V.generate 512 $ f True - fwd = V.generate 512 $ f False - -{-# SPECIALIZE univDamage :: DamageParameters Double -> DamageModel Double #-} -univDamage :: Fractional a => DamageParameters a -> DamageModel a -univDamage DP{..} r l = V.generate l mat - where - mat i = genSubstMat (p1+p2) (q1+q2) - where - (p1, q1) = if r then let lam5 = ssd_lambda ^ (l-i) - lam3 = ssd_kappa ^ (1+i) - lam = lam3 + lam5 - lam3 * lam5 - p = ssd_sigma * lam + ssd_delta * (1-lam) - in (0,p) - else let lam5 = ssd_lambda ^ (1+i) - lam3 = ssd_kappa ^ (l-i) - lam = lam3 + lam5 - lam3 * lam5 - p = ssd_sigma * lam + ssd_delta * (1-lam) - in (p,0) - - p2 = dsd_sigma * lam5_ds + dsd_delta * (1-lam5_ds) - q2 = dsd_sigma * lam3_ds + dsd_delta * (1-lam3_ds) - lam5_ds = dsd_lambda ^ (1+i) - lam3_ds = dsd_lambda ^ (l-i) - diff --git a/src/Bio/Genocall/Metadata.hs b/src/Bio/Genocall/Metadata.hs --- a/src/Bio/Genocall/Metadata.hs +++ b/src/Bio/Genocall/Metadata.hs @@ -2,19 +2,19 @@ -- | Metadata necessary for a sensible genotyping workflow. module Bio.Genocall.Metadata where -import Bio.Genocall.Adna ( DamageParameters(..) ) +import Bio.Adna ( DamageParameters(..) ) import Control.Applicative hiding ( empty ) import Control.Concurrent ( threadDelay ) import Control.Exception ( bracket, onException, handleJust ) import Control.Monad ( forM_ ) -import Data.Text ( Text, pack ) -import Data.HashMap.Strict ( HashMap ) import Data.Aeson -import Data.ByteString.Char8 ( readFile ) -import Data.ByteString.Lazy ( toChunks ) +import Data.Binary +import Data.Binary.Get ( runGetOrFail ) +import Data.Binary.Put ( runPut ) +import Data.ByteString.Lazy ( toChunks, readFile ) import Data.ByteString.Unsafe ( unsafeUseAsCStringLen ) import Data.Monoid -import Data.Vector.Unboxed ( Vector ) +import Data.Text ( Text, pack ) import Foreign.Ptr ( castPtr ) import GHC.IO.Exception ( IOErrorType(..) ) import Prelude hiding ( writeFile, readFile ) @@ -23,13 +23,17 @@ import System.Posix.IO import qualified Data.HashMap.Strict as M +import qualified Data.Vector.Unboxed as U +data DivTable = DivTable !Double !(U.Vector Int) + deriving Show + data Sample = Sample { sample_libraries :: [Library], - sample_avro_files :: HashMap Text Text, -- ^ maps a region to the av file - sample_bcf_files :: HashMap Text Text, -- ^ maps a region to the bcf file - sample_div_tables :: HashMap Text (Double, Vector Int), -- ^ maps a region to the table needed for div. estimation - sample_divergences :: HashMap Text DivEst + sample_avro_files :: M.HashMap Text Text, -- ^ maps a region to the av file + sample_bcf_files :: M.HashMap Text Text, -- ^ maps a region to the bcf file + sample_div_tables :: M.HashMap Text DivTable, -- ^ maps a region to the table needed for div. estimation + sample_divergences :: M.HashMap Text DivEst } deriving Show data Library = Library { @@ -47,12 +51,16 @@ } deriving Show -type Metadata = HashMap Text Sample +type Metadata = M.HashMap Text Sample instance ToJSON DivEst where toJSON DivEst{..} = object $ [ "estimate" .= point_est , "confidence-region" .= conf_region ] +instance Binary DivEst where + put DivEst{..} = put point_est >> put conf_region + get = DivEst <$> get <*> get + instance FromJSON DivEst where parseJSON (Object o) = DivEst <$> o .: "estimate" <*> o .:? "confidence-region" .!= [] parseJSON (Array a) = flip DivEst [] <$> parseJSON (Array a) @@ -67,6 +75,11 @@ , "ds-delta" .= dsd_delta , "ds-lambda" .= dsd_lambda ] +instance Binary float => Binary (DamageParameters float) where + put DP{..} = put ssd_sigma >> put ssd_delta >> put ssd_lambda >> put ssd_kappa >> + put dsd_sigma >> put dsd_delta >> put dsd_lambda + get = DP <$> get <*> get <*> get <*> get <*> get <*> get <*> get + instance FromJSON float => FromJSON (DamageParameters float) where parseJSON = withObject "damage parameters" $ \o -> DP <$> o .: "ss-sigma" @@ -81,6 +94,10 @@ toJSON (Library name files dp) = object ( maybe id ((:) . ("damage" .=)) dp $ [ "name" .= name, "files" .= files ] ) +instance Binary Library where + put Library{..} = put library_name >> put library_files >> put library_damage + get = Library <$> get <*> get <*> get + instance FromJSON Library where parseJSON (String name) = return $ Library name [name <> ".bam"] Nothing parseJSON (Object o) = Library <$> o .: "name" @@ -99,6 +116,12 @@ hashToJson k vs = if M.null vs then id else (:) (k .= vs) listToJson k vs = if null vs then id else (:) (k .= vs) +instance Binary Sample where + put Sample{..} = put sample_libraries >> putObject sample_avro_files >> + putObject sample_bcf_files >> putObject sample_div_tables >> + putObject sample_divergences + get = Sample <$> get <*> getObject <*> getObject <*> getObject <*> getObject + instance FromJSON Sample where parseJSON (String s) = pure $ Sample [Library s [s <> ".bam"] Nothing] M.empty M.empty M.empty M.empty parseJSON (Array ls) = (\ll -> Sample ll M.empty M.empty M.empty M.empty) <$> parseJSON (Array ls) @@ -109,11 +132,51 @@ <*> (M.singleton "" <$> o .: "divergence" <|> o.:? "divergences" .!= M.empty) parseJSON _ = fail $ "String, Array or Object expected for Sample" +instance FromJSON DivTable where + parseJSON x = parseJSON x >>= \[a,b] -> DivTable <$> parseJSON a <*> parseJSON b +instance Binary DivTable where + put (DivTable a b) = put a >> putVector b + get = DivTable <$> get <*> getVector + +instance ToJSON DivTable where + toJSON (DivTable a b) = toJSON [toJSON a, toJSON b] + +putObject :: Binary value => M.HashMap Text value -> Put +putObject m = put (M.size m) >> M.foldrWithKey (\k v a -> put k >> put v >> a) (return ()) m + +getObject :: Binary value => Get (M.HashMap Text value) +getObject = get >>= \l -> get_map M.empty (l::Int) + where + get_map !acc 0 = return acc + get_map !acc n = get >>= \k -> get >>= \v -> get_map (M.insert k v acc) (n-1) + +-- Hm. I'm putting the vector in reverse order, because the +-- accumulation when reading reverses it. +putVector :: (U.Unbox a, Binary a) => U.Vector a -> Put +putVector v = put (U.length v) >> U.mapM_ put (U.reverse v) + +getVector :: (U.Unbox a, Binary a) => Get (U.Vector a) +getVector = get >>= \l -> U.fromListN l <$> get_list [] l + where + get_list acc 0 = return acc + get_list acc n = get >>= \ !x -> get_list (x:acc) (n-1) + + -- | Read the configuration file. Retries, because NFS tends to result -- in 'ResourceVanished' if the file is replaced while we try to read it. +readJsonMetadata :: FilePath -> IO Metadata +readJsonMetadata fn = either error return . eitherDecode =<< go (15::Int) + where + go !n = handleJust -- retry every sec for 15 seconds + (\e -> case ioeGetErrorType e of ResourceVanished | n > 0 -> Just () ; _ -> Nothing) + (\_ -> threadDelay 1000000 >> go (n-1)) + (readFile fn) + +-- | Read the configuration file. Retries, because NFS tends to result +-- in 'ResourceVanished' if the file is replaced while we try to read it. readMetadata :: FilePath -> IO Metadata -readMetadata fn = either error return . eitherDecodeStrict =<< go (15::Int) +readMetadata fn = either (error . show) (\(_,_,r) -> return r) . runGetOrFail getObject =<< go (15::Int) where go !n = handleJust -- retry every sec for 15 seconds (\e -> case ioeGetErrorType e of ResourceVanished | n > 0 -> Just () ; _ -> Nothing) @@ -141,12 +204,19 @@ (openFd fpn WriteOnly (Just 0o666) defaultFileFlags{ exclusive = True }) (closeFd) $ \fd -> (do mdata <- readMetadata fp - forM_ (toChunks . encode . toJSON $ f mdata) $ \ch -> + forM_ (toChunks . runPut . putObject $ f mdata) $ \ch -> unsafeUseAsCStringLen ch $ \(p,l) -> fdWriteBuf fd (castPtr p) (fromIntegral l) rename fpn fp) `onException` removeLink fpn +writeMetadata :: FilePath -> Metadata -> IO () +writeMetadata fp mdata = bracket + (openFd fp WriteOnly (Just 0o666) defaultFileFlags{ exclusive = True }) + (closeFd) $ \fd -> + forM_ (toChunks . runPut . putObject $ mdata) $ \ch -> + unsafeUseAsCStringLen ch $ \(p,l) -> + fdWriteBuf fd (castPtr p) (fromIntegral l) split_sam_rgns :: Metadata -> [String] -> [( String, [Maybe String] )] split_sam_rgns _meta [ ] = [] diff --git a/src/Bio/TwoBit.hs b/src/Bio/TwoBit.hs --- a/src/Bio/TwoBit.hs +++ b/src/Bio/TwoBit.hs @@ -12,6 +12,8 @@ getSubseqAscii, getSubseqMasked, getLazySubseq, + getFragment, + getFwdSubseqV, getSeqnames, lookupSequence, getSeqLength, @@ -36,6 +38,7 @@ import qualified Data.IntMap as I import qualified Data.HashMap.Lazy as M import Data.Maybe +import qualified Data.Vector.Unboxed as U import Numeric import System.IO.Posix.MMap import System.Random @@ -257,4 +260,31 @@ mask2maybe n Soft = Just n mask2maybe _ Hard = Nothing mask2maybe _ Both = Nothing + +-- | Gets a fragment from a 2bit file. The result always has the +-- desired length; if necessary, it is padded with Ns. Be careful about +-- the unconventional encoding: 0..4 == TCAGN +getFragment :: TwoBitFile -> Seqid -> Int -> Int -> U.Vector Word8 +getFragment tbf chr p l = + case lookupSequence tbf chr of + Nothing -> U.replicate l 4 + Just tbs -> getFwdSubseqV tbf tbs p l + +-- Careful about weird encoding: 0..4 == TCAGN +getFwdSubseqV :: TwoBitFile -> TwoBitSequence -> Int -> Int -> U.Vector Word8 +getFwdSubseqV TBF{..} TBS{..} start len = U.unfoldrN len step ini + where + ini = (start, takeOverlap start tbs_n_blocks) + + step (off, nbs) + | off < 0 = Just (4, (succ off, nbs)) + | off >= tbs_dna_size = Just (4, (succ off, nbs)) + | otherwise = case nbs of + [ ] -> Just (y, (succ off, [ ])) + (s,l):nbs' | off < s -> Just (y, (succ off, nbs)) + | off < s+l -> Just (4, (succ off, nbs)) + | otherwise -> Just (y, (succ off, nbs')) + where + x = B.index tbf_raw (tbs_dna_offset + off `shiftR` 2) + y = x `shiftR` (6 - 2 * (off .&. 3)) .&. 3 -- T,C,A,G diff --git a/tools/redeye-dar.hs b/tools/redeye-dar.hs --- a/tools/redeye-dar.hs +++ b/tools/redeye-dar.hs @@ -36,11 +36,11 @@ -- - needs support for multiple input files -- - needs to deal with long (unmerged) reads (by ignoring them?) +import Bio.Adna hiding ( bang ) import Bio.Bam.Header import Bio.Bam.Index import Bio.Bam.Rec import Bio.Base -import Bio.Genocall.Adna import Bio.Genocall.Metadata import Bio.Iteratee import Bio.Util.AD @@ -132,8 +132,8 @@ -- Likelihood given precomputed damage table. We compute the giant -- table ahead of time, which maps length, index and base pair to a -- likelihood. - lk tab_m _ _ (Merged b) = U.ifoldl' (\a i np -> a * tab_m `bang` index' my_bounds (U.length b, i, NP np)) 1 b - lk _ tab_f _ (Mate1st b) = U.ifoldl' (\a i np -> a * tab_f `bang` index' my_bounds (U.length b, i, NP np)) 1 b + lk tab_m _ _ (Merged b) = U.ifoldl' (\a i np -> a * tab_m `bang` index' my_bounds (U.length b, i, NP np)) 1 b + lk _ tab_f _ (Mate1st b) = U.ifoldl' (\a i np -> a * tab_f `bang` index' my_bounds (U.length b, i, NP np)) 1 b lk _ _ tab_s (Mate2nd b) = U.ifoldl' (\a i np -> a * tab_s `bang` index' my_bounds (U.length b, i, NP np)) 1 b index' bnds x | inRange bnds x = index bnds x diff --git a/tools/redeye-div.hs b/tools/redeye-div.hs --- a/tools/redeye-div.hs +++ b/tools/redeye-div.hs @@ -75,8 +75,8 @@ Just smp -> async $ do ests <- forM eff_regions >=> mapM wait $ \rgn -> async $ fmap ((,) rgn) - $ uncurry (estimateSingle conf_verbose) - $ foldl1' (\(a,u) (b,v) -> (a+b, U.zipWith (+) u v)) + $ estimateSingle conf_verbose + $ foldl1' (\(DivTable a u) (DivTable b v) -> DivTable (a+b) (U.zipWith (+) u v)) $ H.elems $ H.filterWithKey (match rgn) $ sample_div_tables smp @@ -102,8 +102,8 @@ -- XXX we should estimate an indel rate, to be appended as the fourth -- result (but that needs different tables) -estimateSingle :: Bool -> Double -> U.Vector Int -> IO DivEst -estimateSingle verbose llk_rr tab = do +estimateSingle :: Bool -> DivTable -> IO DivEst +estimateSingle verbose (DivTable llk_rr tab) = do (fit, res, stats) <- minimize quietParameters 0.0001 (llk tab) (U.fromList [0,0,0]) let xform = map (\x -> exp x / (1 + exp x)) . VS.toList diff --git a/tools/redeye-pileup.hs b/tools/redeye-pileup.hs --- a/tools/redeye-pileup.hs +++ b/tools/redeye-pileup.hs @@ -30,6 +30,7 @@ -- overhead of 150% useless likelihood values for the sex chromosomes -- and maybe estimate heterozygosity where there is none. +import Bio.Adna import Bio.Base import Bio.Bam.Header import Bio.Bam.Index @@ -37,7 +38,6 @@ import Bio.Bam.Reader import Bio.Bam.Rec import Bio.Genocall -import Bio.Genocall.Adna import Bio.Genocall.AvroFile import Bio.Genocall.Metadata import Bio.Iteratee @@ -283,11 +283,11 @@ -- estimation afterwards. Returns the product of the -- parameter-independent parts of the likehoods and the histogram -- indexed by D and H (see @genotyping.pdf@ for details). -tabulateSingle :: (Functor m, MonadIO m) => Iteratee [Calls] m (Double, U.Vector Int) +tabulateSingle :: (Functor m, MonadIO m) => Iteratee [Calls] m DivTable tabulateSingle = do tab <- liftIO $ M.replicate (12 * maxD * maxD) (0 :: Int) - (,) <$> foldStreamM (\acc -> accum tab acc . p_snp_pile) (0 :: Double) - <*> liftIO (U.unsafeFreeze tab) + DivTable <$> foldStreamM (\acc -> accum tab acc . p_snp_pile) (0 :: Double) + <*> liftIO (U.unsafeFreeze tab) where -- We need GL values for the invariant, the three homozygous variant -- and the three single-event heterozygous variant cases. The