Very low percentage of merged reads in DADA2 with 16S V4, 2×150 bp paired-end sequencing

[Sergio]

Question

Hi,

I am processing 16S rRNA amplicon sequencing data using QIIME 2 and DADA2.

My samples were amplified targeting the V4 region of the 16S rRNA gene using the 515F and 806R primers. The libraries were sequenced on an Illumina MiSeq i100 using 2×150 bp paired-end sequencing.

The expected V4 amplicon is approximately 250 bp, so theoretically I should have around 50 bp of overlap between R1 and R2:

150 + 150 − 250 ≈ 50 bp overlap

However, I am obtaining a very low percentage of merged reads after DADA2.

he main loss therefore appears to occur during the merging step, rather than during filtering or denoising.

I would like to understand what could be causing this low merging rate.

Could this be related to:

• the actual amplicon length being different from the expected ~250 bp?

• insufficient overlap after trimming?

• primers still being present?

• the quality of the last bases of R1/R2?

• another issue with the paired-end reads?

I would appreciate any advice on how to diagnose the cause before changing my DADA2 parameters.

I used the following DADA2 command:

qiime dada2 denoise-paired \
  --i-demultiplexed-seqs demux-paired-end.qza \
  --p-trim-left-f 0 \
  --p-trim-left-r 0 \
  --p-trunc-len-f 150 \
  --p-trunc-len-r 150 \
  --o-base-transition-stats base-transition-stats.qza \
  --o-representative-sequences rep-seqs.qza \
  --o-table table.qza \
  --o-denoising-stats denoising-stats.qza \
  --p-n-threads 1 \
  --verbose

I truncated both reads at 150 bp because the quality of R2 drops at position 151. The quality of the preceding bases appears to be good, so I did not want to truncate the reads further. I also did not trim any bases from the beginning of either read.

[SoilRotifer (Mike Robeson)]

Hi @Sergio

@Sergio

he expected V4 amplicon is approximately 250 bp, so theoretically I should have around 50 bp of overlap between R1 and R2:

150 + 150 − 250 ≈ 50 bp overlap

I'd still look into truncating your reads a little, given the amount of overlap you have. That is, there is going to be in increase in the chance of mismatches with longer overlap, even with high-quality data. See this thread.

I also had a similar case in which I also had to trim a little from the 5' end, with similar scores. I assume these might be from the new score/binning methods used? If so, hopefully some one else can provide recommendations / insight.

-Cheers!

[Sergio]

Thank you for your response!

I tried truncating the reads as you suggested, but I am still having trouble finding a good balance.

When I keep the reads at 150 bp, the merging rate is very low (around 7–14%). However, when I reduce the minimum overlap to 8 bp, the merging rate improves, but I get a much higher number of chimeric sequences. If I truncate the reads further, the merging rate decreases substantially again.

So I seem to be getting this trade-off:

• 150/150 bp → very low merging
• min-overlap 8 → better merging, but many more chimeras
• more truncation → merging drops again

This makes me think that the main issue could be related to the amount of overlap available between my reads. I am also confused about the expected amplicon length for 515F–806R, since I have found values ranging from approximately 250 to 290 bp. With 2×150 bp reads, that would mean approximately 50 bp versus only ~10 bp of theoretical overlap.

I really do not understand what could be causing the problem. Could the low merging rate mainly be due to insufficient overlap?

I used the following command for the version with only min-overlap 8:

qiime dada2 denoise-paired \
  --i-demultiplexed-seqs ../raw_fastq_flat/demux-paired-end.qza \
  --p-trim-left-f 0 \
  --p-trim-left-r 0 \
  --p-trunc-len-f 150 \
  --p-trunc-len-r 150 \
  --p-min-overlap 8 \
  --o-representative-sequences rep-seqs.qza \
  --o-table table.qza \
  --o-denoising-stats denoising-stats.qza \
  --o-base-transition-stats base-transition-stats.qza \
  --p-n-threads 2 \
  --verbose

And this one for the version with more aggressive truncation:

qiime dada2 denoise-paired \
  --i-demultiplexed-seqs ../raw_fastq_flat/demux-paired-end.qza \
  --p-trim-left-f 20 \
  --p-trim-left-r 20 \
  --p-trunc-len-f 145 \
  --p-trunc-len-r 145 \
  --p-min-overlap 6 \
  --o-representative-sequences rep-seqs.qza \
  --o-table table.qza \
  --o-denoising-stats denoising-stats.qza \
  --o-base-transition-stats base-transition-stats.qza \
  --p-n-threads 2 \
  --verbose

Thanks again for your help!

[SoilRotifer (Mike Robeson)]

Hi @Sergio,

Thanks for the informaiton.

@Sergio

I am also confused about the expected amplicon length for 515F–806R, since I have found values ranging from approximately 250 to 290 bp. With 2×150 bp reads, that would mean approximately 50 bp versus only ~10 bp of theoretical overlap.

The range is often due to including (~290 bp) or not including (~250 bp) the primers as part of the fragment length. That is, standard sequencing approaches sequence through the PCR primers. Meaning the PCR primers appear in the first ~20 bases of your reads at the 5' end. Often users will run cutadapt to explicitly remove these primers prior to downstream processing (search the forum for more threads on this topic).

If the approach you used does in fact sequence through the primers... then unfortunately 2x150bp is generally not long enough to merge the V4 reads. In which case I'd recommend only using the forward reads.

The only way 2x150 works well for the V4 region is if you are using the EMP 16S rRNA gene sequencing protocol, which essentially does not sequence through the PCR primers, thus allowing you to reach further down towards the 3' end, as you are not wasting precious sequencing on the PCR primer for the 2x150 run. See here, for more details.