Sequence length variation in dada2

[JoBrz]

Hello,
at first, I would like to thank you for dada2. I try to analyze microbiome of soil samples using this software and I have found some problems. Therefore, I would like to ask you for some advice.
The 341F and 785R primers were used to amplified V3-V4 region. It was PE 2x250 bp sequencing on Illumina. I received fastq files after demultiplexing.
So I used:

out <- filterAndTrim(fnFs, filtFs, fnRs, filtRs, truncLen=c(245,245), trimLeft=c(17,21), maxN=0, maxEE=c(2,5), truncQ=2, rm.phix=TRUE, compress=TRUE, multithread=TRUE) # On Windows set multithread=FALSE.

                                                reads.in     reads.out 

BA-C1_S3_L001_R1_001.fastq 181045 157561 (87%)
BA-C2_S4_L001_R1_001.fastq 181708 149812 (82%)
BS_S2_L001_R1_001.fastq 176511 151360 (86%)
C1_S5_L001_R1_001.fastq 303361 273360 (90%)
control_S1_L001_R1_001.fastq 223694 193257 (86%)

I tried truncLen=c(250,240) but the results were worse. Then,

errF <- learnErrors(filtFs,nbases=2e8, multithread=TRUE)
errR <- learnErrors(filtRs,nbases=2e8, multithread=TRUE)
dada-class: object describing DADA2 denoising results
1001 sequence variants were inferred from 46392 input unique sequences.
Key parameters: OMEGA_A = 1e-40, OMEGA_C = 1e-40, BAND_SIZE = 16
mergers <- mergePairs(dadaFs, derepFs, dadaRs, derepRs, verbose=TRUE)
126338 paired-reads (in 3512 unique pairings) successfully merged out of 150445 (in 14585 pairings) input.
116622 paired-reads (in 3120 unique pairings) successfully merged out of 142473 (in 19071 pairings) input.
123671 paired-reads (in 4339 unique pairings) successfully merged out of 144841 (in 16839 pairings) input.
232560 paired-reads (in 6805 unique pairings) successfully merged out of 242312 (in 11466 pairings) input.
147899 paired-reads (in 5076 unique pairings) successfully merged out of 182765 (in 20537 pairings) input.
head(mergers[[1]])
1 TGGGGAATATTGGACAATGGGCGCAAGCCTGATCCAGCCATGCCGCGTGAGTGAAGAAGGCCTTCGGGTTGTAAAGCTCTTTTGTCCGGAAAGAAAAGCTTTCGGTTAATACCCGGAAGTCCTGACGGTACCGGAAGAATAAGCACCGGCTAACTTCGTGCCAGCAGCCGCGGTAATACGAAGGGTGCAAGCGTTACTCGGAATTACTGGGCGTAAAGCGTGCGTAGGTGGTTCGTTAAGTCTGATGTGAAAGCCCTGGGCTCAACCTGGGAATTGCACTGGATACTGGCGGGCTAGAGTGCGGTAGAGGATGGCGGAATTCCCGGTGTAGCAGTGAAATGCGTAGAGATCGGGAGGAACATCCGTGGCGAAGGCGGCCATCTGGACCAGCACTGACACTGAGGCACGAAAGCGTGGGGAGCAAACA
2 TGAGGAATATTGGTCAATGGACGCAAGTCTGAACCAGCCACGTCGCGTGAAGGAAGACGGCCCTACGGGTTGTAAACTTCTTTTGTAAGGGAATAAAGTGCAGTACGTGTACTGTTTTGCATGTACCTTACGAATAAGGATCGGCTAACTCCGTGCCAGCAGCCGCGGTAATACGGAGGATCCAAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGCAGGCGGAAGATTAAGTCGGCGGTGAAATTTTGCAGCTCAACTGTAAAAGTGCCTTCGAAACTGGTTTTCTTGAGTGTGGATGAAGTAGGCGGAATTTGTGGTGTAGCGGTGAAATGCATAGATATCACGAGGAACTCCGATTGCGAAGGCAGCTTACTAAGCCATAACTGACGCTCAGGCACGAAGGCGTGGGGATCAAACA
3 TGAGGAATATTGGTCAATGGGCGCAAGCCTGAACCAGCCATCCCGCGTGAAGGAAGACGGCCCTATGGGTTGTAAACTTCTTTTGTACATGAAGAACGTCTCTTACGTGTAAGGGAGTGACGGTAATGTACGAATAAGCATCGGCTAACTCCGTGCCAGCAGCCGCGGTAATACGGAGGATGCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTACGTAGGTGGCCAGTTAAGTCAGTGGTGAAAACCTGTCGCTCAACGATAGGCGTGCCATTGATACTGATTGGCTTGAGTACAGATGAGGTAAGCGGAATGTGTAGTGTAGCGGTGAAATGCATAGATATTACACAGAACTCCGATTGCGAAGGCAGCTTACTAAACTGACACTGACGCTGAGGTACGAAAGCGTGGGTAGCGAACA
4 TGAGGAATATTGGTCAATGGACGCAAGTCTGAACCAGCCACGTCGCGTGAAGGAAGACGGCCCTACGGGTTGTAAACTTCTTTTGTAAGGGAATAAAGTGCAGTACGAGTACTGTTTTGCATGTACCTTACGAATAAGGATCGGCTAACTCCGTGCCAGCAGCCGCGGTAATACGGAGGATCCAAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGCAGGCGGAAGATTAAGTCGGCGGTGAAATTTTGCAGCTCAACTGTAAAAGTGCCTTCGAAACTGGTTTTCTTGAGTGTGGATGAAGTAGGCGGAATTTGTGGTGTAGCGGTGAAATGCATAGATATCACGAGGAACTCCGATTGCGAAGGCAGCTTACTAAGCCATAACTGACGCTCAGGCACGAAGGCGTGGGGATCAAACA
5 TGAGGAATATTGGTCAATGGGCGCAAGCCTGAACCAGCCACGTCGCGTGAAGGAAGACGGCCCTACGGGTTGTAAACTTCTTTTGTAAGGGAATAAAGTGTAGTACGTGTACTATTTTGCATGTACCTTACGAATAAGGATCGGCTAACTCCGTGCCAGCAGCCGCGGTAATACGGAGGATCCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGCAGGCGGGAGATTAAGTCGGCGGTGAAATTTTGCAGCTCAACTGTAAAAGTGCCTTCGAAACTGGTTTCCTTGAGTGTGGATGAAGTAGGCGGAATTTGTGGTGTAGCGGTGAAATGCATAGATATCACGAGGAACTCCGATTGCGCAGGCAGCTTACTAAGCCATAACTGACGCTCAGGCACGAAGGCGTGGGGATCAAACA
6 TGGGGAATATTGGACAATGGGGGCAACCCTGATCCAGCCATACCGCGTGTGTGAAGAAGGCCTTCGGGTTGTAAAGCACTTTCAGTGGGGAGGAAAAGTTGCCGGCGAATACCCGGCGACCGTGACGCTACCTGCAGAAGAAGCACCGGCAAACTCCGTGCCAGCAGCCGCGGTAATACGGAGGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGCGCGCGTAGGCGGTTTGGTAAGCTGGATGTGAAAGCCCTGGGCTCAACCTGGGAACTGCGTCCAGAACTGCCAAGCTAGAGTATGGTAGAGGGCAGTGGAATTTCCGGTGTAGCGGTGAAATGCGTAGAGATCGGAAGGAACACCAGTGGCGAAGGCGACTGCCTGGACCAATACTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACA
abundance forward reverse nmatch nmismatch nindel prefer accept
1 15093 1 2 25 0 0 1 TRUE
2 9849 2 1 30 0 0 1 TRUE
3 6097 4 4 30 0 0 1 TRUE
4 5605 3 1 30 0 0 2 TRUE
5 3686 5 3 30 0 0 1 TRUE
6 2242 9 6 25 0 0 1 TRUE
seqtab <- makeSequenceTable(mergers)
dim(seqtab)
[1] 5 20976
table(nchar(getSequences(seqtab)))

258 259 260 261 317 320 342 367 373 374 382 383 384
1 4 21 1 2 1 3 2 2 1 1 1 1
386 390 393 396 397 400 401 402 403 404 405 406 407
1 1 1 6 3 2 749 4020 409 164 128 28 401
408 409 410 411 412 413 414 415 416 417 418 419 420
23 81 15 19 4 2 3 10 6 5 3 181 70
421 422 423 424 425 426 427 428 429 430 431 436 437
120 1239 49 118 28 622 10629 1723 64 3 1 2 2

tapply(colSums(seqtab), nchar(colnames(seqtab)), sum)
258 259 260 261 317 320 342 367 373 374 382
3 12 177 5 11 2 41 7 33 27 56
383 384 386 390 393 396 397 400 401 402 403
5 12 64 3 33 27 777 16 15761 164020 32315
404 405 406 407 408 409 410 411 412 413 414
19478 10242 2769 27413 713 2486 1112 933 357 142 211
415 416 417 418 419 420 421 422 423 424 425
305 207 178 131 13227 4691 3218 79321 715 7079 1459
426 427 428 429 430 431 436 437
13813 288531 53348 1343 7 97 59 98
The amplicon should be 444 long, so minus 38 nucleotides from both primers, the peak should be 406 nucs long. Here the highest peak is from 427 (like 444-17=427).
I run BLAST as you recommended to for some sequence from the range 406-408 , 426-428 and all of them belonged to Bacteria. Please, could you give me some advice what I should do?

Thank you,
Joanna

plot.errF_245.pdf
plot.errR_245.pdf
qualplot_filtFs245.pdf
qualplot_filtRs245.pdf
qualplot_fnFs.pdf
qualplot_fnRs.pdf
seqComplexity.pdf

[benjjneb]

The V3V4 region has significant length variability, and the biological length variation is actual bimodal with a peak at ~440 and another ~460 (before removing the primers). So the two peaks you are seeing after removing primers are exactly what is expected.

The much shorter sequences are common as well, and likely represent off-target amplifcation of some sort. You can just remove all ASVs shorter than say 390 to get rid of most of that.