inclusion of DADA2 EM algorithm for sequence correction

[killidude]

Is there the possibility of including a raw data correction algorithm on par with that implemented in DADA2? This algorithm would take the raw R1 and R2 FASTQ files (for computational efficiency could be deduplicated already), and use the error model derived from the raw reads themselves to correct the raw reads and then returns the corrected FASTQ file. For downstream analyses, any clustering step (e.g. unoise in vsearch or swarm) would then cluster biological variation rather than both "noise" caused by PCR and sequencing errors, and biological variation.

I realize this might be a relatively big task, but I believe it would be a very useful feature. Once the raw data have been corrected, then all unique sequences would represent "true" biological variation observed in the sample, and would be appropriate for downstream analyses without the need for additional clustering.

Thank you,

Tomas

[torognes]

Thanks for the suggestion! I'll consider whether this is possible and doable.

Some relevant links:
https://github.com/benjjneb/dada2
https://benjjneb.github.io/dada2/
https://doi.org/10.1038/nmeth.3869
https://doi.org/10.1186/1471-2105-13-283

I am not sure it is possible to estimate the error model from the raw reads without performing some kind of partitioning or clustering, similar to what is done in DADA2, but perhaps some simplification is possible.

[killidude]

Thank you @torognes.

From what I can see of the DADA2 sequence correction algorithm, the sequences are first clustered. The sequences in the cluster are then corrected to minimize the difference between the observed and expected SFS given an average per position Q score. Given that the expected SFS is calculated using a coalescent model, I suspect this is from Watterson's theta assuming intraspecific variation, thus the clusters should be small (2-3% within cluster divergence to reflect within species divergence). Clusters with larger within cluster divergence would result in the retention of more sequencing/PCR errors as true SNPs. Given that PCR errors are effectively estimated from the expected vs. observed SFS, I do not see how this can be done without clustering the data to species-level clusters first, if one uses an intraspecific coalescent model to calculate the expected SFS.

There should be no need for alignment prior to clustering since even in ribosomal DNA sequences there will be very few if any intraspecific indels, and thus all sequences of the same species flanking the same primer in the FASTQ file (coming from an Illumina machine) will already be aligned.

Thanks again,

Tomas

[cjfields]

I'd love to see this implemented, the step is a major bottleneck for us, but unless I'm oversimplifying this there are several complicating factors for any implementation which may be why you haven't seen other implementations yet. So before you jump in it'd be worth getting feedback from @benjjneb on this one, particularly if this is something already being developed.

Some complications off the top of my head:

1. The error estimation step (learnErrors), which assesses the quality information and base calls in your sequence data to estimate error rates that can be used in the following denoising step. This requires a function to model rates, and as currrently implemented this function is written in R. Additionally, the function given can vary between sequencing technologies (some like PacBio are more prone to indel errors than Illumina) and even within a technology depending on how the data are given (e.g. if the quality scores are binned or not).
2. During denoising (the dada function) there is the option of pooling data across samples to more consistently capture low abundance ASVs or optionally include prior sequence information for denoising. These are fairly commonly used, and the pseudo pooling option internally uses priors.
3. The additional options within the dada function like BAND_SIZE, OMEGA_A/P/C, etc. These may not be as big a lift (I think these get passed right to the underlying C-based functions) but they are commonly tweaked depending on the experimental data being assessed.

I'm sure I'm butchering some of the details and leaving out other considerations, but I'd def recommend asking Ben. It's his code