Where to trim reads

[Nicholas_Bokulich]

@ariel

@Nicholas_Bokulich Do you have feedback on this?

Sorry for the delay @ariel! I have been out sick .

@WAS1

Hi, Quick one here though I am yet to take a look at the plots, your 27F primer stops at position 27!! (the primer sequence actually starts at 8F) so does not start at position 27

That's a great point! That had not occurred to me, though now that you mention I recall that this primer is often called 8F...

@WAS1

I am sorry the 534R is most probably (because I don’t know much about this one) another modification of 519R/518R which will mean it stops at 534 and start at 518

I also do not know about this one, you should check on that @ariel — though I will note it probably starts at 550 and ends at 534 (since this is a reverse primer).

@ariel

When I look at my quality plots , the quality drops to even 24 and 22 in several points in the reverse, especially after 228r

Even at 291f and 240r your median quality scores are 29 and 27, respectively, so I would not worry too much.

@ariel

Unfortunately, even if i extend to 291f, and 240 reverse, I still have 12 samples with zero merged reads, and almost 500 with fewer than 5000. Here are the stats from dada2. I’m a bit concerned about how to proceed.

For all but 1 of those 12, it looks like the # of input sequences is very very low, so that is not a merging issue and those samples would be lost anyway. For one of these samples (HMP2_J15247_M_ST_T0_B0_0120_ZWFDEY0-01_GES15-03642_L001), there is a large number of input sequences but all are dropped during denoising! So it is likely that that sample failed somehow — you may want to inspect that one individually to see what the quality profile looks like on it.

So I would say this is not a merging issue since many reads are passing, but clearly it still is since you are losing many reads at the merging step. This probably comes down to a few factors:

1. Are you sure primers are trimmed? I expect so, since we now have reads joining successfully, but you should double-check just to make sure.
2. We may want to reassess this statement:

@Nicholas_Bokulich

there is so little variation in 16S length that it does not really matter all that much compared to something like ITS.

While it is true that 16S has relatively low length heterogeneity, this may not be entirely true about all domains — my experience is mostly with V4, which has a very tight length distribution. V1-V3 may be different, and a quick search shows that V3 variability alone may have a range of ~50 nt

Also, the developer of dada2 has recommended up to 80 nt overlap as a safe overestimate whenever possible:

Quality score and trimming in Dada2

One addition to that: You need to have 20 y's, plus enough additional y's to handle the biological length variation in your amplicon. If you aim for the minimum, you will lose real variants that are a few nts shorter than the average amplicon length, and therefore systematically drop certain taxonomic groups. 20 is the bare minimum, but err on having more overlap. If you can get an overlap of 80 with acceptable quality, I would go for it.

and also that non-target DNA hit by these primers can yield potentially useful results:

Comparison of DADA2 with Deblur

This is a completely valid option (which I'm assuming there is a command to do in Q2?). "Cutting a band in silico" to choose amplicons of the expected length is a good way to clean up certain kinds of residual artefacts that can crop up. The DADA2 pipeline does not do this by default because in some cases those non-target length amplicons are real interesting things. For example, Trichomonas vaginalis, the protozoan cause of trichomoniasis, is detected by the EMP primer set, but is ~290 nts rat…

All of this above advice may or may not fit the actual length distribution of your amplicons: if you have a gel image or bioanalyzer trace or something else along those lines, use that to calculate the total starting amplicon length distribution.

@WAS1

I don’t think so. yes region with this should be 490/491, I would agree, but by truncating R at 226 after removing 17 leaves u with 209 length so overlap is not 26. the R primer is not like the F primer

I am not following this either. How is the R primer different from the F? One way or another the primer length can be subtracted from the total amplicon length if primers have been trimmed, so @ariel's calculations look correct to me — but please elaborate in case I am missing something.

@ariel

Should i at this point try to proceed with the forward reads only?

That is the last resort, and easy to do (just use dada2 denoise-single and the reverse reads will be ignored, no need to re-import)

But there are two other options:

1. as noted above, look at the actual length distribution
2. Brute force: truncate your sequences less and less (i.e., make longer input sequences) to see if you can get more joined (merged) sequences out of dada2 in spite of the quality drop off (your qualities are good, especially on the forward read. Look at the median accuracy scores to determine where to cut)
3. Try merging with q2-vsearch to get a "second opinion" (though note that many sequences may not join due to low quality at the 5' ends).

@ariel

I’m assuming this would also mean that I need to retrain the classifier with just the forward primer

No — you can use the same taxonomy classifier trained to the entire V1-V3 region, no need to train (and training with only the forward primer would not work — you would train with both primers and truncate to the same length as your sequences if you want to be really precise, but that much trimming does not matter much in my experience)

@ariel

It would be great to have explicit instructions in the tutorial on how to calculate overlap, and also great if DADA2 spit out some statistic on the overlap as well.

cc:ing @Mehrbod_Estaki — he is working on a new tutorial and this would be a great topic to cover! Thanks for the suggestion @ariel!

@Mehrbod_Estaki may also have some more advice on this, he has tons of experience wrangling with dada2 and noisy datasets!