Dada2: ITS2 with poor error model + bad chimera detection

[mestaki]

Hello,
I’m working with 2x300 PE, MiSeq ITS2 data from infant feces. I’m following along the ITS DADA2 workflow (great resource by the way!) in R 3.6 with dada 1.12 but have ran into some concerns along the way and am looking for some recommendations. The forwards seem to be in ok quality with the Reverse being expectedly a bit worse (qplot). I have removed the primers with cutadapt as per the tutorial and only modified the error model setting slightly:

errF <- learnErrors(filtFs, multithread = TRUE, nbases=2e8, randomize=TRUE)
errR <- learnErrors(filtRs, multithread = TRUE, nbases=2e8, randomize=TRUE)

The error model in the forward plot doesn't look very good:

The reverse however looks better.. oddly enough.

My first question here would be, what would be the recommendation in this situation and if there are any ways to improve the error model. Would adding trim/trunc parameters help? I avoided these as per the tutorial.

The stats summary also shows a very high level of sequences being filtered out in the initial filtering. I found it to be much higer than what I expect from typical 16S data but I thought maybe it was due to high PhiX (20%) and very short (<50) length reads, and not including any truncating parameters. But it does still seem a bit too high. Is this something to be concerned about?
Anyways, I carried on with denoising and chimera detection with pool=TRUE in both instances and created my feature-table. The sequence length distribution of the final feature-table seems a bit all over the place. I was expecting a bit of this because its ITS2 but is this much normal?

And another odditiy is that, the 73 chimeras that were detected have length distributions that look typical of the non-chimera sequences:

Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   63.0   148.0   184.0   186.3   233.0   385.0

Being a bit concerned, I blasted them all and turns out almost all of them hit very nicely with fungi (and some bacteria..<- another confusing occurance) with almost perfect % identity and alignment (report summary).
So I'm wondering if these are to be trusted or not?

Any thoughts and recommendations with regards to any of these would be greatly appreciated.

[benjjneb]

My first question here would be, what would be the recommendation in this situation and if there are any ways to improve the error model.

Your error model looks great. I agree this is very clean reverse read error model, but that's a good thing not a bad thing! I see no reason from your plotErrors plots to be concerned about the error model.

The stats summary also shows a very high level of sequences being filtered out in the initial filtering. I found it to be much higer than what I expect from typical 16S data but I thought maybe it was due to high PhiX (20%) and very short (<50) length reads, and not including any truncating parameters. But it does still seem a bit too high. Is this something to be concerned about?

Typically no this is not something to be concerned about. Filtering a high fraction of intiial low quality reads is usually more of a benefit than a detriment to the overall analysis. However, I'm a bit surprised as well with how many reads are being filtered out, but maybe I shouldn't be. What does your plotQualityProfile plots look like? And what was your exact filtering command? If your dead set on keeping more reads, you can always loosen the quality filtering parameters.

[benjjneb]

And another odditiy is that, the 73 chimeras that were detected have length distributions that look typical of the non-chimera sequences:

That in itself is not a cause for concern. Chimeras form more readily between closely related sequences, and thus tend to have length distributions that look like the lengths of the real sequences.

Being a bit concerned, I blasted them all and turns out almost all of them hit very nicely with fungi (and some bacteria..<- another confusing occurance) with almost perfect % identity and alignment (report summary).

For the reason above (Chimeras form more readily between closely related sequences) seeing chimeras that are close to other sequences in BLAST is not itself a cause for concern. You are getting a fairly high fraction of exact matches in your chimeras though (100%ID, 100%coverage). It is possible for sequences to be falsly flagged chimeric. We've found the best balance of sensitivity to specificity to come from the "consensus" method that calls chimeras based on a vote across all samples, but that does depend on you having a fair number of samples so that the vote is informative. How many samples are you processing here? And what was the exact chimera removal command you used?

Again, a simple parameter change can make chimera calling more conservative, and may be warranted in some situations (see the section in this preprint on chimera filtering for some discussion). The parameter to play with is removeBimeraDenovo(..., minFoldParentOverAbundance=XX) which you can increase up to 8 or even 16.

[mestaki]

Hi @benjjneb,
Thank you for your input!

I see no reason from your plotErrors plots to be concerned about the error model.

That's good to hear, I was just worried that in the forward reads' model the observed values were deviating away from the estiamtes starting ~ 25. Perhaps I need to reavluate my understanding of those plots a bit.

What does your plotQualityProfile plots look like?

The first link (qplot) I posted above is a qiime2 summary of the reads, if you click on the Interactive quality plot tab at the very top you'll see similar qplots to those in DADA2. I thought that would be a bit more informative as it sumsamples across all samples instead of one by one. But here is a dada2 qplot from a couple of the samples:
Forwards:

Reverse:

And what was your exact filtering command?

out <- filterAndTrim(cutFs, filtFs, cutRs, filtRs, maxN = 0, maxEE = c(2, 2), 
                     truncQ = 2, minLen = 50, rm.phix = TRUE, compress = TRUE, multithread = TRUE)

If your dead set on keeping more reads, you can always loosen the quality filtering parameters.

I'm actually pretty content with the final read counts, especially since I'm not expecting a super diverse community, my concern was more about the error-model fit.

We've found the best balance of sensitivity to specificity to come from the "consensus" method

I could try this again with the "consensus" method but I had used the "pooled" method because pooled was also used in the denoising step and your comment that the ideal situation would be to use same method in both steps. My intention is to use diversity estimation from the Breakaway package which depends on singletons/rare variants, thus the pooling options. Would setting chimera removal to consensus here be a better option?

How many samples are you processing here?

I have 62 samples including one positive and one negative samples. I imagine that's enough for informed sharing? This MiSeq run however actually had about 230 samples in it total, ~ 125 wine samples targeting ITS1, 62 samples ITS2 (infants), and another 40 adult samples ITS2, but I am processing them separately because they went through separate PCRs. Would there be any benefit including some of the other ITS2 samples which I expect to be more diverse?

And what was the exact chimera removal command you used?
seqtab.nochim <- removeBimeraDenovo(seqtab, method="pooled", multithread=T, verbose=TRUE)

Let me know what you think and the minFoldParentOverAbundance=XX seems like a good point to play around with in the mean time. Thanks!

[benjjneb]

I could try this again with the "consensus" method but I had used the "pooled" method because pooled was also used in the denoising step and your comment that the ideal situation would be to use same method in both steps. My intention is to use diversity estimation from the Breakaway package which depends on singletons/rare variants, thus the pooling options. Would setting chimera removal to consensus here be a better option?

Nope you are right on here. I didn't realize you were using the dada(..., pooled=TRUE) mode.

I have 62 samples including one positive and one negative samples. I imagine that's enough for informed sharing? This MiSeq run however actually had about 230 samples in it total, ~ 125 wine samples targeting ITS1, 62 samples ITS2 (infants), and another 40 adult samples ITS2, but I am processing them separately because they went through separate PCRs. Would there be any benefit including some of the other ITS2 samples which I expect to be more diverse?

Nope, adding different samples (either loci or types) won't help the power of pooling because they won't be informative about your focal samples.

Let me know what you think and the minFoldParentOverAbundance=XX seems like a good point to play around with in the mean time. Thanks!

That looks fine, but I'd recommend even more considering increasing minFoldParentOverAbundance now that I realize you are doing pooled processing. 8 might be a good value to try initially. I also might consider triple-checking that the primers really were removed, as this is the other potential reason for artefactual chimera identification (if there are ambig bases in the primers anyway).