Inconsitent inference of ITS ASVs after upgrade to 1.16 (many more ASVs)

[nr0cinu]

Hi! :)

I upgraded our DADA2 pipeline from R 3.6 to R 4.0 which included upgrading Bioconductor and the upgrade of DADA2 from 1.14.1 to 1.16.0. I have (re-)analysed several 16S rRNA gene datasets with the updated R packages, which looked fine and resulted only in very minor differences.

However, when I (re-)analysed ITS datasets, I suddenly observed many more ASVs:

For example, these are the difference in a ITS1 dataset:

DADA2 v1.14: resulted in 400 ASVs, 7 of which were classified as gSaccharomyces scerevisiae (the most dominant taxon in the dataset) DADA2 v1.16: resulted in 7945 ASVs, 7460 of which were classified as gSaccharomyces scerevisiae. Many of these ASVs are only one mismatch different from the previously found ASVs. Most of the new ASVs also only have very little reads.

The final yield also went down from 73% to 58%:

DADA2 v1.14:

DADA2 v1.16:

The difference starts happing at the error modelling:

DADA2 v1.14:

DADA2 v1.16:

This is from read 1, read 2 looks similar.

One difference to 16S that I can think of is that the input FASTQ reads are of variable length, because I trim away the reverse primer using cutadapt before filtering. Otherwise I use the same pipeline for both types of amplicons.

I have changed nothing in the pipeline, only upgraded R and the R packages. I essentially follow https://f1000research.com/articles/5-1492 and use pool=TRUE.

I am not sure what is happening here? Any advise is appreciated!

Thanks!

Best, Bela

[benjjneb]

My guess would be that somehow there is a change in the input reads, i.e. after cutadapt and filtering.

Did you use the same filtered/cutadapted files in both cases? Or were those regenerated via the pipeline from the raw data?

If the second, can you post plotQualityProfile of the same sample from both the 1.14 pipeline and the 1.16 pipeline?

[nr0cinu]

Thank you for the quick answer!

I use the same version of cutadapt and the same commands in both. I call cutadapt from the R script, but there the report/summary is identical, and the number of reads after filtering is also identical in both versions.

As far as I can see, everything is identical before denoising.

Here are the plotQualityProfiles for one of the samples (read 1, read 1 filtered, read 2, read 2 filtered)

1.14:

1.16:

Thanks, Bela

[nr0cinu]

PS: in the above sample, there weren’t many primers trimmed, but there are also samples, were the primer trimming creates reads shorter than the length trim, e.g.:

[benjjneb]

Can you put together a minimal example?

I.e. starting from the same cutadapted/filtered sample (to factor out cutadapt/filtering), run denoising/merging in 1.14.1 and in 1.16, and report the results.

This can actually be done w/in a single R script by using devtools::install_github("benjjneb/dada2", ref="1.14.1") and devtools::install_github("benjjneb/dada2", ref="1.16") before the blocks of code doing version-specific analysis.

[nr0cinu]

Dear Benjamin,

I have narrowed the problem down to dada(…, multithread = TRUE). If multithreading is activated, the error models suddenly look different (worse) which is presumably why I get so many new ASVs. multithread = TRUE did not have this effect in v1.14.1. Please see the linked R report where I run the same commands with multithread = FALSE/TRUE with both DADA2 versions: https://gist.github.com/and3k/1d58564b500236b94a66b86bff58fc90

Some notes:

• This minimal example was run on a small subset of my data. If I use even less data the effect is less strong it seems, and with all samples the effect is stronger.
• The problem seems to be limited to read 1 in this subset. In the full dataset we also include FASTQ files which contain reverse complementary amplicons, and there I see both read 1 and read 2 effected.
• multithread = 1 has the same effect as multithread = TRUE

Thank you!

Best, Bela

[benjjneb]

Thanks for the report. I can reproduce this inconsistency, and have localized it to this commit: https://github.com/benjjneb/dada2/commit/c1adc9403fc5b8a988fde60f84283566f476747c

Perhaps also relevant is the previous (meaningful) commit: https://github.com/benjjneb/dada2/commit/d350746175d300f924fb7e41f59b5209119bec99

Good news is that the inconsistency between single- and multi-threading operation doesn't appear to affect denoising directly, it only affects the enumeration of the errors after denoising -- but this does feed back into denoising indirectly by affecting the estimated error rates.