diff --git a/Setup.hs b/Setup.hs --- a/Setup.hs +++ b/Setup.hs @@ -1,3 +1,4 @@ +import Control.Exception ( try, IOException ) import Distribution.PackageDescription ( PackageDescription(..) ) import Distribution.Simple import Distribution.Simple.InstallDirs ( docdir, mandir, CopyDest (NoCopyDest) ) @@ -6,11 +7,10 @@ import Distribution.Simple.Program.Run ( runProgramInvocation, programInvocation, progInvokeCwd ) import Distribution.Simple.Program.Types ( ConfiguredProgram, simpleProgram ) import Distribution.Simple.Setup ( copyDest, copyVerbosity, fromFlag, installVerbosity, haddockVerbosity ) -import Distribution.Simple.Utils ( installOrdinaryFile, installOrdinaryFiles, notice ) +import Distribution.Simple.Utils import Distribution.Verbosity ( Verbosity, moreVerbose ) import System.Exit ( exitSuccess ) import System.FilePath ( splitDirectories, joinPath, takeExtension, replaceExtension, () ) -import System.Directory ( getCurrentDirectory, setCurrentDirectory, createDirectoryIfMissing, doesFileExist ) main :: IO () main = do @@ -39,11 +39,10 @@ , takeExtension (last p) == ".tex" ] installOrdinaryFiles' :: Verbosity -> FilePath -> [(FilePath, FilePath)] -> IO () -installOrdinaryFiles' verb dest = mapM_ (uncurry go) +installOrdinaryFiles' verb dest = mapM_ go where - go base src = do e <- doesFileExist (base src) - if e then installOrdinaryFile verb (base src) (dest src) - else notice verb $ show (base src) ++ " was not built, can't install." + go :: (FilePath, FilePath) -> IO (Either IOException ()) + go (base,src) = try $ installOrdinaryFile verb (base src) (dest src) withLatex :: LocalBuildInfo -> (ConfiguredProgram -> IO ()) -> IO () withLatex lbi k = maybe (return ()) k $ lookupProgram (simpleProgram "pdflatex") $ withPrograms lbi @@ -51,11 +50,12 @@ runPdflatex :: PackageDescription -> LocalBuildInfo -> Verbosity -> IO () runPdflatex pkg lbi verb = withLatex lbi $ \cmd -> do - cwd <- getCurrentDirectory - createDirectoryIfMissing True (buildDir lbi "latex") + createDirectoryIfMissingVerbose verb True (buildDir lbi "latex") sequence_ [ runProgramInvocation (moreVerbose verb) $ - (programInvocation cmd [ "-interaction=nonstopmode", cwd joinPath ("doc":f) ]) - { progInvokeCwd = Just (buildDir lbi "latex") } + (programInvocation cmd [ "-interaction=nonstopmode", ddir joinPath ("doc":f) ]) + { progInvokeCwd = Just bdir } | ("doc":f@(_:_)) <- map splitDirectories $ extraSrcFiles pkg , takeExtension (last f) == ".tex" ] - + where + bdir = buildDir lbi "latex" + ddir = joinPath (map (const "..") $ splitDirectories bdir) diff --git a/biohazard.cabal b/biohazard.cabal --- a/biohazard.cabal +++ b/biohazard.cabal @@ -1,5 +1,5 @@ Name: biohazard -Version: 0.6.3 +Version: 0.6.5 Synopsis: bioinformatics support library Description: This is a collection of modules I separated from various bioinformatics tools. The hope is to make @@ -16,6 +16,8 @@ Maintainer: udo.stenzel@eva.mpg.de Copyright: (C) 2010-2015 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 @@ -37,7 +39,7 @@ source-repository this type: git location: git://github.com/udo-stenzel/biohazard.git - tag: 0.6.3 + tag: 0.6.5 Library @@ -59,40 +61,45 @@ Bio.Genocall, Bio.Genocall.Adna, Bio.Genocall.AvroFile, - Bio.Glf, + Bio.Genocall.Metadata, Bio.Iteratee, Bio.Iteratee.Bgzf, Bio.Iteratee.Builder, Bio.Iteratee.ZLib, Bio.PriorityQueue, Bio.TwoBit, - Bio.Util, + Bio.Util.AD, + Bio.Util.AD2, + Bio.Util.Numeric, + Bio.Util.Regex, + Data.MiniFloat, Data.Avro Other-modules: Paths_biohazard Build-depends: aeson >= 0.7 && < 0.9, - array >= 0.4 && < 0.6, async == 2.0.*, attoparsec >= 0.10 && < 0.13, - base >= 4.5 && < 4.9, - binary >= 0.7 && < 0.8, + base >= 4.5 && < 4.10, + binary >= 0.7 && < 0.9, bytestring >= 0.10.2 && < 0.11, bytestring-mmap >= 0.2 && < 1.0, containers >= 0.4.1 && < 0.6, + deepseq >= 1.3 && < 1.5, directory >= 1.2 && < 2.0, exceptions >= 0.6 && < 0.9, filepath >= 1.3 && < 2.0, iteratee >= 0.8.9.6 && < 0.8.10, ListLike >= 3.0 && < 5.0, + nonlinear-optimization == 0.3.*, primitive >= 0.5 && < 0.7, random >= 1.0 && < 1.2, scientific == 0.3.*, stm == 2.4.*, template-haskell == 2.*, text >= 1.0 && < 2.0, - transformers >= 0.3 && < 0.5, - unix == 2.*, + transformers >= 0.3 && < 0.6, + unix >= 2.5 && < 2.8, unordered-containers >= 0.2.3 && < 0.3, Vec == 1.*, vector >= 0.9 && < 0.11, @@ -117,8 +124,88 @@ -- 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 + Build-depends: async, + base, + biohazard, + filepath, + unordered-containers, + text, + vector + +Executable redeye-div + Main-is: redeye-div.hs + Ghc-options: -Wall -auto-all -threaded -rtsopts -with-rtsopts=-N + Hs-Source-Dirs: tools + Build-depends: async, + base, + biohazard, + filepath, + hmatrix == 0.16.*, + unordered-containers, + text, + vector + +Executable redeye-pileup + Main-is: redeye-pileup.hs + Ghc-options: -Wall -auto-all -threaded -rtsopts -with-rtsopts=-N + Hs-Source-Dirs: tools + Build-depends: aeson, + base, + biohazard, + bytestring, + containers, + directory, + filepath, + iteratee, + text, + unordered-containers, + Vec, + vector + +Executable redeye-single + Main-is: redeye-single.hs + Ghc-options: -Wall -auto-all -rtsopts + Hs-Source-Dirs: tools + Build-depends: aeson, + base, + biohazard, + bytestring, + containers, + directory, + iteratee, + filepath, + text, + unix, + unordered-containers, + vector + +Executable gt-scan + Main-is: gt-scan.hs + Ghc-options: -Wall -auto-all + Hs-Source-Dirs: tools + Build-depends: aeson, + base, + biohazard, + bytestring, + containers, + iteratee, + nonlinear-optimization, + primitive, + strict == 0.3.*, + unordered-containers, + text, + vector + +-- ------ + Executable afroengineer - Main-Is: afroengineer.hs + Main-is: afroengineer.hs + Ghc-options: -Wall -auto-all Hs-source-dirs: tools -- Ghc-options: -Wall -auto-all -threaded -rtsopts -with-rtsopts=-N Ghc-options: -Wall -auto-all -rtsopts @@ -134,6 +221,7 @@ Executable bam-fixpair Main-is: bam-fixpair.hs + Ghc-options: -Wall -auto-all Hs-Source-Dirs: tools -- Ghc-options: -Wall -auto-all -threaded -rtsopts -with-rtsopts=-N Ghc-options: -Wall -auto-all -rtsopts @@ -142,10 +230,12 @@ biohazard, bytestring, hashable >= 1.0 && < 1.3, - transformers + transformers, + vector Executable bam-meld Main-is: bam-meld.hs + Ghc-options: -Wall -auto-all Hs-Source-Dirs: tools -- Ghc-options: -Wall -auto-all -threaded -rtsopts -with-rtsopts=-N Ghc-options: -Wall -auto-all -rtsopts @@ -156,6 +246,7 @@ Executable bam-resample Main-is: bam-resample.hs + Ghc-options: -Wall -auto-all Hs-Source-Dirs: tools -- Ghc-options: -Wall -auto-all -threaded -rtsopts -with-rtsopts=-N Ghc-options: -Wall -auto-all -rtsopts @@ -166,6 +257,7 @@ Executable bam-rewrap Main-is: bam-rewrap.hs + Ghc-options: -Wall -auto-all Hs-Source-Dirs: tools -- Ghc-options: -Wall -auto-all -threaded -rtsopts -with-rtsopts=-N Ghc-options: -Wall -auto-all -rtsopts @@ -176,6 +268,7 @@ Executable bam-rmdup Main-is: bam-rmdup.hs + Ghc-options: -Wall -auto-all Hs-Source-Dirs: tools -- Ghc-options: -Wall -auto-all -threaded -rtsopts -with-rtsopts=-N Ghc-options: -Wall -auto-all -rtsopts @@ -190,6 +283,7 @@ Executable bam-trim Main-is: bam-trim.hs + Ghc-options: -Wall -auto-all Hs-Source-Dirs: tools -- Ghc-options: -Wall -auto-all -threaded -rtsopts -with-rtsopts=-N Ghc-options: -Wall -auto-all -rtsopts @@ -197,29 +291,9 @@ biohazard, bytestring -Executable count-coverage - Main-is: count-coverage.hs - Hs-Source-Dirs: tools - -- Ghc-options: -Wall -auto-all -threaded -rtsopts -with-rtsopts=-N - Ghc-options: -Wall -auto-all -rtsopts - Build-depends: base, - biohazard, - iteratee - -Executable dmg-est - Main-is: dmg-est.hs - Hs-Source-Dirs: tools - -- Ghc-options: -Wall -auto-all -threaded -rtsopts -with-rtsopts=-N - Ghc-options: -Wall -auto-all -rtsopts - Other-Modules: AD - Build-depends: async, - base, - biohazard, - nonlinear-optimization == 0.3.*, - vector - Executable fastq2bam Main-is: fastq2bam.hs + Ghc-options: -Wall -auto-all Hs-Source-Dirs: tools -- Ghc-options: -Wall -auto-all -threaded -rtsopts -with-rtsopts=-N Ghc-options: -Wall -auto-all -rtsopts @@ -230,31 +304,6 @@ iteratee, vector -Executable glf-consensus - Main-is: glf-consensus.hs - Hs-Source-Dirs: tools - -- Ghc-options: -Wall -auto-all -threaded -rtsopts -with-rtsopts=-N - Ghc-options: -Wall -auto-all -rtsopts - Build-depends: base, - biohazard, - bytestring, - containers, - exceptions, - iteratee - -Executable gt-call - Main-is: gt-call.hs - Hs-Source-Dirs: tools - -- Ghc-options: -Wall -auto-all -threaded -rtsopts -with-rtsopts=-N - Ghc-options: -Wall -auto-all -rtsopts - Build-depends: base, - biohazard, - bytestring, - deepseq, - iteratee, - text, - vector - Executable jivebunny Main-is: jivebunny.hs Hs-Source-Dirs: tools @@ -279,7 +328,8 @@ vector-th-unbox Executable mt-anno - Main-Is: mt-anno.hs + Main-is: mt-anno.hs + Ghc-options: -Wall -auto-all Hs-Source-Dirs: tools -- Ghc-options: -Wall -auto-all -threaded -rtsopts -with-rtsopts=-N Ghc-options: -Wall -auto-all -rtsopts @@ -290,7 +340,8 @@ containers Executable mt-ccheck - Main-Is: mt-ccheck.hs + Main-is: mt-ccheck.hs + Ghc-options: -Wall -auto-all Hs-Source-Dirs: tools -- Ghc-options: -Wall -auto-all -threaded -rtsopts -with-rtsopts=-N Ghc-options: -Wall -auto-all -rtsopts @@ -299,13 +350,5 @@ biohazard, containers, unordered-containers - -executable wiggle-coverage - Main-is: wiggle-coverage.hs - Hs-Source-Dirs: tools - -- Ghc-options: -Wall -auto-all -threaded -rtsopts -with-rtsopts=-N - Ghc-options: -Wall -auto-all -rtsopts - Build-depends: base, - biohazard -- :vim:tw=132: diff --git a/data/index_db.json b/data/index_db.json --- a/data/index_db.json +++ b/data/index_db.json @@ -502,6 +502,8 @@ "710": "CGAGGCTG", "711": "AAGAGGCA", "712": "GTAGAGGA", + "714": "GCTCATGA", + "715": "ATCTCAGG", "716": "ACTCGCTA", "718": "GGAGCTAC", "719": "GCGTAGTA", diff --git a/doc/genotyping.tex b/doc/genotyping.tex --- a/doc/genotyping.tex +++ b/doc/genotyping.tex @@ -7,6 +7,7 @@ \newcommandx{\beware}[2][1=]{\todo[nolist,noline,inline,linecolor=red,backgroundcolor=red!25,bordercolor=red,#1]{#2}} \newcommandx{\idea}[2][1=]{\todo[nolist,noline,inline,linecolor=blue,backgroundcolor=blue!25,bordercolor=blue,#1]{#2}} \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}} \title{Deamination Aware Genotype Calling} \author{Udo Stenzel} @@ -234,14 +235,74 @@ nonsense, since $X_i$ already happened. So we must condition on it and write $P(H_i|X_i,Q_i,G)$.} +\section{Parameter Fitting for Single Sample} + +We seek to fit the Allele Frequency Spectrum to a single sample, or in +other words, to estimate heterozygosity and divergence. Let $X,Y$ be +the genotype at some site and $R$ be the reference allele. Now set: + +\begin{equation*} +P(x,y|r) = \left\{ \begin{array}{cl} + dh/3 & \mbox{if} \quad x \neq y \wedge (x=r \vee y=r) \\ + d(h-1)/3 & \mbox{if} \quad x=y \wedge x \neq r \\ + 1 - d & \mbox{if} \quad x=y \wedge x=r \\ + 0 & \mbox{otherwise} + \end{array} \right. +\end{equation*} + +Here we assume that all heterozygous mutations happen at the same +uniform rate, and all homozygous mutations at a different uniform rate. +We declare heterozygous genotypes with two alternative alleles to be +impossibru, getting us out of the need to fit a third parameter, which +would have little impact anyway. We get a divergence of $d$ and a +heterozygosity of $dh$. Labelling the genotype likelihoods as $G_{XY}$ +and assuming the reference allele is $A$ for convenience, we can compute +the likelihood per site: + +\begin{align*} +L &= G_{AA} \left( 1-d \right) + \frac{d(1-h)}{3} \left( G_{CC} + G_{GG} + G_{TT} \right) + + \frac{dh}{3} \left( G_{AC} + G_{AG} + G_{AT} \right) \\ + &= G_{AA} \left( 1-d + d(1-h) \frac{G_{CC} + G_{GG} + G_{TT}}{3 G_{AA}} + + dh \frac{G_{AC} + G_{AG} + G_{AT}}{3 G_{AA}} \right) +\end{align*} + +We now define three quantities of interest, which we sort by magnitude. +The smallest one is accumulated directly, the differences to the other +two are discretized and tabulated. Just six small tables are needed. + +\begin{equation*} +R_i := \ln 3 G_{AA}, \quad +D_i := \ln G_{CC} + G_{GG} + G_{TT}, \quad +H_i := \ln G_{AC} + G_{AG} + G_{AT} +\end{equation*} + +We further set $\delta = \ln \frac{d}{1-d}$ and $\eta = \ln +\frac{h}{1-h}$, because the log-odds-ratios pose fewer numerical +problems and provide cleaner confidence intervals. Thus we recover the +log-likelihood as follows: + +\begin{align*} +L_l &= - \sum_{i} R_i - \sum_{i} \ln \left( \frac{1}{1+e^\delta} + + \frac{e^\delta}{1+e^\delta} \frac{1}{1+e^\eta} e^{D_i} + + \frac{e^\delta}{1+e^\delta} \frac{e^\eta}{1+e^\eta} e^{H_i} \right) +\end{align*} + +\idea{Need to do something similar for indels.} + +\result{Hand-crafting partial derivatives of this equation didn't seem to +yield anything that simplifies, so I applied Automatic +Differentiation and handed it over to Hager-Zhang.} + \section{Testing Method} \subsection{Handcrafted Data} -To test for egregious bugs, we can write a couple of SMA or BAM files by +To test for egregious bugs, we can write a couple of SAM or BAM files by hand. This shouldn't really be called a test; it's ordinary, boring debugging. +\result{Nothing to see here. Code runs.} + \subsection{Simulated Data} For all of the simulations, the genome used does not matter at all. For @@ -260,21 +321,27 @@ into the genotype calling, which makes testing harder while providing no insight at all.} +\beware{It is not actually possible to test the precision of a genotype +caller by counting miscalls. In any such tests, a genotype callers that +doesn't try to call heterozygotes ``wins''. Since that's undesirable, +``simply counting errors'' is a terrible testing method.} + \paragraph{Simulated Modern Data} Starting from a genome with known divergence and heterozygosity, we simulate plain reads with some sequencing error and suitable quality -scores, then genotype call. +scores, then call genotypes. -Called genotypes can be compared to the correct genome, but more -importantly, parameters (divergence, heterozygosity) should be fitted -and compared to their true values, particularly at low (roughly one or -twofold) coverage. +Parameters (divergence, heterozygosity) should be fitted and compared to +their true values, particularly at low (roughly one or twofold) +coverage. \beware{There is no point in simulating fancy sequencing error. Here, we assume the simple model is correct and show that maximum likelihood estimation of parameters works in this setting.} +\result{We accidentally skipped this part.} + \paragraph{Simulated Ancient Data} This is the same idea, this time including damage with known parameters. @@ -295,20 +362,34 @@ use an empirical distribution of overhang lengths, if that could be obtained.} -\idea{Since damage should not correlate with genotypes, estimating -damage in a separate first pass might work and would be lot cheaper, +\result{Since damage should not correlate with genotypes, estimating +damage in a separate first pass works and is a lot cheaper, both conceptually and operationally, than co-estimating damage with -heterozygosity. This is a good time to try it.} +heterozygosity.} +\result{Homa simulated two diploid genomes with divergence and aligned, +damaged reads from these. I estimated damage from the reads, and the +estimate is nearly spot on. On a dataset with 20\% divergence, the +estimated divergence and heterozygosity are unaffected by damage, +probably because damage is a minor effect compared to divergence here. +On a 0.1\% divergent dataset, naive genotype calling overestimates +heterozygosity by a factor of 15, but gets homozygous divergence right +(estimated $D=1.56\%, H/D=0.93$, but should be $D=0.3\%, H/D=0.67$. +Running with the estimated damage parameters gives a practically perfect +result.} + + \subsection{Real Sequencing Data} \paragraph{Clean, high-coverage, modern, haploid data} We need actual sequencing data from a haploid region at sensible coverage. The goal is to test the two available error models in a -setting without confounding factors, especially heterozygosity. -This should be used to select the better error model and to fit the -$\theta$ parameter if the \texttt{Maq} model is selected. +setting without confounding factors, especially heterozygosity. This +should be used to select the better error model and to fit the $\theta$ +parameter by the maximum likelihood method, if the \texttt{Maq} model is +selected. It's a good idea to check whether $\theta$ varies between +sample sor sequencing runs. The haploid region of choice might be the mitochondrion, which is haploid, but the data will be somewhat contaminated with nuMT sequences. @@ -316,6 +397,12 @@ male specimen would work, here the difficulty is to find that unique region. +\idea{A couple of variations on the error model need to be investigated: +handle in bases in ascending, descending, random or input order of +quality; treat errors on different strands as dependent or independent; +various values of $\theta$. For $\theta = 1$, all of these should match +the naïve error model.} + \paragraph{Clean, high-coverage, modern data, two mixed haploids} Just like the previous test, but this time with two haploid samples @@ -324,6 +411,10 @@ test which error model is better and assess the correctness of the calls. +\beware{Counting miscalls suffers from the usual problem that it's +entirely unclear how to count errors. The plan is therefore to estimate +heterozygosity.} + \paragraph{Clean, high-coverage, diploid modern data} We test the two error models and select the better one. This must be an @@ -338,10 +429,10 @@ \paragraph{Clean, low-coverage, modern data} Assuming we fixed the error model, assuming we can reliably estimate -difficult parameters like heterozygosity, he we investigate the bahviour -at low coverage. The sample can be a high coverage sample suitably -downsampled. In this case, we have an expectation for the estimated -parameters. +difficult parameters like heterozygosity, here we investigate the +behaviour at low coverage. The sample can be a high coverage sample +suitably downsampled. In this case, we have a reasonable expectation +for the estimated parameters. \paragraph{Ancient data, one mitochondrion} @@ -358,7 +449,7 @@ To investigate interaction of heterozygosity, deamination, error model in a setting where true heterozygous genotypes are known. Data should be clean and ideally from the same run (we don't want to deal with -additional contamination and different error profiles). The assumtion +additional contamination and different error profiles). The assumption is that we can correctly call either sample on its own. In principle, if we haven't encountered difficulties so far, this should @@ -435,29 +526,72 @@ \subsection{Covariance-Matrix as Prior} When co-calling individuals from multiple populations, the correct prior -for the genotypes would be based on a covariance matrix. Estimating -that matrix allows Treemix, Patterson's~D and Pruefer's Divergence. +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) 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 +straight forward, but it would require integrating over the space of possible combinations of allele frequencies, which sounds impractical, -and becomes more impractical the more samples are considered. +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. -Instead, we can estimate the joint probability (genotype(s), allele -frequency) and maximize that, which is much easier\footnote{Effectively, -we estimate the allele frequency at every position for every sample. -Which is impossible, but the aggregate makes sense for populations.}. -Only summation over the possible genotypes is necesssary, which is just -10 per individual, and individuals are independent; allele frequencies -and covariance matrix are co-estimated using something resembling the EM -algorithm. (One idea would be to not store the aforementioned 600GB of -likelihoods, but only 6GB or thereabout of allele frequency data. The -likelihoods can be generated from the smaller BAM files on the fly.) +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.}. +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 +\emph{one} allele frequency at a time, holding the others +constant\todo{Is this correct? At any rate, even if not, multivariate +maximization should do it.}. Maximization is finding a covariance +matrix from allele frequencies, and this is again easy. The obvious +moment based estimator should work. +\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.} + +\idea{Dealing with contamination becomes obvious in this framework. For +a contaminated sample, every position has two allele frequencies (or two +sets of three allele frequencies). We need to assign a contaminant +probability to each read, which should derive from the alignment score +given the currently estimated allele frequencies.} + +\subsection{Principal Component Analysis} + +\beware{I personally think PCA is not an analysis, it barely(!) serves +as visualisation method. Nonetheless, it is frequently requested.} + +``Modern'' PCA starts with a matrix where rows correspond to $m$ +individuals and columns to $n$ markers, such as allele frequencies. +Means are subtracted from each column, then each column is divided by +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.) + + + \appendix \section{Random Base vs. Random Error} @@ -534,6 +668,11 @@ \emph{mapDamage: testing for damage patterns in ancient DNA sequences}. Bioinformatics 27 (15), 2153--2155 (2011). +\bibitem{seqem} + E. R. Martin et. al., + \emph{SeqEM: an adaptive genotype-calling approach for next-generation + sequencing studies}. + Bioinformatics 26 (22), 2803-2810 (2010). \end{thebibliography} \end{document} diff --git a/src/Bio/Align.hs b/src/Bio/Align.hs --- a/src/Bio/Align.hs +++ b/src/Bio/Align.hs @@ -42,9 +42,9 @@ -- -- The algorithm is the O(nd) algorithm by Myers, implemented in C. A -- gap and a mismatch score the same. The strings are supposed to code --- for DNA, the code understands IUPAC ambiguity codes. Two characters --- match iff there is at least one nucleotide both can code for. Note --- that N is a wildcard, while X matches nothing. +-- for DNA, the code understands IUPAC-IUB ambiguity codes. Two +-- characters match iff there is at least one nucleotide both can code +-- for. Note that N is a wildcard, while X matches nothing. myersAlign :: Int -> S.ByteString -> Mode -> S.ByteString -> (Int, S.ByteString, S.ByteString) myersAlign maxd seqA mode seqB = diff --git a/src/Bio/Bam.hs b/src/Bio/Bam.hs --- a/src/Bio/Bam.hs +++ b/src/Bio/Bam.hs @@ -1,14 +1,24 @@ -module Bio.Bam ( module X ) where +module Bio.Bam ( + module Bio.Bam.Fastq, + module Bio.Bam.Filter, + module Bio.Bam.Header, + module Bio.Bam.Index, + module Bio.Bam.Reader, + module Bio.Bam.Rec, + module Bio.Bam.Trim, + module Bio.Bam.Writer, + module Bio.Iteratee + ) where -import Bio.Bam.Fastq as X -import Bio.Bam.Filter as X -import Bio.Bam.Header as X -import Bio.Bam.Index as X -import Bio.Bam.Reader as X -import Bio.Bam.Rec as X -import Bio.Bam.Trim as X -import Bio.Bam.Writer as X -import Bio.Iteratee as X +import Bio.Bam.Fastq +import Bio.Bam.Filter +import Bio.Bam.Header +import Bio.Bam.Index +import Bio.Bam.Reader +import Bio.Bam.Rec +import Bio.Bam.Trim +import Bio.Bam.Writer +import Bio.Iteratee -- ^ Umbrella module for most of what's under 'Bio.Bam'. diff --git a/src/Bio/Bam/Fastq.hs b/src/Bio/Bam/Fastq.hs --- a/src/Bio/Bam/Fastq.hs +++ b/src/Bio/Bam/Fastq.hs @@ -58,12 +58,13 @@ -- start with @\>@ or @\@@, we treat both equally. The first word of -- the header becomes the read name, the remainder of the header is -- ignored. The sequence can be split across multiple lines; --- whitespace, dashes and dots are ignored, IUPAC ambiguity codes are --- accepted as bases, anything else causes an error. The sequence ends --- at a line that is either a header or starts with @\+@, in the latter --- case, that line is ignored and must be followed by quality scores. --- There must be exactly as many Q-scores as there are bases, followed --- immediately by a header or end-of-file. Whitespace is ignored. +-- whitespace, dashes and dots are ignored, IUPAC-IUB ambiguity codes +-- are accepted as bases, anything else causes an error. The sequence +-- ends at a line that is either a header or starts with @\+@, in the +-- latter case, that line is ignored and must be followed by quality +-- scores. There must be exactly as many Q-scores as there are bases, +-- followed immediately by a header or end-of-file. Whitespace is +-- ignored. {-# WARNING parseFastq "parseFastq no longer removes syntactic warts!" #-} parseFastq :: Monad m => Enumeratee S.ByteString [ BamRec ] m a diff --git a/src/Bio/Bam/Header.hs b/src/Bio/Bam/Header.hs --- a/src/Bio/Bam/Header.hs +++ b/src/Bio/Bam/Header.hs @@ -54,7 +54,6 @@ import Data.ByteString.Builder import Data.Ix import Data.List ( (\\), foldl' ) -import Data.Monoid import Data.Sequence ( (><), (|>) ) import Data.String import Data.Version ( Version, showVersion ) @@ -299,8 +298,8 @@ getRef :: Refs -> Refseq -> BamSQ getRef refs (Refseq i) - | 0 <= i && fromIntegral i <= Z.length refs = Z.index refs (fromIntegral i) - | otherwise = BamSQ "*" 0 [] + | 0 <= i && fromIntegral i < Z.length refs = Z.index refs (fromIntegral i) + | otherwise = BamSQ "*" 0 [] flagPaired, flagProperlyPaired, flagUnmapped, flagMateUnmapped, flagReversed, flagMateReversed, flagFirstMate, flagSecondMate, diff --git a/src/Bio/Bam/Pileup.hs b/src/Bio/Bam/Pileup.hs --- a/src/Bio/Bam/Pileup.hs +++ b/src/Bio/Bam/Pileup.hs @@ -2,22 +2,17 @@ {-# OPTIONS_GHC -funbox-strict-fields #-} module Bio.Bam.Pileup where --- import Text.Printf - import Bio.Base import Bio.Bam.Header import Bio.Bam.Rec -import Bio.Genocall.Adna import Bio.Iteratee -import Control.Arrow ( (&&&) ) import Control.Applicative import Control.Monad hiding ( mapM_ ) import Control.Monad.Fix ( fix ) import Data.Foldable hiding ( sum, product ) -import Data.Monoid import Data.Ord -import Data.Vec.Packed ( Mat44D, packMat ) +import Data.Vec.Packed ( Mat44D ) import qualified Data.ByteString as B import qualified Data.Vector.Generic as V @@ -79,7 +74,6 @@ -- *TODO* -- -- * A whole lot of testing. --- * Actual genotype calling. -- * ML fitting and evaluation of parameters for different possible -- error and damage models. -- * Maybe specialize to ploidy one and two. @@ -89,18 +83,19 @@ -- length of a deleted sequence. The logic is that we look at a base -- followed by some indel, and all those indels are combined into a -- single insertion and a single deletion. -data PrimChunks = Seek !Int !PrimBase -- ^ skip to position (at start or after N operation) - | Indel !Int [DamagedBase] !PrimBase -- ^ observed deletion and insertion between two bases +data PrimChunks = Seek Int PrimBase -- ^ skip to position (at start or after N operation) + | Indel [Nucleotides] [DamagedBase] PrimBase -- ^ observed deletion and insertion between two bases | EndOfRead -- ^ nothing anymore deriving Show -data PrimBase = Base { _pb_wait :: !Int -- ^ number of bases to wait due to a deletion - , _pb_likes :: !DamagedBase -- ^ four likelihoods - , _pb_mapq :: !Qual -- ^ map quality - , _pb_rev :: !Bool -- ^ reverse strand? +data PrimBase = Base { _pb_wait :: Int -- ^ number of bases to wait due to a deletion + , _pb_likes :: DamagedBase -- ^ four likelihoods + , _pb_mapq :: Qual -- ^ map quality + , _pb_rev :: Bool -- ^ reverse strand? , _pb_chunks :: PrimChunks } -- ^ more chunks deriving Show +type PosPrimChunks = (Refseq, Int, PrimChunks) -- | Represents our knowledge about a certain base, which consists of -- the base itself (A,C,G,T, encoded as 0..3; no Ns), the quality score @@ -111,33 +106,33 @@ -- Unfortunately, none of this can be rolled into something more simple, -- because damage and sequencing error behave so differently. -data DamagedBase = DB { db_call :: !Nucleotide - , db_qual :: !Qual - , db_dmg :: !Mat44D } +data DamagedBase = DB { db_call :: {-# UNPACK #-} !Nucleotide -- ^ called base + , db_qual :: {-# UNPACK #-} !Qual -- ^ quality of called base + , db_ref :: {-# UNPACK #-} !Nucleotides -- ^ reference base from MD field + , db_dmg :: {-# UNPACK #-} !Mat44D } -- ^ damage matrix instance Show DamagedBase where - showsPrec _ (DB n q _) = shows n . (:) '@' . shows q + showsPrec _ (DB n q r _) + | nucToNucs n == r = shows n . (:) '@' . shows q + | otherwise = shows n . (:) '/' . shows r . (:) '@' . shows q -- | Decomposes a BAM record into chunks suitable for piling up. We --- pick apart the CIGAR field, and combine it with sequence and quality --- as appropriate. We ignore the @MD@ field, even if it is present. --- Clipped bases are removed/skipped as appropriate. We also ignore the --- reference allele, in fact, we don't even know it, which nicely avoids --- any possible reference bias by construction. But we do apply a --- substitution matrix to each base, which must be supplied along with --- the read. +-- pick apart the CIGAR and MD fields, and combine them with sequence +-- and quality as appropriate. Clipped bases are removed/skipped as +-- appropriate. We also do apply a substitution matrix to each base, +-- which must be supplied along with the read. -decompose :: BamRaw -> [Mat44D] -> PrimChunks -decompose br matrices - | isUnmapped b || b_rname == invalidRefseq = EndOfRead - | otherwise = firstBase b_pos 0 0 matrices +decompose :: [Mat44D] -> BamRaw -> Maybe PosPrimChunks +decompose matrices br = + if isUnmapped b || isDuplicate b || not (isValidRefseq b_rname) + then Nothing else Just (b_rname, b_pos, pchunks) where b@BamRec{..} = unpackBam br + pchunks = firstBase b_pos 0 0 (maybe [] id $ getMd b) matrices !max_cig = V.length b_cigar !max_seq = V.length b_seq - -- !mapq = br_mapq br !baq = extAsString "BQ" b -- This will compute the effective quality. As far as I can see @@ -145,7 +140,7 @@ -- and BAQ. If QUAL is invalid, we replace it (arbitrarily) with -- 23 (assuming a rather conservative error rate of ~0.5%), BAQ is -- added to QUAL, and MAPQ is an upper limit for effective quality. - get_seq :: Int -> Mat44D -> DamagedBase + get_seq :: Int -> Nucleotides -> Mat44D -> DamagedBase get_seq i = case b_seq V.! i of -- nucleotide n | n == nucsA -> DB nucA qe | n == nucsC -> DB nucC qe @@ -158,30 +153,65 @@ | otherwise = Q (unQ q + (B.index baq i - 64)) -- else correct for BAQ !qe = min q' b_mapq -- use MAPQ as upper limit + get_seq' :: Int -> Mat44D -> DamagedBase + get_seq' i = case b_seq V.! i of -- nucleotide + n | n == nucsA -> DB nucA qe nucsA + | n == nucsC -> DB nucC qe nucsC + | n == nucsG -> DB nucG qe nucsG + | n == nucsT -> DB nucT qe nucsT + | otherwise -> DB nucA (Q 0) n + where + !q = case b_qual V.! i of Q 0xff -> Q 30 ; x -> x -- quality; invalid (0xff) becomes 30 + !q' | i >= B.length baq = q -- no BAQ available + | otherwise = Q (unQ q + (B.index baq i - 64)) -- else correct for BAQ + !qe = min q' b_mapq -- use MAPQ as upper limit + -- Look for first base following the read's start or a gap (CIGAR -- code N). Indels are skipped, since these are either bugs in the -- aligner or the aligner getting rid of essentially unalignable -- bases. - firstBase :: Int -> Int -> Int -> [Mat44D] -> PrimChunks - firstBase !_ !_ !_ [ ] = EndOfRead - firstBase !pos !is !ic mms@(m:ms) + firstBase :: Int -> Int -> Int -> [MdOp] -> [Mat44D] -> PrimChunks + firstBase !_ !_ !_ _ [ ] = EndOfRead + firstBase !pos !is !ic mds mms@(m:ms) | is >= max_seq || ic >= max_cig = EndOfRead | otherwise = case b_cigar V.! ic of - Ins :* cl -> firstBase pos (cl+is) (ic+1) mms - SMa :* cl -> firstBase pos (cl+is) (ic+1) mms - Del :* cl -> firstBase (pos+cl) is (ic+1) mms - Nop :* cl -> firstBase (pos+cl) is (ic+1) mms - HMa :* _ -> firstBase pos is (ic+1) mms - Pad :* _ -> firstBase pos is (ic+1) mms - Mat :* 0 -> firstBase pos is (ic+1) mms - Mat :* _ -> Seek pos $ nextBase 0 pos is ic 0 m ms + Ins :* cl -> firstBase pos (cl+is) (ic+1) mds mms + SMa :* cl -> firstBase pos (cl+is) (ic+1) mds mms + Del :* cl -> firstBase (pos+cl) is (ic+1) (drop_del cl mds) mms + Nop :* cl -> firstBase (pos+cl) is (ic+1) mds mms + HMa :* _ -> firstBase pos is (ic+1) mds mms + Pad :* _ -> firstBase pos is (ic+1) mds mms + Mat :* 0 -> firstBase pos is (ic+1) mds mms + Mat :* _ -> Seek pos $ nextBase 0 pos is ic 0 mds m ms + where + -- We have to treat (MdNum 0), because samtools actually + -- generates(!) it all over the place and if not handled as a + -- special case, it looks like an incinsistend MD field. + drop_del n (MdDel ns : mds') + | n < length ns = MdDel (drop n ns) : mds' + | n > length ns = drop_del (n - length ns) mds' + | otherwise = mds' + drop_del n (MdNum 0 : mds') = drop_del n mds' + drop_del _ mds' = mds' -- Generate likelihoods for the next base. When this gets called, -- we are looking at an M CIGAR operation and all the subindices are -- valid. - nextBase :: Int -> Int -> Int -> Int -> Int -> Mat44D -> [Mat44D] -> PrimBase - nextBase !wt !pos !is !ic !io m ms = Base wt (get_seq is m) b_mapq (isReversed b) - $ nextIndel [] 0 (pos+1) (is+1) ic (io+1) ms + -- I don't think we can ever get (MdDel []), but then again, who + -- knows what crazy shit samtools decides to generate. There is + -- little harm in special-casing it. + nextBase :: Int -> Int -> Int -> Int -> Int -> [MdOp] -> Mat44D -> [Mat44D] -> PrimBase + nextBase !wt !pos !is !ic !io mds m ms = case mds of + MdNum 0 : mds' -> nextBase wt pos is ic io mds' m ms + MdDel [] : mds' -> nextBase wt pos is ic io mds' m ms + MdNum 1 : mds' -> nextBase' (get_seq' is m) mds' + MdNum n : mds' -> nextBase' (get_seq' is m) (MdNum (n-1) : mds') + MdRep ref : mds' -> nextBase' (get_seq is ref m) mds' + MdDel _ : _ -> nextBase' (get_seq is nucsN m) mds + [ ] -> nextBase' (get_seq is nucsN m) [ ] + where + nextBase' ref mds' = Base wt ref b_mapq (isReversed b) + $ nextIndel [] [] (pos+1) (is+1) ic (io+1) mds' ms -- Look for the next indel after a base. We collect all indels (I -- and D codes) into one combined operation. If we hit N or the @@ -190,33 +220,46 @@ -- isn't valid in the middle of a read (H and S), but then what -- would we do about it anyway? Just ignoring it is much easier and -- arguably at least as correct. - nextIndel :: [[DamagedBase]] -> Int -> Int -> Int -> Int -> Int -> [Mat44D] -> PrimChunks - nextIndel _ _ !_ !_ !_ !_ [ ] = EndOfRead - nextIndel ins del !pos !is !ic !io mms@(m:ms) + nextIndel :: [[DamagedBase]] -> [Nucleotides] -> Int -> Int -> Int -> Int -> [MdOp] -> [Mat44D] -> PrimChunks + nextIndel _ _ !_ !_ !_ !_ _ [ ] = EndOfRead + nextIndel ins del !pos !is !ic !io mds mms@(m:ms) | is >= max_seq || ic >= max_cig = EndOfRead | otherwise = case b_cigar V.! ic of - Ins :* cl -> nextIndel (isq cl) del pos (cl+is) (ic+1) 0 (drop cl mms) - SMa :* cl -> nextIndel ins del pos (cl+is) (ic+1) 0 (drop cl mms) - Del :* cl -> nextIndel ins (cl+del) (pos+cl) is (ic+1) 0 mms - Pad :* _ -> nextIndel ins del pos is (ic+1) 0 mms - HMa :* _ -> nextIndel ins del pos is (ic+1) 0 mms - Mat :* cl | io == cl -> nextIndel ins del pos is (ic+1) 0 mms - | otherwise -> Indel del out $ nextBase del pos is ic io m ms -- ends up generating a 'Base' - Nop :* cl -> firstBase (pos+cl) is (ic+1) mms -- ends up generating a 'Seek' + Ins :* cl -> nextIndel (isq cl) del pos (cl+is) (ic+1) 0 mds (drop cl mms) + SMa :* cl -> nextIndel ins del pos (cl+is) (ic+1) 0 mds (drop cl mms) + Del :* cl -> nextIndel ins (del++dsq) (pos+cl) is (ic+1) 0 mds' mms + where (dsq,mds') = split_del cl mds + Pad :* _ -> nextIndel ins del pos is (ic+1) 0 mds mms + HMa :* _ -> nextIndel ins del pos is (ic+1) 0 mds mms + Nop :* cl -> firstBase (pos+cl) is (ic+1) mds mms -- ends up generating a 'Seek' + Mat :* cl | io == cl -> nextIndel ins del pos is (ic+1) 0 mds mms + | otherwise -> indel del out $ nextBase (length del) pos is ic io mds m ms -- ends up generating a 'Base' where + indel d o k = rlist o `seq` Indel d o k out = concat $ reverse ins - isq cl = zipWith ($) [ get_seq i | i <- [is..is+cl-1] ] (take cl mms) : ins + isq cl = zipWith ($) [ get_seq i gap | i <- [is..is+cl-1] ] (take cl mms) : ins + rlist [] = () + rlist (a:as) = a `seq` rlist as + -- We have to treat (MdNum 0), because samtools actually + -- generates(!) it all over the place and if not handled as a + -- special case, it looks like an incinsistend MD field. + split_del n (MdDel ns : mds') + | n < length ns = (take n ns, MdDel (drop n ns) : mds') + | n > length ns = let (ns', mds'') = split_del (n - length ns) mds' in (ns++ns', mds'') + | otherwise = (ns, mds') + split_del n (MdNum 0 : mds') = split_del n mds' + split_del n mds' = (replicate n nucsN, mds') -- | Statistics about a genotype call. Probably only useful for -- fitlering (so not very useful), but we keep them because it's easy to -- track them. -data CallStats = CallStats { read_depth :: !Int -- number of contributing reads - , reads_mapq0 :: !Int -- number of (non-)contributing reads with MAPQ==0 - , sum_mapq :: !Int -- sum of map qualities of contributing reads - , sum_mapq_squared :: !Int } -- sum of squared map qualities of contributing reads - deriving Show +data CallStats = CallStats { read_depth :: {-# UNPACK #-} !Int -- number of contributing reads + , reads_mapq0 :: {-# UNPACK #-} !Int -- number of (non-)contributing reads with MAPQ==0 + , sum_mapq :: {-# UNPACK #-} !Int -- sum of map qualities of contributing reads + , sum_mapq_squared :: {-# UNPACK #-} !Int } -- sum of squared map qualities of contributing reads + deriving (Show, Eq) instance Monoid CallStats where mempty = CallStats { read_depth = 0 @@ -244,32 +287,38 @@ type GL = U.Vector Prob -newtype V_Nuc = V_Nuc (U.Vector Nucleotide) deriving (Eq, Ord, Show) -data IndelVariant = IndelVariant { deleted_bases :: !Int - , inserted_bases :: !V_Nuc } +newtype V_Nuc = V_Nuc (U.Vector Nucleotide) deriving (Eq, Ord, Show) +newtype V_Nucs = V_Nucs (U.Vector Nucleotides) deriving (Eq, Ord, Show) + +data IndelVariant = IndelVariant { deleted_bases :: {-# UNPACK #-} !V_Nucs + , inserted_bases :: {-# UNPACK #-} !V_Nuc } deriving (Eq, Ord, Show) --- Both types of piles carry along the map quality. We'll only need it --- in the case of Indels. -type BasePile = [( Qual, DamagedBase )] -- a list of encountered bases -type IndelPile = [( Qual, (Int, [DamagedBase]) )] -- a list of indel variants +-- | Map quality and a list of encountered bases, with damage +-- information and reference base if known. +type BasePile = [( Qual, DamagedBase )] + +-- | Map quality and a list of encountered indel variants. The deletion +-- has the reference sequence, if known, an insertion has the inserted +-- sequence with damage information. +type IndelPile = [( Qual, ([Nucleotides], [DamagedBase]) )] -- a list of indel variants + -- | Running pileup results in a series of piles. A 'Pile' has the -- basic statistics of a 'VarCall', but no GL values and a pristine list -- of variants instead of a proper call. We emit one pile with two -- 'BasePile's (one for each strand) and one 'IndelPile' (the one -- immediately following) at a time. -data Pile' a b = Pile { p_refseq :: !Refseq - , p_pos :: !Int - , p_snp_stat :: !CallStats +data Pile' a b = Pile { p_refseq :: {-# UNPACK #-} !Refseq + , p_pos :: {-# UNPACK #-} !Int + , p_snp_stat :: {-# UNPACK #-} !CallStats , p_snp_pile :: a - , p_indel_stat :: !CallStats + , p_indel_stat :: {-# UNPACK #-} !CallStats , p_indel_pile :: b } deriving Show type Pile = Pile' (BasePile, BasePile) IndelPile -type Calls = Pile' GL (GL, [IndelVariant]) -- | The pileup enumeratee takes 'BamRaw's, decomposes them, interleaves -- the pieces appropriately, and generates 'Pile's. The output will @@ -281,14 +330,11 @@ -- Processing stops when the first read with invalid 'br_rname' is -- encountered or a t end of file. -pileup :: Monad m => DamageModel Double -> Enumeratee [BamRaw] [Pile] m a -pileup dm = takeWhileE (isValidRefseq . b_rname . unpackBam) ><> filterStream useable ><> - eneeCheckIfDonePass (icont . runPileM pileup' finish (Refseq 0) 0 [] Empty dm) +pileup :: Monad m => Enumeratee [PosPrimChunks] [Pile] m a +pileup = eneeCheckIfDonePass (icont . runPileM pileup' finish (Refseq 0) 0 [] Empty) where - useable = not . (\b -> isUnmapped b || isDuplicate b) . unpackBam - - finish () _r _p [] Empty _dm out inp = idone (liftI out) inp - finish () _ _ _ _ _ _ _ = error "logic error: leftovers after pileup" + finish () _r _p [] Empty out inp = idone (liftI out) inp + finish () _ _ _ _ _ _ = error "logic error: leftovers after pileup" -- | The pileup logic keeps a current coordinate (just two integers) and @@ -315,99 +361,130 @@ type PileF m r = Refseq -> Int -> -- current position [PrimBase] -> -- active queue Heap -> -- waiting queue - DamageModel Double -> (Stream [Pile] -> Iteratee [Pile] m r) -> -- output function - Stream [BamRaw] -> -- pending input - Iteratee [BamRaw] m (Iteratee [Pile] m r) + Stream [PosPrimChunks] -> -- pending input + Iteratee [PosPrimChunks] m (Iteratee [Pile] m r) instance Functor (PileM m) where + {-# INLINE fmap #-} fmap f (PileM m) = PileM $ \k -> m (k . f) instance Applicative (PileM m) where + {-# INLINE pure #-} pure a = PileM $ \k -> k a + {-# INLINE (<*>) #-} u <*> v = PileM $ \k -> runPileM u (\a -> runPileM v (k . a)) instance Monad (PileM m) where + {-# INLINE return #-} return a = PileM $ \k -> k a + {-# INLINE (>>=) #-} m >>= k = PileM $ \k' -> runPileM m (\a -> runPileM (k a) k') -instance MonadIO m => MonadIO (PileM m) where - liftIO m = PileM $ \k r p a w d o i -> liftIO m >>= \x -> k x r p a w d o i - +{-# INLINE get_refseq #-} get_refseq :: PileM m Refseq get_refseq = PileM $ \k r -> k r r +{-# INLINE get_pos #-} get_pos :: PileM m Int get_pos = PileM $ \k r p -> k p r p +{-# INLINE upd_pos #-} upd_pos :: (Int -> Int) -> PileM m () upd_pos f = PileM $ \k r p -> k () r $! f p +{-# INLINE set_pos #-} set_pos :: (Refseq, Int) -> PileM m () set_pos (!r,!p) = PileM $ \k _ _ -> k () r p +{-# INLINE get_active #-} get_active :: PileM m [PrimBase] get_active = PileM $ \k r p a -> k a r p a +{-# INLINE upd_active #-} upd_active :: ([PrimBase] -> [PrimBase]) -> PileM m () upd_active f = PileM $ \k r p a -> k () r p $! f a +{-# INLINE add_active #-} +add_active :: PrimBase -> PileM m () +add_active !pb = PileM $ \k r p a -> k () r p (pb:a) + +{-# INLINE clr_active #-} +clr_active :: PileM m [PrimBase] +clr_active = PileM $ \k r p a -> k a r p [] + +{-# INLINE ins_waiting #-} +ins_waiting :: Int -> PrimBase -> PileM m () +ins_waiting !q !v = PileM $ \ k r p a w -> k () r p a $! Node q v Empty Empty `union` w + +{-# INLINE get_waiting #-} get_waiting :: PileM m Heap get_waiting = PileM $ \k r p a w -> k w r p a w -upd_waiting :: (Heap -> Heap) -> PileM m () -upd_waiting f = PileM $ \k r p a w -> k () r p a $! f w - -get_damage_model :: PileM m (DamageModel Double) -get_damage_model = PileM $ \k r p a w d -> k d r p a w d +{-# INLINE set_waiting #-} +set_waiting :: Heap -> PileM m () +set_waiting !w = PileM $ \k r p a _ -> k () r p a w +-- | Sends one piece of output downstream. You are not expected to +-- understand how this works, but at last it doesn't leak memory. +{-# INLINE yield #-} yield :: Monad m => Pile -> PileM m () -yield x = PileM $ \k r p a w d out inp -> - eneeCheckIfDone (\out' -> k () r p a w d out' inp) . out $ Chunk [x] +yield x = PileM $ \ !kont !r !p !a !w !out !inp -> Iteratee $ \od oc -> + let loop = kont () r p a w + onDone y s = od (idone y s) inp + onCont k Nothing = runIter (loop k inp) od oc + onCont k (Just e) = runIter (throwRecoverableErr e (loop k . (<>) inp)) od oc + in runIter (out (Chunk [x])) onDone onCont -- | Inspect next input element, if any. Returns @Just b@ if @b@ is the -- next input element, @Nothing@ if no such element exists. Waits for -- more input if nothing is available immediately. -peek :: PileM m (Maybe BamRaw) -peek = PileM $ \k r p a w d out inp -> case inp of - EOF _ -> k Nothing r p a w d out inp - Chunk [ ] -> liftI $ runPileM peek k r p a w d out - Chunk (b:_) -> k (Just b) r p a w d out inp +{-# INLINE peek #-} +peek :: PileM m (Maybe PosPrimChunks) +peek = PileM $ \k r p a w out inp -> case inp of + EOF _ -> k Nothing r p a w out inp + Chunk [ ] -> liftI $ runPileM peek k r p a w out + Chunk (b:_) -> k (Just b) r p a w out inp -- | Discard next input element, if any. Does nothing if input has -- already ended. Waits for input to discard if nothing is available -- immediately. +{-# INLINE bump #-} bump :: PileM m () -bump = PileM $ \k r p a w d out inp -> case inp of - EOF _ -> k () r p a w d out inp - Chunk [ ] -> liftI $ runPileM bump k r p a w d out - Chunk (_:x) -> k () r p a w d out (Chunk x) +bump = PileM $ \k r p a w out inp -> case inp of + EOF _ -> k () r p a w out inp + Chunk [ ] -> liftI $ runPileM bump k r p a w out + Chunk (_:x) -> k () r p a w out (Chunk x) +{-# INLINE consume_active #-} consume_active :: a -> (a -> PrimBase -> PileM m a) -> PileM m a consume_active nil cons = do ac <- get_active upd_active (const []) foldM cons nil ac --- | The actual pileup algorithm. +-- | The actual pileup algorithm. If /active/ contains something, +-- continue here. Else find the coordinate to continue from, which is +-- the minimum of the next /waiting/ coordinate and the next coordinate +-- in input; if found, continue there, else we're all done. pileup' :: Monad m => PileM m () -pileup' = do - refseq <- get_refseq - active <- get_active - next_waiting <- fmap ((,) refseq) . getMinKey <$> get_waiting - next_input <- fmap ((b_rname &&& b_pos) . unpackBam) <$> peek +pileup' = PileM $ \ !k !refseq !pos !active !waiting !out !inp -> - -- If /active/ contains something, continue here. Else find the coordinate - -- to continue from, which is the minimum of the next /waiting/ coordinate - -- and the next coordinate in input; if found, continue there, else we're - -- all done. - case (active, next_waiting, next_input) of - ( (_:_), _, _ ) -> pileup'' - ( [ ], Just nw, Nothing ) -> set_pos nw >> pileup'' - ( [ ], Nothing, Just ni ) -> set_pos ni >> pileup'' - ( [ ], Just nw, Just ni ) -> set_pos (min nw ni) >> pileup'' - ( [ ], Nothing, Nothing ) -> return () + let recurse = runPileM pileup' k refseq pos active waiting out + cont2 rs po = runPileM pileup'' k rs po active waiting out inp + leave = k () refseq pos active waiting out inp + in case (active, getMinKey waiting, inp) of + ( _:_, _, _ ) -> cont2 refseq pos + ( [ ], Just nw, EOF _ ) -> cont2 refseq nw + ( [ ], Nothing, EOF _ ) -> leave + ( _, _, Chunk [ ] ) -> liftI recurse + ( [ ], Nothing, Chunk ((r,p,_):_) ) -> cont2 r p + ( [ ], Just nw, Chunk ((r,p,_):_) ) + | (refseq,nw) <= (r,p) -> cont2 refseq nw + | otherwise -> cont2 r p + + pileup'' :: Monad m => PileM m () pileup'' = do -- Input is still 'BamRaw', since these can be relied on to be @@ -415,62 +492,19 @@ -- if so, decompose it and add it to the appropriate queue. rs <- get_refseq po <- get_pos - dm <- get_damage_model - -- liftIO $ printf "pileup' @%d:%d, %d active, %d waiting\n" - -- (unRefseq rs) po (-1::Int) (-1::Int) - - -- feed in input as long as it starts at the current position - fix $ \loop -> peek >>= mapM_ (\br -> - let b = unpackBam br - in when (b_rname b == rs && b_pos b == po) $ do - bump - case decompose br $ map packMat $ toList $ dm (isReversed b) (V.length (b_seq b)) of - Seek p pb -> upd_waiting (insert p pb) - Indel _ _ pb -> upd_active (pb:) - EndOfRead -> return () - loop) - - - -- Check /waiting/ queue. If there is anything waiting for the - -- current position, move it to /active/ queue. - fix $ \loop -> (viewMin <$> get_waiting) >>= mapM_ (\(mk,pb,w') -> - when (mk == po) $ do upd_active (pb:) - upd_waiting (const w') - loop) - - -- Scan /active/ queue and make a 'BasePile'. Also see what's next in the - -- 'PrimChunks': 'Indel's contribute to an 'IndelPile', 'Seek's and - -- deletions are pushed back to the /waiting/ queue, 'EndOfRead's are - -- removed, and everything else is added to the fresh /active/ queue. - ((fin_bs, fin_bp), (fin_is, fin_ip)) <- consume_active (mempty, mempty) $ - \(!bpile, !ipile) (Base wt qs mq str pchunks) -> - let put (Q q) x (!st,!vs) = ( st { read_depth = read_depth st + 1 - , reads_mapq0 = reads_mapq0 st + (if q == 0 then 1 else 0) - , sum_mapq = sum_mapq st + fromIntegral q - , sum_mapq_squared = sum_mapq_squared st + fromIntegral q * fromIntegral q } - , (Q q, x) : vs ) - b' = Base (wt-1) qs mq str pchunks - put' = put mq (if str then Left qs else Right qs) - in case pchunks of - _ | wt > 0 -> do upd_active (b' :) ; return ( bpile, ipile ) - Seek p' pb' -> do upd_waiting (insert p' pb') ; return ( put' bpile, ipile ) - Indel del ins pb' -> do upd_active (pb' :) ; return ( put' bpile, put mq (del,ins) ipile ) - EndOfRead -> do return ( put' bpile, ipile ) - - -- We just reversed /active/ inplicitly, which is no desaster, but may come - -- as a surprise downstream. So reverse it back. - upd_active reverse + p'feed_input + p'check_waiting + ((fin_bs, fin_bp), (fin_is, fin_ip)) <- p'scan_active -- Output, but don't bother emitting empty piles. Note that a plain -- basecall still yields an entry in the 'IndelPile'. This is necessary, -- because actual indel calling will want to know how many reads /did not/ -- show the variant. However, if no reads show any variant, and here is the -- first place where we notice that, the pile is useless. - let uninteresting (_,(d,i)) = d == 0 && null i - - unless (null fin_bp && all uninteresting fin_ip) - $ yield $ Pile rs po fin_bs (partitionPairEithers fin_bp) fin_is fin_ip + let uninteresting (_,(d,i)) = null d && null i + unless (null fin_bp && all uninteresting fin_ip) . yield $ + Pile rs po fin_bs (partitionPairEithers fin_bp) fin_is fin_ip -- Bump coordinate and loop. (Note that the bump to the next -- reference /sequence/ is done implicitly, because we will run out of @@ -478,6 +512,54 @@ upd_pos succ pileup' +-- | Feeds input as long as it starts at the current position +p'feed_input :: PileM m () +p'feed_input = do + rs <- get_refseq + po <- get_pos + + fix $ \loop -> peek >>= mapM_ (\(rs', po', prim) -> + when (rs == rs' && po == po') $ do + bump + case prim of Seek !p !pb -> ins_waiting p pb + Indel _ _ !pb -> add_active pb + EndOfRead -> return () + loop) + +-- | Checks /waiting/ queue. If there is anything waiting for the +-- current position, moves it to /active/ queue. +p'check_waiting :: PileM m () +p'check_waiting = do + po <- get_pos + fix $ \loop -> (viewMin <$> get_waiting) >>= mapM_ (\(!mk,!pb,w') -> + when (mk == po) $ do add_active pb + set_waiting w' + loop) + +-- | Scans /active/ queue and makes a 'BasePile'. Also sees what's next +-- in the 'PrimChunks': 'Indel's contribute to an 'IndelPile', 'Seek's +-- and deletions are pushed back to the /waiting/ queue, 'EndOfRead's +-- are removed, and everything else is added to the fresh /active/ +-- queue. +p'scan_active :: PileM m (( CallStats, [( Qual, Either DamagedBase DamagedBase )] ), + ( CallStats, [( Qual, ([Nucleotides], [DamagedBase]) )] )) +p'scan_active = + consume_active (mempty, mempty) $ + \(!bpile, !ipile) (Base wt qs mq str pchunks) -> + let put (Q !q) !x (!st,!vs) = ( st { read_depth = read_depth st + 1 + , reads_mapq0 = reads_mapq0 st + (if q == 0 then 1 else 0) + , sum_mapq = sum_mapq st + fromIntegral q + , sum_mapq_squared = sum_mapq_squared st + fromIntegral q * fromIntegral q } + , (Q q, x) : vs ) + b' = Base (wt-1) qs mq str pchunks + put' = put mq (if str then Left qs else Right qs) + in case pchunks of + _ | wt > 0 -> do add_active b' ; return ( bpile, ipile ) + Seek p' pb' -> do ins_waiting p' pb' ; return ( put' bpile, ipile ) + Indel del ins pb' -> do add_active pb' ; return ( put' bpile, put mq (del,ins) ipile ) + EndOfRead -> do return ( put' bpile, ipile ) + + partitionPairEithers :: [(a, Either b c)] -> ([(a,b)], [(a,c)]) partitionPairEithers = foldr either' ([],[]) where @@ -497,9 +579,6 @@ t1@(Node k1 x1 l1 r1) `union` t2@(Node k2 x2 l2 r2) | k1 <= k2 = Node k1 x1 (t2 `union` r1) l1 | otherwise = Node k2 x2 (t1 `union` r2) l2 - -insert :: Int -> PrimBase -> Heap -> Heap -insert k v heap = Node k v Empty Empty `union` heap getMinKey :: Heap -> Maybe Int getMinKey Empty = Nothing diff --git a/src/Bio/Bam/Reader.hs b/src/Bio/Bam/Reader.hs --- a/src/Bio/Bam/Reader.hs +++ b/src/Bio/Bam/Reader.hs @@ -45,7 +45,6 @@ import Control.Monad import Data.Attoparsec.ByteString ( anyWord8 ) import Data.Char ( digitToInt ) -import Data.Monoid import Data.Sequence ( (|>) ) import Data.String ( fromString ) import System.Environment ( getArgs ) diff --git a/src/Bio/Bam/Rec.hs b/src/Bio/Bam/Rec.hs --- a/src/Bio/Bam/Rec.hs +++ b/src/Bio/Bam/Rec.hs @@ -1,7 +1,5 @@ -{-# LANGUAGE OverloadedStrings, PatternGuards, BangPatterns #-} -{-# LANGUAGE NoMonomorphismRestriction, FlexibleContexts, FlexibleInstances #-} -{-# LANGUAGE RecordWildCards, TypeFamilies, MultiParamTypeClasses #-} -{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE RecordWildCards, BangPatterns, TypeFamilies, FlexibleContexts #-} +{-# LANGUAGE OverloadedStrings, FlexibleInstances, MultiParamTypeClasses #-} -- | Parsers and Printers for BAM and SAM. We employ an @Iteratee@ -- interface, and we strive to support everything possible in BAM. So @@ -55,14 +53,13 @@ isMerged, type_mask, - progressPos, + progressBam, Word32 ) where import Bio.Base import Bio.Bam.Header import Bio.Iteratee -import Bio.Util ( showNum ) import Control.Monad import Control.Monad.Primitive ( unsafePrimToPrim, unsafeInlineIO ) @@ -384,15 +381,6 @@ s' = if c `S.elem` s then s else c `S.cons` s -- | A simple progress indicator that prints sequence id and position. -progressPos :: MonadIO m => String -> (String -> IO ()) -> Refs -> Enumeratee [BamRaw] [BamRaw] m a -progressPos msg put refs = eneeCheckIfDonePass (icont . go 0) - where - go !_ k (EOF mx) = idone (liftI k) (EOF mx) - go !n k (Chunk [ ]) = liftI $ go n k - go !n k (Chunk as@(a:_)) = do let !n' = n + length as - when (n `div` 65536 /= n' `div` 65536) $ liftIO $ do - let BamRec{..} = unpackBam a - nm = unpackSeqid (sq_name (getRef refs b_rname)) ++ ":" - put $ "\27[K" ++ msg ++ nm ++ showNum b_pos ++ "\r" - eneeCheckIfDonePass (icont . go n') . k $ Chunk as +progressBam :: MonadIO m => String -> (String -> IO ()) -> Refs -> Enumeratee [BamRaw] [BamRaw] m a +progressBam = progressPos (\br -> case unpackBam br of b -> (b_rname b, b_pos b)) diff --git a/src/Bio/Bam/Trim.hs b/src/Bio/Bam/Trim.hs --- a/src/Bio/Bam/Trim.hs +++ b/src/Bio/Bam/Trim.hs @@ -90,4 +90,17 @@ trim_low_quality q = const $ all (< q) - +-- | Overlap-merging of read pairs. We shall compute the likelihood +-- for every possible overlap, the select the most likely one (unless it +-- looks completely random), compute a quality from the second best +-- merge, then merge and clamp the quality accordingly. Output should +-- be the pair *and* the merged representation, suitably flagged. +-- We might try looking for chimaera after completing the merge, if only +-- we knew which ones to expect. +-- +-- Likelihoods are straight forward; for adapters we have to assume a +-- realistic error rate (Q30 sounds good). To make it robust, we reduce +-- the adapters to their invariant prefix (about 20nt long), and try to +-- trim all adapters known to us (should be only two or three). +-- +-- Single-end reads are treated as pairs with an empty second read. diff --git a/src/Bio/Bam/Writer.hs b/src/Bio/Bam/Writer.hs --- a/src/Bio/Bam/Writer.hs +++ b/src/Bio/Bam/Writer.hs @@ -15,12 +15,10 @@ import Bio.Iteratee import Bio.Iteratee.Builder -import Control.Applicative import Data.ByteString.Builder ( toLazyByteString ) import Data.Bits import Data.Char ( ord, chr ) -import Data.Foldable ( foldMap ) -import Data.Monoid +import Data.Foldable ( foldMap ) import Foreign.Marshal.Alloc ( alloca ) import Foreign.Storable ( pokeByteOff, peek ) import System.IO diff --git a/src/Bio/Base.hs b/src/Bio/Base.hs --- a/src/Bio/Base.hs +++ b/src/Bio/Base.hs @@ -1,5 +1,5 @@ -{-# LANGUAGE GeneralizedNewtypeDeriving, TypeFamilies, FlexibleInstances #-} -{-# LANGUAGE MultiParamTypeClasses, BangPatterns, TemplateHaskell #-} +{-# LANGUAGE GeneralizedNewtypeDeriving, TypeFamilies, FlexibleInstances, CPP #-} +{-# LANGUAGE MultiParamTypeClasses, BangPatterns, TemplateHaskell, RankNTypes #-} -- | Common data types used everywhere. This module is a collection of -- very basic "bioinformatics" data types that are simple, but don't -- make sense to define over and over. @@ -7,7 +7,7 @@ module Bio.Base( Nucleotide(..), Nucleotides(..), Qual(..), toQual, fromQual, fromQualRaised, probToQual, - Prob(..), toProb, fromProb, qualToProb, pow, + Prob'(..), Prob, toProb, fromProb, qualToProb, pow, Word8, nucA, nucC, nucG, nucT, @@ -19,7 +19,6 @@ isProperBase, properBases, compl, compls, - everything, Seqid, unpackSeqid, @@ -39,18 +38,18 @@ w2c, c2w, - findAuxFile + findAuxFile, + module Data.Monoid ) where -import Bio.Util ( log1p ) -import Data.Array.Unboxed +import Bio.Util.Numeric ( log1p ) import Data.Bits import Data.ByteString.Internal ( c2w, w2c ) import Data.Char ( isAlpha, isSpace, ord, toUpper ) +import Data.Ix ( Ix ) +import Data.Monoid import Data.Word ( Word8 ) -import Data.Vector.Generic ( Vector(..) ) -import Data.Vector.Generic.Mutable ( MVector(..) ) -import Data.Vector.Unboxed.Deriving +import Data.Vector.Unboxed.Deriving ( derivingUnbox ) import Foreign.Storable ( Storable(..) ) import Numeric ( showFFloat ) import System.Directory ( doesFileExist ) @@ -58,10 +57,15 @@ import System.Environment ( getEnvironment ) import qualified Data.ByteString.Char8 as S - +import qualified Data.Vector.Unboxed as U --- | A nucleotide base. We only represent A,C,G,T. +#if __GLASGOW_HASKELL__ == 704 +import Data.Vector.Generic ( Vector(..) ) +import Data.Vector.Generic.Mutable ( MVector(..) ) +#endif +-- | A nucleotide base. We only represent A,C,G,T. The contained +-- 'Word8' ist guaranteed to be 0..3. newtype Nucleotide = N { unN :: Word8 } deriving ( Eq, Ord, Enum, Ix, Storable ) derivingUnbox "Nucleotide" [t| Nucleotide -> Word8 |] [| unN |] [| N |] @@ -70,17 +74,15 @@ minBound = N 0 maxBound = N 3 -everything :: (Bounded a, Ix a) => [a] -everything = range (minBound, maxBound) - -- | A nucleotide base in an alignment. -- Experience says we're dealing with Ns and gaps all the type, so -- purity be damned, they are included as if they were real bases. -- --- To allow @Nucleotides@s to be unpacked and incorparated into +-- To allow @Nucleotides@s to be unpacked and incorporated into -- containers, we choose to represent them the same way as the BAM file -- format: as a 4 bit wide field. Gaps are encoded as 0 where they --- make sense, N is 15. +-- make sense, N is 15. The contained 'Word8' is guaranteed to be +-- 0..15. newtype Nucleotides = Ns { unNs :: Word8 } deriving ( Eq, Ord, Enum, Ix, Storable ) @@ -98,6 +100,7 @@ -- directly on the \"Phred\" value, as the name suggests. The same goes -- for the 'Ord' instance: greater quality means higher \"Phred\" -- score, meand lower error probability. + newtype Qual = Q { unQ :: Word8 } deriving ( Eq, Ord, Storable, Bounded ) derivingUnbox "Qual" [t| Qual -> Word8 |] [| unQ |] [| Q |] @@ -114,18 +117,21 @@ fromQualRaised :: Double -> Qual -> Double fromQualRaised k (Q q) = 10 ** (- k * fromIntegral q / 10) --- | A positive 'Double' value stored in log domain. We store the +-- | A positive floating point value stored in log domain. We store the -- natural logarithm (makes computation easier), but allow conversions -- to the familiar \"Phred\" scale used for 'Qual' values. -newtype Prob = Pr { unPr :: Double } deriving ( Eq, Ord, Storable ) +newtype Prob' a = Pr { unPr :: a } deriving ( Eq, Ord, Storable ) -derivingUnbox "Prob" [t| Prob -> Double |] [| unPr |] [| Pr |] +-- | Common way of using 'Prob''. +type Prob = Prob' Double -instance Show Prob where +derivingUnbox "Prob'" [t| forall a . U.Unbox a => Prob' a -> a |] [| unPr |] [| Pr |] + +instance RealFloat a => Show (Prob' a) where showsPrec _ (Pr p) = (:) 'q' . showFFloat (Just 1) q where q = - 10 * p / log 10 -instance Num Prob where +instance (Floating a, Ord a) => Num (Prob' a) where fromInteger a = Pr (log (fromInteger a)) Pr x + Pr y = Pr $ if x >= y then x + log1p ( exp (y-x)) else y + log1p (exp (x-y)) Pr x - Pr y = Pr $ if x >= y then x + log1p (- exp (y-x)) else error "no negative error probabilities" @@ -134,26 +140,26 @@ abs x = x signum _ = Pr 0 -instance Fractional Prob where +instance (Floating a, Fractional a, Ord a) => Fractional (Prob' a) where fromRational a = Pr (log (fromRational a)) Pr a / Pr b = Pr (a - b) recip (Pr a) = Pr (negate a) infixr 8 `pow` -pow :: Prob -> Double -> Prob -pow (Pr a) e = Pr (a*e) +pow :: Num a => Prob' a -> a -> Prob' a +pow (Pr a) e = Pr $ a * e -toProb :: Double -> Prob +toProb :: Floating a => a -> Prob' a toProb p = Pr (log p) -fromProb :: Prob -> Double +fromProb :: Floating a => Prob' a -> a fromProb (Pr q) = exp q -qualToProb :: Qual -> Prob +qualToProb :: Floating a => Qual -> Prob' a qualToProb (Q q) = Pr (- log 10 * fromIntegral q / 10) -probToQual :: Prob -> Qual +probToQual :: (Floating a, RealFrac a) => Prob' a -> Qual probToQual (Pr p) = Q (round (- 10 * p / log 10)) nucA, nucC, nucG, nucT :: Nucleotide @@ -173,9 +179,8 @@ -- | Sequence identifiers are ASCII strings -- Since we tend to store them for a while, we use strict byte strings. --- If you get a lazy bytestring from somewhere, use 'shelve' to convert --- it for storage. Use @unpackSeqid@ and @packSeqid@ to avoid the --- import of @Data.ByteString@. +-- Use @unpackSeqid@ and @packSeqid@ to avoid the qualified import of +-- @Data.ByteString@. type Seqid = S.ByteString -- | Unpacks a @Seqid@ into a @String@ @@ -222,10 +227,9 @@ -- The usual codes for A,C,G,T and U are understood, '-' and '.' become -- gaps and everything else is an N. toNucleotide :: Char -> Nucleotide -toNucleotide c = if inRange (bounds arr) (ord c) then N (arr ! ord c) else N 0 +toNucleotide c = if ord c < 128 then N (arr `U.unsafeIndex` ord c) else N 0 where - arr :: UArray Int Word8 - arr = listArray (0,127) (repeat 0) // + arr = U.replicate 128 0 U.// ( [ (ord x, n) | (x, N n) <- pairs ] ++ [ (ord (toUpper x), n) | (x, N n) <- pairs ] ) @@ -236,10 +240,9 @@ -- The usual codes for A,C,G,T and U are understood, '-' and '.' become -- gaps and everything else is an N. toNucleotides :: Char -> Nucleotides -toNucleotides c = if inRange (bounds arr) (ord c) then Ns (arr ! ord c) else nucsN +toNucleotides c = if ord c < 128 then Ns (arr `U.unsafeIndex` ord c) else nucsN where - arr :: UArray Int Word8 - arr = listArray (0,127) (repeat (unNs nucsN)) // + arr = U.replicate 128 (unNs nucsN) U.// ( [ (ord x, n) | (x, Ns n) <- pairs ] ++ [ (ord (toUpper x), n) | (x, Ns n) <- pairs ] ) @@ -329,10 +332,9 @@ -- | Complements a Nucleotides. {-# INLINE compls #-} compls :: Nucleotides -> Nucleotides -compls (Ns x) = Ns $ arr ! (x .&. 15) +compls (Ns x) = Ns $ arr `U.unsafeIndex` fromIntegral (x .&. 15) where - arr :: UArray Word8 Word8 - !arr = listArray (0,15) [ 0, 8, 4, 12, 2, 10, 6, 14, 1, 9, 5, 13, 3, 11, 7, 15 ] + !arr = U.fromListN 16 [ 0, 8, 4, 12, 2, 10, 6, 14, 1, 9, 5, 13, 3, 11, 7, 15 ] -- | Moves a @Position@. The position is moved forward according to the diff --git a/src/Bio/Genocall.hs b/src/Bio/Genocall.hs --- a/src/Bio/Genocall.hs +++ b/src/Bio/Genocall.hs @@ -1,8 +1,6 @@ {-# LANGUAGE BangPatterns #-} module Bio.Genocall where -import Debug.Trace - import Bio.Bam.Pileup import Bio.Base import Bio.Genocall.Adna @@ -20,8 +18,8 @@ -- | Simple indel calling. We don't bother with it too much, so here's -- the gist: We collect variants (simply different variants, details --- don't matter), so @n@ variants give rise to (n+1)*n/2 GL values. --- (That's two out of @(n+1)@, the reference allele, represented here as +-- don't matter), so \(n\) variants give rise to \((n+1)*n/2\) GL values. +-- (That's two out of \((n+1)\), the reference allele, represented here as -- no deletion and no insertion, is there, too.) To assign these, we -- need a likelihood for an observed variant given an assumed genotype. -- @@ -31,36 +29,44 @@ -- though the real sequence is a different variant. For variants of -- different length, the likelihood is the map quality. This -- corresponds to the assumption that indel errors in sequencing are --- much less likely than mapping errors. Since this hardly our --- priority, the approximations are declared good enough. +-- much less likely than mapping errors. Since this is hardly our +-- priority, the approximations are hereby declared good enough. simple_indel_call :: Int -> IndelPile -> (GL, [IndelVariant]) -simple_indel_call ploidy vars = (simple_call ploidy mkpls vars, vars') +simple_indel_call _ [ ] = ( V.empty, [] ) +simple_indel_call _ [_] = ( V.empty, [] ) +simple_indel_call ploidy vars = ( simple_call ploidy mkpls vars, vars' ) where - vars' = Set.toList $ Set.fromList - [ IndelVariant d (V_Nuc $ V.fromList $ map db_call i) | (_q,(d,i)) <- vars ] + vars' = IndelVariant (V_Nucs V.empty) (V_Nuc V.empty) : + (Set.toList . Set.fromList) + [ IndelVariant (V_Nucs $ V.fromList d) + (V_Nuc $ V.fromList $ map db_call i) + | (_q,(d,i)) <- vars + , not (null d) || not (null i) ] - match = zipWith $ \(DB b q m) n -> let p = m ! n :-> b - p' = fromQual q - in toProb $ p + p' - p * p' + match = zipWith $ \(DB b q _ m) n -> let p = m ! n :-> b + p' = fromQual q + in toProb $ p + p' - p * p' - mkpls (q,(d,i)) = let !q' = qualToProb q - in [ if d /= dr || length i /= V.length ir - then q' else q' + product (match i $ V.toList ir) - | IndelVariant dr (V_Nuc ir) <- vars' ] + mkpls :: (Qual, ([Nucleotides], [DamagedBase])) -> [Prob] + mkpls (q,(d,i)) = [ qualToProb q + + if length d /= V.length dr || length i /= V.length ir + then 0 else product (match i $ V.toList ir) + | IndelVariant (V_Nucs dr) (V_Nuc ir) <- vars' ] -- | Naive SNP call; essentially the GATK model. We create a function -- that computes a likelihood for a given base, then hand over to simple -- call. Since everything is so straight forward, this works even in -- the face of damage. -simple_snp_call :: Int -> BasePile -> GL -simple_snp_call ploidy vars = simple_call ploidy mkpls vars +simple_snp_call :: (Qual -> Double) -> Int -> BasePile -> Snp_GLs +simple_snp_call from_qual ploidy vars = snp_gls (simple_call ploidy mkpls vars) ref where - mkpls (q, DB b qq m) = [ toProb $ x + pe*(s-x) | n <- [0..3], let x = m ! N n :-> b ] + 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 ] where - !p1 = fromQual q - !p2 = fromQual qq + !p1 = from_qual q + !p2 = from_qual qq !pe = p1 + p2 - p1*p2 !s = sum [ m ! N n :-> b | n <- [0..3] ] / 4 @@ -72,8 +78,9 @@ -- getting the current read, for every variant assuming that variant was -- sampled. -- --- NOTE, this may warrant specialization to diploidy and four alleles --- (common SNPs) and diploidy and two alleles (common indels). +-- XXX This eats up ~40% of total runtime; it *screams out* for +-- specialization to diploidy and four alleles (common SNPs) and maybe +-- diploidy and two alleles (common indels). simple_call :: Int -> (a -> [Prob]) -> [a] -> GL simple_call ploidy pls = foldl1' (V.zipWith (*)) . map step @@ -81,7 +88,7 @@ foldl1' _ [ ] = V.singleton 1 foldl1' f (a:as) = foldl' f a as - !mag = toProb (fromIntegral ploidy) `pow` (-1) + !mag = recip $ toProb (fromIntegral ploidy) -- XXX This could probably be simplified given the mk_pls function -- below. @@ -134,32 +141,36 @@ -- | SNP call according to maq/samtools/bsnp model. The matrix k counts -- how many errors we made, approximately. -maq_snp_call :: Int -> Double -> BasePile -> GL -maq_snp_call ploidy theta bases = V.fromList $ map l $ mk_snp_gts ploidy +maq_snp_call :: Int -> Double -> BasePile -> Snp_GLs +maq_snp_call ploidy theta bases = snp_gls (V.fromList $ map l $ mk_snp_gts ploidy) ref where -- Bases with effective qualities in order of decreasing(!) quality. -- A vector based algorithm may fit here. bases' = sortBy (flip $ comparing db_qual) [ db { db_qual = mq `min` db_qual db } | (mq,db) <- bases ] + ref = case bases of (_, DB _ _ r _) : _ -> r ; _ -> nucsN + everynuc :: Vec.Vec4 Nucleotide everynuc = nucA :. nucC :. nucG :. nucT :. () -- L(G) l gt = l' gt (toProb 1) (0 :: Mat44D) bases' - l' _ !acc _ [ ] = acc + l' !_ !acc !_ [ ] = acc l' !gt !acc !k (!x:xs) = 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 + + -- 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 -- P(H|X), again a vector of four p_x__q = dot p_x__q_h dg p_h__x = Vec.zipWith (\p p_h -> p / p_x__q * p_h) p_x__q_h dg - dg = (db_dmg x `multmv` gt) + dg = db_dmg x `multmv` gt kk = Vec.getElem (fromIntegral . unN $ db_call x) k + pack p_h__x k' = Vec.setElem (fromIntegral . unN $ db_call x) kk k @@ -201,22 +212,21 @@ show_indel (d, ins) = shows ins $ '-' : show d -} -{- showCall :: (a -> ShowS) -> VarCall (GL,a) -> ShowS -showCall f vc = shows (vc_refseq vc) . (:) ':' . - shows (vc_pos vc) . (:) '\t' . - f (snd $ vc_vars vc) . (++) "\tDP=" . - shows (vc_depth vc) . (++) ":MQ0=" . - shows (vc_mapq0 vc) . (++) ":MAPQ=" . - shows mapq . (:) '\t' . - show_pl (fst $ vc_vars vc) - where - show_pl :: Vector Prob -> ShowS - show_pl = (++) . intercalate "," . map show . V.toList +-- Error model with dependency parameter. Since both strands are +-- supposed to still be independent, we feed in only one pile, and +-- later combine both calls. XXX What's that doing HERE?! - mapq = vc_sum_mapq vc `div` vc_depth vc -} +type Calls = Pile' Snp_GLs (GL, [IndelVariant]) +-- | This pairs up GL values and the reference allele. When +-- constructing it, we make sure the GL values are in the correct order +-- if the reference allele is listed first. +data Snp_GLs = Snp_GLs !GL !Nucleotides + deriving Show --- | Error model with dependency parameter. Since both strands are --- supposed to still be independent, we feed in only one pile, and --- later combine both calls. XXX What's that doing HERE?! +snp_gls :: GL -> Nucleotides -> Snp_GLs +snp_gls pls ref | ref == nucsT = Snp_GLs (pls `V.backpermute` V.fromList [9,6,0,7,1,2,8,3,4,5]) ref + | ref == nucsG = Snp_GLs (pls `V.backpermute` V.fromList [5,3,0,4,1,2,8,6,7,9]) ref + | ref == nucsC = Snp_GLs (pls `V.backpermute` V.fromList [2,1,0,4,3,5,7,6,8,9]) ref + | otherwise = Snp_GLs pls ref diff --git a/src/Bio/Genocall/Adna.hs b/src/Bio/Genocall/Adna.hs --- a/src/Bio/Genocall/Adna.hs +++ b/src/Bio/Genocall/Adna.hs @@ -84,50 +84,15 @@ -- 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 ) - --- Forward strand first, C->T only; reverse strand next, G->A instead - -{- -{-# SPECIALIZE ssDamage :: SsDamageParameters Double -> DamageModel Double #-} -ssDamage :: Fractional a => SsDamageParameters a -> DamageModel a -ssDamage SSD{..} r l = V.generate l $ if r then ssd_rev else ssd_fwd where - ssd_fwd i = genSubstMat p 0 - where - !lam5 = ssd_lambda ^ (1+i) - !lam3 = ssd_kappa ^ (l-i) - !lam = lam3 + lam5 - lam3 * lam5 - !p = ssd_sigma * lam + ssd_delta * (1-lam) - - ssd_rev i = genSubstMat 0 p - where - !lam5 = ssd_lambda ^ (l-i) - !lam3 = ssd_kappa ^ (1+i) - !lam = lam3 + lam5 - lam3 * lam5 - !p = ssd_sigma * lam + ssd_delta * (1-lam) - - - -{-# SPECIALIZE dsDamage :: DsDamageParameters Double -> DamageModel Double #-} -dsDamage :: Fractional a => DsDamageParameters a -> DamageModel a -dsDamage DSD{..} _ l = V.generate l mat - where - mat i = genSubstMat p q - where - p = dsd_sigma * lam5 + dsd_delta * (1-lam5) - q = dsd_sigma * lam3 + dsd_delta * (1-lam3) - lam5 = dsd_lambda ^ (1+i) - lam3 = dsd_lambda ^ (l-i) --} - -{-# INLINE vec4 #-} -vec4 :: a -> a -> a -> a -> Vec4 a -vec4 a b c d = a :. b :. c :. d :. () + 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 diff --git a/src/Bio/Genocall/AvroFile.hs b/src/Bio/Genocall/AvroFile.hs --- a/src/Bio/Genocall/AvroFile.hs +++ b/src/Bio/Genocall/AvroFile.hs @@ -1,17 +1,25 @@ -{-# LANGUAGE TemplateHaskell, OverloadedStrings #-} +{-# LANGUAGE TemplateHaskell, OverloadedStrings, PatternGuards #-} module Bio.Genocall.AvroFile where import Bio.Base +import Bio.Bam.Header import Bio.Bam.Pileup +import Control.Applicative import Data.Aeson -import Data.Avro hiding ((.=)) +import Data.Avro import Data.Binary.Builder import Data.Binary.Get -import Data.Monoid +import Data.List ( intersperse ) +import Data.MiniFloat +import Data.Scientific ( toBoundedInteger ) +import Data.Text.Encoding ( encodeUtf8 ) import qualified Data.ByteString as B +import qualified Data.HashMap.Strict as H import qualified Data.Text as T +import qualified Data.Vector as V import qualified Data.Vector.Unboxed as U +import qualified Data.Sequence as Z -- ^ File format for genotype calls. @@ -23,23 +31,87 @@ -- the current one is getting too large. data GenoCallBlock = GenoCallBlock - { reference_name :: T.Text - , start_position :: Int + { reference_name :: {-# UNPACK #-} !Refseq + , start_position :: {-# UNPACK #-} !Int , called_sites :: [ GenoCallSite ] } + deriving (Show, Eq) data GenoCallSite = GenoCallSite - { snp_stats :: CallStats - , snp_likelihoods :: [ Int ] -- B.ByteString - , indel_stats :: CallStats + { snp_stats :: {-# UNPACK #-} !CallStats + -- snp likelihoods appear in the same order as in VCF, the reference + -- allele goes first if it is A, C, G or T. Else A goes first---not + -- my problem how to express that in VCF. + , snp_likelihoods :: {-# UNPACK #-} !(U.Vector Mini) -- B.ByteString? + , ref_allele :: {-# UNPACK #-} !Nucleotides + , indel_stats :: {-# UNPACK #-} !CallStats , indel_variants :: [ IndelVariant ] - , indel_likelihoods :: [ Int ] -- B.ByteString - } + , indel_likelihoods :: {-# UNPACK #-} !(U.Vector Mini) } -- B.ByteString? + deriving (Show, Eq) -$( deriveAvros [ ''GenoCallBlock, ''GenoCallSite, ''CallStats, ''IndelVariant ] ) +-- | Storing likelihoods: we take the natural logarithm (GL values are +-- already in a log scale) and convert to minifloat 0.4.4 +-- representation. Range and precision should be plenty. +compact_likelihoods :: U.Vector Prob -> U.Vector Mini -- B.ByteString +compact_likelihoods = U.map $ float2mini . negate . unPr +-- compact_likelihoods = map fromIntegral {- B.pack -} . U.toList . U.map (float2mini . negate . unPr) + +deriveAvros [ ''GenoCallBlock, ''GenoCallSite, ''CallStats, ''IndelVariant ] + instance Avro V_Nuc where - toSchema _ = return $ object [ "type" .= String "bytes", "doc" .= String "A,C,G,T" ] + toSchema _ = return $ object [ "type" .= String "bytes", "doc" .= String doc ] + where doc = T.pack $ intersperse ',' $ show $ [minBound .. maxBound :: Nucleotide] toBin (V_Nuc v) = encodeIntBase128 (U.length v) <> U.foldr ((<>) . singleton . unN) mempty v - fromBin = decodeIntBase128 >>= fmap (V_Nuc . U.fromList . map N . B.unpack) . getByteString + fromBin = decodeIntBase128 >>= \l -> V_Nuc . U.fromListN l . map N . B.unpack <$> getByteString l toAvron (V_Nuc v) = String . T.pack . map w2c . U.toList $ U.map unN v + +instance Avro V_Nucs where + toSchema _ = return $ object [ "type" .= String "bytes", "doc" .= String doc ] + where doc = T.pack $ intersperse ',' $ show $ [minBound .. maxBound :: Nucleotides] + toBin (V_Nucs v) = encodeIntBase128 (U.length v) <> U.foldr ((<>) . singleton . unNs) mempty v + fromBin = decodeIntBase128 >>= \l -> V_Nucs . U.fromListN l . map Ns . B.unpack <$> getByteString l + toAvron (V_Nucs v) = String . T.pack . map w2c . U.toList $ U.map unNs v + +instance Avro Nucleotides where + toSchema _ = return $ String "int" + toBin = encodeIntBase128 . unNs + fromBin = Ns <$> decodeIntBase128 + toAvron = Number . fromIntegral . unNs + +instance Avro Mini where + toSchema _ = return $ String "int" + toBin = encodeIntBase128 . unMini + fromBin = Mini <$> decodeIntBase128 + toAvron = Number . fromIntegral . unMini + +-- | We encode the Refseq as an Avro enum, which serves as a kind of +-- symbol table. To make this work, the environment of the 'MkSchema' +-- monad has to be prepopulated with a suitable schema. +instance Avro Refseq where + toSchema _ = getNamedSchema "Refseq" + toBin = encodeIntBase128 . unRefseq + fromBin = Refseq <$> decodeIntBase128 + + -- This is cheating, we should use the enum names, but they are not + -- available. Doesn't matter, this is mostly for debugging anyway. + toAvron = Number . fromIntegral . unRefseq + + +-- | Reconstructs the list of reference sequences from Avro metadata. +-- If a type named @Refseq@ is defined in the schema and is an enum, it +-- defines the symbol table, otherwise an empty list is returned. If +-- @biohazard.refseq_length@ exists, and is an array, it's elements are +-- interpreted as the lengths in order, otherwise the lengths are set to +-- zero. +getRefseqs :: AvroMeta -> Refs +getRefseqs meta + | Object o <- findSchema "Refseq" meta + , Just (String "enum") <- H.lookup "type" o + , Just (Array syms) <- H.lookup "symbols" o + = Z.fromList [ BamSQ (encodeUtf8 nm) ln [] | (String nm, ln) <- V.toList syms `zip` lengths ] + | otherwise = Z.empty + where + lengths = case decodeStrict =<< H.lookup "biohazard.refseq_length" meta of + Just (Array lns) -> [ case l of Number n -> maybe 0 id $ toBoundedInteger n ; _ -> 0 | l <- V.toList lns ] + _ -> repeat 0 diff --git a/src/Bio/Genocall/Metadata.hs b/src/Bio/Genocall/Metadata.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Genocall/Metadata.hs @@ -0,0 +1,154 @@ +{-# LANGUAGE OverloadedStrings, RecordWildCards, FlexibleContexts, BangPatterns #-} +-- | Metadata necessary for a sensible genotyping workflow. +module Bio.Genocall.Metadata where + +import Bio.Genocall.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.ByteString.Unsafe ( unsafeUseAsCStringLen ) +import Data.Monoid +import Data.Vector.Unboxed ( Vector ) +import Foreign.Ptr ( castPtr ) +import GHC.IO.Exception ( IOErrorType(..) ) +import Prelude hiding ( writeFile, readFile ) +import System.IO.Error ( isAlreadyExistsErrorType, ioeGetErrorType ) +import System.Posix.Files ( rename, removeLink ) +import System.Posix.IO + +import qualified Data.HashMap.Strict as M + +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 + } deriving Show + +data Library = Library { + library_name :: Text, + library_files :: [Text], + library_damage :: Maybe (DamageParameters Double) + } deriving Show + +-- | Divergence estimate. Lists contain three or four floats, these are +-- divergence, heterozygosity at W sites, heterozygosity at S sites, and +-- optionally gappiness in this order. +data DivEst = DivEst { + point_est :: [Double], + conf_region :: [( [Double], [Double] )] + } deriving Show + + +type Metadata = HashMap Text Sample + +instance ToJSON DivEst where + toJSON DivEst{..} = object $ [ "estimate" .= point_est + , "confidence-region" .= conf_region ] + +instance FromJSON DivEst where + parseJSON (Object o) = DivEst <$> o .: "estimate" <*> o .:? "confidence-region" .!= [] + parseJSON (Array a) = flip DivEst [] <$> parseJSON (Array a) + parseJSON _ = fail $ "divergence estimate should be an array or an object" + +instance ToJSON float => ToJSON (DamageParameters float) where + toJSON DP{..} = object [ "ss-sigma" .= ssd_sigma + , "ss-delta" .= ssd_delta + , "ss-lambda" .= ssd_lambda + , "ss-kappa" .= ssd_kappa + , "ds-sigma" .= dsd_sigma + , "ds-delta" .= dsd_delta + , "ds-lambda" .= dsd_lambda ] + +instance FromJSON float => FromJSON (DamageParameters float) where + parseJSON = withObject "damage parameters" $ \o -> + DP <$> o .: "ss-sigma" + <*> o .: "ss-delta" + <*> o .: "ss-lambda" + <*> o .: "ss-kappa" + <*> o .: "ds-sigma" + <*> o .: "ds-delta" + <*> o .: "ds-lambda" + +instance ToJSON Library where + toJSON (Library name files dp) = object ( maybe id ((:) . ("damage" .=)) dp + $ [ "name" .= name, "files" .= files ] ) + +instance FromJSON Library where + parseJSON (String name) = return $ Library name [name <> ".bam"] Nothing + parseJSON (Object o) = Library <$> o .: "name" + <*> (maybe id (:) <$> o .:? "file" + <*> o .:? "files" .!= []) + <*> o .:? "damage" + parseJSON _ = fail "String or Object expected for library" + +instance ToJSON Sample where + toJSON (Sample ls avfs bcfs dts ds) = object $ hashToJson "divergences" ds $ + listToJson "libraries" ls $ + hashToJson "avro-files" avfs $ + hashToJson "bcf-files" bcfs $ + hashToJson "div-tables" dts [] + where + hashToJson k vs = if M.null vs then id else (:) (k .= vs) + listToJson k vs = if null vs then id else (:) (k .= vs) + +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) + parseJSON (Object o) = Sample <$> o .: "libraries" + <*> (M.singleton "" <$> o .: "avro-file" <|> o .:? "avro-files" .!= M.empty) + <*> (M.singleton "" <$> o .: "bcf-file" <|> o .:? "bcf-files" .!= M.empty) + <*> o .:? "div-tables" .!= M.empty + <*> (M.singleton "" <$> o .: "divergence" <|> o.:? "divergences" .!= M.empty) + parseJSON _ = fail $ "String, Array or Object expected for Sample" + + +-- | 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) + 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) + +-- | Update the configuration file. Open a new file (fn++"~new") in +-- exclusive mode. Then read the old file, write the update to the new +-- file, rename it atomically, then close it. Use of O_EXCL should +-- ensure that nobody interferes. This is atomic even on NFS, provided +-- NFS and kernel are new enough. For older NFS, I cannot be bothered. +-- +-- (The first idea was to base this on the supposed fact that link(2) is +-- atomic and fails if the new filename exists. This approach does seem +-- to contain a race condition, though.) +updateMetadata :: (Metadata -> Metadata) -> FilePath -> IO () +updateMetadata f fp = go (360::Int) -- retry every 5 secs for 30 minutes + where + fpn = fp <> "~new" + + go !n = handleJust + (\e -> if isAlreadyExistsErrorType (ioeGetErrorType e) && n > 0 then Just () else Nothing) + (\_ -> threadDelay 5000000 >> go (n-1)) $ do + bracket + (openFd fpn WriteOnly (Just 0o666) defaultFileFlags{ exclusive = True }) + (closeFd) $ \fd -> + (do mdata <- readMetadata fp + forM_ (toChunks . encode . toJSON $ f mdata) $ \ch -> + unsafeUseAsCStringLen ch $ \(p,l) -> + fdWriteBuf fd (castPtr p) (fromIntegral l) + rename fpn fp) + `onException` removeLink fpn + + +split_sam_rgns :: Metadata -> [String] -> [( String, [Maybe String] )] +split_sam_rgns _meta [ ] = [] +split_sam_rgns meta (s:ss) = (s, if null rgns then [Nothing] else map Just rgns) : split_sam_rgns meta rest + where (rgns, rest) = break (\x -> pack x `M.member` meta) ss diff --git a/src/Bio/Glf.hs b/src/Bio/Glf.hs deleted file mode 100644 --- a/src/Bio/Glf.hs +++ /dev/null @@ -1,133 +0,0 @@ -{-# LANGUAGE FlexibleContexts #-} -module Bio.Glf ( - GlfSeq(..), - GlfRec(..), - enee_glf_file, - enum_glf_file, - enum_glf_handle - ) where - -import Bio.Iteratee -import Bio.Iteratee.Bgzf -import Control.Monad -import Data.Bits -import System.IO - -import qualified Data.ByteString.Char8 as S -import qualified Data.Iteratee.ListLike as I - - -data GlfRec = SNP { glf_refbase :: {-# UNPACK #-} !Char - , glf_offset :: {-# UNPACK #-} !Int - , glf_depth :: {-# UNPACK #-} !Int - , glf_min_lk :: {-# UNPACK #-} !Int - , glf_mapq :: {-# UNPACK #-} !Int - , glf_lk :: [Int] } - | Indel { glf_refbase :: {-# UNPACK #-} !Char - , glf_offset :: {-# UNPACK #-} !Int - , glf_depth :: {-# UNPACK #-} !Int - , glf_min_lk :: {-# UNPACK #-} !Int - , glf_mapq :: {-# UNPACK #-} !Int - , glf_lk_hom1 :: {-# UNPACK #-} !Int - , glf_lk_hom2 :: {-# UNPACK #-} !Int - , glf_lk_het :: {-# UNPACK #-} !Int - , glf_is_ins1 :: !Bool - , glf_is_ins2 :: !Bool - , glf_seq1 :: {-# UNPACK #-} !S.ByteString - , glf_seq2 :: {-# UNPACK #-} !S.ByteString } - deriving Show - -data GlfSeq = GlfSeq { glf_seqname :: {-# UNPACK #-} !S.ByteString - , glf_seqlen :: {-# UNPACK #-} !Int } - deriving Show - - -enee_glf_recs :: Monad m => Enumeratee S.ByteString [GlfRec] m b -enee_glf_recs = eneeCheckIfDone step - where - step oit' = I.isFinished >>= step' oit' - - step' oit' True = return $ liftI oit' - step' oit' False = do - type_ref <- I.head - let refbase = "XACMGRSVTWYHKDBN" !! fromIntegral (type_ref .&. 0xf) - case type_ref `shiftR` 4 of - 0 -> return $ oit' $ EOF Nothing - 1 -> do r <- get_snp $ get_common (SNP refbase) - eneeCheckIfDone step . oit' $ Chunk [r] - 2 -> do r <- get_indel $ get_common (Indel refbase) - eneeCheckIfDone step . oit' $ Chunk [r] - x -> fail $ "unknown GLF record #" ++ show x - - get_common f = return f - `ap` (fromIntegral `liftM` endianRead4 LSB) - `ap` (fromIntegral `liftM` endianRead3 LSB) - `ap` (fromIntegral `liftM` I.head) - `ap` (fromIntegral `liftM` I.head) - - get_snp f = f `ap` get_lk_arr - get_lk_arr = replicateM 10 (fromIntegral `liftM` I.head) - - get_indel f = do - f' <- f `ap` (fromIntegral `liftM` I.head) - `ap` (fromIntegral `liftM` I.head) - `ap` (fromIntegral `liftM` I.head) - l1 <- getInt16le - l2 <- getInt16le - liftM2 (f' (l1 >= 0) (l2 >= 0)) (iGetString (abs l1)) (iGetString (abs l2)) - - getInt16le = do i <- endianRead2 LSB - return $ if i > 0x7fff then fromIntegral i - 0x10000 - else fromIntegral i - -enee_glf_seq :: Monad m => (GlfSeq -> Enumeratee [GlfRec] a m b) -> Enumeratee S.ByteString a m b -enee_glf_seq per_seq oit = do l <- endianRead4 LSB - s <- liftM2 GlfSeq (S.init `liftM` iGetString (fromIntegral l)) - (fromIntegral `liftM` endianRead4 LSB) - enee_glf_recs ><> per_seq s $ oit - --- | Iterates over a GLF file. In @get_glf_file per_seq per_file@, the --- enumerator @per_file genome_name@, where @genome_name@ is the name --- stored in the GLF header, is run once, then the enumeratee @per_seq --- glfseq@ is iterated over the records in each sequence. -enee_glf_file :: Monad m => (GlfSeq -> Enumeratee [GlfRec] a m b) - -> (S.ByteString -> Enumerator a m b) - -> Enumeratee S.ByteString a m b -enee_glf_file per_seq per_file oit = do - matched <- I.heads (S.pack "GLF\003") - when (matched /= 4) (fail "GLF signature not found") - nm <- endianRead4 LSB >>= iGetString . fromIntegral - lift (per_file nm oit) >>= loop - where - -- loop :: Monad m => Iteratee a m b -> Iteratee S.ByteString m (Iteratee a m b) - loop it = I.isFinished >>= loop' it - loop' it True = return it - loop' it False = loop =<< enee_glf_seq per_seq it - - --- | Enumerate the contents of a GLF file, apply suitable Enumeratees to --- both sequences and records, resulting in an Enumerator of /whatever/, --- typically output Strings or records... --- --- This type is positively weird and I'm not entirely sure this is the --- right way to go about it. -enum_glf_file :: (MonadIO m, MonadMask m) - => FilePath - -> (GlfSeq -> Enumeratee [GlfRec] a m b) - -> (S.ByteString -> Enumerator a m b) - -> Enumerator a m b -enum_glf_file fp per_seq per_file output = - enumFile defaultBufSize fp >=> run $ - joinI $ decompressBgzf $ - enee_glf_file per_seq per_file output - -enum_glf_handle :: (MonadIO m, MonadMask m) - => Handle - -> (GlfSeq -> Enumeratee [GlfRec] a m b) - -> (S.ByteString -> Enumerator a m b) - -> Enumerator a m b -enum_glf_handle hdl per_seq per_file output = - enumHandle defaultBufSize hdl >=> run $ - joinI $ decompressBgzf $ - enee_glf_file per_seq per_file output - diff --git a/src/Bio/Iteratee.hs b/src/Bio/Iteratee.hs --- a/src/Bio/Iteratee.hs +++ b/src/Bio/Iteratee.hs @@ -7,6 +7,8 @@ groupStreamBy, groupStreamOn, iGetString, + iterGet, + iterLoop, iLookAhead, headStream, peekStream, @@ -25,6 +27,7 @@ mapMaybeStream, parMapChunksIO, progressNum, + progressPos, I.mapStream, I.takeWhileE, @@ -67,19 +70,26 @@ Fd, withFileFd, - module X ) where + module Data.Iteratee.Binary, + module Data.Iteratee.Char, + module Data.Iteratee.IO, + module Data.Iteratee.Iteratee + ) where import Bio.Base ( findAuxFile ) -import Bio.Util ( showNum ) +import Bio.Bam.Header +import Bio.Util.Numeric ( showNum ) import Control.Concurrent.Async ( Async, async, wait, cancel ) import Control.Monad import Control.Monad.Catch import Control.Monad.IO.Class import Control.Monad.Trans.Class -import Data.Iteratee.Binary as X -import Data.Iteratee.Char as X -import Data.Iteratee.IO as X hiding ( defaultBufSize ) -import Data.Iteratee.Iteratee as X hiding ( identity ) +import Data.Binary.Get +import Data.Bits ( shiftR ) +import Data.Iteratee.Binary +import Data.Iteratee.Char +import Data.Iteratee.IO hiding ( defaultBufSize ) +import Data.Iteratee.Iteratee hiding ( identity ) import Data.ListLike ( ListLike ) import Data.Monoid import Data.Typeable @@ -89,7 +99,7 @@ import System.Posix ( Fd, openFd, closeFd, OpenMode(..), defaultFileFlags ) import qualified Data.Attoparsec.ByteString as A -import qualified Data.ByteString as S +import qualified Data.ByteString.Char8 as S import qualified Data.Iteratee as I import qualified Data.ListLike as LL import qualified Data.Vector.Generic as VG @@ -204,6 +214,28 @@ in idone r (Chunk $ S.drop (n-l) c) | otherwise = liftI $ step (c:acc) (l + S.length c) +-- | Repeatedly apply an 'Iteratee' to a value until end of stream. +-- Returns the final value. +iterLoop :: (Nullable s, Monad m) => (a -> Iteratee s m a) -> a -> Iteratee s m a +iterLoop it a = do e <- I.isFinished + if e then return a + else it a >>= iterLoop it + +-- | Convert a 'Get' into an 'Iteratee'. The 'Get' is applied once, the +-- decoded data is returned, unneded input remains in the stream. +iterGet :: Monad m => Get a -> Iteratee S.ByteString m a +iterGet = go . runGetIncremental + where + go (Fail _ _ err) = throwErr (iterStrExc err) + go (Done rest _ a) = idone a (Chunk rest) + go (Partial dec) = liftI $ \ck -> case ck of + Chunk s -> go (dec $ Just s) + EOF mx -> case dec Nothing of + Fail _ _ err -> throwErr (iterStrExc err) + Partial _ -> throwErr (iterStrExc "") + Done rest _ a | S.null rest -> idone a (EOF mx) + | otherwise -> idone a (Chunk rest) + {-# INLINE mBind #-} -- | Lifts a monadic action and combines it with a continuation. -- @mBind m f@ is the same as @lift m >>= f@, but does not require a @@ -327,8 +359,8 @@ foldStream :: (Monad m, Nullable s, ListLike s a) => (b -> a -> b) -> b -> Iteratee s m b foldStream f = foldChunksM (\b s -> return $! LL.foldl' f b s) - -zipStreams :: (Monad m, Nullable s, ListLike s e) +-- | Apply two 'Iteratee's to the same stream. +zipStreams :: (Nullable s, ListLike s el, Monad m) => Iteratee s m a -> Iteratee s m b -> Iteratee s m (a, b) zipStreams = I.zip @@ -420,6 +452,19 @@ put $ "\27[K" ++ msg ++ showNum n ++ "\r" eneeCheckIfDonePass (icont . go n') . k $ Chunk as +-- | A simple progress indicator that prints a position. +progressPos :: (MonadIO m, ListLike s a, NullPoint s) + => (a -> (Refseq, Int)) -> String -> (String -> IO ()) -> Refs -> Enumeratee s s m b +progressPos f msg put refs = eneeCheckIfDonePass (icont . go invalidRefseq 0) + where + go !_ !_ k (EOF mx) = idone (liftI k) (EOF mx) + go !rs0 !po0 k c@(Chunk as) + | LL.null as = liftI $ go rs0 po0 k + | otherwise = when (rs1 /= rs0 || po1 `shiftR` 19 /= po0 `shiftR` 19) + (let nm = S.unpack (sq_name (getRef refs rs1)) ++ ":" + in put $ "\27[K" ++ msg ++ nm ++ showNum po1 ++ "\r") + `ioBind_` eneeCheckIfDonePass (icont . go rs1 po1) (k c) + where (!rs1, !po1) = f (LL.head as) -- A very simple queue data type. -- Invariants: q = QQ l f b --> l == length f + length b diff --git a/src/Bio/Iteratee/Builder.hs b/src/Bio/Iteratee/Builder.hs --- a/src/Bio/Iteratee/Builder.hs +++ b/src/Bio/Iteratee/Builder.hs @@ -12,26 +12,25 @@ module Bio.Iteratee.Builder where -import Control.Monad -import Control.Monad.IO.Class +import Bio.Iteratee +import Bio.Iteratee.Bgzf import Data.Bits import Data.Monoid import Data.Primitive.Addr import Data.Primitive.ByteArray +import Data.Word ( Word8, Word16, Word32 ) +import Foreign.Marshal.Alloc +import Foreign.Marshal.Utils +import Foreign.Ptr +import Foreign.Storable ( peek, poke ) import GHC.Exts -import GHC.Word ( Word8, Word16, Word32 ) +import System.IO.Unsafe ( unsafePerformIO ) import qualified Data.ByteString as B import qualified Data.ByteString.Unsafe as B import qualified Data.ByteString.Builder as B ( Builder, toLazyByteString ) import qualified Data.ByteString.Lazy as B ( foldrChunks ) -import Bio.Iteratee -import Bio.Iteratee.Bgzf - -import Foreign.Marshal.Utils -import Foreign.Ptr - -- | The 'MutableByteArray' is garbage collected, so we don't get leaks. -- Once it has grown to a practical size (and the initial 128k should be -- very practical), we don't get fragmentation either. We also avoid @@ -39,7 +38,8 @@ -- lazy or strict have to be allocated. data BB = BB { buffer :: {-# UNPACK #-} !(MutableByteArray RealWorld) , len :: {-# UNPACK #-} !Int - , mark :: {-# UNPACK #-} !Int } + , mark :: {-# UNPACK #-} !Int + , mark2 :: {-# UNPACK #-} !Int } -- This still seems to have considerable overhead. Don't know if this -- can be improved by effectively inlining IO and turning the BB into an @@ -58,7 +58,7 @@ -- | Creates a buffer with initial capacity of ~128k. newBuffer :: IO BB -newBuffer = newPinnedByteArray 128000 >>= \arr -> return $ BB arr 0 0 +newBuffer = newPinnedByteArray 128000 >>= \arr -> return $ BB arr 0 0 0 -- | Ensures a given free space in the buffer by doubling its capacity -- if necessary. @@ -118,6 +118,17 @@ pushByteString :: B.ByteString -> Push pushByteString bs = ensureBuffer (B.length bs) <> unsafePushByteString bs +{-# INLINE unsafePushFloat #-} +unsafePushFloat :: Float -> Push +unsafePushFloat f = unsafePushWord32 i + where + i :: Word32 + i = unsafePerformIO $ alloca $ \b -> poke (castPtr b) f >> peek b + +{-# INLINE pushFloat #-} +pushFloat :: Float -> Push +pushFloat f = ensureBuffer 4 <> unsafePushFloat f + {-# INLINE pushBuilder #-} pushBuilder :: B.Builder -> Push pushBuilder = B.foldrChunks ((<>) . pushByteString) mempty . B.toLazyByteString @@ -139,6 +150,36 @@ endRecord :: Push endRecord = Push $ \b -> do let !l = len b - mark b - 4 + writeByteArray (buffer b) (mark b + 0) (fromIntegral $ shiftR l 0 :: Word8) + writeByteArray (buffer b) (mark b + 1) (fromIntegral $ shiftR l 8 :: Word8) + writeByteArray (buffer b) (mark b + 2) (fromIntegral $ shiftR l 16 :: Word8) + writeByteArray (buffer b) (mark b + 3) (fromIntegral $ shiftR l 24 :: Word8) + return b + +-- | Ends the first part of a record. The length is filled in *before* +-- the mark, which is specifically done to support the *two* length +-- fields in BCF. It also remembers the current position. Horrible +-- things happen if this isn't preceeded by *two* succesive invocations +-- of 'setMark'. +{-# INLINE endRecordPart1 #-} +endRecordPart1 :: Push +endRecordPart1 = Push $ \b -> do + let !l = len b - mark b - 4 + writeByteArray (buffer b) (mark b - 4) (fromIntegral $ shiftR l 0 :: Word8) + writeByteArray (buffer b) (mark b - 3) (fromIntegral $ shiftR l 8 :: Word8) + writeByteArray (buffer b) (mark b - 2) (fromIntegral $ shiftR l 16 :: Word8) + writeByteArray (buffer b) (mark b - 1) (fromIntegral $ shiftR l 24 :: Word8) + return $ b { mark2 = len b } + +-- | Ends the second part of a record. The length is filled in at the +-- mark, but computed from the sencond mark only. This is specifically +-- done to support the *two* length fields in BCF. Horrible things +-- happen if this isn't preceeded by *two* succesive invocations of +-- 'setMark' and one of 'endRecordPart1'. +{-# INLINE endRecordPart2 #-} +endRecordPart2 :: Push +endRecordPart2 = Push $ \b -> do + let !l = len b - mark2 b writeByteArray (buffer b) (mark b + 0) (fromIntegral $ shiftR l 0 :: Word8) writeByteArray (buffer b) (mark b + 1) (fromIntegral $ shiftR l 8 :: Word8) writeByteArray (buffer b) (mark b + 2) (fromIntegral $ shiftR l 16 :: Word8) diff --git a/src/Bio/TwoBit.hs b/src/Bio/TwoBit.hs --- a/src/Bio/TwoBit.hs +++ b/src/Bio/TwoBit.hs @@ -1,20 +1,25 @@ -{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE BangPatterns, NamedFieldPuns, RecordWildCards #-} module Bio.TwoBit ( module Bio.Base, - TwoBitFile, + TwoBitFile(..), + TwoBitSequence(..), openTwoBit, + getFwdSubseqWith, getSubseq, getSubseqWith, getSubseqAscii, getSubseqMasked, + getLazySubseq, getSeqnames, - hasSequence, + lookupSequence, getSeqLength, clampPosition, getRandomSeq, + takeOverlap, + mergeBlocks, Mask(..) ) where @@ -26,8 +31,10 @@ import qualified Data.ByteString as B import qualified Data.ByteString.Lazy as L import Data.Char (toLower) +-- can't use Data.IntMap.Strict to remain compatible with +-- containers-0.4.1 and therefore ghc 7.4 :-( import qualified Data.IntMap as I -import qualified Data.Map as M +import qualified Data.HashMap.Lazy as M import Data.Maybe import Numeric import System.IO.Posix.MMap @@ -44,17 +51,18 @@ -- genome), which can be interpreted in whatever way fits. And that's why -- we have 'Mask' and 'getSubseqWith'. -- --- TODO: use binary search for the Int->Int mappings? +-- TODO: use binary search for the Int->Int mappings on the raw data? data TwoBitFile = TBF { tbf_raw :: B.ByteString, - tbf_seqs :: !(M.Map Seqid TwoBitSequence) + -- This map is intentionally lazy. May or may not be important. + tbf_seqs :: !(M.HashMap Seqid TwoBitSequence) } -data TwoBitSequence = Indexed { tbs_n_blocks :: !(I.IntMap Int) - , tbs_m_blocks :: !(I.IntMap Int) - , tbs_dna_offset :: {-# UNPACK #-} !Int - , tbs_dna_size :: {-# UNPACK #-} !Int } +data TwoBitSequence = TBS { tbs_n_blocks :: !(I.IntMap Int) + , tbs_m_blocks :: !(I.IntMap Int) + , tbs_dna_offset :: {-# UNPACK #-} !Int + , tbs_dna_size :: {-# UNPACK #-} !Int } -- | Brings a 2bit file into memory. The file is mmap'ed, so it will -- not work on streams that are not actual files. It's also unsafe if @@ -87,11 +95,12 @@ nb <- readBlockList mb <- readBlockList len <- getWord32 >> bytesRead - return $! Indexed (I.fromList nb) (I.fromList mb) (ofs + fromIntegral len) ds + return $! TBS (I.fromList nb) (I.fromList mb) (ofs + fromIntegral len) ds readBlockList = getWord32 >>= \n -> liftM2 zip (repM n getWord32) (repM n getWord32) --- | Repeat monadic action 'n' times. Returns result in reverse(!) order. +-- | Repeat monadic action 'n' times. Returns result in reverse(!) +-- order, but doesn't build a huge list of thunks in memory. repM :: Monad m => Int -> m a -> m [a] repM n0 m = go [] n0 where @@ -109,16 +118,14 @@ data Mask = None | Soft | Hard | Both deriving (Eq, Ord, Enum, Show) -getFwdSubseqWith :: B.ByteString -> Int -- raw data, dna offset - -> I.IntMap Int -> I.IntMap Int -- N blocks, M blocks +getFwdSubseqWith :: TwoBitFile -> TwoBitSequence -- raw data, sequence -> (Word8 -> Mask -> a) -- mask function - -> Int -> Int -> [a] -- start, len, result -getFwdSubseqWith raw ofs n_blocks m_blocks nt start len = - do_mask (takeOverlap start n_blocks `mergeblocks` takeOverlap start m_blocks) start . - take len . drop (start .&. 3) . + -> Int -> [a] -- start, lazy result +getFwdSubseqWith TBF{..} TBS{..} nt start = + do_mask (takeOverlap start tbs_n_blocks `mergeBlocks` takeOverlap start tbs_m_blocks) start . + drop (start .&. 3) . B.foldr toDNA [] . - B.take (len `shiftR` 2 + 2) . -- needed?! - B.drop (fromIntegral $ ofs + (start `shiftR` 2)) $ raw + B.drop (fromIntegral $ tbs_dna_offset + (start `shiftR` 2)) $ tbf_raw where toDNA b = (++) [ 3 .&. (b `shiftR` x) | x <- [6,4,2,0] ] @@ -130,19 +137,19 @@ -- | Merge blocks of Ns and blocks of Ms into single list of blocks with -- masking annotation. Gaps remain. Used internally only. -mergeblocks :: [(Int,Int)] -> [(Int,Int)] -> [(Int,Int,Mask)] -mergeblocks ((_,0):nbs) mbs = mergeblocks nbs mbs -mergeblocks nbs ((_,0):mbs) = mergeblocks nbs mbs +mergeBlocks :: [(Int,Int)] -> [(Int,Int)] -> [(Int,Int,Mask)] +mergeBlocks ((_,0):nbs) mbs = mergeBlocks nbs mbs +mergeBlocks nbs ((_,0):mbs) = mergeBlocks nbs mbs -mergeblocks ((ns,nl):nbs) ((ms,ml):mbs) - | ns < ms = let l = min (ms-ns) nl in (ns,l, Hard) : mergeblocks ((ns+l,nl-l):nbs) ((ms,ml):mbs) - | ms < ns = let l = min (ns-ms) ml in (ms,l, Soft) : mergeblocks ((ns,nl):nbs) ((ms+l,ml-l):mbs) - | otherwise = let l = min nl ml in (ns,l, Both) : mergeblocks ((ns+l,nl-l):nbs) ((ms+l,ml-l):mbs) +mergeBlocks ((ns,nl):nbs) ((ms,ml):mbs) + | ns < ms = let l = min (ms-ns) nl in (ns,l, Hard) : mergeBlocks ((ns+l,nl-l):nbs) ((ms,ml):mbs) + | ms < ns = let l = min (ns-ms) ml in (ms,l, Soft) : mergeBlocks ((ns,nl):nbs) ((ms+l,ml-l):mbs) + | otherwise = let l = min nl ml in (ns,l, Both) : mergeBlocks ((ns+l,nl-l):nbs) ((ms+l,ml-l):mbs) -mergeblocks ((ns,nl):nbs) [] = (ns,nl, Hard) : mergeblocks nbs [] -mergeblocks [] ((ms,ml):mbs) = (ms,ml, Soft) : mergeblocks [] mbs +mergeBlocks ((ns,nl):nbs) [] = (ns,nl, Hard) : mergeBlocks nbs [] +mergeBlocks [] ((ms,ml):mbs) = (ms,ml, Soft) : mergeBlocks [] mbs -mergeblocks [ ] [ ] = [] +mergeBlocks [ ] [ ] = [] -- | Extract a subsequence and apply masking. TwoBit file can represent @@ -153,15 +160,27 @@ getSubseqWith :: (Nucleotide -> Mask -> a) -> TwoBitFile -> Range -> [a] getSubseqWith maskf tbf (Range { r_pos = Pos { p_seq = chr, p_start = start }, r_length = len }) = do let sq1 = maybe (error $ unpackSeqid chr ++ " doesn't exist") id $ M.lookup chr (tbf_seqs tbf) - let go = getFwdSubseqWith (tbf_raw tbf) (tbs_dna_offset sq1) (tbs_n_blocks sq1) (tbs_m_blocks sq1) + let go = getFwdSubseqWith tbf sq1 if start < 0 - then reverse $ go (maskf . cmp_nt) (-start-len) len - else go (maskf . fwd_nt) start len + then reverse $ take len $ go (maskf . cmp_nt) (-start-len) + else take len $ go (maskf . fwd_nt) start where fwd_nt = (!!) [nucT, nucC, nucA, nucG] . fromIntegral cmp_nt = (!!) [nucA, nucG, nucT, nucC] . fromIntegral +-- | Works only in forward direction. +getLazySubseq :: TwoBitFile -> Position -> [Nucleotide] +getLazySubseq tbf (Pos { p_seq = chr, p_start = start }) = do + let sq1 = maybe (error $ unpackSeqid chr ++ " doesn't exist") id $ M.lookup chr (tbf_seqs tbf) + let go = getFwdSubseqWith tbf sq1 + if start < 0 + then error "sorry, can't go backwards" + -- then reverse $ take len $ go (maskf . cmp_nt) (-start-len) + else go fwd_nt start + where + fwd_nt n _ = [nucT, nucC, nucA, nucG] !! fromIntegral n + -- | Extract a subsequence without masking. getSubseq :: TwoBitFile -> Range -> [Nucleotide] getSubseq = getSubseqWith const @@ -190,8 +209,8 @@ getSeqnames :: TwoBitFile -> [Seqid] getSeqnames = M.keys . tbf_seqs -hasSequence :: TwoBitFile -> Seqid -> Bool -hasSequence tbf sq = isJust . M.lookup sq . tbf_seqs $ tbf +lookupSequence :: TwoBitFile -> Seqid -> Maybe TwoBitSequence +lookupSequence tbf sq = M.lookup sq . tbf_seqs $ tbf getSeqLength :: TwoBitFile -> Seqid -> Int getSeqLength tbf chr = diff --git a/src/Bio/Util.hs b/src/Bio/Util.hs deleted file mode 100644 --- a/src/Bio/Util.hs +++ /dev/null @@ -1,226 +0,0 @@ -module Bio.Util ( - wilson, invnormcdf, choose, - estimateComplexity, showNum, showOOM, - float2mini, mini2float, log1p, expm1, - phredplus, phredminus, phredsum, (<#>), phredconverse - ) where - -import Data.Bits -import Data.Char ( intToDigit ) -import Data.List ( foldl' ) -import Data.Word ( Word8 ) - --- ^ Random useful stuff I didn't know where to put. - --- | calculates the Wilson Score interval. --- If @(l,m,h) = wilson c x n@, then @m@ is the binary proportion and --- @(l,h)@ it's @c@-confidence interval for @x@ positive examples out of --- @n@ observations. @c@ is typically something like 0.05. - -wilson :: Double -> Int -> Int -> (Double, Double, Double) -wilson c x n = ( (m - h) / d, p, (m + h) / d ) - where - nn = fromIntegral n - p = fromIntegral x / nn - - z = invnormcdf (1-c*0.5) - h = z * sqrt (( p * (1-p) + 0.25*z*z / nn ) / nn) - m = p + 0.5 * z * z / nn - d = 1 + z * z / nn - -showNum :: Show a => a -> String -showNum = triplets [] . reverse . show - where - triplets acc [] = acc - triplets acc (a:[]) = a:acc - triplets acc (a:b:[]) = b:a:acc - triplets acc (a:b:c:[]) = c:b:a:acc - triplets acc (a:b:c:s) = triplets (',':c:b:a:acc) s - -showOOM :: Double -> String -showOOM x | x < 0 = '-' : showOOM (negate x) - | otherwise = findSuffix (x*10) ".kMGTPEZY" - where - findSuffix _ [] = "many" - findSuffix y (s:ss) | y < 100 = intToDigit (round y `div` 10) : case (round y `mod` 10, s) of - (0,'.') -> [] ; (0,_) -> [s] ; (d,_) -> [s, intToDigit d] - | y < 1000 = intToDigit (round y `div` 100) : intToDigit ((round y `mod` 100) `div` 10) : - if s == '.' then [] else [s] - | y < 10000 = intToDigit (round y `div` 1000) : intToDigit ((round y `mod` 1000) `div` 100) : - '0' : if s == '.' then [] else [s] - | otherwise = findSuffix (y*0.001) ss - --- Stolen from Lennart Augustsson's erf package, who in turn took it rom --- http://home.online.no/~pjacklam/notes/invnorm/ Accurate to about 1e-9. -invnormcdf :: (Ord a, Floating a) => a -> a -invnormcdf p = - let a1 = -3.969683028665376e+01 - a2 = 2.209460984245205e+02 - a3 = -2.759285104469687e+02 - a4 = 1.383577518672690e+02 - a5 = -3.066479806614716e+01 - a6 = 2.506628277459239e+00 - - b1 = -5.447609879822406e+01 - b2 = 1.615858368580409e+02 - b3 = -1.556989798598866e+02 - b4 = 6.680131188771972e+01 - b5 = -1.328068155288572e+01 - - c1 = -7.784894002430293e-03 - c2 = -3.223964580411365e-01 - c3 = -2.400758277161838e+00 - c4 = -2.549732539343734e+00 - c5 = 4.374664141464968e+00 - c6 = 2.938163982698783e+00 - - d1 = 7.784695709041462e-03 - d2 = 3.224671290700398e-01 - d3 = 2.445134137142996e+00 - d4 = 3.754408661907416e+00 - - pLow = 0.02425 - - nan = 0/0 - - in if p < 0 then - nan - else if p == 0 then - -1/0 - else if p < pLow then - let q = sqrt(-2 * log p) - in (((((c1*q+c2)*q+c3)*q+c4)*q+c5)*q+c6) / - ((((d1*q+d2)*q+d3)*q+d4)*q+1) - else if p < 1 - pLow then - let q = p - 0.5 - r = q*q - in (((((a1*r+a2)*r+a3)*r+a4)*r+a5)*r+a6)*q / - (((((b1*r+b2)*r+b3)*r+b4)*r+b5)*r+1) - else if p <= 1 then - - invnormcdf (1 - p) - else - nan - - --- | Try to estimate complexity of a whole from a sample. Suppose we --- sampled @total@ things and among those @singles@ occured only once. --- How many different things are there? --- --- Let the total number be @m@. The copy number follows a Poisson --- distribution with paramter @\lambda@. Let @z := e^{\lambda}@, then --- we have: --- --- P( 0 ) = e^{-\lambda} = 1/z --- P( 1 ) = \lambda e^{-\lambda} = ln z / z --- P(>=1) = 1 - e^{-\lambda} = 1 - 1/z --- --- singles = m ln z / z --- total = m (1 - 1/z) --- --- D := total/singles = (1 - 1/z) * z / ln z --- f := z - 1 - D ln z = 0 --- --- To get @z@, we solve using Newton iteration and then substitute to --- get @m@: --- --- df/dz = 1 - D/z --- z' := z - z (z - 1 - D ln z) / (z - D) --- m = singles * z /log z --- --- It converges as long as the initial @z@ is large enough, and @10D@ --- (in the line for @zz@ below) appears to work well. - -estimateComplexity :: (Integral a, Floating b, Ord b) => a -> a -> Maybe b -estimateComplexity total singles | total <= singles = Nothing - | singles <= 0 = Nothing - | otherwise = Just m - where - d = fromIntegral total / fromIntegral singles - step z = z * (z - 1 - d * log z) / (z - d) - iter z = case step z of zd | abs zd < 1e-12 -> z - | otherwise -> iter $! z-zd - zz = iter $! 10*d - m = fromIntegral singles * zz / log zz - - --- | Computes @-10 * log_10 (10 ** (-x\/10) + 10 ** (-y\/10))@ without --- losing precision. Used to add numbers on "the Phred scale", --- otherwise known as (deci-)bans. -{-# INLINE phredplus #-} -phredplus :: Double -> Double -> Double -phredplus x y = if x < y then pp x y else pp y x where - pp u v = u - 10 / log 10 * log1p (exp ((u-v) * log 10 / 10)) - --- | Computes @-10 * log_10 (10 ** (-x\/10) - 10 ** (-y\/10))@ without --- losing precision. Used to subtract numbers on "the Phred scale", --- otherwise known as (deci-)bans. -{-# INLINE phredminus #-} -phredminus :: Double -> Double -> Double -phredminus x y = if x < y then pm x y else pm y x where - pm u v = u - 10 / log 10 * log1p (- exp ((u-v) * log 10 / 10)) - --- | Computes @-10 * log_10 (sum [10 ** (-x\/10) | x <- xs])@ without losing --- precision. -{-# INLINE phredsum #-} -phredsum :: [Double] -> Double -phredsum = foldl' (<#>) (1/0) - -infixl 3 <#>, `phredminus`, `phredplus` -{-# INLINE (<#>) #-} -(<#>) :: Double -> Double -> Double -(<#>) = phredplus - --- | Computes @1-p@ without leaving the "Phred scale" -phredconverse :: Double -> Double -phredconverse v = - 10 / log 10 * log1p (- exp ((-v) * log 10 / 10)) - --- | Computes @log (1+x)@ to a relative precision of @10^-8@ even for --- very small @x@. Stolen from http://www.johndcook.com/cpp_log_one_plus_x.html -{-# INLINE log1p #-} -log1p :: (Floating a, Ord a) => a -> a -log1p x | x < -1 = error "log1p: argument must be greater than -1" - -- x is large enough that the obvious evaluation is OK: - | x > 0.0001 || x < -0.0001 = log $ 1 + x - -- Use Taylor approx. log(1 + x) = x - x^2/2 with error roughly x^3/3 - -- Since |x| < 10^-4, |x|^3 < 10^-12, relative error less than 10^-8: - | otherwise = (1 - 0.5*x) * x - - --- | Computes @exp x - 1@ to a relative precision of @10^-10@ even for --- very small @x@. Stolen from http://www.johndcook.com/cpp_expm1.html -expm1 :: (Floating a, Ord a) => a -> a -expm1 x | x > -0.00001 && x < 0.00001 = (1 + 0.5 * x) * x -- Taylor approx - | otherwise = exp x - 1 -- direct eval - - --- | Binomial coefficient: @n `choose` k == n! / ((n-k)! k!)@ -{-# INLINE choose #-} -choose :: Integral a => a -> a -> a -n `choose` k = product [n-k+1 .. n] `div` product [2..k] - - --- | Conversion to 0.4.4 format minifloat: This minifloat fits into a --- byte. It has no sign, four bits of precision, and the range is from --- 0 to 63488, initially in steps of 1/8. Nice to store quality scores --- with reasonable precision and range. -float2mini :: RealFloat a => a -> Word8 -float2mini f | f' < 0 = error "no negative minifloats" -- negative zero is fine! - | f < 2 = f' - | e >= 17 = 0xff - | s < 16 = error $ "oops: " ++ show (e,s) - | s < 32 = (e-1) `shiftL` 4 .|. (s .&. 0xf) - | s == 32 = e `shiftL` 4 - | otherwise = error $ "oops: " ++ show (e,s) - where - f' = round (8*f) - e = fromIntegral $ exponent f - s = round $ 32 * significand f - --- | Conversion from 0.4.4 format minifloat, see 'float2mini'. -mini2float :: Fractional a => Word8 -> a -mini2float w | e == 0 = fromIntegral w / 8.0 - | otherwise = 2^e * fromIntegral m / 16.0 - where - m = (w .&. 0xF) .|. 0x10 - e = w `shiftR` 4 - diff --git a/src/Bio/Util/AD.hs b/src/Bio/Util/AD.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Util/AD.hs @@ -0,0 +1,133 @@ +{-# LANGUAGE BangPatterns #-} +module Bio.Util.AD + ( AD(..), paramVector, minimize + , module Numeric.Optimization.Algorithms.HagerZhang05 + , debugParameters, quietParameters + ) where + +import Numeric.Optimization.Algorithms.HagerZhang05 +import qualified Data.Vector.Unboxed as U +import qualified Data.Vector.Storable as V + +-- | Simple forward-mode AD to get a scalar valued function with gradient. +data AD = C !Double | D !Double !(U.Vector Double) deriving Show + +instance Eq AD where + C x == C y = x == y + C x == D y _ = x == y + D x _ == C y = x == y + D x _ == D y _ = x == y + +instance Ord AD where + C x `compare` C y = x `compare` y + C x `compare` D y _ = x `compare` y + D x _ `compare` C y = x `compare` y + D x _ `compare` D y _ = x `compare` y + +instance Num AD where + {-# INLINE (+) #-} + C x + C y = C (x+y) + C x + D y v = D (x+y) v + D x u + C y = D (x+y) u + D x u + D y v = D (x+y) (U.zipWith (+) u v) + + {-# INLINE (-) #-} + C x - C y = C (x-y) + C x - D y v = D (x-y) (U.map negate v) + D x u - C y = D (x-y) u + D x u - D y v = D (x-y) (U.zipWith (-) u v) + + {-# INLINE (*) #-} + C x * C y = C (x*y) + C x * D y v = D (x*y) (U.map (x*) v) + D x u * C y = D (x*y) (U.map (y*) u) + D x u * D y v = D (x*y) (U.zipWith (+) (U.map (x*) v) (U.map (y*) u)) + + {-# INLINE negate #-} + negate (C x) = C (negate x) + negate (D x u) = D (negate x) (U.map negate u) + + {-# INLINE fromInteger #-} + fromInteger = C . fromInteger + + {-# INLINE abs #-} + abs (C x) = C (abs x) + abs (D x u) | x < 0 = D (negate x) (U.map negate u) + | otherwise = D x u + + {-# INLINE signum #-} + signum (C x) = C (signum x) + signum (D x _) = C (signum x) + + +instance Fractional AD where + {-# INLINE (/) #-} + C x / C y = C (x/y) + D x u / C y = D (x*z) (U.map (z*) u) where z = recip y + C x / D y v = D (x/y) (U.map (w*) v) where w = negate $ x * z * z ; z = recip y + D x u / D y v = D (x/y) (U.zipWith (-) (U.map (z*) u) (U.map (w*) v)) + where z = recip y ; w = x * z * z + + {-# INLINE recip #-} + recip = liftF recip (\x -> - recip (x*x)) + + {-# INLINE fromRational #-} + fromRational = C . fromRational + + +instance Floating AD where + {-# INLINE pi #-} + pi = C pi + + {-# INLINE exp #-} + exp = liftF exp exp + + {-# INLINE sqrt #-} + sqrt = liftF sqrt $ \x -> recip (2 * sqrt x) + + {-# INLINE log #-} + log = liftF log recip + + sin = liftF sin cos + cos = liftF cos (negate . sin) + sinh = liftF sinh cosh + cosh = liftF cosh sinh + + tan = liftF tan $ \x -> recip (cos x * cos x) + tanh = liftF tanh $ \x -> recip (cosh x * cosh x) + asin = liftF asin $ \x -> recip (sqrt (1 - x * x)) + acos = liftF acos $ \x -> - recip (sqrt (1 - x * x)) + atan = liftF atan $ \x -> recip (1 + x * x) + asinh = liftF asinh $ \x -> recip (sqrt (x * x + 1)) + acosh = liftF acosh $ \x -> - recip (sqrt (x * x - 1)) + atanh = liftF atanh $ \x -> recip (1 - x * x) + + +{-# INLINE liftF #-} +liftF :: (Double -> Double) -> (Double -> Double) -> AD -> AD +liftF f _ (C x) = C (f x) +liftF f g (D x u) = D (f x) (U.map (* g x) u) + +{-# INLINE paramVector #-} +paramVector :: [Double] -> [AD] +paramVector xs = [ D x (U.generate l (\j -> if i == j then 1 else 0)) | (i,x) <- zip [0..] xs ] + where l = length xs + +{-# INLINE minimize #-} +minimize :: Parameters -> Double -> ([AD] -> AD) -> U.Vector Double -> IO (V.Vector Double, Result, Statistics) +minimize params eps func v0 = + optimize params eps v0 (VFunction $ fst . combofn) + (VGradient $ snd . combofn) + (Just . VCombined $ combofn) + where + combofn parms = case func $ paramVector $ U.toList parms of + D x g -> ( x, g ) + C x -> ( x, U.replicate (U.length parms) 0 ) + + +quietParameters :: Parameters +quietParameters = defaultParameters { printFinal = False, verbose = Quiet, maxItersFac = 123 } + +debugParameters :: Parameters +debugParameters = defaultParameters { verbose = Verbose } + diff --git a/src/Bio/Util/AD2.hs b/src/Bio/Util/AD2.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Util/AD2.hs @@ -0,0 +1,132 @@ +{-# LANGUAGE BangPatterns #-} +module Bio.Util.AD2 ( AD2(..), paramVector2 ) where + +import qualified Data.Vector.Unboxed as U + +-- | Simple forward-mode AD to get a scalar valued function +-- with gradient and Hessian. +data AD2 = C2 !Double | D2 !Double !(U.Vector Double) !(U.Vector Double) + +instance Show AD2 where + show (C2 x) = show x + show (D2 x y z) = show x ++ " " ++ show (U.toList y) ++ " " + ++ show [ U.toList (U.slice i d z) | i <- [0, d .. d*d-1] ] + where d = U.length y + +instance Eq AD2 where + C2 x == C2 y = x == y + C2 x == D2 y _ _ = x == y + D2 x _ _ == C2 y = x == y + D2 x _ _ == D2 y _ _ = x == y + +instance Ord AD2 where + C2 x `compare` C2 y = x `compare` y + C2 x `compare` D2 y _ _ = x `compare` y + D2 x _ _ `compare` C2 y = x `compare` y + D2 x _ _ `compare` D2 y _ _ = x `compare` y + +instance Num AD2 where + {-# INLINE (+) #-} + C2 x + C2 y = C2 (x+y) + C2 x + D2 y v h = D2 (x+y) v h + D2 x u g + C2 y = D2 (x+y) u g + D2 x u g + D2 y v h = D2 (x+y) (U.zipWith (+) u v) (U.zipWith (+) g h) + + {-# INLINE (-) #-} + C2 x - C2 y = C2 (x-y) + C2 x - D2 y v h = D2 (x-y) (U.map negate v) (U.map negate h) + D2 x u g - C2 y = D2 (x-y) u g + D2 x u g - D2 y v h = D2 (x-y) (U.zipWith (-) u v) (U.zipWith (-) g h) + + {-# INLINE (*) #-} + C2 x * C2 y = C2 (x*y) + C2 x * D2 y v h = D2 (x*y) (U.map (x*) v) (U.map (x*) h) + D2 x u g * C2 y = D2 (x*y) (U.map (y*) u) (U.map (y*) g) + D2 x u g * D2 y v h = D2 (x*y) grad hess + where grad = U.zipWith (+) (U.map (x*) v) (U.map (y*) u) + hess = U.zipWith (+) + (U.zipWith (+) (U.map (x*) h) (U.map (y*) g)) + (U.zipWith (+) (cross u v) (cross v u)) + + {-# INLINE negate #-} + negate (C2 x) = C2 (negate x) + negate (D2 x u g) = D2 (negate x) (U.map negate u) (U.map negate g) + + {-# INLINE fromInteger #-} + fromInteger = C2 . fromInteger + + {-# INLINE abs #-} + abs (C2 x) = C2 (abs x) + abs (D2 x u g) | x < 0 = D2 (negate x) (U.map negate u) (U.map negate g) + | otherwise = D2 x u g + + {-# INLINE signum #-} + signum (C2 x) = C2 (signum x) + signum (D2 x _ _) = C2 (signum x) + + +instance Fractional AD2 where + {-# INLINE (/) #-} + C2 x / C2 y = C2 (x/y) + D2 x u g / C2 y = D2 (x*z) (U.map (z*) u) (U.map (z*) g) where z = recip y + x / y = x * recip y + + {-# INLINE recip #-} + recip = liftF recip (\x -> - recip (sqr x)) (\x -> 2 * recip (cube x)) + + {-# INLINE fromRational #-} + fromRational = C2 . fromRational + +instance Floating AD2 where + {-# INLINE pi #-} + pi = C2 pi + + {-# INLINE exp #-} + exp = liftF exp exp exp + + {-# INLINE sqrt #-} + sqrt = liftF sqrt (\x -> recip $ 2 * sqrt x) (\x -> - recip (sqrt (cube x))) + + {-# INLINE log #-} + log = liftF log recip (\x -> - recip (sqr x)) + + sin = liftF sin cos (negate . sin) + cos = liftF cos (negate . sin) (negate . cos) + sinh = liftF sinh cosh sinh + cosh = liftF cosh sinh cosh + + tan = liftF tan (\x -> recip (sqr (cos x))) (\x -> 2 * tan x / sqr (cos x)) + tanh = liftF tanh (\x -> recip (sqr (cosh x))) (\x -> -2 * tanh x / sqr (cosh x)) + + asin = liftF asin (\x -> recip (sqrt (1 - sqr x))) (\x -> x / sqrt (cube (1 - sqr x))) + acos = liftF acos (\x -> - recip (sqrt (1 - sqr x))) (\x -> -x / sqrt (cube (1 - sqr x))) + asinh = liftF asinh (\x -> recip (sqrt (sqr x + 1))) (\x -> -x / sqrt (cube (sqr x + 1))) + acosh = liftF acosh (\x -> - recip (sqrt (sqr x - 1))) (\x -> x / sqrt (cube (sqr x - 1))) + atan = liftF atan (\x -> recip (1 + sqr x)) (\x -> -2 * x / sqr (1 + sqr x)) + atanh = liftF atanh (\x -> recip (1 - sqr x)) (\x -> 2 * x / sqr (1 - sqr x)) + +{-# INLINE sqr #-} +sqr :: Double -> Double +sqr x = x * x + +{-# INLINE cube #-} +cube :: Double -> Double +cube x = x * x * x + +{-# INLINE liftF #-} +liftF :: (Double -> Double) -> (Double -> Double) -> (Double -> Double) -> AD2 -> AD2 +liftF f _ _ (C2 x) = C2 (f x) +liftF f f' f'' (D2 x v g) = D2 (f x) (U.map (* f' x) v) hess + where + hess = U.zipWith (+) (U.map (* f' x) g) (U.map (* f'' x) (cross v v)) + +{-# INLINE cross #-} +cross :: U.Vector Double -> U.Vector Double -> U.Vector Double +cross u v = U.concatMap (\dy -> U.map (dy*) u) v + +{-# INLINE paramVector2 #-} +paramVector2 :: [Double] -> [AD2] +paramVector2 xs = [ D2 x (U.generate l (\j -> if i == j then 1 else 0)) nil + | (i,x) <- zip [0..] xs ] + where l = length xs ; nil = U.replicate (l*l) 0 + diff --git a/src/Bio/Util/Numeric.hs b/src/Bio/Util/Numeric.hs new file mode 100644 --- /dev/null +++ b/src/Bio/Util/Numeric.hs @@ -0,0 +1,201 @@ +module Bio.Util.Numeric ( + wilson, invnormcdf, choose, + estimateComplexity, showNum, showOOM, + log1p, expm1, (<#>), + lsum, llerp, + sigmoid2, isigmoid2 + ) where + +import Data.List ( foldl1' ) +import Data.Char ( intToDigit ) + +-- ^ Random useful stuff I didn't know where to put. + +-- | calculates the Wilson Score interval. +-- If @(l,m,h) = wilson c x n@, then @m@ is the binary proportion and +-- @(l,h)@ it's @c@-confidence interval for @x@ positive examples out of +-- @n@ observations. @c@ is typically something like 0.05. + +wilson :: Double -> Int -> Int -> (Double, Double, Double) +wilson c x n = ( (m - h) / d, p, (m + h) / d ) + where + nn = fromIntegral n + p = fromIntegral x / nn + + z = invnormcdf (1-c*0.5) + h = z * sqrt (( p * (1-p) + 0.25*z*z / nn ) / nn) + m = p + 0.5 * z * z / nn + d = 1 + z * z / nn + +showNum :: Show a => a -> String +showNum = triplets [] . reverse . show + where + triplets acc [] = acc + triplets acc (a:[]) = a:acc + triplets acc (a:b:[]) = b:a:acc + triplets acc (a:b:c:[]) = c:b:a:acc + triplets acc (a:b:c:s) = triplets (',':c:b:a:acc) s + +showOOM :: Double -> String +showOOM x | x < 0 = '-' : showOOM (negate x) + | otherwise = findSuffix (x*10) ".kMGTPEZY" + where + findSuffix _ [] = "many" + findSuffix y (s:ss) | y < 100 = intToDigit (round y `div` 10) : case (round y `mod` 10, s) of + (0,'.') -> [] ; (0,_) -> [s] ; (d,_) -> [s, intToDigit d] + | y < 1000 = intToDigit (round y `div` 100) : intToDigit ((round y `mod` 100) `div` 10) : + if s == '.' then [] else [s] + | y < 10000 = intToDigit (round y `div` 1000) : intToDigit ((round y `mod` 1000) `div` 100) : + '0' : if s == '.' then [] else [s] + | otherwise = findSuffix (y*0.001) ss + +-- Stolen from Lennart Augustsson's erf package, who in turn took it rom +-- http://home.online.no/~pjacklam/notes/invnorm/ Accurate to about 1e-9. +invnormcdf :: (Ord a, Floating a) => a -> a +invnormcdf p = + let a1 = -3.969683028665376e+01 + a2 = 2.209460984245205e+02 + a3 = -2.759285104469687e+02 + a4 = 1.383577518672690e+02 + a5 = -3.066479806614716e+01 + a6 = 2.506628277459239e+00 + + b1 = -5.447609879822406e+01 + b2 = 1.615858368580409e+02 + b3 = -1.556989798598866e+02 + b4 = 6.680131188771972e+01 + b5 = -1.328068155288572e+01 + + c1 = -7.784894002430293e-03 + c2 = -3.223964580411365e-01 + c3 = -2.400758277161838e+00 + c4 = -2.549732539343734e+00 + c5 = 4.374664141464968e+00 + c6 = 2.938163982698783e+00 + + d1 = 7.784695709041462e-03 + d2 = 3.224671290700398e-01 + d3 = 2.445134137142996e+00 + d4 = 3.754408661907416e+00 + + pLow = 0.02425 + + nan = 0/0 + + in if p < 0 then + nan + else if p == 0 then + -1/0 + else if p < pLow then + let q = sqrt(-2 * log p) + in (((((c1*q+c2)*q+c3)*q+c4)*q+c5)*q+c6) / + ((((d1*q+d2)*q+d3)*q+d4)*q+1) + else if p < 1 - pLow then + let q = p - 0.5 + r = q*q + in (((((a1*r+a2)*r+a3)*r+a4)*r+a5)*r+a6)*q / + (((((b1*r+b2)*r+b3)*r+b4)*r+b5)*r+1) + else if p <= 1 then + - invnormcdf (1 - p) + else + nan + + +-- | Try to estimate complexity of a whole from a sample. Suppose we +-- sampled @total@ things and among those @singles@ occured only once. +-- How many different things are there? +-- +-- Let the total number be @m@. The copy number follows a Poisson +-- distribution with paramter @\lambda@. Let @z := e^{\lambda}@, then +-- we have: +-- +-- P( 0 ) = e^{-\lambda} = 1/z +-- P( 1 ) = \lambda e^{-\lambda} = ln z / z +-- P(>=1) = 1 - e^{-\lambda} = 1 - 1/z +-- +-- singles = m ln z / z +-- total = m (1 - 1/z) +-- +-- D := total/singles = (1 - 1/z) * z / ln z +-- f := z - 1 - D ln z = 0 +-- +-- To get @z@, we solve using Newton iteration and then substitute to +-- get @m@: +-- +-- df/dz = 1 - D/z +-- z' := z - z (z - 1 - D ln z) / (z - D) +-- m = singles * z /log z +-- +-- It converges as long as the initial @z@ is large enough, and @10D@ +-- (in the line for @zz@ below) appears to work well. + +estimateComplexity :: (Integral a, Floating b, Ord b) => a -> a -> Maybe b +estimateComplexity total singles | total <= singles = Nothing + | singles <= 0 = Nothing + | otherwise = Just m + where + d = fromIntegral total / fromIntegral singles + step z = z * (z - 1 - d * log z) / (z - d) + iter z = case step z of zd | abs zd < 1e-12 -> z + | otherwise -> iter $! z-zd + zz = iter $! 10*d + m = fromIntegral singles * zz / log zz + + +-- | Computes @log (exp x + exp y)@ without leaving the log domain and +-- hence without losing precision. +infixl 5 <#> +{-# INLINE (<#>) #-} +(<#>) :: (Floating a, Ord a) => a -> a -> a +x <#> y = if x >= y then x + log1p (exp (y-x)) else y + log1p (exp (x-y)) + +-- | Computes @log (1+x)@ to a relative precision of @10^-8@ even for +-- very small @x@. Stolen from http://www.johndcook.com/cpp_log_one_plus_x.html +{-# INLINE log1p #-} +log1p :: (Floating a, Ord a) => a -> a +log1p x | x < -1 = error "log1p: argument must be greater than -1" + -- x is large enough that the obvious evaluation is OK: + | x > 0.0001 || x < -0.0001 = log $ 1 + x + -- Use Taylor approx. log(1 + x) = x - x^2/2 with error roughly x^3/3 + -- Since |x| < 10^-4, |x|^3 < 10^-12, relative error less than 10^-8: + | otherwise = (1 - 0.5*x) * x + + +-- | Computes @exp x - 1@ to a relative precision of @10^-10@ even for +-- very small @x@. Stolen from http://www.johndcook.com/cpp_expm1.html +expm1 :: (Floating a, Ord a) => a -> a +expm1 x | x > -0.00001 && x < 0.00001 = (1 + 0.5 * x) * x -- Taylor approx + | otherwise = exp x - 1 -- direct eval + +-- | Computes \( \log ( \sum_i e^{x_i} ) \) sensibly. The list must be +-- sorted in descending(!) order. +{-# INLINE lsum #-} +lsum :: (Floating a, Ord a) => [a] -> a +lsum xs = foldl1' (\x y -> if x >= y then x + log1p (exp (y-x)) else err) xs + where err = error $ "lsum: argument list must be in descending order" + +-- | Computes \( \log \left( c e^x + (1-c) e^y \right) \). +{-# INLINE llerp #-} +llerp :: (Floating a, Ord a) => a -> a -> a -> a +llerp c x y | c <= 0.0 = y + | c >= 1.0 = x + | x >= y = log c + x + log1p ( (1-c)/c * exp (y-x) ) + | otherwise = log1p (-c) + y + log1p ( c/(1-c) * exp (x-y) ) + +-- | Binomial coefficient: @n `choose` k == n! / ((n-k)! k!)@ +{-# INLINE choose #-} +choose :: Integral a => a -> a -> a +n `choose` k = product [n-k+1 .. n] `div` product [2..k] + + +-- | Kind-of sigmoid function that maps the reals to the interval +-- @[0,1)@. Good to compute a probability without introducing boundary +-- conditions. +sigmoid2 :: (Num a, Fractional a, Floating a) => a -> a +sigmoid2 x = y*y where y = (exp x - 1) / (exp x + 1) + +-- | Inverse of 'sigmoid2'. +isigmoid2 :: (Num a, Fractional a, Floating a) => a -> a +isigmoid2 y = log $ (1 + sqrt y) / (1 - sqrt y) + + diff --git a/src/Bio/Util/Regex.hsc b/src/Bio/Util/Regex.hsc new file mode 100644 --- /dev/null +++ b/src/Bio/Util/Regex.hsc @@ -0,0 +1,44 @@ +{-# LANGUAGE CPP, ForeignFunctionInterface #-} +-- | The absolute minimum necessary for regex matching using POSIX regexec. +module Bio.Util.Regex ( Regex, regComp, regMatch )where + +#include +#include + +import Control.Applicative +import Control.Monad +import Foreign.Ptr +import Foreign.ForeignPtr +import Foreign.Marshal.Alloc +import Foreign.C.String +import Foreign.C.Types +import System.IO.Unsafe + +newtype Regex = Regex (ForeignPtr Regex) + +regComp :: String -> Regex +regComp re = unsafePerformIO $ do + fp <- mallocForeignPtrBytes #{size regex_t} + withForeignPtr fp $ \p -> do + withCString re $ \pre -> do + ec <- regcomp p pre (#{const REG_EXTENDED} + #{const REG_NOSUB}) + when (ec /= 0) $ do + sz <- regerror ec p nullPtr 0 + allocaBytes (fromIntegral sz) $ \err -> do + _ <- regerror ec p err sz + peekCString err >>= error . (++) "regexec: " + addForeignPtrFinalizer regfree fp + return $ Regex fp + +regMatch :: Regex -> String -> Bool +regMatch (Regex fp) str = + unsafePerformIO $ + withForeignPtr fp $ \p -> + withCString str $ \s -> + (==) 0 <$> regexec p s 0 nullPtr 0 + + +foreign import ccall unsafe regcomp :: Ptr Regex -> CString -> CInt -> IO CInt +foreign import ccall unsafe regexec :: Ptr Regex -> CString -> CSize -> Ptr () -> CInt -> IO CInt +foreign import ccall unsafe regerror :: CInt -> Ptr Regex -> CString -> CSize -> IO CSize +foreign import ccall unsafe "®free" regfree :: FunPtr (Ptr Regex -> IO ()) diff --git a/src/Data/Avro.hs b/src/Data/Avro.hs --- a/src/Data/Avro.hs +++ b/src/Data/Avro.hs @@ -1,15 +1,12 @@ {-# LANGUAGE OverloadedStrings, FlexibleInstances, TemplateHaskell #-} {-# LANGUAGE RecordWildCards, BangPatterns, FlexibleContexts #-} +{-# LANGUAGE PatternGuards #-} module Data.Avro where import Bio.Iteratee import Control.Applicative import Control.Monad -import Control.Monad.ST ( runST, ST ) -import Data.Aeson hiding ((.=)) -import Data.Array.MArray -import Data.Array.ST ( STUArray ) -import Data.Array.Unsafe ( castSTUArray ) +import Data.Aeson import Data.Binary.Get import Data.Bits import Data.Binary.Builder @@ -19,18 +16,20 @@ import Data.Monoid import Data.Scientific import Data.Text.Encoding -import Data.Word ( Word32, Word64 ) -import Foreign.Storable ( Storable, sizeOf ) +import Data.Word ( Word8, Word32, Word64 ) +import Foreign.Marshal.Alloc ( alloca ) +import Foreign.Storable ( Storable, sizeOf, peek, pokeByteOff ) import Language.Haskell.TH import System.Random +import System.IO.Unsafe ( unsafeDupablePerformIO ) -import qualified Data.ByteString as B -import qualified Data.ByteString.Lazy as BL -import qualified Data.HashMap.Strict as H -import qualified Data.ListLike as LL -import qualified Data.Text as T -import qualified Data.Vector as V -import qualified Data.Vector.Unboxed as U +import qualified Data.ByteString as B +import qualified Data.ByteString.Lazy as BL +import qualified Data.HashMap.Strict as H +import qualified Data.ListLike as LL +import qualified Data.Text as T +import qualified Data.Vector as V +import qualified Data.Vector.Unboxed as U -- ^ Support for Avro. -- Current status is that we can generate schemas for certain Haskell @@ -43,12 +42,6 @@ -- product uses record syntax and the top level is a plain record. -- The obvious primitives are supported. -(.=) :: ToJSON a => String -> a -> (T.Text, Value) -k .= v = (T.pack k, toJSON v) - -string :: String -> Value -string = String . T.pack - -- | This is the class of types we can embed into the Avro -- infrastructure. Right now, we can derive a schema, encode to -- the Avro binary format, and encode to the Avro JSON encoding. @@ -75,6 +68,7 @@ toAvron :: a -> Value +-- | Making schemas requires a memo table of type definitions. newtype MkSchema a = MkSchema { mkSchema :: (a -> H.HashMap T.Text Value -> Value) -> H.HashMap T.Text Value -> Value } @@ -93,8 +87,17 @@ Just obj' | obj == obj' -> k (String nm') h | otherwise -> error $ "same type name, different schema: " ++ nm -runMkSchema :: MkSchema Value -> Value -runMkSchema x = mkSchema x postproc H.empty +getNamedSchema :: String -> MkSchema Value +getNamedSchema nm = MkSchema $ \k h -> + let nm' = T.pack nm + in case H.lookup nm' h of + Nothing -> error $ "Schema for " ++ nm ++ " not provided." + -- Use the provided schema now, use only the name next time. + Just obj -> k obj $! H.insert nm' (String nm') h + + +runMkSchema :: MkSchema Value -> H.HashMap T.Text Value -> Value +runMkSchema x = mkSchema x postproc where -- Objects are fine as is. postproc (Object o) _ = Object o @@ -126,6 +129,12 @@ fromBin = decodeIntBase128 toAvron = Number . fromIntegral +instance Avro Word8 where + toSchema _ = return $ String "long" + toBin = encodeIntBase128 + fromBin = decodeIntBase128 + toAvron = Number . fromIntegral + instance Avro Int64 where toSchema _ = return $ String "long" toBin = encodeIntBase128 @@ -156,28 +165,31 @@ fromBin = decodeUtf8 <$> fromBin toAvron = String - -- Integer<->Float conversions, stolen from cereal. {-# INLINE wordToFloat #-} wordToFloat :: Word32 -> Float -wordToFloat x = runST (cast x) +wordToFloat x = cast x {-# INLINE wordToDouble #-} wordToDouble :: Word64 -> Double -wordToDouble x = runST (cast x) +wordToDouble x = cast x {-# INLINE floatToWord #-} floatToWord :: Float -> Word32 -floatToWord x = runST (cast x) +floatToWord x = cast x {-# INLINE doubleToWord #-} doubleToWord :: Double -> Word64 -doubleToWord x = runST (cast x) +doubleToWord x = cast x {-# INLINE cast #-} -cast :: ( MArray (STUArray s) b (ST s), MArray (STUArray s) a (ST s) ) => a -> ST s b -cast x = (newArray (0 :: Int, 0) x >>= castSTUArray >>= flip readArray 0) +cast :: ( Storable a, Storable b ) => a -> b +cast x | sizeOf x == sizeOf y = y + | otherwise = error "cannot cast: size mismatch" + where + y = unsafeDupablePerformIO $ alloca $ \buf -> + pokeByteOff buf 0 x >> peek buf -- | Implements Zig-Zag-Coding like in Protocol Buffers and Avro. zig :: (Storable a, Bits a) => a -> a @@ -199,10 +211,8 @@ decodeWordBase128 = go 0 0 where go acc sc = do x <- getWord8 - let !acc' = acc .|. fromIntegral x `shiftL` sc - if x .&. 0x80 == 0 - then return acc' - else go acc' (sc+7) + let !acc' = acc .|. (fromIntegral x .&. 0x7f) `shiftL` sc + if x .&. 0x80 == 0 then return acc' else go acc' (sc+7) -- | Encodes an int of any size by combining the zig-zag coding with the -- base 128 encoding. @@ -218,7 +228,7 @@ zigInt = encodeIntBase128 zagInt :: Get Int -zagInt = decodeWordBase128 +zagInt = decodeIntBase128 -- Complex Types @@ -289,6 +299,7 @@ where get_blocks !acc = zagInt >>= \l -> if l == 0 then return acc else get_block acc l >>= get_blocks + get_block !acc l = if l == 0 then return acc else fromBin >>= \k -> fromBin >>= \v -> get_block (H.insert k v acc) (l-1) @@ -344,8 +355,8 @@ mk_enum_inst :: [Name] -> Q [Dec] mk_enum_inst nms = [d| instance Avro $(conT nm) where - toSchema _ = return $ object [ "type" .= string "enum" - , "name" .= string $(tolit nm) + toSchema _ = return $ object [ "type" .= String "enum" + , "name" .= String $(tolit nm) , "symbols" .= $(tolitlist nms) ] toBin x = $( return $ CaseE (VarE 'x) @@ -364,7 +375,7 @@ toAvron x = $( return $ CaseE (VarE 'x) [ Match (ConP nm1 []) - (NormalB (AppE (VarE 'string) + (NormalB (AppE (ConE 'String) (LitE (StringL (nameBase nm1))))) [] | nm1 <- nms ] ) |] @@ -411,7 +422,7 @@ mk_product_schema nm1 tps = [| $( fieldlist tps ) >>= \flds -> memoObject $( tolit nm1 ) - [ "type" .= string "record" + [ "type" .= String "record" , "fields" .= Array (V.fromList flds) ] |] fieldlist = foldr go [| return [] |] @@ -419,7 +430,7 @@ go (nm1,_,tp) k = [| do sch <- toSchema $(sigE (varE 'undefined) (return tp)) obs <- $k - return $ object [ "name" .= string $(tolit nm1) + return $ object [ "name" .= String $(tolit nm1) , "type" .= sch ] : obs |] @@ -439,7 +450,9 @@ data ContainerOpts = ContainerOpts { objects_per_block :: Int - , filetype_label :: B.ByteString } + , filetype_label :: B.ByteString + , initial_schemas :: H.HashMap T.Text Value + , meta_info :: H.HashMap T.Text B.ByteString } -- Writing a container file. This is an 'Enumeratee', we read a list of -- suitable types, we write a header containing the generated schema, @@ -450,12 +463,11 @@ ma <- peekStream sync_marker <- liftIO $ B.pack <$> replicateM 16 randomIO - let schema = encode . runMkSchema . toSchema . fromJust $ ma + let schema = encode $ runMkSchema (toSchema $ fromJust ma) initial_schemas - meta :: H.HashMap T.Text B.ByteString - meta = H.fromList [( "avro.schema", B.concat $ BL.toChunks schema ) - ,( "avro.codec", "null" ) - ,( "biohazard.filetype", filetype_label )] + meta = H.insert "avro.schema" (B.concat $ BL.toChunks schema) $ + H.insert "avro.codec" "null" $ + H.insert "biohazard.filetype" filetype_label $ meta_info hdr = fromByteString "Obj\1" <> toBin meta <> fromByteString sync_marker @@ -463,46 +475,67 @@ foldStream (\(!n,c) o -> (n+1, c <> toBin o)) (0::Int,mempty) let code1 = toLazyByteString code - block = toBin num <> toBin (BL.length code1) <> + block = zigInt num <> toBin (BL.length code1) <> fromLazyByteString code1 <> fromByteString sync_marker lift (enumList (BL.toChunks $ toLazyByteString block) out') lift (enumList (BL.toChunks $ toLazyByteString hdr) out) >>= enc_blocks +-- | Avro Meta Data is currently unprocessed. Contains the codec, the +-- schema, a version number. +type AvroMeta = H.HashMap T.Text B.ByteString + +-- | Decodes an AVRO container file into a list. Meta data is passed +-- on. Note that if this blows up, it's usually due to it being applied +-- at the wrong type. Be sure to correctly count the brackets... +-- -- XXX Possible codecs: null, zlib, snappy, lzma; all missing -- XXX Should check schema on reading. -readAvroContainer :: (Monad m, ListLike s a, Avro a) => Enumeratee B.ByteString s m r +readAvroContainer :: (Monad m, Avro a) => Enumeratee' AvroMeta B.ByteString [a] m r readAvroContainer out = do 4 <- heads "Obj\1" -- enough magic? - meta <- iterGet (fromBin :: Get (H.HashMap T.Text B.ByteString)) + meta <- iterGet fromBin sync_marker <- iGetString 16 - flip iterLoop out $ \o -> do num <- iterGet zagInt - sz <- iterGet fromBin - o' <- joinI $ takeStream sz $ -- codec goes here - convStream (LL.singleton `liftM` iterGet fromBin) o - 16 <- heads sync_marker - return o' + flip iterLoop (out meta) $ \o -> do + _num <- iterGet zagInt + sz <- iterGet zagInt + -- liftIO $ hPutStrLn stderr $ "got block: " ++ showNum num + -- ++ " things in " ++ showNum sz ++ " bytes." + o' <- joinI $ takeStream sz -- codec goes here + $ convStream (LL.singleton `liftM` iterGet fromBin) o + 16 <- heads sync_marker + -- liftIO $ hPutStrLn stderr "got good sync" + return o' --- | Repeatedly apply an 'Iteratee' to a value until end of stream. --- Returns the final value. -iterLoop :: (Nullable s, Monad m) => (a -> Iteratee s m a) -> a -> Iteratee s m a -iterLoop it a = do e <- isFinished - if e then return a - else it a >>= iterLoop it +-- | Finds a names schema from the meta data of an Avro container. +findSchema :: T.Text -> AvroMeta -> Value +findSchema nm meta = maybe Null go $ decodeStrict =<< H.lookup "avro.schema" meta + where + go :: Value -> Value + go (Object obj) + | Just (String k) <- H.lookup "name" obj, k == nm -- found it + = Object obj + | Just (String "record") <- H.lookup "type" obj -- record, descend into "fields" + = maybe Null go_struct $ H.lookup "fields" obj + | Just (String "array") <- H.lookup "type" obj -- array, descend into "items" + = maybe Null go_union $ H.lookup "items" obj + | Just (String "map") <- H.lookup "type" obj -- map, descend into "values" + = maybe Null go $ H.lookup "values" obj + go _ = Null -iterGet :: Monad m => Get a -> Iteratee B.ByteString m a -iterGet = go . runGetIncremental - where - go (Fail _ _ err) = throwErr (iterStrExc err) - go (Done rest _ a) = idone a (Chunk rest) - go (Partial dec) = liftI $ \ck -> case ck of - Chunk s -> go (dec $ Just s) - EOF mx -> case dec Nothing of - Fail _ _ err -> throwErr (iterStrExc err) - Partial _ -> throwErr (iterStrExc "") - Done rest _ a | B.null rest -> idone a (EOF mx) - | otherwise -> idone a (Chunk rest) + go_struct (Array arr) = V.foldr (try_next . go') Null arr -- struct fields, recurse into "type" subfield + go_struct _ = Null + + go' (Object o) | Just o' <- H.lookup "type" o = go o' + go' _ = Null + + go_union (Array arr) = V.foldr (try_next . go) Null arr -- union arms, recurse + go_union _ = Null + + try_next Null b = b + try_next a _ = a + diff --git a/src/Data/MiniFloat.hs b/src/Data/MiniFloat.hs new file mode 100644 --- /dev/null +++ b/src/Data/MiniFloat.hs @@ -0,0 +1,44 @@ +{-# LANGUAGE TypeFamilies, FlexibleInstances, CPP #-} +{-# LANGUAGE MultiParamTypeClasses, TemplateHaskell #-} +module Data.MiniFloat ( Mini(..), float2mini, mini2float ) where + +import Data.Bits +import Data.Ix +import Data.Word ( Word8 ) +import Data.Vector.Unboxed.Deriving ( derivingUnbox ) + +#if __GLASGOW_HASKELL__ == 704 +import Data.Vector.Generic ( Vector(..) ) +import Data.Vector.Generic.Mutable ( MVector(..) ) +#endif + +data Mini = Mini { unMini :: Word8 } deriving ( Eq, Ord, Show, Ix, Bounded ) + +derivingUnbox "Mini" [t| Mini -> Word8 |] [| unMini |] [| Mini |] + +-- | Conversion to 0.4.4 format minifloat: This minifloat fits into a +-- byte. It has no sign, four bits of precision, and the range is from +-- 0 to 63488, initially in steps of 1/8. Nice to store quality scores +-- with reasonable precision and range. +float2mini :: RealFloat a => a -> Mini +float2mini f | f' < 0 = error "no negative minifloats" -- negative zero is fine! + | f < 2 = Mini f' + | e >= 17 = Mini 0xff + | s < 16 = error $ "oops: " ++ show (e,s) + | s < 32 = Mini $ (e-1) `shiftL` 4 .|. (s .&. 0xf) + | s == 32 = Mini $ e `shiftL` 4 + | otherwise = error $ "oops: " ++ show (e,s) + where + f' = round (8*f) + e = fromIntegral $ exponent f + s = round $ 32 * significand f + +-- | Conversion from 0.4.4 format minifloat, see 'float2mini'. +mini2float :: Fractional a => Mini -> a +mini2float (Mini w) | e == 0 = fromIntegral w / 8.0 + | otherwise = 2^e * fromIntegral m / 16.0 + where + m = (w .&. 0xF) .|. 0x10 + e = w `shiftR` 4 + + diff --git a/src/cbits/myers_align.h b/src/cbits/myers_align.h --- a/src/cbits/myers_align.h +++ b/src/cbits/myers_align.h @@ -9,12 +9,12 @@ //! \brief aligns two sequences in O(nd) time //! This alignment algorithm following Eugene W. Myers: "An O(ND) //! Difference Algorithm and Its Variations". -//! Both input sequences are ASCIIZ-encoded with IUPAC ambiguity codes. -//! By definition, if ambiguity codes overlap, that's a match, else a -//! mismatch. Mismatches and gaps count a unit penalty. If mode is -//! myers_align_globally, both sequences must align completely. If mode -//! is myers_align_is_prefix, seq_a must align completely as prefix of -//! seq_b. If mode is myers_align_has_prefix, seq_b must align +//! Both input sequences are ASCIIZ-encoded with IUPAC-IUB ambiguity +//! codes. By definition, if ambiguity codes overlap, that's a match, +//! else a mismatch. Mismatches and gaps count a unit penalty. If mode +//! is myers_align_globally, both sequences must align completely. If +//! mode is myers_align_is_prefix, seq_a must align completely as prefix +//! of seq_b. If mode is myers_align_has_prefix, seq_b must align //! completely as prefix of seq_a. //! //! Note that the calculation time is O(nd) where n is the length of the @@ -37,9 +37,8 @@ const char* seq_b, int len_b, int maxd, char *bt_a, char *bt_b ) ; -//! \brief converts an IUPAC ambiguity code to a bitmap -//! Each base is represented by a bit, makes checking for matches -//! easier. +//! \brief converts an IUPAC-IUB ambiguity code to a bitmap Each base is +//! represented by a bit, makes checking for matches easier. inline int char_to_bitmap( char x ) { switch( x & ~32 ) diff --git a/tools/AD.hs b/tools/AD.hs deleted file mode 100644 --- a/tools/AD.hs +++ /dev/null @@ -1,99 +0,0 @@ -{-# LANGUAGE BangPatterns #-} -module AD where - -import qualified Data.Vector.Unboxed as U - --- Simple forward-mode AD to get a scalar valued function and a --- gradient. - -data AD = C !Double | D !Double !(U.Vector Double) - deriving Show - -instance Num AD where - {-# INLINE (+) #-} - C x + C y = C (x+y) - C x + D y v = D (x+y) v - D x u + C y = D (x+y) u - D x u + D y v = D (x+y) (U.zipWith (+) u v) - - {-# INLINE (-) #-} - C x - C y = C (x-y) - C x - D y v = D (x-y) (U.map negate v) - D x u - C y = D (x-y) u - D x u - D y v = D (x-y) (U.zipWith (-) u v) - - {-# INLINE (*) #-} - C x * C y = C (x*y) - C x * D y v = D (x*y) (U.map (x*) v) - D x u * C y = D (x*y) (U.map (y*) u) - D x u * D y v = D (x*y) (U.zipWith (+) (U.map (x*) v) (U.map (y*) u)) - - {-# INLINE negate #-} - negate (C x) = C (negate x) - negate (D x u) = D (negate x) (U.map negate u) - - {-# INLINE fromInteger #-} - fromInteger = C . fromInteger - - {-# INLINE abs #-} - abs (C x) = C (abs x) - abs (D x u) | x < 0 = D (negate x) (U.map negate u) - | otherwise = D x u - - {-# INLINE signum #-} - signum (C x) = C (signum x) - signum (D x _) = C (signum x) - - -instance Fractional AD where - {-# INLINE (/) #-} - C x / C y = C (x/y) - D x u / C y = D (x*z) (U.map (z*) u) where z = recip y - C x / D y v = D (x/y) (U.map (w*) v) where w = negate $ x * z * z ; z = recip y - D x u / D y v = D (x/y) (U.zipWith (-) (U.map (z*) u) (U.map (w*) v)) - where z = recip y ; w = x * z * z - - {-# INLINE recip #-} - recip (C x) = C (recip x) - recip (D x u) = D (recip x) (U.map (y*) u) where y = negate $ recip $ x*x - - {-# INLINE fromRational #-} - fromRational = C . fromRational - - -instance Floating AD where - {-# INLINE pi #-} - pi = C pi - - {-# INLINE exp #-} - exp (C x) = C (exp x) - exp (D x u) = D (exp x) (U.map (* exp x) u) - - {-# INLINE sqrt #-} - sqrt (C x) = C (sqrt x) - sqrt (D x u) = D (sqrt x) (U.map (*w) u) where w = recip $ 2 * sqrt x - - {-# INLINE log #-} - log (C x) = C (log x) - log (D x u) = D (log x) (U.map (*w) u) where w = recip x - - {- (**) = undefined -- :: a -> a -> a - logBase = undefined -- :: a -> a -> a - sin = undefined -- :: a -> a - tan = undefined -- :: a -> a - cos = undefined -- :: a -> a - asin = undefined -- :: a -> a - atan = undefined -- :: a -> a - acos = undefined -- :: a -> a - sinh = undefined -- :: a -> a - tanh = undefined -- :: a -> a - cosh = undefined -- :: a -> a - asinh = undefined -- :: a -> a - atanh = undefined -- :: a -> a - acosh = undefined -- :: a -> a -} - - -paramVector :: [Double] -> [AD] -paramVector xs = [ D x (U.generate l (\j -> if i == j then 1 else 0)) | (i,x) <- zip [0..] xs ] - where l = length xs - diff --git a/tools/Index.hs b/tools/Index.hs --- a/tools/Index.hs +++ b/tools/Index.hs @@ -1,4 +1,5 @@ -{-# LANGUAGE TemplateHaskell, GeneralizedNewtypeDeriving, MultiParamTypeClasses, TypeFamilies #-} +{-# LANGUAGE TemplateHaskell, GeneralizedNewtypeDeriving #-} +{-# LANGUAGE MultiParamTypeClasses, TypeFamilies, CPP #-} module Index where -- ^ This tiny module defines the 'Index' type and derives the 'Unbox' @@ -11,8 +12,11 @@ import Data.Vector.Unboxed.Deriving import Data.Word ( Word64 ) import Foreign.Storable ( Storable ) + +#if __GLASGOW_HASKELL__ == 704 import Data.Vector.Generic ( Vector(..) ) import Data.Vector.Generic.Mutable ( MVector(..) ) +#endif -- | An index sequence must have at most eight bases. We represent a -- base and its quality score in a single byte: the top three bits are diff --git a/tools/afroengineer.hs b/tools/afroengineer.hs --- a/tools/afroengineer.hs +++ b/tools/afroengineer.hs @@ -22,7 +22,6 @@ import Data.Bits import Data.Char import Data.List ( isSuffixOf ) -import Data.Monoid import Numeric import Prelude hiding ( round ) import System.Console.GetOpt @@ -34,7 +33,6 @@ import qualified Bio.Iteratee.ZLib as ZLib import qualified Data.ByteString.Char8 as S import qualified Data.ByteString.Lazy.Char8 as L -import qualified Data.Foldable as F import qualified Data.Iteratee as I import qualified Data.Sequence as Z import qualified Data.Vector.Generic as V diff --git a/tools/bam-fixpair.hs b/tools/bam-fixpair.hs --- a/tools/bam-fixpair.hs +++ b/tools/bam-fixpair.hs @@ -32,10 +32,11 @@ import Bio.Bam.Header import Bio.Bam.Reader hiding ( mergeInputs, combineCoordinates ) import Bio.Bam.Rec +import Bio.Bam.Trim import Bio.Bam.Writer import Bio.Iteratee import Bio.PriorityQueue -import Bio.Util ( showNum ) +import Bio.Util.Numeric ( showNum ) import Control.Arrow ( (&&&) ) import Control.Applicative import Control.Monad @@ -53,6 +54,7 @@ import Text.Printf import qualified Data.ByteString as S +import qualified Data.Vector.Generic as V data Verbosity = Silent | Errors | Warnings | Notices deriving (Eq, Ord) data KillMode = KillNone | KillUu | KillAll deriving (Eq, Ord) @@ -65,10 +67,11 @@ , report_ixs :: !Bool , verbosity :: Verbosity , killmode :: KillMode - , output :: BamMeta -> Iteratee [BamRec] IO () } + , output :: BamMeta -> Iteratee [BamRec] IO () + , fixsven :: Maybe Int } config0 :: IO Config -config0 = return $ CF True True False True False True Errors KillNone (protectTerm . pipeBamOutput) +config0 = return $ CF True True False True False True Errors KillNone (protectTerm . pipeBamOutput) Nothing options :: [OptDescr (Config -> IO Config)] options = [ @@ -96,13 +99,15 @@ Option "" ["no-report-fflag"] (NoArg (\c -> return $ c { report_fflag = False })) "Do not report commonly inconsistent flags", Option "" ["no-report-fflag"] (NoArg (\c -> return $ c { report_ixs = False })) "Do not report mismatched index fields", + Option "" ["fix-sven"] (ReqArg set_fixsven "QUAL") "Trim 3' ends of avg qual lower than QUAL", + Option "h?" ["help","usage"] (NoArg usage) "Print this helpful message and exit", Option "V" ["version"] (NoArg vrsn) "Print version number and exit" ] where usage _ = do pn <- getProgName - let blah = "Usage: " ++ pn ++ " [OPTION...] [FILE...]\n\ - \Merge BAM files, rearrange them to move mate pairs together, \ - \output a file with consistent mate pair information." + let blah = "Usage: " ++ pn ++ " [OPTION...] [FILE...]\n" ++ + "Merge BAM files, rearrange them to move mate pairs together, " ++ + "output a file with consistent mate pair information." hPutStrLn stderr $ usageInfo blah options exitSuccess @@ -113,6 +118,7 @@ set_output "-" c = return $ c { output = pipeBamOutput } set_output f c = return $ c { output = writeBamFile f } set_validate c = return $ c { output = \_ -> skipToEof } + set_fixsven a c = readIO a >>= \q -> return $ c { fixsven = Just q } -- XXX placeholder... @@ -127,6 +133,7 @@ withQueues $ \queues -> mergeInputs files >=> run $ \hdr -> re_pair queues config (meta_refs hdr) =$ + mapChunks (maybe id do_trim (fixsven config)) =$ (output config) (add_pg hdr) @@ -591,4 +598,30 @@ bp_pos (Singleton u) = b_pos $ unpackBam u bp_pos (Pair u _) = b_pos $ unpackBam u bp_pos (LoneMate u) = b_pos $ unpackBam u + + +do_trim :: Int -> [BamRec] -> [BamRec] +do_trim q = scan_empties . map trim1 + where + trim1 b = case [ l | l <- [0 .. V.length (b_qual b) -1], avquallow (V.drop l qs) ] of + [ ] -> b + l:_ -> trim_3 l b + where + qs | isReversed b = V.reverse (b_qual b) + | otherwise = b_qual b + + scan_empties (x:y:z) + | b_qname x == b_qname y + = if V.null (b_qual x) || V.null (b_qual y) + then scan_empties z + else x : y : scan_empties z + + scan_empties (x:z) + = if V.null (b_qual x) + then scan_empties z + else x : scan_empties z + + scan_empties [] = [] + + avquallow vec = V.sum (V.map (fromIntegral . unQ) vec) <= q * V.length vec diff --git a/tools/bam-meld.hs b/tools/bam-meld.hs --- a/tools/bam-meld.hs +++ b/tools/bam-meld.hs @@ -18,7 +18,6 @@ import Bio.Iteratee import Control.Monad ( unless, foldM ) import Data.List ( sortBy ) -import Data.Monoid import Data.String ( fromString ) import Data.Version ( showVersion ) import Paths_biohazard ( version ) diff --git a/tools/bam-rmdup.hs b/tools/bam-rmdup.hs --- a/tools/bam-rmdup.hs +++ b/tools/bam-rmdup.hs @@ -2,14 +2,13 @@ import Bio.Bam import Bio.Bam.Rmdup import Bio.Base -import Bio.Util ( showNum, showOOM, estimateComplexity ) +import Bio.Util.Numeric ( showNum, showOOM, estimateComplexity ) import Control.Monad import Control.Monad.ST ( runST ) import Data.Bits import Data.Foldable ( toList ) import Data.List ( intercalate ) import Data.Maybe -import Data.Monoid ( mempty ) import Data.Ord ( comparing ) import Data.Vector.Algorithms.Intro ( sortBy ) import Data.Version ( showVersion ) @@ -46,7 +45,7 @@ circulars :: Refs -> IO (IM.IntMap (Seqid,Int), Refs) } -- | Which reference sequences to scan -data Which = All | Some Refseq Refseq | Unaln deriving Show +data Which = Allrefs | Some Refseq Refseq | Unaln deriving Show defaults :: Conf defaults = Conf { output = Nothing @@ -62,7 +61,7 @@ , get_label = get_library , putResult = putStr , debug = \_ -> return () - , which = All + , which = Allrefs , circulars = \rs -> return (IM.empty, rs) } options :: [OptDescr (Conf -> IO Conf)] @@ -106,7 +105,7 @@ set_multi c = return $ c { clean_multimap = clean_multi_flags } set_range a c - | a == "A" || a == "a" = return $ c { which = All } + | a == "A" || a == "a" = return $ c { which = Allrefs } | a == "U" || a == "u" = return $ c { which = Unaln } | otherwise = case reads a of [ (x,"") ] -> return $ c { which = Some (Refseq $ x-1) (Refseq $ x-1) } @@ -202,7 +201,7 @@ debug "mapping of read groups to libraries:\n" mapM_ debug [ unpackSeqid k ++ " --> " ++ unpackSeqid v ++ "\n" | (k,v) <- M.toList tbl ] - let filters = progressPos "Rmdup at " debug refs' ><> + let filters = progressBam "Rmdup at " debug refs' ><> mapChunks (mapMaybe (transform . unpackBam)) ><> mapChunksM (mapMM clean_multimap) ><> filterStream (\br -> (keep_unaligned || is_aligned br) && @@ -305,11 +304,11 @@ mergeInputRanges :: (MonadIO m, MonadMask m) => Which -> [FilePath] -> Enumerator' BamMeta [BamRaw] m a -mergeInputRanges All fps = mergeInputs combineCoordinates fps +mergeInputRanges Allrefs fps = mergeInputs combineCoordinates fps mergeInputRanges _ [ ] = \k -> return $ k mempty mergeInputRanges rng (fp0:fps0) = go fp0 fps0 where - enum1 fp k1 = case rng of All -> decodeAnyBamFile fp k1 + enum1 fp k1 = case rng of Allrefs -> decodeAnyBamFile fp k1 Some x y -> decodeBamFileRange x y fp k1 Unaln -> decodeWithIndex eneeBamUnaligned fp k1 diff --git a/tools/count-coverage.hs b/tools/count-coverage.hs deleted file mode 100644 --- a/tools/count-coverage.hs +++ /dev/null @@ -1,61 +0,0 @@ -{-# LANGUAGE BangPatterns, NoMonomorphismRestriction, FlexibleContexts #-} -import Bio.Bam.Header -import Bio.Bam.Reader -import Bio.Bam.Rec -import Bio.Base -import Bio.Iteratee -import Data.Version ( showVersion ) -import Paths_biohazard ( version ) -import System.Environment -import System.Exit -import System.IO ( hPutStr ) - -main :: IO () -main = do - mq <- getArgs >>= \args -> case (args, reads (head args)) of - ([ ], _) -> return (Q 0) - ([_], [(x,[])]) -> return (Q x) - _ -> do pn <- getProgName - hPutStr stderr $ pn ++ ", version " ++ showVersion version - ++ "\nUsage: " ++ pn ++ "[]\n" - exitFailure - - let putLine nm cv = putStr $ nm ++ '\t' : shows cv "\n" - - printOne :: Refs -> (Refseq, Int) -> IO () - printOne refs (r,c) = putLine (unpackSeqid (sq_name (getRef refs r))) c - - do_count :: Monad m => Iteratee [(a,Int)] m Int - do_count = foldStream (\a -> (+) a . snd) 0 - - (total,()) <- enumHandle defaultBufSize stdin >=> run $ - joinI $ decodeAnyBam $ \hdr -> - joinI $ mapMaybeStream ( \br -> case unpackBam br of - b | not (isUnmapped b) && b_mapq b >= mq - -> Just $! P (b_rname b) (b_pos b) (alignedLength (b_cigar b)) - _ -> Nothing ) $ - joinI $ groupStreamOn ref count_cov $ - zipStreams do_count (mapStreamM_ $ printOne $ meta_refs hdr) - - putLine "total" total - -data P = P { ref :: !Refseq, pos :: !Int, alen :: !Int } - -count_cov :: Monad m => a -> m (Iteratee [P] m Int) -count_cov _ = return $ liftI $ step 0 - where - step !a (EOF ex) = idone a (EOF ex) - step !a (Chunk [ ]) = liftI $ step a - step !a (Chunk (r:rs)) = extend a (pos r) (pos r + alen r) (Chunk rs) - - extend !a !u !v (EOF ex) = idone (a+v-u) (EOF ex) - extend !a !u !v (Chunk [ ]) = liftI $ extend a u v - extend !a !u !v (Chunk (r:rs)) - | pos r <= v = extend a u (max v (pos r + alen r)) (Chunk rs) - | otherwise = step (a+v-u) (Chunk (r:rs)) - - - - - - diff --git a/tools/dmg-est.hs b/tools/dmg-est.hs deleted file mode 100644 --- a/tools/dmg-est.hs +++ /dev/null @@ -1,369 +0,0 @@ -{-# LANGUAGE RecordWildCards, NamedFieldPuns, BangPatterns, TypeFamilies #-} --- Estimates aDNA damage. Crude first version. --- --- - Read or subsample a BAM file, make compact representation of the reads. --- - Compute likelihood of each read under simple model of --- damage, error/divergence, contamination. --- --- For the fitting, we simplify radically: ignore sequencing error, --- assume damage and simple, symmetric substitutions which subsume error --- and divergence. --- --- Trying to compute symbolically is too much, the high power terms get --- out of hand quickly, and we get mixed powers of \lambda and \kappa. --- The fastest version so far uses the cheap implementation of automatic --- differentiation in AD.hs together with the Hager-Zhang method from --- package nonlinear-optimization. BFGS from hmatrix-gsl takes longer --- to converge. Didn't try an actual Newton iteration (yet?), AD from --- package ad appears slower. --- --- If I include parameters, whose true value is zero, the transformation --- to the log-odds-ratio doesn't work, because then the maximum doesn't --- exist anymore. For many parameters, zero makes sense, but one --- doesn't. A different transformation ('sigmoid2'/'isigmoid2' --- below) allows for an actual zero (but not an actual one), while --- avoiding ugly boundary conditions. That appears to work well. --- --- The current hack assumes all molecules have an overhang at both ends, --- then each base gets deaminated with a position dependent probability --- following a geometric distribution. If we try to model a fraction of --- undeaminated molecules (a contaminant) in addition, this fails. To --- rescue the idea, I guess we must really decide if the molecule has an --- overhang at all (probability 1/2) at each end, then deaminate it. --- --- TODO --- - needs better packaging, better output --- - needs support for multiple input files(?) --- - needs read group awareness(?) --- - needs to deal with long (unmerged) reads (by ignoring them?) - -import Bio.Bam.Header -import Bio.Bam.Index -import Bio.Bam.Rec -import Bio.Base -import Bio.Genocall.Adna -import Bio.Iteratee -import Control.Concurrent.Async -import Data.Bits -import Data.Foldable -import Data.Ix -import Data.Maybe -import Numeric.Optimization.Algorithms.HagerZhang05 -import System.Environment - -import qualified Data.Vector as V -import qualified Data.Vector.Generic as G -import qualified Data.Vector.Unboxed as U - -import AD -import Prelude hiding ( sequence_, mapM, mapM_, concatMap, sum, minimum, foldr1 ) - --- | Roughly @Maybe (Nucleotide, Nucleotide)@, encoded compactly -newtype NP = NP { unNP :: Word8 } deriving (Eq, Ord, Ix) -data Seq = Merged { unSeq :: U.Vector Word8 } - | First { unSeq :: U.Vector Word8 } - | Second { unSeq :: U.Vector Word8 } - -instance Show NP where - show (NP w) - | w == 16 = "NN" - | w > 16 = "XX" - | otherwise = [ "ACGT" !! fromIntegral (w `shiftR` 2) - , "ACGT" !! fromIntegral (w .&. 3) ] - - -sigmoid2, isigmoid2 :: (Num a, Fractional a, Floating a) => a -> a -sigmoid2 x = y*y where y = (exp x - 1) / (exp x + 1) -isigmoid2 y = log $ (1 + sqrt y) / (1 - sqrt y) - -{-# INLINE lk_fun1 #-} -lk_fun1 :: (Num a, Show a, Fractional a, Floating a, Memorable a) => Int -> [a] -> V.Vector Seq -> a -lk_fun1 lmax parms = case length parms of - 1 -> V.foldl' (\a b -> a - log (lk tab00 tab00 tab00 b)) 0 . guardV -- undamaged case - where - !tab00 = fromListN (rangeSize my_bounds) [ l_epq p_subst 0 0 x - | (_,_,x) <- range my_bounds ] - - 4 -> V.foldl' (\a b -> a - log (lk tabDS tabDS1 tabDS1 b)) 0 . guardV -- double strand case - where - !tabDS = fromListN (rangeSize my_bounds) [ l_epq p_subst p_d p_e x - | (l,i,x) <- range my_bounds - , let p_d = mu $ lambda ^^ (1+i) - , let p_e = mu $ lambda ^^ (l-i) ] - - !tabDS1 = fromListN (rangeSize my_bounds) [ l_epq p_subst p_d 0 x - | (_,i,x) <- range my_bounds - , let p_d = mu $ lambda ^^ (1+i) ] - - 5 -> V.foldl' (\a b -> a - log (lk tabSS tabSS1 tabSS2 b)) 0 . guardV -- single strand case - where - !tabSS = fromListN (rangeSize my_bounds) [ l_epq p_subst p_d 0 x - | (l,i,x) <- range my_bounds - , let lam5 = lambda ^^ (1+i) ; lam3 = kappa ^^ (l-i) - , let p_d = mu $ lam3 + lam5 - lam3 * lam5 ] - - !tabSS1 = fromListN (rangeSize my_bounds) [ l_epq p_subst p_d 0 x - | (_,i,x) <- range my_bounds - , let p_d = mu $ lambda ^^ (1+i) ] - - !tabSS2 = fromListN (rangeSize my_bounds) [ l_epq p_subst 0 p_d x - | (_,i,x) <- range my_bounds - , let p_d = mu $ lambda ^^ (1+i) ] - - _ -> error "Not supposed to happen: unexpected number of model parameters." - where - ~(l_subst : ~(l_sigma : ~(l_delta : ~(l_lam : ~(l_kap : _))))) = parms - - p_subst = 0.33333 * sigmoid2 l_subst - sigma = sigmoid2 l_sigma - delta = sigmoid2 l_delta - lambda = sigmoid2 l_lam - kappa = sigmoid2 l_kap - - guardV = V.filter (\u -> U.length (unSeq u) >= lmin && U.length (unSeq u) <= lmax) - - -- 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 _ (First b) = U.ifoldl' (\a i np -> a * tab_f `bang` index' my_bounds (U.length b, i, NP np)) 1 b - lk _ _ tab_s (Second 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 - | otherwise = error $ "Huh? " ++ show x ++ " \\nin " ++ show bnds - - my_bounds = ((lmin,0,NP 0),(lmax,lmax,NP 16)) - mu p = sigma * p + delta * (1-p) - - --- Likelihood for a certain pair of bases given error rate, C-T-rate --- and G-A rate. -l_epq :: (Num a, Fractional a, Floating a) => a -> a -> a -> NP -> a -l_epq e p q (NP x) = case x of { - 0 -> s ; 1 -> e ; 2 -> e ; 3 -> e ; - 4 -> e ; 5 -> s-p+4*e*p ; 6 -> e ; 7 -> e+p-4*e*p ; - 8 -> e+q-4*e*q ; 9 -> e ; 10 -> s-q+4*e*q ; 11 -> e ; - 12 -> e ; 13 -> e ; 14 -> e ; 15 -> s ; - _ -> 1 } where s = 1 - 3 * e - - -lkfun :: Int -> V.Vector Seq -> U.Vector Double -> Double -lkfun lmax brs parms = lk_fun1 lmax (U.toList parms) brs - -combofn :: Int -> V.Vector Seq -> U.Vector Double -> (Double, U.Vector Double) -combofn lmax brs parms = (x,g) - where D x g = lk_fun1 lmax (paramVector $ U.toList parms) brs - -params :: Parameters -params = defaultParameters { printFinal = False, verbose = Quiet, maxItersFac = 20 } - -lmin :: Int -lmin = 25 - -main :: IO () -main = do - [fp] <- getArgs - brs <- subsampleBam fp >=> run $ \_ -> - joinI $ filterStream (\b -> not (isUnmapped (unpackBam b)) && G.length (b_seq (unpackBam b)) >= lmin) $ - joinI $ takeStream 100000 $ - joinI $ mapStream pack_record $ - joinI $ filterStream (\u -> U.length (U.filter (<16) (unSeq u)) * 10 >= 9 * U.length (unSeq u)) $ - stream2vectorN 30000 - - let lmax = V.maximum $ V.map (U.length . unSeq) brs - v0 = crude_estimate brs - opt v = optimize params 0.0001 v - (VFunction $ lkfun lmax brs) - (VGradient $ snd . combofn lmax brs) - (Just . VCombined $ combofn lmax brs) - - results <- mapConcurrently opt [ v0, U.take 4 v0, U.take 1 v0 ] - - let mlk = minimum [ finalValue st | (_,_,st) <- results ] - tot = sum [ exp $ mlk - finalValue st | (_,_,st) <- results ] - p l = exp (mlk - l) / tot - - [ (p_ss, [ _, ssd_sigma_, ssd_delta_, ssd_lambda, ssd_kappa ]), - (p_ds, [ _, dsd_sigma_, dsd_delta_, dsd_lambda ]), - (_ , [ _ ]) ] = [ (p (finalValue st), map sigmoid2 $ G.toList xs) | (xs,_,st) <- results ] - - ssd_sigma = p_ss * ssd_sigma_ - ssd_delta = p_ss * ssd_delta_ - dsd_sigma = p_ds * dsd_sigma_ - dsd_delta = p_ds * dsd_delta_ - - print DP{..} - --- We'll require the MD field to be present. Then we cook each read --- into a list of paired bases. Deleted bases are dropped, inserted --- bases replaced with an escape code. --- --- XXX This is annoying... almost, but not quite the same as the code --- in the "Pileup" module. This also relies on MD and doesn't offer the --- alternative of accessing a reference genome. (The latter may not be --- worth the trouble.) It also resembles the 'ECig' logic from --- "Bio.Bam.Rmdup". - -pack_record :: BamRaw -> Seq -pack_record br = if isReversed b then k (revcom u1) else k u1 - where - b@BamRec{..} = unpackBam br - - k | isMerged b = Merged - | isTrimmed b = Merged - | isSecondMate b = Second - | otherwise = First - - revcom = U.reverse . U.map (\x -> if x > 15 then x else xor x 15) - u1 = U.fromList . map unNP $ go (G.toList b_cigar) (G.toList b_seq) (fromMaybe [] $ getMd b) - - go :: [Cigar] -> [Nucleotides] -> [MdOp] -> [NP] - - go (_:*0 :cs) ns mds = go cs ns mds - go cs ns (MdNum 0:mds) = go cs ns mds - go cs ns (MdDel []:mds) = go cs ns mds - go _ [] _ = [] - go [] _ _ = [] - - go (Mat:*nm :cs) (n:ns) (MdNum mm:mds) = mk_pair n n : go (Mat:*(nm-1):cs) ns (MdNum (mm-1):mds) - go (Mat:*nm :cs) (n:ns) (MdRep n':mds) = mk_pair n n' : go (Mat:*(nm-1):cs) ns mds - go (Mat:*nm :cs) ns (MdDel _ :mds) = go (Mat:* nm :cs) ns mds - - go (Ins:*nm :cs) ns mds = replicate nm esc ++ go cs (drop nm ns) mds - go (SMa:*nm :cs) ns mds = replicate nm esc ++ go cs (drop nm ns) mds - go (Del:*nm :cs) ns (MdDel (_:ds):mds) = go (Del:*(nm-1):cs) ns (MdDel ds:mds) - go (Del:*nm :cs) ns ( _:mds) = go (Del:* nm :cs) ns mds - - go (_:cs) nd mds = go cs nd mds - - -esc :: NP -esc = NP 16 - -mk_pair :: Nucleotides -> Nucleotides -> NP -mk_pair (Ns a) = case a of 1 -> mk_pair' 0 - 2 -> mk_pair' 1 - 4 -> mk_pair' 2 - 8 -> mk_pair' 3 - _ -> const esc - where - mk_pair' u (Ns b) = case b of 1 -> NP $ u .|. 0 - 2 -> NP $ u .|. 4 - 4 -> NP $ u .|. 8 - 8 -> NP $ u .|. 12 - _ -> esc - - -infix 7 /%/ -(/%/) :: Integral a => a -> a -> Double -0 /%/ 0 = 0 -a /%/ b = fromIntegral a / fromIntegral b - --- Crude estimate. Need two overhang lengths, two deamination rates, --- undamaged fraction, SS/DS, substitution rate. --- --- DS or SS: look whether CT or GA is greater at 3' terminal position √ --- Left overhang length: ratio of damage at second position to first √ --- Right overang length: ratio of CT at last to snd-to-last posn √ --- + ratio of GA at last to snd-to-last posn √ --- SS rate: condition on damage on one end, compute rate at other √ --- DS rate: condition on damage, compute rate in interior √ --- substitution rate: count all substitutions not due to damage √ --- undamaged fraction: see below √ --- --- Contaminant fraction: let f5 (f3, f1) be the fraction of reads --- showing damage at the 5' end (3' end, both ends). Let a (b) be --- the probability of an endogenous reads to show damage at the 5' --- end (3' end). Let e be the fraction of endogenous reads. Then --- we have: --- --- f5 = e * a --- f3 = e * b --- f1 = e * a * b --- --- f5 * f3 / f1 = e --- --- Straight forward and easy to understand, but in practice, this method --- produces ridiculous overestimates, ridiculous underestimates, --- negative contamination rates, and general grief. It's actually --- better to start from a constant number. - - -crude_estimate :: V.Vector Seq -> U.Vector Double -crude_estimate seqs0 = U.fromList [ l_subst, l_sigma, l_delta, l_lam, l_kap ] - where - seqs = V.filter ((>= 10) . U.length) $ V.map unSeq seqs0 - - total_equals = V.sum (V.map (U.length . U.filter isNotSubst) seqs) - total_substs = V.sum (V.map (U.length . U.filter isOrdinarySubst) seqs) * 6 `div` 5 - l_subst = isigmoid2 $ max 0.001 $ total_substs /%/ (total_equals + total_substs) - - c_to_t, g_to_a, c_to_c :: Word8 - c_to_t = 7 - g_to_a = 8 - c_to_c = 5 - - isNotSubst x = x < 16 && x `shiftR` 2 == x .&. 3 - isOrdinarySubst x = x < 16 && x `shiftR` 2 /= x .&. 3 && - x /= c_to_t && x /= g_to_a - - ct_at_alpha = V.length $ V.filter (\v -> v U.! 0 == c_to_t && dmg_omega v) seqs - cc_at_alpha = V.length $ V.filter (\v -> v U.! 0 == c_to_c && dmg_omega v) seqs - ct_at_beta = V.length $ V.filter (\v -> v U.! 1 == c_to_t && dmg_omega v) seqs - cc_at_beta = V.length $ V.filter (\v -> v U.! 1 == c_to_c && dmg_omega v) seqs - - dmg_omega v = v U.! (l-1) == c_to_t || v U.! (l-1) == g_to_a - || v U.! (l-2) == c_to_t || v U.! (l-2) == g_to_a - || v U.! (l-3) == c_to_t || v U.! (l-3) == g_to_a - where l = U.length v - - l_lam = isigmoid2 lambda - lambda = min 0.9 $ max 0.1 $ - (ct_at_beta * (cc_at_alpha + ct_at_alpha)) /%/ - ((cc_at_beta + ct_at_beta) * ct_at_alpha) - - ct_at_omega = V.length $ V.filter (\v -> v U.! (U.length v -1) == c_to_t && dmg_alpha v) seqs - cc_at_omega = V.length $ V.filter (\v -> v U.! (U.length v -1) == c_to_c && dmg_alpha v) seqs - ct_at_psi = V.length $ V.filter (\v -> v U.! (U.length v -2) == c_to_t && dmg_alpha v) seqs - cc_at_psi = V.length $ V.filter (\v -> v U.! (U.length v -2) == c_to_c && dmg_alpha v) seqs - - dmg_alpha v = v U.! 0 == c_to_t || v U.! 1 == c_to_t || v U.! 2 == c_to_t - - l_kap = isigmoid2 $ min 0.9 $ max 0.1 $ - (ct_at_psi * (cc_at_omega+ct_at_omega)) /%/ - ((cc_at_psi+ct_at_psi) * ct_at_omega) - - total_inner_CCs = V.sum $ V.map (U.length . U.filter (== c_to_c) . takeInner) seqs - total_inner_CTs = V.sum $ V.map (U.length . U.filter (== c_to_t) . takeInner) seqs - takeInner v = U.slice 5 (U.length v - 10) v - - delta = (total_inner_CTs /%/ (total_inner_CTs+total_inner_CCs)) - raw_rate = ct_at_alpha /%/ (ct_at_alpha + cc_at_alpha) - - -- clamping is necessary if f_endo ends up wrong - l_delta = isigmoid2 $ min 0.99 delta - l_sigma = isigmoid2 . min 0.99 $ raw_rate / lambda - - -class Memorable a where - type Memo a :: * - - fromListN :: Int -> [a] -> Memo a - bang :: Memo a -> Int -> a - -instance Memorable Double where - type Memo Double = U.Vector Double - - fromListN = U.fromListN - bang = (U.!) - -instance Memorable AD where - type Memo AD = (Int, U.Vector Double) - - fromListN n xs@(D _ v:_) = (1+d, U.fromListN (n * (1+d)) $ concatMap unpack xs) - where - !d = U.length v - unpack (C a) = a : replicate d 0 - unpack (D a da) = a : U.toList da - - bang (d, v) i = D (v U.! (d*i+0)) (U.slice (d*i+1) (d-1) v) diff --git a/tools/fastq2bam.hs b/tools/fastq2bam.hs --- a/tools/fastq2bam.hs +++ b/tools/fastq2bam.hs @@ -5,7 +5,6 @@ import Bio.Iteratee.ZLib import Control.Monad import Data.Bits -import Data.Monoid ( mempty ) import System.Console.GetOpt import System.Environment import System.Exit diff --git a/tools/glf-consensus.hs b/tools/glf-consensus.hs deleted file mode 100644 --- a/tools/glf-consensus.hs +++ /dev/null @@ -1,205 +0,0 @@ -{-# LANGUAGE BangPatterns #-} -import Control.Applicative ( (<$>) ) -import Control.Monad -import Control.Monad.Catch -import Data.Char ( isSpace, toLower, chr ) -import Data.List ( intercalate, sort ) -import Data.Version ( showVersion ) -import Paths_biohazard ( version ) -import System.Console.GetOpt -import System.IO -import System.Environment ( getArgs, getProgName ) -import System.Exit - -import qualified Data.ByteString.Char8 as S -import qualified Data.ByteString.Lazy.Char8 as L -import qualified Data.Map as M - -import qualified Data.Iteratee.ListLike as I - -import Bio.Base -import Bio.Glf -import Bio.Iteratee - -data Config = Config { - conf_min_qual :: Int, - conf_call :: [Int] -> [(Int, Char)], - conf_output :: Iteratee String IO (), - conf_input :: GlfInput, - conf_conv :: Formatter, - conf_mkname :: S.ByteString -> String } - -type GlfInput = (GlfSeq -> Enumeratee [GlfRec] String IO ()) - -> (S.ByteString -> Enumerator String IO ()) - -> Enumerator String IO () - -options :: [ OptDescr (Config -> IO Config) ] -options = [ - Option "1" ["haploid"] - (NoArg (\c -> return $ c { conf_call = haploid_call })) - "Force haploid consensus", - Option "2" ["diploid"] - (NoArg (\c -> return $ c { conf_call = diploid_call })) - "Allow diploid consensus", - Option "m" ["min-qual"] - (ReqArg (\a c -> readIO a >>= \m -> return $ c { conf_min_qual = m }) "Q") - "Require minimum quality of Q", - Option "o" ["output"] - (ReqArg (\fp c -> return $ c { conf_output = iterToFile fp }) "FILE") - "Write output to FILE instead of stdout", - Option "q" ["fastq"] - (NoArg (\c -> return $ c { conf_conv = print_fastq })) - "Write FastQ instead of FastA", - Option "I" ["identifier"] - (ReqArg (\n c -> return $ c { conf_mkname = subst_name n }) "ID") - "Use ID as identifier for consensus", - Option "if" ["input"] - (ReqArg (\fp c -> return $ c { conf_input = enum_glf_file fp }) "FILE") - "Read input from FILE instead of stdin", - Option "h?" ["help", "usage"] - (NoArg (usage exitSuccess)) - "Print this help", - Option "V" ["version"] - (NoArg vrsn) - "Print version number and exit" ] - -vrsn :: Config -> IO a -vrsn _ = do pn <- getProgName - hPutStrLn stderr $ pn ++ ", version " ++ showVersion version - exitSuccess - -usage :: IO a -> Config -> IO a -usage e _ = getProgName >>= \p -> putStrLn (usageInfo (blurb p) options) >> e - where blurb prg = - "Usage: " ++ prg ++ " [Option...] [FastA-File...]\n" ++ - "Reads GLF from stdin and prints the contained consensus sequence in\n" ++ - "FastA/FastQ format. Gaps are filled with a reference sequence if known\n" ++ - "from the FastA files on the command line, otherwise with Ns." - -iterToFile :: FilePath -> Iteratee String IO () -iterToFile fp = bracket (lift $ openFile fp WriteMode) - (lift . hClose) - (mapChunksM_ . hPutStr) - -defaultConfig :: Config -defaultConfig = Config 0 diploid_call (mapChunksM_ putStr) (enum_glf_handle stdin) print_fasta S.unpack - -main :: IO () -main = do (opts, files, errors) <- getOpt Permute options <$> getArgs - unless (null errors) $ mapM_ (hPutStrLn stderr) errors >> exitFailure - Config min_qual call output input conv mkname <- foldl (>>=) (return defaultConfig) opts - refs <- M.fromList . concatMap readFasta <$> mapM L.readFile files - - hPutStrLn stderr $ - "known reference sequences: [" ++ intercalate ", " - [ show (L.unpack k) ++ " (" ++ show (L.length v) ++ ")" | (k,v) <- M.toList refs ] - ++ "]" - - let per_file :: Seqid -> Enumerator String IO () - per_file _genome_name = return - - per_seq :: GlfSeq -> Enumeratee [GlfRec] String IO () - per_seq glfseq = extract1consensus (mkRef refs glfseq) call min_qual - ><> conv (mkname $ glf_seqname glfseq) - - input per_seq per_file output >>= run - --- get the "most likely" consensus, defined as: --- - as many reference bases or else Ns as were skipped from the previous record, then --- - if there's an insert, the most likely insert sequence (may be empty) --- - if there's a deletion, skip the most likely number of bases (may be zero) --- - else the most likely base - -mkRef :: M.Map L.ByteString L.ByteString -> GlfSeq -> Int -> Int -> QSeq -mkRef refs glfseq = case M.lookup (L.fromChunks [glf_seqname glfseq]) refs of - Nothing -> \o l -> replicate (min l (glf_seqlen glfseq - o)) ('N',2) - Just s -> \o l -> let l' = fromIntegral $ min l (glf_seqlen glfseq - o) - in [ (toLower b,30) | b <- L.unpack $ L.take l' $ L.drop (fromIntegral o) s ] - -type QSeq = [(Char,Int)] -- sequence w/ quality - -extract1consensus :: Monad m - => (Int -> Int -> QSeq) - -> ([Int] -> [(Int,Char)]) -- call function - -> Int -- minimum quality - -> Enumeratee [GlfRec] QSeq m r -- eats records, emits calls -extract1consensus ref call min_qual oit = liftI $ scan oit 0 0 - where - -- rec_pos: position of last record - -- ref_pos: first position in reference we haven't handled - scan k !_ !ref_pos (EOF x) = lift $ enumPure1Chunk (ref ref_pos maxBound) >=> enumChunk (EOF x) $ k - scan k !rec_pos_ !ref_pos (Chunk [ ]) = liftI $ scan k rec_pos_ ref_pos - scan k !rec_pos_ !ref_pos (Chunk (r:rs)) = - case r of SNP {} -> let (_,!base) : (!qual,_) : _ = sort $ call (glf_lk r) - in ( if qual >= min_qual - then lift $ enumPure1Chunk (ref ref_pos (rec_pos - ref_pos)) k - >>= enumPure1Chunk [(base,qual)] - else lift $ enumPure1Chunk (ref ref_pos (1 + rec_pos - ref_pos)) k ) - >>= \k' -> scan k' rec_pos (1+rec_pos) (Chunk rs) - - Indel {} | ins && iqual >= min_qual -> lift (enumPure1Chunk (ref ref_pos (rec_pos + 1 - ref_pos)) k >>= - enumPure1Chunk [ (b,iqual) | b <- S.unpack sq ]) >>= \k'' -> - scan k'' rec_pos ref_pos (Chunk rs) - | not ins && iqual >= min_qual -> lift (enumPure1Chunk (ref ref_pos (rec_pos - ref_pos)) k) >>= \k' -> - scan k' rec_pos (ref_pos + S.length sq) (Chunk (drop (S.length sq) rs)) - | otherwise -> lift (enumPure1Chunk (ref ref_pos (rec_pos - ref_pos)) k) >>= \k' -> - scan k' rec_pos ref_pos (Chunk rs) - where - !rec_pos = rec_pos_ + glf_offset r - (ins,sq) = if glf_lk_hom1 r > glf_lk_hom2 r - then (glf_is_ins2 r, glf_seq2 r) else (glf_is_ins1 r, glf_seq1 r) - iqual = abs $ glf_lk_hom1 r - glf_lk_hom2 r - - -diploid_call, haploid_call :: [Int] -> [(Int, Char)] -diploid_call lks = zip lks "AMRWCSYGKT" -haploid_call lks = zip (map (lks !!) [0,4,7,9]) "ACGT" - - -type Formatter = String -> Enumeratee QSeq String IO () - -print_fasta :: Formatter -print_fasta name = eneeCheckIfDone (\k -> mapStream fst ><> toLines 60 $ k $ Chunk ('>' : name ++ "\n")) - -print_fastq :: Formatter -print_fastq name = eneeCheckIfDone p'header - where - p'header k = p'seq . k $ Chunk ('@' : name ++ "\n") - p'seq it = I.zip ((mapStream fst ><> toLines 60) it) (liftI $ coll []) - >>= \(it', qs) -> eneeCheckIfDone (p'sep qs) it' - p'sep qs k = lift $ (enumList (map S.unpack qs) >=> run) (toLines 60 . k $ Chunk "+\n") - - mkqual = chr . max 33 . min 126 . (+) 33 . fromIntegral - coll !acc (EOF x) = lift (putStrLn $ show $ length acc) >> idone (reverse acc) (EOF x) - coll !acc (Chunk []) = liftI $ coll acc - coll !acc (Chunk c) = liftI . coll $! norm (S.pack (map (mkqual . snd) c)) acc - - -- ensure that we don't build many small ByteStrings - norm !x [] = [x] - norm !x (y:ys) | S.length x > S.length y = norm (y `S.append` x) ys - | otherwise = x:y:ys - - -toLines :: Monad m => Int -> Enumeratee String String m r -toLines n = eneeCheckIfDone (\k -> I.isFinished >>= go k) - where - go k True = return $ liftI k - go k False = do s <- I.take n I.stream2list >>= lift . run - eneeCheckIfDone (\k1 -> toLines n . k1 $ Chunk "\n") . k $ Chunk s - - -readFasta :: L.ByteString -> [(L.ByteString, L.ByteString)] -readFasta = rd . dropWhile (not . isHeader) . L.lines - where - isHeader l = not (L.null l) && L.head l == '>' - rd [] = [] - rd (l:ls) = let name = L.takeWhile (not . isSpace) $ L.drop 1 l - (sqs,rest) = break isHeader ls - in (name, L.filter (`elem` "ACGTBDHVSWMKYRNU") $ L.concat sqs) : rd rest - -subst_name :: String -> S.ByteString -> String -subst_name [] _ = [] -subst_name ('%':'s':t) s = S.unpack s ++ subst_name t s -subst_name ('%':'%':t) s = '%' : subst_name t s -subst_name (t:ts) s = t : subst_name ts s - diff --git a/tools/gt-call.hs b/tools/gt-call.hs deleted file mode 100644 --- a/tools/gt-call.hs +++ /dev/null @@ -1,392 +0,0 @@ -{-# LANGUAGE RecordWildCards, BangPatterns, OverloadedStrings #-} -{-# LANGUAGE TemplateHaskell, FlexibleContexts #-} --- Command line driver for simple genotype calling. - -import Bio.Base -import Bio.Bam.Header -import Bio.Bam.Reader -import Bio.Bam.Rec -import Bio.Bam.Pileup -import Bio.Genocall -import Bio.Genocall.Adna -import Bio.Genocall.AvroFile -import Bio.Iteratee -import Bio.Util ( float2mini ) -import Control.Applicative -import Control.DeepSeq -import Control.Monad -import Data.Avro -import Data.Function -import System.Console.GetOpt -import System.Environment -import System.Exit -import System.IO - --- import qualified Data.ByteString as B -import qualified Data.ByteString.Char8 as S -import qualified Data.Iteratee as I --- import qualified Data.Text as T -import qualified Data.Text.Encoding as T -import qualified Data.Vector.Unboxed as V - --- import Debug.Trace - --- Ultimately, we might produce a VCF file looking somewhat like this: --- --- ##FORMAT= --- ##FORMAT= --- ##FORMAT= --- ##FORMAT= --- (we should count bases on both strands for this) --- --- ##FORMAT= --- ##INFO= --- ##INFO= --- (basic statistics. we keep these) --- --- ##FORMAT= --- ##FORMAT= --- ##INFO= --- (this is bullshit) --- --- ##FORMAT= --- ##FORMAT= --- ##FORMAT= --- (these are straight forward to compute?) --- --- ##INFO= --- ##INFO= --- ##INFO= --- ##INFO= --- ##INFO= --- ##INFO= --- ##INFO= --- ##INFO= --- ##INFO= --- ##INFO= --- ##INFO= --- ##INFO= --- ##INFO= --- ##INFO= --- ##INFO= --- ##INFO= --- ##INFO= --- ##INFO= --- ##INFO= --- ##INFO= --- (this is external, will not be generated) --- --- ##INFO= --- ##INFO= --- (this is computable, isn't it?!) --- --- ##INFO= --- (this is from VarScan 2, a program that uses fixed cutoffs. It --- is not clear that this has any use at all.) --- --- ##INFO= --- ##INFO= --- ##INFO= --- ##INFO= --- ##INFO= --- ##INFO= --- ##INFO= --- ##INFO= --- ##INFO= --- ##INFO= --- ##INFO= --- (WTF?) - --- parameters used for the Unified Genotyper: --- downsample_to_coverage=250 --- heterozygosity=0.001 --- pcr_error_rate=1.0E-4 --- indel_heterozygosity=1.25E-4 - - --- auxilliary files (from Martin's option parser): --- --- ancestor_path EMF /mnt/expressions/martin/sequence_db/epo/epo_6_primate_v64/split/ --- G1000 VCF /mnt/expressions/martin/sequence_db/snps/20110521_G1000_release/phase1_intergrated_calls.20101123.snps_indels_svs.sites.vcf.gz --- bscores TSV1i /mnt/454/Altaiensis/users/martin/HighCoverage/additional_information/bscores/liftover/human.tsv.gz --- mammalscores TSV2f /mnt/454/Altaiensis/users/martin/HighCoverage/additional_information/mammal_conservation/liftover/human.tsv.gz --- primatescores TSV2f /mnt/454/Altaiensis/users/martin/HighCoverage/additional_information/primate_conservation/liftover/human.tsv.gz --- gerpscores TSV2f /mnt/454/Altaiensis/users/fernando/sequencedb/GERP/liftover/human.tsv.gz --- mapability TSV2i /mnt/454/Altaiensis/users/martin/HighCoverage/additional_information/mapability/liftover/human.tsv.gz --- uregions TSV1 /mnt/454/Altaiensis/users/martin/HighCoverage/additional_information/EL_control_regions/liftover/human.tsv.gz --- syserrors TSV1 /mnt/454/Altaiensis/users/martin/HighCoverage/additional_information/sys_errors/liftover/human.tsv.gz --- syserrorsHCB TSV1 /mnt/454/Altaiensis/users/fernando/sequencedb/SysErrHCB/human.tsv.gz - --- TSV1: chr start end score --- TSV2: chr pos score - --- About damage parameters: We effectively have three different models --- (SS, DS, no damage) and it may not be possible to choose one a --- priori. To manage this cleanly, we should have one universal model, --- but the three we have are not generalizations of each other. --- --- So we treat the choice of model as another parameter. We feed --- parameters for all three in, together with probabilities for each. --- Said probabilities are derived from the likelihoods obtained when --- fitting the parameters individually. Genotype calling then involves --- calling once under each model and summing them (effectively --- marginalizing on the choice of model). - -data Conf = Conf { - conf_output :: Maybe Output, - conf_sample :: S.ByteString, - conf_ploidy :: S.ByteString -> Int, - conf_damage :: Maybe (DamageParameters Double), - conf_loverhang :: Maybe Double, - conf_roverhang :: Maybe Double, - conf_ds_deam :: Double, - conf_ss_deam :: Double, - conf_theta :: Maybe Double, - conf_report :: String -> IO (), - conf_prior_het :: Prob, - conf_prior_indel :: Prob } - -defaultConf :: Conf -defaultConf = Conf Nothing "John_Doe" (const 2) Nothing Nothing Nothing - 0.02 0.45 Nothing (\_ -> return ()) - (qualToProb $ Q 30) (qualToProb $ Q 45) - -options :: [OptDescr (Conf -> IO Conf)] -options = [ - Option "o" ["output", "avro-output"] (ReqArg set_avro_out "FILE") "Write AVRO output to FILE", - Option [ ] ["fasta-output"] (ReqArg set_fa_output "FILE") "Write FA output to FILE", - Option "N" ["name","sample-name"] (ReqArg set_sample "NAME") "Set sample name to NAME", - Option "1" ["haploid-chromosomes"] (ReqArg set_haploid "PRF") "Targets starting with PRF are haploid", - Option "2" ["diploid-chromosomes"] (ReqArg set_diploid "PRF") "Targets starting with PRF are diploid", - Option "D" ["damage"] (ReqArg set_damage "PARMS") "Set universal damage parameters", - Option "l" ["overhang-param","left-overhang-param"] - (ReqArg set_loverhang "PROB") "Parameter for 5' overhang length is PROB", - Option "r" ["right-overhang-param"] (ReqArg set_roverhang "PROB") "Parameter for 3' overhang length is PROB, assume single-strand prep", - Option "d" ["deamination-rate","ds-deamination-rate","double-strand-deamination-rate"] - (ReqArg set_ds_deam "FRAC") "Deamination rate in double stranded section is FRAC", - Option "s" ["ss-deamination-rate","single-strand-deamination-rate"] - (ReqArg set_ss_deam "FRAC") "Deamination rate in single stranded section is FRAC", - Option "t" ["theta","dependency-coefficient"] - (ReqArg set_theta "FRAC") "Set dependency coefficient to FRAC (\"N\" to turn off)", - Option "H" ["prior-heterozygous", "heterozygosity"] - (ReqArg set_phet "PROB") "Set prior for a heterozygous variant to PROB", - -- Removed this, because it needs access to a reference. - -- But maybe we can derive this from a suitable BAM file? - -- Or move it to another tool? - -- Option "S" ["prior-snp","snp-rate","divergence"] - -- (ReqArg set_pdiv "PROB") "Set prior for an indel variant to PROB", - Option "I" ["prior-indel","indel-rate"] (ReqArg set_pindel "PROB") "Set prior for an indel variant to PROB", - Option "v" ["verbose"] (NoArg be_verbose) "Print more diagnostics", - Option "h?" ["help","usage"] (NoArg disp_usage) "Display this message" ] - where - disp_usage _ = do pn <- getProgName - let blah = "Usage: " ++ pn ++ " [OPTION...] [BAM-FILE...]" - putStrLn $ usageInfo blah options - exitFailure - - be_verbose c = return $ c { conf_report = hPutStrLn stderr } - - set_fa_output fn = add_output $ output_fasta fn - set_avro_out fn = add_output $ output_avro fn - - add_output ofn cf = - return $ cf { conf_output = Just $ \k -> - ofn $ \oit1 -> maybe (k oit1) ($ \oit2 -> k (\c r -> () <$ I.zip (oit1 c r) (oit2 c r))) (conf_output cf) } - - set_sample nm c = return $ c { conf_sample = S.pack nm } - - set_haploid arg c = return $ c { conf_ploidy = \chr -> if S.pack arg `S.isPrefixOf` chr then 1 else conf_ploidy c chr } - set_diploid arg c = return $ c { conf_ploidy = \chr -> if S.pack arg `S.isPrefixOf` chr then 2 else conf_ploidy c chr } - - set_theta "N" c = return $ c { conf_theta = Nothing } - set_theta a c = (\t -> c { conf_theta = Just t }) <$> readIO a - set_loverhang a c = (\l -> c { conf_loverhang = Just l }) <$> readIO a - set_roverhang a c = (\l -> c { conf_roverhang = Just l }) <$> readIO a - set_ss_deam a c = (\r -> c { conf_ss_deam = r }) <$> readIO a - set_ds_deam a c = (\r -> c { conf_ds_deam = r }) <$> readIO a - set_phet a c = (\r -> c { conf_prior_het = toProb r }) <$> readIO a - set_pindel a c = (\r -> c { conf_prior_indel = toProb r }) <$> readIO a - set_damage a c = (\u -> c { conf_damage = Just u }) <$> readIO a - -main :: IO () -main = do - (opts, files, errs) <- getOpt Permute options <$> getArgs - unless (null errs) $ mapM_ (hPutStrLn stderr) errs >> exitFailure - conf@Conf{..} <- foldl (>>=) (return defaultConf) opts - - let no_damage = conf_report "using no damage model" >> return noDamage - ss_damage p = conf_report ("using single strand damage model with " ++ show p) >> return (univDamage p) - ds_damage p = conf_report ("using double strand damage model with " ++ show p) >> return (univDamage p) - u_damage p = conf_report ("using universal damage parameters " ++ show p) >> return (univDamage p) - - dmg_model <- case (conf_damage, conf_loverhang, conf_roverhang) of - (Just u, _, _) -> u_damage u - (_, Nothing, Nothing) -> no_damage - (_, Just pl, Nothing) -> ds_damage $ DP 0 0 0 0 conf_ss_deam conf_ds_deam pl - (_, Nothing, Just pr) -> ss_damage $ DP conf_ss_deam conf_ds_deam pr pr 0 0 0 - (_, Just pl, Just pr) -> ss_damage $ DP conf_ss_deam conf_ds_deam pl pr 0 0 0 - - maybe (output_fasta "-") id conf_output $ \oiter -> - mergeInputs combineCoordinates files >=> run $ \hdr -> - filterStream ((\b -> not (isUnmapped b) && isValidRefseq (b_rname b)) . unpackBam) =$ - progressPos "GT call at " conf_report (meta_refs hdr) =$ - by_groups ((==) `on` b_rname . unpackBam) (\br out -> do - let sname = sq_name $ getRef (meta_refs hdr) $ b_rname $ unpackBam br - pl = conf_ploidy sname - liftIO $ conf_report $ S.unpack sname ++ ["",": haploid call",": diploid call"] !! pl - pileup dmg_model =$ mapStream (calls conf_theta pl) out) =$ - oiter conf (meta_refs hdr) - - -type OIter = Conf -> Refs -> Iteratee [Calls] IO () -type Output = (OIter -> IO ()) -> IO () - -output_fasta :: FilePath -> (OIter -> IO r) -> IO r -output_fasta fn k = if fn == "-" then k (fa_out stdout) - else withFile fn WriteMode $ k . fa_out - where - fa_out :: Handle -> Conf -> Refs -> Iteratee [Calls] IO () - fa_out hdl Conf{..} refs = - by_groups ((==) `on` p_refseq) (\cs out -> do - let sname = sq_name $ getRef refs $ p_refseq cs - out' <- lift $ enumPure1Chunk [S.concat [">", conf_sample, "--", sname]] out - convStream (do callz <- headStream - let s1 = format_snp_call conf_prior_het callz - S.append s1 <$> format_indel_call conf_prior_indel callz) - =$ collect_lines out') =$ - mapStreamM_ (S.hPut hdl . (flip S.snoc '\n')) - - --- | We do calls of any ploidy, but the FastA output code will fail if --- the ploidy isn't 1 or 2. For indel calls, the FastA output will also --- cheat and pretend it was a haploid call. --- --- XXX For the time being, forward and reverse piles get concatenated. --- For the naive call, this doesn't matter. For the MAQ call, it feels --- more correct to treat them separately and multiply (add?) the results. - -calls :: Maybe Double -> Int -> Pile -> Calls -calls Nothing pl pile = pile { p_snp_pile = s, p_indel_pile = i } - where - !s = simple_snp_call pl $ uncurry (++) $ p_snp_pile pile - !i = force $ simple_indel_call pl $ p_indel_pile pile - -calls (Just theta) pl pile = pile { p_snp_pile = s, p_indel_pile = i } - where - !s = maq_snp_call pl theta $ uncurry (++) $ p_snp_pile pile -- XXX - !i = force $ simple_indel_call pl $ p_indel_pile pile - -instance NFData IndelVariant where - rnf (IndelVariant d (V_Nuc i)) = rnf d `seq` rnf i `seq` () - - --- | Formatting a SNP call. If this was a haplopid call (four GL --- values), we pick the most likely base and pass it on. If it was --- diploid, we pick the most likely dinucleotide and pass it on. - -format_snp_call :: Prob -> Calls -> S.ByteString -format_snp_call p cs - | V.length gl == 4 = S.take 1 $ S.drop (maxQualIndex gl) hapbases - | V.length gl == 10 = S.take 1 $ S.drop (maxQualIndex $ V.zipWith (*) ps gl) dipbases - | otherwise = error "Thou shalt not try to format_snp_call unless thou madeth a haploid or diploid call!" - where - gl = p_snp_pile cs - ps = V.fromListN 10 [p,1,p,1,1,p,1,1,1,p] - dipbases = "NAMCRSGWYKT" - hapbases = "NACGT" - --- | Formatting an Indel call. We pick the most likely variant and --- pass its sequence on. Then we drop incoming calls that should be --- deleted according to the chosen variant. Note that this will blow up --- unless the call was done assuming a haploid genome (which is --- guaranteeed /in this program/)! - -format_indel_call :: Monad m => Prob -> Calls -> Iteratee [Calls] m S.ByteString -format_indel_call p cs - | V.length gl0 == nv = go gl0 - | V.length gl0 == nv * (nv+1) `div` 2 = go homs - | otherwise = error "Thou shalt not try to format_indel_call unless thou madeth a haploid or diploid call!" - where - (gl0,vars) = p_indel_pile cs - !nv = length vars - !homs = V.fromListN nv [ gl0 V.! (i*(i+1) `div` 2 -1) | i <- [1..nv] ] - - go gl = I.dropWhile skip >> return (S.pack $ show $ V.toList ins) - where - eff_gl = V.fromList $ zipWith adjust (V.toList gl) vars - adjust q (IndelVariant ds (V_Nuc is)) = if ds == 0 && V.null is then q else p * q - - IndelVariant del (V_Nuc ins) = ( IndelVariant 0 (V_Nuc V.empty) : vars ) !! maxQualIndex eff_gl - skip ocs = p_refseq ocs == p_refseq cs && p_pos ocs < p_pos cs + del - -maxQualIndex :: V.Vector Prob -> Int -maxQualIndex vec = case V.ifoldl' step (0, 0, 0) vec of - (!i, !m, !m2) -> if m / m2 > 2 then i else 0 - where - step (!i,!m,!m2) j v = if v >= m then (j+1,v,m) else (i,m,m2) - -collect_lines :: Monad m => Enumeratee S.ByteString [S.ByteString] m r -collect_lines = eneeCheckIfDone (liftI . go S.empty) - where - go acc k (EOF mx) = idone (k $ Chunk [acc]) $ EOF mx - go acc k (Chunk s) = case S.splitAt 60 (acc `S.append` s) of - (left, right) | S.null right -> liftI $ go left k - | otherwise -> eneeCheckIfDone (liftI . go right) . k $ Chunk [left] - -by_groups :: ( Monad m, ListLike s a, Nullable s ) - => (a -> a -> Bool) -> (a -> Enumeratee s b m r) -> Enumeratee s b m r -by_groups pr k out = do - mhd <- peekStream - case mhd of - Nothing -> return out - Just hd -> takeWhileE (pr hd) =$ k hd out >>= by_groups pr k - - -output_avro :: FilePath -> (OIter -> IO r) -> IO r -output_avro fn k = if fn == "-" then k (av_out stdout) - else withFile fn WriteMode $ k . av_out - where - av_out :: Handle -> Conf -> Refs -> Iteratee [Calls] IO () - av_out hdl _cfg refs = compileBlocks refs =$ - writeAvroContainer ContainerOpts{..} =$ - mapChunksM_ (S.hPut hdl) - - objects_per_block = 16 - filetype_label = "Genotype Likelihoods V0.1" - - --- Serialize the results from genotype calling in a sensible way. We --- write an Avro file, but we add another blocking layer on top so we --- don't need to endlessly repeat coordinates. - -compileBlocks :: Monad m => Refs -> Enumeratee [Calls] [GenoCallBlock] m a -compileBlocks refs = convStream $ do - c1 <- headStream - tailBlock (p_refseq c1) (p_pos c1) (p_pos c1) (16*1024 :: Int) [pack c1] - where - tailBlock !rs !p0 !po !n acc = do - mc <- peekStream - case mc of - Just c1 | rs == p_refseq c1 && po+1 == p_pos c1 && n > 0 -> do - _ <- headStream - tailBlock rs p0 (po+1) (n-1) $ pack c1 : acc - - _ -> return [ GenoCallBlock - { reference_name = T.decodeLatin1 $ sq_name $ getRef refs rs - , start_position = p0 - , called_sites = reverse acc } ] - - pack c1 = GenoCallSite{..} - where - snp_stats = p_snp_stat c1 - indel_stats = p_indel_stat c1 - snp_likelihoods = compact_likelihoods $ p_snp_pile c1 - indel_likelihoods = compact_likelihoods $ fst $ p_indel_pile c1 - indel_variants = snd $ p_indel_pile c1 - --- | Storing likelihoods: we take the natural logarithm (GL values are --- already in a log scale) and convert to minifloat 0.4.4 --- representation. Range and precision should be plenty. -compact_likelihoods :: V.Vector Prob -> [Int] -- B.ByteString -compact_likelihoods = map fromIntegral {- B.pack -} . V.toList . V.map (float2mini . negate . unPr) - diff --git a/tools/gt-scan.hs b/tools/gt-scan.hs new file mode 100644 --- /dev/null +++ b/tools/gt-scan.hs @@ -0,0 +1,140 @@ +{-# LANGUAGE OverloadedStrings, BangPatterns, RecordWildCards, FlexibleContexts, TypeFamilies #-} +-- Scan file with GT likelihoods, fit something... +-- +-- First iteration: Mitochondrion only. We don't need to fit +-- anything. So far, the likelihoods behave strangely in that smaller +-- \theta is always better, as long as it doesn't become zero. +-- +-- Second iteration: Mitochondrion only, but with a divergence +-- parameter. Needs to be scanned in parallel with a TwoBit file. + +import Bio.Base +import Bio.Bam.Header +import Bio.Genocall.AvroFile +import Bio.Iteratee +import Bio.TwoBit +import Bio.Util.AD +import Bio.Util.Numeric +import Data.Avro +import Data.List ( intercalate ) +import Data.MiniFloat ( mini2float ) +import Data.Strict.Tuple ( Pair((:!:)) ) +import Numeric ( showFFloat ) + +import qualified Data.Vector.Storable as S +import qualified Data.Vector.Unboxed as U +import qualified Data.Vector.Unboxed.Mutable as UM + +main :: IO () +main = do hg19 <- openTwoBit "/mnt/datengrab/hg19.2bit" + mtbl <- UM.replicate (max_lk-min_lk+1) 0 + + let all_lk tbl (p1 :!: p2) ref site = (lk0 p1 site :!:) `fmap` lk1 tbl p2 ref site + + p0 :!: pe <- enumDefaultInputs >=> run $ + joinI $ readAvroContainer $ \meta -> do + foldStreamM (lk_block (getRefseqs meta) (all_lk mtbl) hg19) (1 :!: 1) + + tbl <- U.unsafeFreeze mtbl + + -- optimize llk1 vs. d argument. + let plainfn :: U.Vector Double -> Double + plainfn args = llk1 tbl (unPr pe) $ args U.! 0 + + combofn :: U.Vector Double -> (Double, U.Vector Double) + combofn args = case llk1 tbl (C (unPr pe)) $ D (args U.! 0) (U.singleton 1) of + (D x dx) -> ( x, dx ) + + params = defaultParameters { printFinal = False, verbose = {- Verbose -} Quiet, maxItersFac = 20 } + + (x,q,s) <- optimize params 0.0001 (U.singleton 0.01) + (VFunction plainfn) + (VGradient $ snd . combofn) + (Just $ VCombined combofn) + + -- print $ llk1 tbl (C (unPr pe)) (D 0.001 (U.singleton 1)) + putStrLn $ intercalate "\t" + [ showNum (round $ unPr p0 :: Int), showFFloat (Just 5) (sigmoid2 $ x S.! 0) [] + , show q, showNum (round $ finalValue s :: Int), show s ] + -- print (map sigmoid2 $ S.toList x, q, s) + + +-- | Scans block together with reference sequence. Folds a monadic +-- action over the called sites. +lk_block :: Monad m => Refs -> (b -> Nucleotide -> GenoCallSite -> m b) -> TwoBitFile -> b -> GenoCallBlock -> m b +lk_block refs f tbf b GenoCallBlock{..} = foldM3f b start_position refseq called_sites + where + refseq = getLazySubseq tbf $ Pos (sq_name $ getRef refs reference_name) start_position + + foldM2 acc (x:xs) (y:ys) = do !acc' <- f acc x y ; foldM2 acc' xs ys + foldM2 acc [ ] _ = return acc + foldM2 acc _ [ ] = return acc + + + -- XXX terrible hack to deal with PhiX! Remove this as soon as + -- sensible! + foldM3f acc n (x:xs) (y:ys) + | n `elem` bad = foldM3f acc (succ n) xs ys + | otherwise = do !acc' <- f acc x y ; foldM3f acc' (succ n) xs ys + foldM3f acc _ [ ] _ = return acc + foldM3f acc _ _ [ ] = return acc + + bad = [1400,1643] + -- bad = [586,832,1649,2810,4517] + +{- p_block tbf GenoCallBlock{..} = do + printf "Block %s:%d-%d\n" (show reference_name) start_position + (start_position + length called_sites) + zipWithM_ (curry print) refseq (map snp_likelihoods called_sites) + where + refseq = getLazySubseq tbf (Pos (encodeUtf8 reference_name) start_position) -} + + + +-- | Likelihood with flat prior (no parameters). +lk0 :: Prob -> GenoCallSite -> Prob +lk0 !pp GenoCallSite{..} | U.length snp_likelihoods == 4 = + pp * 0.25 * U.sum (U.map (Pr . negate . mini2float) snp_likelihoods) + | otherwise = pp + +-- | Likelihood precomputation. Total likelihood computes as product +-- over sites @i@ with reference alleles @X_i@: +-- @ +-- L(d) = \prod_i ( (1-d) * GL(X_i) + 1/3 * d * \sum_{Y/=X_i} GL(Y) ) +-- = \prod_i GL(X_i) * \prod_i ( 1 - d + 1/3 * d * \sum_{Y/=X_i} GL(Y)/GL(X) ) +-- @ +-- +-- We compute the first term on the first pass and tabulate a quantized +-- form of the second term: @round (log \sum_{Y/=X_i} GL(Y)/GL(X))@. +-- (Maybe add a scaling factor, though the plain natural log seems +-- pretty good.) + +type LkTableM = UM.IOVector Int +type LkTable = U.Vector Int + +min_lk, max_lk :: Int +min_lk = -256 +max_lk = 255 + +-- | Likelihood with one parameter, the divergence. Computes one +-- part directly, bins the variable part into a mutable table. +lk1 :: LkTableM -> Prob -> Nucleotide -> GenoCallSite -> IO Prob +lk1 tbl !pp ref GenoCallSite{..} | U.length snp_likelihoods == 4 = do + let lx = Pr . negate . mini2float $ snp_likelihoods U.! fromEnum ref + odds = U.ifoldl' (\a i v -> if i == fromEnum ref then a else a + Pr (- mini2float v)) 0 snp_likelihoods / lx + qq = round (unPr odds) `min` max_lk `max` min_lk - min_lk + UM.write tbl qq . succ =<< UM.read tbl qq + return $! pp * lx + | otherwise = return pp + +-- | Actual negative log-likelihood. Gets a table and a divergence +-- value. Returns likelihoods and first two derivatives with respect to +-- the divergence value. +llk1 :: (Ord a, Floating a) => LkTable -> a -> a -> a +llk1 tbl p d = U.ifoldl' step (-p) tbl + where + !d1 = log1p (- sigmoid2 d) + !d3 = log (sigmoid2 d) - log 3 + + step acc qq num = acc - fromIntegral num * (d1 <#> d3 + fromIntegral (min_lk + qq)) + diff --git a/tools/jivebunny.hs b/tools/jivebunny.hs --- a/tools/jivebunny.hs +++ b/tools/jivebunny.hs @@ -23,7 +23,7 @@ -- assignment rates (Done.) import Bio.Bam -import Bio.Util ( showNum ) +import Bio.Util.Numeric ( showNum ) import Control.Applicative import Control.Arrow ( (&&&) ) import Control.Monad ( when, unless, forM_, foldM ) @@ -308,7 +308,8 @@ cf_loudness :: Loudness, cf_single :: Bool, cf_samplesize :: Int, - cf_readgroups :: [FilePath] } + cf_readgroups :: [FilePath], + cf_implied :: [T.Text] } defaultConf :: IO Conf defaultConf = do ixdb <- getDataFileName "index_db.json" @@ -321,7 +322,8 @@ cf_loudness = Normal, cf_single = False, cf_samplesize = 50000, - cf_readgroups = [] } + cf_readgroups = [], + cf_implied = [default_rgs] } options :: [OptDescr (Conf -> IO Conf)] options = [ @@ -332,6 +334,7 @@ Option [ ] ["threshold"] (ReqArg set_thresh "FRAC") "Iterate till uncertainty is below FRAC", Option [ ] ["sample"] (ReqArg set_sample "NUM") "Sample NUM reads for mixture estimation", Option [ ] ["components"] (ReqArg set_compo "NUM") "Print NUM components of the mixture", + Option [ ] ["nocontrol"] (NoArg set_no_control) "Suppress implied read groups for controls", Option "v" ["verbose"] (NoArg set_loud) "Print more diagnostic messages", Option "q" ["quiet"] (NoArg set_quiet) "Print fewer diagnostic messages", Option "h?" ["help", "usage"] (NoArg (const usage)) "Print this message and exit", @@ -344,6 +347,7 @@ set_loud c = return $ c { cf_loudness = Loud } set_quiet c = return $ c { cf_loudness = Quiet } set_single c = return $ c { cf_single = True } + set_no_control c = return $ c { cf_implied = [] } set_thresh a c = readIO a >>= \x -> return $ c { cf_threshold = x } set_sample a c = readIO a >>= \x -> return $ c { cf_samplesize = x } set_compo a c = readIO a >>= \x -> return $ c { cf_num_stats = const x } @@ -381,7 +385,7 @@ Just x | cf_single -> return $ x { p5is = single_placeholder } | otherwise -> return x - rgdefs <- concatMap (readRGdefns (alias_names p7is) (alias_names p5is)) . (:) default_rgs <$> mapM T.readFile cf_readgroups + rgdefs <- concatMap (readRGdefns (alias_names p7is) (alias_names p5is)) . (++) cf_implied <$> mapM T.readFile cf_readgroups notice $ "Got " ++ showNum (U.length (unique_indices p7is)) ++ " unique P7 indices and " ++ showNum (U.length (unique_indices p5is)) ++ " unique P5 indices.\n" notice $ "Declared " ++ showNum (length rgdefs) ++ " read groups.\n" @@ -506,8 +510,13 @@ (take num $ U.foldr fmt_one [] v') inspect' :: HM.HashMap (Int,Int) (B.ByteString, t) -> V.Vector T.Text -> V.Vector T.Text -> Handle -> Int -> Mix -> IO () -inspect' rgs n7 n5 hdl num mix = do - v <- U.unsafeThaw $ U.fromListN (VS.length mix) $ zip [0..] $ VS.toList mix +inspect' rgs n7 n5 hdl num_ mix = do + -- Due to padding, we get invalid indices here. Better filter them + -- out, because we sure can't print them later. + let num = min num_ $ V.length n5 * V.length n7 + v <- U.unsafeThaw $ U.fromListN (V.length n5 * V.length n7) $ + filter (\(i,_) -> i `rem` stride' (V.length n5) < V.length n5) $ + zip [0..] $ VS.toList mix V.partialSortBy (\(_,a) (_,b) -> compare b a) v num -- meh. v' <- U.unsafeFreeze v diff --git a/tools/mt-ccheck.hs b/tools/mt-ccheck.hs --- a/tools/mt-ccheck.hs +++ b/tools/mt-ccheck.hs @@ -50,7 +50,6 @@ import Control.Applicative import Control.Monad import Data.Bits -import Data.Monoid import Data.List import Numeric import System.Console.GetOpt diff --git a/tools/redeye-dar.hs b/tools/redeye-dar.hs new file mode 100644 --- /dev/null +++ b/tools/redeye-dar.hs @@ -0,0 +1,453 @@ +{-# LANGUAGE RecordWildCards, NamedFieldPuns, BangPatterns, TypeFamilies #-} +-- Estimates aDNA damage. Crude first version. +-- +-- - Read or subsample a BAM file, make compact representation of the reads. +-- - Compute likelihood of each read under simple model of +-- damage, error/divergence, contamination. +-- +-- For the fitting, we simplify radically: ignore sequencing error, +-- assume damage and simple, symmetric substitutions which subsume error +-- and divergence. +-- +-- Trying to compute symbolically is too much, the high power terms get +-- out of hand quickly, and we get mixed powers of \lambda and \kappa. +-- The fastest version so far uses the cheap implementation of automatic +-- differentiation in AD.hs together with the Hager-Zhang method from +-- package nonlinear-optimization. BFGS from hmatrix-gsl takes longer +-- to converge. Didn't try an actual Newton iteration (yet?), AD from +-- package ad appears slower. +-- +-- If I include parameters, whose true value is zero, the transformation +-- to the log-odds-ratio doesn't work, because then the maximum doesn't +-- exist anymore. For many parameters, zero makes sense, but one +-- doesn't. A different transformation ('sigmoid2'/'isigmoid2' +-- below) allows for an actual zero (but not an actual one), while +-- avoiding ugly boundary conditions. That appears to work well. +-- +-- The current hack assumes all molecules have an overhang at both ends, +-- then each base gets deaminated with a position dependent probability +-- following a geometric distribution. If we try to model a fraction of +-- undeaminated molecules (a contaminant) in addition, this fails. To +-- rescue the idea, I guess we must really decide if the molecule has an +-- overhang at all (probability 1/2) at each end, then deaminate it. +-- +-- TODO +-- - needs better output +-- - needs support for multiple input files +-- - needs to deal with long (unmerged) reads (by ignoring them?) + +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 +import Bio.Util.AD2 +import Bio.Util.Numeric +import Control.Applicative +import Control.Concurrent.Async +import Control.Monad ( unless ) +import Data.Bits +import Data.Foldable +import Data.Ix +import Data.Maybe +import Data.String ( fromString ) +import Data.Text ( unpack ) +import System.Console.GetOpt +import System.Environment +import System.Exit +import System.FilePath +import System.IO ( hPutStrLn ) + +import qualified Data.HashMap.Strict as M +import qualified Data.Vector as V +import qualified Data.Vector.Generic as G +import qualified Data.Vector.Unboxed as U + +import Prelude hiding ( sequence_, mapM, mapM_, concatMap, sum, minimum, foldr1, foldl ) + +-- | Roughly @Maybe (Nucleotide, Nucleotide)@, encoded compactly +newtype NP = NP { unNP :: Word8 } deriving (Eq, Ord, Ix) +data Seq = Merged { unSeq :: U.Vector Word8 } + | Mate1st { unSeq :: U.Vector Word8 } + | Mate2nd { unSeq :: U.Vector Word8 } + +instance Show NP where + show (NP w) + | w == 16 = "NN" + | w > 16 = "XX" + | otherwise = [ "ACGT" !! fromIntegral (w `shiftR` 2) + , "ACGT" !! fromIntegral (w .&. 3) ] + + +{-# INLINE lk_fun1 #-} +lk_fun1 :: (Num a, Show a, Fractional a, Floating a, Memorable a) + => Int -> Int -> [a] -> V.Vector Seq -> a +lk_fun1 lmin lmax parms = case length parms of + 1 -> V.foldl' (\a b -> a - log (lk tab00 tab00 tab00 b)) 0 . guardV -- undamaged case + where + !tab00 = fromListN (rangeSize my_bounds) [ l_epq p_subst 0 0 x + | (_,_,x) <- range my_bounds ] + + 4 -> V.foldl' (\a b -> a - log (lk tabDS tabDS1 tabDS1 b)) 0 . guardV -- double strand case + where + !tabDS = fromListN (rangeSize my_bounds) [ l_epq p_subst p_d p_e x + | (l,i,x) <- range my_bounds + , let p_d = mu $ lambda ^^ (1+i) + , let p_e = mu $ lambda ^^ (l-i) ] + + !tabDS1 = fromListN (rangeSize my_bounds) [ l_epq p_subst p_d 0 x + | (_,i,x) <- range my_bounds + , let p_d = mu $ lambda ^^ (1+i) ] + + 5 -> V.foldl' (\a b -> a - log (lk tabSS tabSS1 tabSS2 b)) 0 . guardV -- single strand case + where + !tabSS = fromListN (rangeSize my_bounds) [ l_epq p_subst p_d 0 x + | (l,i,x) <- range my_bounds + , let lam5 = lambda ^^ (1+i) ; lam3 = kappa ^^ (l-i) + , let p_d = mu $ lam3 + lam5 - lam3 * lam5 ] + + !tabSS1 = fromListN (rangeSize my_bounds) [ l_epq p_subst p_d 0 x + | (_,i,x) <- range my_bounds + , let p_d = mu $ lambda ^^ (1+i) ] + + !tabSS2 = fromListN (rangeSize my_bounds) [ l_epq p_subst 0 p_d x + | (_,i,x) <- range my_bounds + , let p_d = mu $ lambda ^^ (1+i) ] + + _ -> error "Not supposed to happen: unexpected number of model parameters." + where + ~(l_subst : ~(l_sigma : ~(l_delta : ~(l_lam : ~(l_kap : _))))) = parms + + p_subst = 0.33333 * sigmoid2 l_subst + sigma = sigmoid2 l_sigma + delta = sigmoid2 l_delta + lambda = sigmoid2 l_lam + kappa = sigmoid2 l_kap + + guardV = V.filter (\u -> U.length (unSeq u) >= lmin && U.length (unSeq u) <= lmax) + + -- 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_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 + | otherwise = error $ "Huh? " ++ show x ++ " \\nin " ++ show bnds + + my_bounds = ((lmin,0,NP 0),(lmax,lmax,NP 16)) + mu p = sigma * p + delta * (1-p) + + +-- Likelihood for a certain pair of bases given error rate, C-T-rate +-- and G-A rate. +l_epq :: (Num a, Fractional a, Floating a) => a -> a -> a -> NP -> a +l_epq e p q (NP x) = case x of { + 0 -> s ; 1 -> e ; 2 -> e ; 3 -> e ; + 4 -> e ; 5 -> s-p+4*e*p ; 6 -> e ; 7 -> e+p-4*e*p ; + 8 -> e+q-4*e*q ; 9 -> e ; 10 -> s-q+4*e*q ; 11 -> e ; + 12 -> e ; 13 -> e ; 14 -> e ; 15 -> s ; + _ -> 1 } where s = 1 - 3 * e + + +lkfun :: Int -> Int -> V.Vector Seq -> U.Vector Double -> Double +lkfun lmin lmax brs parms = lk_fun1 lmin lmax (U.toList parms) brs + +lkfun' :: Int -> Int -> V.Vector Seq -> [Double] -> AD +lkfun' lmin lmax brs parms = lk_fun1 lmin lmax (paramVector parms) brs + +lkfun'' :: Int -> Int -> V.Vector Seq -> [Double] -> AD2 +lkfun'' lmin lmax brs parms = lk_fun1 lmin lmax (paramVector2 parms) brs + +combofn :: Int -> Int -> V.Vector Seq -> U.Vector Double -> (Double, U.Vector Double) +combofn lmin lmax brs parms = (x,g) + where D x g = lk_fun1 lmin lmax (paramVector $ U.toList parms) brs + + +data Conf = Conf { + conf_lmin :: Int, + conf_metadata :: FilePath, + conf_report :: String -> IO (), + conf_params :: Parameters } + +defaultConf :: Conf +defaultConf = Conf 25 (error "no config file specified") (\_ -> return ()) quietParameters + +options :: [OptDescr (Conf -> IO Conf)] +options = [ + Option "m" ["min-length"] (ReqArg set_lmin "LEN") "Set minimum length to LEN (25)", + Option "c" ["config"] (ReqArg set_conf "FILE") "Configuiration is stored in FILE", + Option "v" ["verbose"] (NoArg set_verbose) "Print progress reports", + Option "h?" ["help","usage"] (NoArg disp_usage) "Print this message and exit" ] + where + set_lmin a c = readIO a >>= \l -> return $ c { conf_lmin = l } + set_conf f c = return $ c { conf_metadata = f } + set_verbose c = return $ c { conf_report = hPutStrLn stderr, conf_params = debugParameters } + + disp_usage _ = do pn <- getProgName + let blah = "Usage: " ++ pn ++ " [OPTION...] [LIBRARY-NAME...]" + putStrLn $ usageInfo blah options + exitSuccess + +main :: IO () +main = do + (opts, lnames, errors) <- getOpt Permute options <$> getArgs + unless (null errors) $ mapM_ (hPutStrLn stderr) errors >> exitFailure + conf <- foldl (>>=) (return defaultConf) opts + mapM_ (main' conf) lnames + +main' :: Conf -> String -> IO () +main' Conf{..} lname = do + [Library _ fs _] <- return . filter ((fromString lname ==) . library_name) . concatMap sample_libraries . M.elems + =<< readMetadata conf_metadata + + -- XXX meh. subsampling from multiple files is not yet supported :( + brs <- subsampleBam (takeDirectory conf_metadata unpack (head fs)) >=> run $ \_ -> + joinI $ filterStream (\b -> not (isUnmapped (unpackBam b)) && G.length (b_seq (unpackBam b)) >= conf_lmin) $ + joinI $ takeStream 100000 $ + joinI $ mapStream pack_record $ + joinI $ filterStream (\u -> U.length (U.filter (<16) (unSeq u)) * 10 >= 9 * U.length (unSeq u)) $ + stream2vectorN 30000 + + let lmax = V.maximum $ V.map (U.length . unSeq) brs + v0 = crude_estimate brs + opt v = optimize conf_params 0.0001 v + (VFunction $ lkfun conf_lmin lmax brs) + (VGradient $ snd . combofn conf_lmin lmax brs) + (Just . VCombined $ combofn conf_lmin lmax brs) + + results <- mapConcurrently opt [ v0, U.take 4 v0, U.take 1 v0 ] + + let mlk = minimum [ finalValue st | (_,_,st) <- results ] + tot = sum [ exp $ mlk - finalValue st | (_,_,st) <- results ] + p l = exp (mlk - l) / tot + + [ (p_ss, [ _, ssd_sigma_, ssd_delta_, ssd_lambda, ssd_kappa ]), + (p_ds, [ _, dsd_sigma_, dsd_delta_, dsd_lambda ]), + (_ , [ _ ]) ] = [ (p (finalValue st), map sigmoid2 $ G.toList xs) | (xs,_,st) <- results ] + + ssd_sigma = p_ss * ssd_sigma_ + ssd_delta = p_ss * ssd_delta_ + dsd_sigma = p_ds * dsd_sigma_ + dsd_delta = p_ds * dsd_delta_ + + putStrLn $ "p_{ss} = " ++ show p_ss ++ ", p_{ds} = " ++ show p_ds + putStrLn $ show DP{..} + updateMetadata (store_dp lname DP{..}) conf_metadata + + -- Trying to get confidence intervals. Right now, just get the + -- gradient and Hessian at the ML point. Gradient should be nearly + -- zero, Hessian should be symmetric and positive definite. + -- (Remember, we minimized.) + mapM_ print [ (r,s) | (_,r,s) <- results ] + putStrLn "" + mapM_ print [ lkfun' conf_lmin lmax brs (G.toList xs) | (xs,_,_) <- results ] + putStrLn "" + mapM_ print [ lkfun'' conf_lmin lmax brs (G.toList xs) | (xs,_,_) <- results ] + +-- We'll require the MD field to be present. Then we cook each read +-- into a list of paired bases. Deleted bases are dropped, inserted +-- bases replaced with an escape code. +-- +-- XXX This is annoying... almost, but not quite the same as the code +-- in the "Pileup" module. This also relies on MD and doesn't offer the +-- alternative of accessing a reference genome. (The latter may not be +-- worth the trouble.) It also resembles the 'ECig' logic from +-- "Bio.Bam.Rmdup". + +pack_record :: BamRaw -> Seq +pack_record br = if isReversed b then k (revcom u1) else k u1 + where + b@BamRec{..} = unpackBam br + + k | isMerged b = Merged + | isTrimmed b = Merged + | isSecondMate b = Mate2nd + | otherwise = Mate1st + + revcom = U.reverse . U.map (\x -> if x > 15 then x else xor x 15) + u1 = U.fromList . map unNP $ go (G.toList b_cigar) (G.toList b_seq) (fromMaybe [] $ getMd b) + + go :: [Cigar] -> [Nucleotides] -> [MdOp] -> [NP] + + go (_:*0 :cs) ns mds = go cs ns mds + go cs ns (MdNum 0:mds) = go cs ns mds + go cs ns (MdDel []:mds) = go cs ns mds + go _ [] _ = [] + go [] _ _ = [] + + go (Mat:*nm :cs) (n:ns) (MdNum mm:mds) = mk_pair n n : go (Mat:*(nm-1):cs) ns (MdNum (mm-1):mds) + go (Mat:*nm :cs) (n:ns) (MdRep n':mds) = mk_pair n n' : go (Mat:*(nm-1):cs) ns mds + go (Mat:*nm :cs) ns (MdDel _ :mds) = go (Mat:* nm :cs) ns mds + + go (Ins:*nm :cs) ns mds = replicate nm esc ++ go cs (drop nm ns) mds + go (SMa:*nm :cs) ns mds = replicate nm esc ++ go cs (drop nm ns) mds + go (Del:*nm :cs) ns (MdDel (_:ds):mds) = go (Del:*(nm-1):cs) ns (MdDel ds:mds) + go (Del:*nm :cs) ns ( _:mds) = go (Del:* nm :cs) ns mds + + go (_:cs) nd mds = go cs nd mds + + +esc :: NP +esc = NP 16 + +mk_pair :: Nucleotides -> Nucleotides -> NP +mk_pair (Ns a) = case a of 1 -> mk_pair' 0 + 2 -> mk_pair' 1 + 4 -> mk_pair' 2 + 8 -> mk_pair' 3 + _ -> const esc + where + mk_pair' u (Ns b) = case b of 1 -> NP $ u .|. 0 + 2 -> NP $ u .|. 4 + 4 -> NP $ u .|. 8 + 8 -> NP $ u .|. 12 + _ -> esc + + +infix 7 /%/ +(/%/) :: Integral a => a -> a -> Double +0 /%/ 0 = 0 +a /%/ b = fromIntegral a / fromIntegral b + +-- Crude estimate. Need two overhang lengths, two deamination rates, +-- undamaged fraction, SS/DS, substitution rate. +-- +-- DS or SS: look whether CT or GA is greater at 3' terminal position √ +-- Left overhang length: ratio of damage at second position to first √ +-- Right overang length: ratio of CT at last to snd-to-last posn √ +-- + ratio of GA at last to snd-to-last posn √ +-- SS rate: condition on damage on one end, compute rate at other √ +-- DS rate: condition on damage, compute rate in interior √ +-- substitution rate: count all substitutions not due to damage √ +-- undamaged fraction: see below √ +-- +-- Contaminant fraction: let f5 (f3, f1) be the fraction of reads +-- showing damage at the 5' end (3' end, both ends). Let a (b) be +-- the probability of an endogenous reads to show damage at the 5' +-- end (3' end). Let e be the fraction of endogenous reads. Then +-- we have: +-- +-- f5 = e * a +-- f3 = e * b +-- f1 = e * a * b +-- +-- f5 * f3 / f1 = e +-- +-- Straight forward and easy to understand, but in practice, this method +-- produces ridiculous overestimates, ridiculous underestimates, +-- negative contamination rates, and general grief. It's actually +-- better to start from a constant number. + + +crude_estimate :: V.Vector Seq -> U.Vector Double +crude_estimate seqs0 = U.fromList [ l_subst, l_sigma, l_delta, l_lam, l_kap ] + where + seqs = V.filter ((>= 10) . U.length) $ V.map unSeq seqs0 + + total_equals = V.sum (V.map (U.length . U.filter isNotSubst) seqs) + total_substs = V.sum (V.map (U.length . U.filter isOrdinarySubst) seqs) * 6 `div` 5 + l_subst = isigmoid2 $ max 0.001 $ total_substs /%/ (total_equals + total_substs) + + c_to_t, g_to_a, c_to_c :: Word8 + c_to_t = 7 + g_to_a = 8 + c_to_c = 5 + + isNotSubst x = x < 16 && x `shiftR` 2 == x .&. 3 + isOrdinarySubst x = x < 16 && x `shiftR` 2 /= x .&. 3 && + x /= c_to_t && x /= g_to_a + + ct_at_alpha = V.length $ V.filter (\v -> v U.! 0 == c_to_t && dmg_omega v) seqs + cc_at_alpha = V.length $ V.filter (\v -> v U.! 0 == c_to_c && dmg_omega v) seqs + ct_at_beta = V.length $ V.filter (\v -> v U.! 1 == c_to_t && dmg_omega v) seqs + cc_at_beta = V.length $ V.filter (\v -> v U.! 1 == c_to_c && dmg_omega v) seqs + + dmg_omega v = v U.! (l-1) == c_to_t || v U.! (l-1) == g_to_a + || v U.! (l-2) == c_to_t || v U.! (l-2) == g_to_a + || v U.! (l-3) == c_to_t || v U.! (l-3) == g_to_a + where l = U.length v + + l_lam = isigmoid2 lambda + lambda = min 0.9 $ max 0.1 $ + (ct_at_beta * (cc_at_alpha + ct_at_alpha)) /%/ + ((cc_at_beta + ct_at_beta) * ct_at_alpha) + + ct_at_omega = V.length $ V.filter (\v -> v U.! (U.length v -1) == c_to_t && dmg_alpha v) seqs + cc_at_omega = V.length $ V.filter (\v -> v U.! (U.length v -1) == c_to_c && dmg_alpha v) seqs + ct_at_psi = V.length $ V.filter (\v -> v U.! (U.length v -2) == c_to_t && dmg_alpha v) seqs + cc_at_psi = V.length $ V.filter (\v -> v U.! (U.length v -2) == c_to_c && dmg_alpha v) seqs + + dmg_alpha v = v U.! 0 == c_to_t || v U.! 1 == c_to_t || v U.! 2 == c_to_t + + l_kap = isigmoid2 $ min 0.9 $ max 0.1 $ + (ct_at_psi * (cc_at_omega+ct_at_omega)) /%/ + ((cc_at_psi+ct_at_psi) * ct_at_omega) + + total_inner_CCs = V.sum $ V.map (U.length . U.filter (== c_to_c) . takeInner) seqs + total_inner_CTs = V.sum $ V.map (U.length . U.filter (== c_to_t) . takeInner) seqs + takeInner v = U.slice 5 (U.length v - 10) v + + delta = (total_inner_CTs /%/ (total_inner_CTs+total_inner_CCs)) + raw_rate = ct_at_alpha /%/ (ct_at_alpha + cc_at_alpha) + + -- clamping is necessary if f_endo ends up wrong + l_delta = isigmoid2 $ min 0.99 delta + l_sigma = isigmoid2 . min 0.99 $ raw_rate / lambda + + +class Memorable a where + type Memo a :: * + + fromListN :: Int -> [a] -> Memo a + bang :: Memo a -> Int -> a + +instance Memorable Double where + type Memo Double = U.Vector Double + + fromListN = U.fromListN + bang = (U.!) + +instance Memorable AD where + type Memo AD = (Int, U.Vector Double) + + fromListN _ [ ] = error "unexpected: tried to memorize an empty list" + fromListN _ (C _ :_) = error "unexpected: tried to memorize a value without derivatives" + fromListN n xs@(D _ v:_) = (1+d, U.fromListN (n * (1+d)) $ concatMap unp xs) + where + !d = U.length v + unp (C a) = a : replicate d 0 + unp (D a da) = a : U.toList da + + bang (d, v) i = D (v U.! (d*i+0)) (U.slice (d*i+1) (d-1) v) + +instance Memorable AD2 where + type Memo AD2 = (Int, U.Vector Double) + + fromListN _ [ ] = error "unexpected: tried to memorize an empty list" + fromListN _ (C2 _ : _) = error "unexpected: tried to memorize a value without derivatives" + fromListN n xs@(D2 _ v _ : _) = (d, U.fromListN (n * (1+d+d*d)) $ concatMap unp xs) + where + !d = U.length v + unp (C2 a) = a : replicate (d+d*d) 0 + unp (D2 a da dda) = a : U.toList da ++ U.toList dda + + bang (d, v) i = D2 (v U.! (stride*i)) + (U.slice (stride*i+1) d v) + (U.slice (stride*i+1+d) (d*d) v) + where + stride = 1 + d + d*d + + +store_dp :: String -> DamageParameters Double -> Metadata -> Metadata +store_dp lname dp = M.map go1 + where + go1 (Sample ls af bf ts dv) = Sample (map go2 ls) af bf ts dv + go2 (Library nm fs dmg) + | nm == fromString lname = Library nm fs (Just dp) + | otherwise = Library nm fs dmg + diff --git a/tools/redeye-div.hs b/tools/redeye-div.hs new file mode 100644 --- /dev/null +++ b/tools/redeye-div.hs @@ -0,0 +1,162 @@ +{-# LANGUAGE BangPatterns, RecordWildCards, OverloadedStrings #-} +-- Here we read the tables from sample_div_tables, add them up as +-- necessary, estimate divergence and heterozygosity from them, and +-- store the result back. The estimate can be done for regions, which +-- are defined by regular expressions. + +import Bio.Base +import Bio.Genocall.Metadata +import Bio.Util.AD +import Bio.Util.AD2 +import Bio.Util.Numeric ( log1p ) +import Bio.Util.Regex ( regComp, regMatch ) +import Control.Concurrent.Async ( async, wait ) +import Control.Monad ( when, unless, forM, (>=>) ) +import Data.Foldable ( foldMap ) +import Data.List ( foldl1' ) +import Data.String ( fromString ) +import Data.Text ( Text, unpack ) +import Numeric ( showFFloat ) +import Numeric.LinearAlgebra.HMatrix ( eigSH', (><), toRows, scale ) +import System.Console.GetOpt +import System.Environment ( getArgs, getProgName ) +import System.Exit ( exitSuccess, exitFailure ) +import System.IO ( hPutStrLn, stderr ) + +import qualified Data.HashMap.Strict as H +import qualified Data.Vector.Storable as VS +import qualified Data.Vector.Unboxed as U + +data Conf = Conf { conf_metadata :: FilePath + , conf_regions :: [Text] + , conf_purge :: Bool + , conf_verbose :: Bool } + +defaultConf :: Conf +defaultConf = Conf (error "no metadata file specified") [] False False + +options :: [OptDescr (Conf -> IO Conf)] +options = [ + Option "c" ["config"] (ReqArg set_conf "FILE") "Set name of json config file to FILE", + Option "r" ["region"] (ReqArg add_region "REGEX") "What matches REGEX becomes a region", + Option "p" ["purge"] (NoArg do_purge) "Purge tables after use", + Option "H" ["human"] (NoArg set_human) "Use regions for a human genome", + Option "v" ["verbose"] (NoArg be_verbose) "Print more diagnostics", + Option "h?" ["help","usage"] (NoArg disp_usage) "Display this message" ] + where + be_verbose c = return $ c { conf_verbose = True } + do_purge c = return $ c { conf_purge = True } + set_conf fn c = return $ c { conf_metadata = fn } + add_region re c = return $ c { conf_regions = fromString re : conf_regions c } + set_human c = return $ c { conf_regions = [ "^(chr)?[0-9]+$", "^(chr)?X$", "^(chr)?Y$" ] } + + disp_usage _ = do pn <- getProgName + let blah = "Usage: " ++ pn ++ " [OPTION...] [SAMPLE...]" + putStrLn $ usageInfo blah options + exitSuccess + + +main :: IO () +main = do + (opts, samples, errs) <- getOpt Permute options `fmap` getArgs + Conf{..} <- foldl (>>=) (return defaultConf) opts + unless (null errs) $ mapM_ (hPutStrLn stderr) errs >> exitFailure + + meta0 <- readMetadata conf_metadata + + let eff_samples = if null samples then H.keys meta0 else map fromString samples + eff_regions = if null conf_regions then [""] else conf_regions + + updates <- forM eff_samples >=> mapM wait $ \sample -> case H.lookup sample meta0 of + + Nothing -> do hPutStrLn stderr $ "unknown sample " ++ show sample + async $ return id + + 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)) + $ H.elems + $ H.filterWithKey (match rgn) + $ sample_div_tables smp + + let app_purge = if conf_purge then appEndo (foldMap purge eff_regions) else id + upd_smp smp' = smp' { sample_divergences = ins_many ests $ sample_divergences smp' + , sample_div_tables = app_purge $ sample_div_tables smp' } + + when conf_verbose $ putStrLn $ "Estimate done for " ++ show sample ++ "." + return $ H.adjust upd_smp sample + updateMetadata (foldr (.) id updates) conf_metadata + + where + match :: Text -> Text -> a -> Bool + match rgn = const . regMatch (regComp $ unpack rgn) . unpack + + purge :: Text -> Endo (H.HashMap Text a) + purge rgn = Endo $ H.filterWithKey ((.) not . match rgn) + + ins_many :: [(Text,v)] -> H.HashMap Text v -> H.HashMap Text v + ins_many = flip $ foldr (uncurry H.insert) + + +-- 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 + (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 + + let showRes [dv,h1,h2] = + "D = " ++ showFFloat (Just 3) dv ", " ++ + "H1 = " ++ showFFloat (Just 3) h1 ", " ++ + "H2 = " ++ showFFloat (Just 3) h2 "" + showRes _ = error "Wtf? (1)" + + -- Confidence interval: PCA on Hessian matrix, then for each + -- eigenvalue λ add/subtract 1.96 / sqrt λ times the corresponding + -- eigenvalue to the estimate. Should describe a nice spheroid. + let D2 _val grd hss = llk2 tab (paramVector2 $ VS.toList fit) + d = U.length grd + (evals, evecs) = eigSH' $ (d >< d) (U.toList hss) + intervs = [ (xform (fit + scale lam evec), xform (fit + scale (-lam) evec)) + | (eval, evec) <- zip (VS.toList evals) (toRows evecs), let lam = 1.96 / sqrt eval ] + + when verbose $ putStrLn $ unlines $ + (:) (show res ++ ", " ++ show stats { finalValue = finalValue stats - llk_rr }) $ + (:) (showRes $ xform fit) $ + map (\(u,v) -> "[ " ++ showRes u ++ " .. " ++ showRes v ++ " ]") intervs + + return $! DivEst (xform fit) intervs + +llk :: U.Vector Int -> [AD] -> AD +llk tab [delta,eta,eta2] = llk' tab 0 delta eta + llk' tab 6 delta eta2 +llk _ _ = error "Wtf? (3)" + +llk2 :: U.Vector Int -> [AD2] -> AD2 +llk2 tab [delta,eta,eta2] = llk' tab 0 delta eta + llk' tab 6 delta eta2 +llk2 _ _ = error "Wtf? (4)" + +{-# INLINE llk' #-} +llk' :: (Ord a, Floating a) => U.Vector Int -> Int -> a -> a -> a +llk' tab base delta eta = block (base+0) g_RR g_RA g_AA + + block (base+1) g_RR g_AA g_RA + + block (base+2) g_RA g_RR g_AA + + block (base+3) g_RA g_AA g_RR + + block (base+4) g_AA g_RR g_RA + + block (base+5) g_AA g_RA g_RR + where + !maxD2 = U.length tab `div` 12 + !maxD = round (sqrt (fromIntegral maxD2) :: Double) + + !g_RR = 1 / Pr (log1p (exp delta)) + !g_AA = Pr delta / Pr (log1p (exp delta)) * 1 / Pr (log1p (exp eta)) + !g_RA = Pr delta / Pr (log1p (exp delta)) * Pr eta / Pr (log1p (exp eta)) + + block ix g1 g2 g3 = U.ifoldl' step 0 $ U.slice (ix * maxD2) maxD2 tab + where + step !acc !i !num = acc - fromIntegral num * unPr p + where + (!d1,!d2) = i `quotRem` maxD + p = g1 + Pr (- fromIntegral d1) * g2 + Pr (- fromIntegral (d1+d2)) * g3 + diff --git a/tools/redeye-pileup.hs b/tools/redeye-pileup.hs new file mode 100644 --- /dev/null +++ b/tools/redeye-pileup.hs @@ -0,0 +1,325 @@ +{-# LANGUAGE RecordWildCards, BangPatterns, OverloadedStrings, FlexibleContexts #-} +-- Command line driver for simple genotype calling. We have three +-- separate steps: Pileup from a BAM file (or multiple merged files) to +-- produce likelihoods (and some auxillary statistics). These are +-- written into an Avro container. Next we need to estimate parameters, +-- in the simplest case divergence and heterozygosity. We can save some +-- time by fusing this with the first step. The final step is calling +-- bases by scnaning the Avro container and applying some model, and +-- again, in the simplest case that's just divergence and +-- heterozygosity. We keep that separate, because different models will +-- require different programs. So here we produce likelihoods and +-- a simple model fit. + +-- The likelihoods depend on damage parameters and an error model, +-- otherwise they are 'eternal'. (For the time being, it's probably +-- wise to go with the naïve error model.) Technically, they also +-- depend on ploidy, but since only diploid organisms are interesting +-- right now, we fix that to two. We pay some overhead on the sex +-- chromosomes, but the simplification is worth it. + +-- About damage parameters: We effectively have three different models +-- (SS, DS, no damage) and it may not be possible to choose one a +-- priori. To manage this cleanly, we should have one universal model, +-- but the three we have are not generalizations of each other. +-- However, all can be generalized into one model with slightly more +-- parameters. See tools/dmg-est.hs for how we fit the model. + +-- Calling is always diploid, for maximum flexibility. We don't really +-- support higher ploidies, so the worst damage is that we output an +-- overhead of 150% useless likelihood values for the sex chromosomes +-- and maybe estimate heterozygosity where there is none. + +import Bio.Base +import Bio.Bam.Header +import Bio.Bam.Index +import Bio.Bam.Pileup +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 +import Control.Applicative +import Control.Monad +import Data.Aeson +import Data.Avro +import Data.String ( fromString ) +import Data.Vec.Packed ( packMat ) +import System.Console.GetOpt +import System.Directory ( renameFile ) +import System.Environment +import System.Exit +import System.FilePath +import System.IO + +import qualified Data.ByteString.Char8 as S +import qualified Data.ByteString.Lazy as BL +import qualified Data.Foldable as F +import qualified Data.HashMap.Strict as H +import qualified Data.Text as T +import qualified Data.Text.Encoding as T +import qualified Data.Vector.Storable as VS +import qualified Data.Vector as V +import qualified Data.Vector.Unboxed as U +import qualified Data.Vector.Unboxed.Mutable as M +import qualified Data.Sequence as Z + +data Conf = Conf { + -- Generator for output file name. Receives sample name and + -- (optional) region as arguments. + conf_output :: String -> Maybe String -> FilePath, + conf_metadata :: FilePath, + conf_theta :: Maybe Double, + conf_report :: String -> IO () } + +defaultConf :: Conf +defaultConf = Conf default_out (error "no metadata file specified") Nothing (\_ -> return ()) + where + default_out smp Nothing = smp <> ".av" + default_out smp (Just rgn) = smp <> "-" <> rgn <> ".av" + +options :: [OptDescr (Conf -> IO Conf)] +options = [ + Option "c" ["config"] (ReqArg set_conf "FILE") "Set name of json config file to FILE", + Option "o" ["output"] (ReqArg set_output "FILE") "Set out file schema to FILE", + Option "t" dep_param (ReqArg set_theta "FRAC") "Set dependency coefficient to FRAC (\"N\" to turn off)", + Option "v" ["verbose"] (NoArg be_verbose) "Print more diagnostics", + Option "h?" ["help","usage"] (NoArg disp_usage) "Display this message" ] + where + dep_param = ["theta","dependency-coefficient"] + + disp_usage _ = do pn <- getProgName + let blah = "Usage: " ++ pn ++ " [OPTION...] [SAMPLE [REGION...] ...]" + putStrLn $ usageInfo blah options + exitSuccess + + be_verbose c = return $ c { conf_report = hPutStrLn stderr } + set_conf fn c = return $ c { conf_metadata = fn } + + set_theta "N" c = return $ c { conf_theta = Nothing } + set_theta a c = (\t -> c { conf_theta = Just t }) <$> readIO a + + set_output fn c = return $ c { conf_output = mkoutput fn } + +mkoutput :: FilePath -> String -> Maybe String -> FilePath +mkoutput str smp rgn = go str + where + go ('%':'s':s) = smp ++ go s + go ('%':'r':s) = maybe id (++) rgn $ go s + go ('%':'%':s) = '%' : go s + go ('%': c :s) = c : go s + go ( c :s) = c : go s + go [ ] = [ ] + +main :: IO () +main = do + (opts, samprgns, errs) <- getOpt Permute options <$> getArgs + Conf{..} <- foldl (>>=) (return defaultConf) opts + unless (null errs) $ mapM_ (hPutStrLn stderr) errs >> exitFailure + + -- samprgns contains samples and regions. We define anything found in the metadata as sample, + -- anything else as region to apply to the previous sample. Name a sample "chr1" and you get + -- what you deserve. + + samples <- flip split_sam_rgns samprgns <$> readMetadata conf_metadata + when (null samples) $ hPutStrLn stderr "need (at least) one sample name" >> exitFailure + + forM_ samples $ \(sample, rgns) -> do + meta <- readMetadata conf_metadata + + case H.lookup (fromString sample) meta of + Nothing -> hPutStrLn stderr $ "unknown sample " ++ show sample + + Just smp -> forM_ rgns $ \rgn -> do + let outstem = conf_output sample rgn + outfile = takeDirectory conf_metadata outstem + tmpfile = outfile ++ ".#" + (tab,()) <- withFile tmpfile WriteMode $ \ohdl -> + mergeLibraries conf_report conf_metadata (sample_libraries smp) rgn >=> run $ \hdr -> + progressPos (\(rs, p, _) -> (rs, p)) "GT call at " conf_report (meta_refs hdr) =$ + pileup =$ + mapStream (calls conf_theta) =$ + zipStreams tabulateSingle (output_avro ohdl $ meta_refs hdr) + + let upd_sample s = s { sample_div_tables = H.insert (maybe T.empty T.pack rgn) tab (sample_div_tables s) + , sample_avro_files = H.insert (maybe T.empty T.pack rgn) (fromString outstem) (sample_avro_files s) } + + updateMetadata (H.adjust upd_sample (fromString sample)) conf_metadata + renameFile tmpfile outfile + +mergeLibraries :: (MonadIO m, MonadMask m) + => (String -> IO ()) -> FilePath + -> [Library] -> Maybe String -> Enumerator' BamMeta [PosPrimChunks] m b +mergeLibraries report cfg [ l ] mrgn = enumLibrary report cfg l mrgn +mergeLibraries report cfg (l:ls) mrgn = mergeEnums' (mergeLibraries report cfg ls mrgn) (enumLibrary report cfg l mrgn) mm + where + mm _ = mergeSortStreams $ \(rs1, p1, _) (rs2, p2, _) -> if (rs1, p1) < (rs2, p2) then Less else NotLess + +enumLibrary :: (MonadIO m, MonadMask m) + => (String -> IO ()) -> FilePath + -> Library -> Maybe String -> Enumerator' BamMeta [PosPrimChunks] m b +enumLibrary report cfg (Library nm fs mdp) mrgn output = do + let (msg, dmg) = case mdp of Nothing -> ("no damage model", noDamage) + Just dp -> ("universal damage parameters" ++ show dp, univDamage dp) + + liftIO . report $ "using " ++ msg ++ " for " ++ T.unpack nm + mergeInputRgns mrgn combineCoordinates (map (() (takeDirectory cfg) . T.unpack) fs) + $== takeWhileE (isValidRefseq . b_rname . unpackBam) + $== mapMaybeStream (\br -> + let b = unpackBam br + m = dmg (isReversed b) (VS.length (b_qual b)) + in decompose (map packMat $ V.toList m) br) + $ output + +mergeInputRgns :: (MonadIO m, MonadMask m) + => Maybe String + -> (BamMeta -> Enumeratee [BamRaw] [BamRaw] (Iteratee [BamRaw] m) a) + -> [FilePath] -> Enumerator' BamMeta [BamRaw] m a +mergeInputRgns _ _ [ ] = \k -> return (k mempty) +mergeInputRgns Nothing (?) fps = mergeInputs (?) fps +mergeInputRgns (Just rs) (?) (fp0:fps0) = go fp0 fps0 + where + enum1 fp k1 = do idx <- liftIO $ readBamIndex fp + enumFileRandom defaultBufSize fp >=> run >=> run $ + decodeAnyBam $ \hdr -> + let Just ri = Z.findIndexL ((==) rs . unpackSeqid . sq_name) (meta_refs hdr) + in eneeBamRefseq idx (Refseq $ fromIntegral ri) $ k1 hdr + + go fp [ ] = enum1 fp + go fp (fp1:fps) = mergeEnums' (go fp1 fps) (enum1 fp) (?) + + +-- | Ploidy is hardcoded as two here. Can be changed if the need +-- arises. +-- +-- XXX For the time being, forward and reverse piles get concatenated. +-- For the naive call, this doesn't matter. For the MAQ call, it feels +-- more correct to treat them separately and multiply (add?) the results. + +calls :: Maybe Double -> Pile -> Calls +calls Nothing pile = pile { p_snp_pile = s, p_indel_pile = i } + where + !s = simple_snp_call fq 2 $ uncurry (++) $ p_snp_pile pile + !i = simple_indel_call 2 $ p_indel_pile pile + -- XXX this should be a cmdline option + -- fq = min 1 . (*) 1.333 . fromQual + fq = fromQual + +calls (Just theta) pile = pile { p_snp_pile = s, p_indel_pile = i } + where + !i = simple_indel_call 2 $ p_indel_pile pile + + -- This lumps the two strands together + -- !s = maq_snp_call 2 theta $ uncurry (++) $ p_snp_pile pile -- XXX + + -- This treats them separately + !s | r == r' = Snp_GLs (U.zipWith (*) x y) r -- same ref base (normal case): multiply + | r == nucsN = Snp_GLs y r' -- forward ref is N, use backward call + | otherwise = Snp_GLs x r -- else use forward call (even if this is incorrect, + where -- there is nothing else we can do here) + Snp_GLs x r = maq_snp_call 2 theta $ fst $ p_snp_pile pile + Snp_GLs y r' = maq_snp_call 2 theta $ snd $ p_snp_pile pile + + +-- | Serialize the results from genotype calling in a sensible way. We +-- write an Avro file, but we add another blocking layer on top so we +-- don't need to endlessly repeat coordinates. + +compileBlocks :: Monad m => Enumeratee [Calls] [GenoCallBlock] m a +compileBlocks = convStream $ do + c1 <- headStream + tailBlock (p_refseq c1) (p_pos c1) (p_pos c1) . (:[]) $! pack c1 + where + tailBlock !rs !p0 !po acc = do + mc <- peekStream + case mc of + Just c1 | rs == p_refseq c1 && po+1 == p_pos c1 && po - p0 < 65536 -> do + _ <- headStream + tailBlock rs p0 (po+1) . (:acc) $! pack c1 + + _ -> return [ GenoCallBlock + { reference_name = rs + , start_position = p0 + , called_sites = reverse acc } ] + + pack c1 = rlist indel_variants `seq` GenoCallSite{..} + where + Snp_GLs snp_pls !ref_allele = p_snp_pile c1 + + !snp_stats = p_snp_stat c1 + !indel_stats = p_indel_stat c1 + !snp_likelihoods = compact_likelihoods snp_pls + !indel_likelihoods = compact_likelihoods $ fst $ p_indel_pile c1 + !indel_variants = snd $ p_indel_pile c1 + + rlist [] = () + rlist (x:xs) = x `seq` rlist xs + + +output_avro :: Handle -> Refs -> Iteratee [Calls] IO () +output_avro hdl refs = compileBlocks =$ + writeAvroContainer ContainerOpts{..} =$ + mapChunksM_ (S.hPut hdl) + where + objects_per_block = 16 + filetype_label = "Genotype Likelihoods V0.1" + initial_schemas = H.singleton "Refseq" $ + object [ "type" .= String "enum" + , "name" .= String "Refseq" + , "symbols" .= Array + (V.fromList . map (String . T.decodeUtf8 . sq_name) $ F.toList refs) ] + meta_info = H.singleton "biohazard.refseq_length" $ + S.concat $ BL.toChunks $ encode $ Array $ V.fromList + [ Number (fromIntegral (sq_length s)) | s <- F.toList refs ] + + +maxD :: Int +maxD = 64 + +-- | Parameter estimation for a single sample. The parameters are +-- divergence and heterozygosity. We tabulate the data here and do the +-- 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 = do + tab <- liftIO $ M.replicate (12 * maxD * maxD) (0 :: Int) + (,) <$> 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 + -- ordering is like in BCF, with the reference first. + -- Ref ~ A ==> PL ~ AA, AC, CC, AG, CG, GG, AT, CT, GT, TT + {-# INLINE accum #-} + accum !tab !acc (Snp_GLs !gls !ref) + | U.length gls /= 10 = error "Ten GL values expected for SNP!" -- should not happen + | ref `elem` [nucsC,nucsG] = accum' 0 tab acc gls + | ref `elem` [nucsA,nucsT] = accum' 6 tab acc gls + | otherwise = return acc -- unknown reference + + -- The simple 2D table didn't work, it lacked resolution in some + -- cases. We make six separate tables instead so we can store two + -- differences with good resolution in every case. + {-# INLINE accum' #-} + accum' refix !tab !acc !gls + | g_RR >= g_RA && g_RA >= g_AA = store 0 g_RR g_RA g_AA + | g_RR >= g_AA && g_AA >= g_RA = store 1 g_RR g_AA g_RA + | g_RA >= g_RR && g_RR >= g_AA = store 2 g_RA g_RR g_AA + | g_RA >= g_AA && g_AA >= g_RR = store 3 g_RA g_AA g_RR + | g_RR >= g_RA = store 4 g_AA g_RR g_RA + | otherwise = store 5 g_AA g_RA g_RR + + where + g_RR = unPr $ U.unsafeIndex gls 0 + g_RA = unPr $ (U.unsafeIndex gls 1 + U.unsafeIndex gls 3 + U.unsafeIndex gls 6) / 3 + g_AA = unPr $ (U.unsafeIndex gls 2 + U.unsafeIndex gls 5 + U.unsafeIndex gls 9) / 3 + + store t a b c = do let d1 = min (maxD-1) . round $ a - b + d2 = min (maxD-1) . round $ b - c + ix = (t + refix) * maxD * maxD + d1 * maxD + d2 + liftIO $ M.read tab ix >>= M.write tab ix . succ + return $! acc + a + diff --git a/tools/redeye-single.hs b/tools/redeye-single.hs new file mode 100644 --- /dev/null +++ b/tools/redeye-single.hs @@ -0,0 +1,287 @@ +{-# LANGUAGE BangPatterns, RecordWildCards, OverloadedStrings, FlexibleContexts #-} +-- Genotype calling for a single individual. The only parameters needed +-- are the (prior) probabilities for a heterozygous or homozygous variant. +-- +-- (This can be extended easily into a caller for a homogenous +-- population where individuals are assumed to be randomly related (i.e. +-- not family). In this case, the prior is the allele frequency +-- spectrum, the call would be the set(!) of genotypes that has maximum +-- posterior probability. Computation is possible in quadratic time and +-- linear space using a DP scheme; see Heng Li's paper for details.) +-- +-- What's the output format? Fasta or Fastq could be useful in limited +-- circumstances, else BCF (not VCF) would be canonical. Or maybe BCF +-- restricted to variant sites. Or BCF restricted to sites not known to +-- be reference. People will come up with filters for sure... + +import Bio.Base +import Bio.Bam +import Bio.Bam.Pileup +import Bio.Genocall.AvroFile +import Bio.Genocall.Metadata +import Bio.Iteratee.Builder +import Bio.Util.Regex ( Regex, regComp, regMatch ) +import Control.Applicative +import Control.Exception ( bracket ) +import Control.Monad +import Data.Avro +import Data.Bits +import Data.Foldable ( toList, foldMap ) +import Data.MiniFloat +import Data.String +import Data.Text ( Text, unpack ) +import Foreign.Ptr ( castPtr ) +import System.Console.GetOpt +import System.Directory ( renameFile ) +import System.Environment +import System.Exit +import System.FilePath +import System.Posix.IO + +import qualified Data.ByteString.Char8 as S +import qualified Data.ByteString.Unsafe as S +import qualified Data.HashMap.Strict as H +import qualified Data.Vector.Unboxed as U +import qualified System.IO as IO + +data Conf = Conf { + conf_metadata :: FilePath, + -- | Generator for output file name. Receives sample name and + -- (optional) region as arguments. + conf_output :: String -> Text -> FilePath, + conf_regions :: Regex, + conf_ploidy :: String -> Int, + conf_report :: String -> IO () } + +defaultConf :: Conf +defaultConf = Conf (error "no metadata file specified") default_out (regComp "") (const 2) (\_ -> return ()) + where + default_out smp "" = smp <> ".bcf" + default_out smp rgn = smp <> "-" <> unpack rgn <> ".bcf" + +options :: [OptDescr (Conf -> IO Conf)] +options = [ + Option "c" ["config"] (ReqArg set_conf "FILE") "Set name of json config file to FILE", + Option "o" ["output"] (ReqArg set_output "FILE") "Set output file pattern to FILE", + Option "r" ["regions"] (ReqArg set_regions "REGEX") "Process only regions matching REGEX", + Option "1" ["haploid-chromosomes"] (ReqArg set_hap "REGEX") "Targets matching REGEX are haploid", + Option "2" ["diploid-chromosomes"] (ReqArg set_dip "REGEX") "Targets matching REGEX are diploid", + Option "v" ["verbose"] (NoArg be_verbose) "Print more diagnostics", + Option "h?" ["help","usage"] (NoArg disp_usage) "Display this message" ] + where + disp_usage _ = do pn <- getProgName + let blah = "Usage: " ++ pn ++ " [OPTION...] [SAMPLE [REGION...] ...]" + putStrLn $ usageInfo blah options + exitFailure + + be_verbose c = return $ c { conf_report = IO.hPutStrLn stderr } + set_conf fn c = return $ c { conf_metadata = fn } + + set_hap a c = return $ c { conf_ploidy = \chr -> if regMatch (regComp a) chr then 1 else conf_ploidy c chr } + set_dip a c = return $ c { conf_ploidy = \chr -> if regMatch (regComp a) chr then 2 else conf_ploidy c chr } + set_regions a c = return $ c { conf_regions = regComp $ "^" ++ a ++ "$" } + + set_output fn c = return $ c { conf_output = mkoutput fn } + +mkoutput :: FilePath -> String -> Text -> FilePath +mkoutput str smp rgn = go str + where + go ('%':'s':s) = smp ++ go s + go ('%':'r':s) = unpack rgn ++ go s + go ('%':'%':s) = '%' : go s + go ('%': c :s) = c : go s + go ( c :s) = c : go s + go [ ] = [ ] + +main :: IO () +main = do + (opts, samples, errs) <- getOpt Permute options <$> getArgs + unless (null errs) $ mapM_ (IO.hPutStrLn stderr) errs >> exitFailure + + conf <- foldl (>>=) (return defaultConf) opts + when (null samples) $ IO.hPutStrLn stderr "need (at least) one sample name" >> exitFailure + + forM_ samples $ \sample -> do + meta <- readMetadata (conf_metadata conf) + + case H.lookup (fromString sample) meta of + Nothing -> IO.hPutStrLn stderr $ "unknown sample " ++ show sample + Just smp -> main' conf sample smp (conf_regions conf) + +-- | Call for a given sample and a set of regions defined by a regex. +-- Input are the av files whose keys match the region regex, output is +-- generated schematically from the keys so that we get one bcf output +-- for every av input. Divergence parameters for each av file are the +-- first set whose key interpreted as a regex matches the key for the av +-- file. +main' :: Conf -> String -> Sample -> Regex -> IO () +main' Conf{..} sample_name smp rgnex = + forM_ (filter (regMatch rgnex . unpack . fst) . H.toList $ sample_avro_files smp) $ \(rgn, avfile) -> + case fmap point_est $ H.foldrWithKey (ifMatch rgn) Nothing (sample_divergences smp) of + Just (prior_div : prior_het : _prior_het2 : more) -> do + liftIO $ conf_report $ "Calling " ++ sample_name ++ "/" ++ unpack rgn ++ "." + let prior_indel = case more of [] -> prior_div * 0.1 ; p : _ -> p + infile = takeDirectory conf_metadata unpack avfile + outfile = takeDirectory conf_metadata conf_output sample_name rgn + tmpfile = outfile ++ ".#" + + bracket (openFd tmpfile WriteOnly (Just 0o666) defaultFileFlags) closeFd $ \ofd -> + enumFile defaultBufSize infile >=> run $ + joinI $ readAvroContainer $ \av_meta -> + joinI $ progressPos getpos "calling at " conf_report (getRefseqs av_meta) $ + bcf_to_fd ofd (getRefseqs av_meta) [fromString sample_name] + (call (prior_div/3) prior_het, call prior_indel prior_het) + + let upd_bcf_files f s = s { sample_bcf_files = f $ sample_bcf_files s } + ins_bcf_file = upd_bcf_files $ H.insert rgn (fromString outfile) + updateMetadata (H.adjust ins_bcf_file (fromString sample_name)) conf_metadata + renameFile tmpfile outfile + + _ -> fail $ sample_name ++ "/" ++ unpack rgn ++ " is missing divergence information" + where + call :: Double -> Double -> U.Vector (Prob' Float) -> Int + call prior prior_h lks = U.maxIndex . U.zipWith (*) lks $ + U.replicate (U.length lks) (toProb . realToFrac $ prior_h * prior) + U.// [ (0, toProb . realToFrac $ 1-prior) ] + U.// [ (i, toProb . realToFrac $ (1-prior_h) * prior) + | i <- takeWhile (< U.length lks) (scanl (+) 2 [3..]) ] + + getpos :: GenoCallBlock -> (Refseq, Int) + getpos b = (reference_name b, start_position b) + + ifMatch :: Text -> Text -> a -> Maybe a -> Maybe a + ifMatch r k v a = if regMatch (regComp (unpack k)) (unpack r) then Just v else a + + +-- | Generates BCF and writes it to a 'Handle'. For the necessary VCF +-- header, we get the /names/ of the reference sequences from the Avro +-- schema and the /lengths/ from the biohazard.refseq_length entry in +-- the meta data. +bcf_to_fd :: MonadIO m => Fd -> Refs -> [S.ByteString] -> CallFuncs -> Iteratee [GenoCallBlock] m () +bcf_to_fd hdl refs name callz = + toBcf refs name callz ><> encodeBgzfWith 9 =$ + mapChunksM_ (\s -> liftIO $ S.unsafeUseAsCStringLen s $ \(p,l) -> + fdWriteBuf hdl (castPtr p) (fromIntegral l)) + + + +type CallFunc = U.Vector (Prob' Float) -> Int +type CallFuncs = (CallFunc, CallFunc) + +vcf_header :: Refs -> [S.ByteString] -> Push +vcf_header refs smps = foldr (\a b -> pushByteString a <> pushByte 10 <> b) mempty $ + [ "##fileformat=VCFv4.2" + , "##INFO=" + , "##INFO=" + , "##FORMAT=" + , "##FORMAT=" + , "##FORMAT=" + , "##FORMAT=" ] ++ + [ S.concat [ "##contig=" ] | s <- toList refs ] ++ + [ S.intercalate "\t" $ "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT" : smps ] + + +-- XXX Ploidy is being ignored. +toBcf :: Monad m => Refs -> [S.ByteString] -> CallFuncs -> Enumeratee [GenoCallBlock] Push m r +toBcf refs smps (snp_call, indel_call) = eneeCheckIfDone go + where + go k = mapChunks (foldMap encode) . k $ Chunk hdr + + hdr = pushByteString "BCF\2\2" <> setMark <> + vcf_header refs smps <> pushByte 0 <> endRecord + + encode :: GenoCallBlock -> Push + encode GenoCallBlock{..} = mconcat $ zipWith (encode1 reference_name) [start_position..] called_sites + + encode1 :: Refseq -> Int -> GenoCallSite -> Push + encode1 ref pos site = + encodeSNP site ref pos snp_call <> + case indel_variants site of + [ ] -> mempty + [_] -> mempty + v:_ | U.null d && U.null i -> encodeIndel site ref pos indel_call + | otherwise -> error "First indel variant should always be the reference." + where + IndelVariant (V_Nucs d) (V_Nuc i) = v + + +encodeSNP :: GenoCallSite -> Refseq -> Int -> CallFunc -> Push +encodeSNP site = encodeVar (map S.singleton alleles) (snp_likelihoods site) (snp_stats site) + where + -- Permuting the reference allele to the front sucks. Since + -- there are only four possibilities, I'm not going to bother + -- with an algorithm and just open-code it. + alleles | ref_allele site == nucsT = "TACG" + | ref_allele site == nucsG = "GACT" + | ref_allele site == nucsC = "CAGT" + | otherwise = "ACGT" + +encodeIndel :: GenoCallSite -> Refseq -> Int -> CallFunc -> Push +encodeIndel site = encodeVar alleles (indel_likelihoods site) (indel_stats site) + where + -- We're looking at the indel /after/ the current position. + -- That's sweet, because we can just prepend the current + -- reference base and make bcftools and friends happy. Longest + -- reported deletion becomes the reference allele. Others may + -- need padding. + rallele = snd $ maximum [ (U.length r, r) | IndelVariant (V_Nucs r) _ <- indel_variants site ] + alleles = [ S.pack $ showNucleotides (ref_allele site) : show (U.toList a) ++ show (U.toList $ U.drop (U.length r) rallele) + | IndelVariant (V_Nucs r) (V_Nuc a) <- indel_variants site ] + +encodeVar :: [S.ByteString] -> U.Vector Mini -> CallStats -> Refseq -> Int -> CallFunc -> Push +encodeVar alleles likelihoods CallStats{..} ref pos do_call = + setMark <> setMark <> -- remember space for two marks + b_share <> endRecordPart1 <> -- store 1st length and 2nd mark + b_indiv <> endRecordPart2 -- store 2nd length + where + b_share = pushWord32 (unRefseq ref) <> + pushWord32 (fromIntegral pos) <> + pushWord32 0 <> -- rlen?! WTF?! + pushFloat gq <> -- QUAL + pushWord16 2 <> -- ninfo + pushWord16 (fromIntegral $ length alleles) <> -- n_allele + pushWord32 0x04000001 <> -- n_fmt, n_sample + pushByte 0x07 <> -- variant name (empty) + foldMap typed_string alleles <> -- alleles + pushByte 0x01 <> -- FILTER (an empty vector) + + pushByte 0x11 <> pushByte 0x01 <> -- INFO key 0 (MQ) + pushByte 0x11 <> pushByte rms_mapq <> -- MQ, typed word8 + pushByte 0x11 <> pushByte 0x02 <> -- INFO key 1 (MQ0) + pushByte 0x12 <> pushWord16 (fromIntegral reads_mapq0) -- MQ0 + + b_indiv = pushByte 0x01 <> pushByte 0x03 <> -- FORMAT key 2 (GT) + pushByte 0x21 <> -- two uint8s for GT + pushByte (2 + 2 * fromIntegral g) <> -- actual GT + pushByte (2 + 2 * fromIntegral h) <> + + pushByte 0x01 <> pushByte 0x04 <> -- FORMAT key 3 (DP) + pushByte 0x12 <> -- one uint16 for DP + pushWord16 (fromIntegral read_depth) <> -- depth + + pushByte 0x01 <> pushByte 0x05 <> -- FORMAT key 4 (PL) + ( let l = U.length lks in if l < 15 + then pushByte (fromIntegral l `shiftL` 4 .|. 2) + else pushWord16 0x02F2 <> pushWord16 (fromIntegral l) ) <> + pl_vals <> -- vector of uint16s for PLs + + pushByte 0x01 <> pushByte 0x06 <> -- FORMAT key 5 (GQ) + pushByte 0x11 <> pushByte gq' -- uint8, genotype + + rms_mapq = round $ sqrt (fromIntegral sum_mapq_squared / fromIntegral read_depth :: Double) + typed_string s | S.length s < 15 = pushByte (fromIntegral $ S.length s `shiftL` 4 .|. 0x7) <> pushByteString s + | otherwise = pushByte 0xF7 <> pushByte 0x03 <> pushWord32 (fromIntegral $ S.length s) <> pushByteString s + + pl_vals = U.foldr ((<>) . pushWord16 . round . max 0 . min 0x7fff . (*) (-10/log 10) . unPr . (/ lks U.! maxidx)) mempty lks + + lks = U.map (Pr . negate . mini2float) likelihoods :: U.Vector (Prob' Float) + maxidx = U.maxIndex lks + + gq = -10 * unPr (U.sum (U.ifilter (\i _ -> i /= maxidx) lks) / U.sum lks) / log 10 + gq' = round . max 0 . min 127 $ gq + + callidx = do_call lks + h = length $ takeWhile (<= callidx) $ scanl (+) 1 [2..] + g = callidx - h * (h+1) `div` 2 + + diff --git a/tools/wiggle-coverage.hs b/tools/wiggle-coverage.hs deleted file mode 100644 --- a/tools/wiggle-coverage.hs +++ /dev/null @@ -1,38 +0,0 @@ -{-# LANGUAGE BangPatterns #-} -import Bio.Bam.Header -import Bio.Bam.Reader -import Bio.Bam.Rec -import Bio.Base -import Bio.Iteratee - -main :: IO () -main = mergeDefaultInputs combineCoordinates >=> run $ \hdr -> - joinI $ filterStream (not . isUnmapped . unpackBam) $ - joinI $ groupStreamOn (b_rname . unpackBam) (cov_to_wiggle hdr) $ - skipToEof - -cov_to_wiggle :: MonadIO m => BamMeta -> Refseq -> m (Iteratee [BamRaw] m ()) -cov_to_wiggle hdr rname = return $ liftI step - where - step (EOF mx) = idone () (EOF mx) - step (Chunk [ ]) = liftI step - step (Chunk (x:xs)) = do - let sid = unpackSeqid . sq_name $ meta_refs hdr `getRef` rname - liftIO $ putStr $ "chrom=" ++ sid ++ " start=" ++ shows (b_pos $ unpackBam x) " step=1\n" - step' (0::Int) [] (b_pos $ unpackBam x) (Chunk (x:xs)) - - step' !cov (e:ends) p str | e == p = step' (cov-1) ends p str - - step' !cov ends p (Chunk [ ]) = liftI (step' cov ends p) - step' !cov ends p (Chunk (x:xs)) | b_pos y == p = let !e' = b_pos y + alignedLength (b_cigar y) - in step' (cov+1) (ins e' ends) p (Chunk xs) - where y = unpackBam x - - step' _ [ ] _ str = step str - step' !cov ends p str = do liftIO $ putStrLn $ show cov - step' cov ends (p+1) str - - ins a [] = [a] - ins a (b:bs) | a <= b = a : b : bs - | otherwise = b : ins a bs -