第4章 序列去重与ASV/OTU聚类：DADA2、VSEARCH的参数选择原则

序列去重和聚类是物种丰度检测中承上启下的关键一步。这一步做不好，前面质控白费，后面注释也容易跑偏。我接触过不少项目，很多人觉得“直接用默认参数跑就完事了”——结果一检查，有的把真实变异当成测序错误给合并了，有的反而把噪音扩增成假ASV。说实话，参数调对了，效果能差出一倍。

这一章我们聚焦两个主流工具：DADA2（ASV级别）和VSEARCH（OTU级别）。我会把每个关键参数背后的逻辑讲清楚，再配上我自己的踩坑经验。

4.1 核心概念：去重、ASV、OTU的区别

先理清几个概念，不然参数调起来容易乱。

• 去重（Dereplication）：把完全相同的序列合并成一条，记录出现频次。这一步能极大压缩数据量。
• ASV（Amplicon Sequence Variant）：DADA2 的做法，区分单碱基差异，相当于把每个真实生物序列当作一个“物种”。分辨率极高，但计算量也大。
• OTU（Operational Taxonomic Unit）：VSEARCH 的做法，按相似度（通常 97%）把序列抱团。牺牲一点分辨率换来稳定性和计算速度。

我个人的习惯：如果是高变异区域（比如16S V3-V4），我倾向用 ASV；如果是全长ITS或者比较保守的片段，OTU 反而更稳——因为同源基因差异本来就不大，强行分 ASV 容易把 PCR 引入的假变异也卷进去。

4.2 DADA2 参数选择：从误差模型到合并

DADA2 的核心逻辑是“用数据自身学习测序错误模型”，然后推断哪些序列是真实的生物变异。它的参数集中在几个函数里：learnErrors、derepFastq、dada、mergePairs、makeSequenceTable。我们一个一个说。

4.2.1 误差模型学习：learnErrors

这个函数负责从序列中统计错配概率。关键参数就是 nbases——用多少碱基来构建模型。我一般这样设：

err <- learnErrors(filtFs, nbases=1e8, randomize=TRUE)

nbases 默认是 1e8，也就是一亿个碱基。但要注意：如果你的样本数少，或者测序深度特别浅（比如只有几万条序列），用默认值会导致模型估计偏。我遇到过一种情况：一个临床样本只有3万条read，用1e8时模型错得离谱，ASV数量暴增。后来我把 nbases 降到 1e6，效果才正常。经验公式：nbases 大约等于你所有样本总碱基数的 30%~50%。

另一个参数 randomize=TRUE 我建议打开。默认是 FALSE，但随机化能让模型不偏向序列出现的顺序——尤其是当你的样本是按实验条件分批上机时，不随机化可能导致批次偏差被学到模型里。

4.2.2 去重与样本内推断：derepFastq 和 dada

derepFastq 就是把完全相同的序列合并。有一个参数容易被忽略：qualityType。默认是 "FastqQuality"，但如果你用的是 Sanger 编码的 Phred+33，这个是对的。如果是 Illumina 1.3+ 的老数据（Phred+64），就得改成 "FastqQualityPhred64"。不然质量值全偏，后面误差模型直接崩。

接着是 dada 函数，它负责做核心的ASV推断。关键参数是 HOMOPOLYMER_GAP_PENALTY 和 BAND_SIZE：

dadaFs <- dada(derepFs, err=err,
               HOMOPOLYMER_GAP_PENALTY=-1,
               BAND_SIZE=32)

• HOMOPOLYMER_GAP_PENALTY：同聚物区域的空位罚分。默认是 -1（罚分）。我建议如果你研究的是高GC区域（比如放线菌16S），可以降到 -2，因为同聚物在该区域更容易出现——罚分重一点能减少误分。
• BAND_SIZE：动态规划时考虑的比对窗口宽度。默认 16，我一般升到 32。为什么呢？因为在长读长（比如 PacBio CCS）或者高变区序列里，真实差异可能跨越几十个碱基，窗口小了可能把大的插入缺失当成多个独立错误，导致分裂出假ASV。但注意：BAND_SIZE 越大，计算时间线性增长。数据量大的样本（比如超百万序列）建议保持 16，否则内存吃不住。

4.2.3 双端合并：mergePairs

这一步把双端 read 拼成完整的扩增子。参数 minOverlap 和 maxMismatch 最重要。

mergers <- mergePairs(dadaFs, derepFs, dadaRs, derepRs,
                      minOverlap=12, maxMismatch=2)

minOverlap 默认是 12 个碱基。如果重叠区太短，拼接后的序列可信度低。我一般设 20，但要注意：如果读长很长（比如 300bp×2），重叠区本身就大，设 20 没问题；如果是短读（150bp×2），重叠区可能才 30 多个碱基，设 20 就太苛刻了，会导致大量序列无法合并。此时我改成 10~12。

maxMismatch 控制重叠区允许的错配数。默认 0，我通常设 2。为什么？因为 PCR 扩增时可能在同一位点引入不同碱基（嵌合体或随机错误），设 0 会把这些合理变异也砍掉。当然，如果你的扩增子高度保守，可以保持 0。

4.3 VSEARCH 参数选择：从去重到聚类

VSEARCH 是 usearch 的开源替代，速度快、内存小。它做 OTU 聚类时，最常用的命令是 --cluster_fast 和 --cluster_size。

4.3.1 去重：--derep_fulllength

VSEARCH 去重很直接：

vsearch --derep_fulllength seqs.fasta \
        --output derep.fasta \
        --sizeout \
        --minseqlength 100

• --sizeout：把每条唯一序列的出现次数写入序列标题里（方便后续聚类算法依丰度排序）。
• --minseqlength：过滤太短的序列。我一般设 100 bp（取决于你的扩增子最小长度），太短的序列大概率是 PCR 引物二聚体或降解片段。

4.3.2 聚类：--cluster_fast 与 --cluster_size

两个常用模式：

• --cluster_fast：按序列长度从大到小排序，以最长的序列为种子，把相似度≥阈值的序列归为一类。速度快，适合大样本。
• --cluster_size：按序列丰度排序，优先把高丰度序列当作种子。如果样本里有些真实物种丰度低但很重要，这个模式容易把它们漏掉——因为低丰度序列会先被高丰度的“吞掉”。

我个人在环境样本（土壤、肠道）里常用 --cluster_fast，因为真实多样性很高，用长度排序更保守。在临床样本（比如血液/脑脊液，多样性低）里，我有时改用 --cluster_size，能更快去除噪音。

# 常用聚类命令（97%相似度）
vsearch --cluster_fast derep.fasta \
        --id 0.97 \
        --centroids centroids.fasta \
        --uc cluster.uc \
        --sizein \
        --sizeout

参数 --id 是核心。默认 0.97，但不同区域应该不同：

4.4 章节知识框架

4.5 实战步骤：我的配置流程

下面是我在项目中实际使用的流程，不绕弯子：

1. 预览数据量：先看每个样本有多少条序列，读长分布。如果总碱基数少于 1e7，就把 nbases 设成 1e6 或更少。
2. 跑 DADA2 默认参数：用 learnErrors 默认参数，出图看误差模型是否平滑。如果不平滑，调 nbases 或 randomize=TRUE。
3. 检查 ASV 数量：如果 ASV 数量异常多（比如比预期多 5 倍），缩小 BAND_SIZE 试试，或者调高 HOMOPOLYMER_GAP_PENALTY。
4. 合并后统计：跑 mergePairs 后看合并成功率。如果低于 60%，降低 minOverlap 或增加 maxMismatch。
5. 做 VSEARCH 对比：对同一个样本用 --id 0.97 跑一遍，看 OTU 数量与 ASV 数量的比值。如果 ASV:OTU 超过 5:1，说明 ASV 分裂太严重，需要收紧参数。
6. 最终确定参数：选一个合理的参数组合，然后用 --nonchimeras 做嵌合体过滤（DADA2 自带 removeBimeraDenovo），再跑一次。

4.6 一个真实案例的教训

我曾经处理过一个海洋沉积物样本，16S V4 区，双向 250bp 测序。用 DADA2 默认参数跑出来 7600 个 ASV，这显然太多了——正常海洋沉积物 OTU 应该 1000~2000 左右。我一开始怀疑是测序错误太多，但看了误差图发现模型拟合良好。后来检查 mergePairs 的合并成功率只有 45%，说明很多双端 read 没拼上。

我试着把 minOverlap 从 12 降到 8，合并率升到 68%，ASV 数量骤降到 4100。再把 maxMismatch 从 0 改成 2，合并率到 82%，ASV 数量落到 2800。最后用 VSEARCH 以 0.97 聚类，最终 OTU 只有 1350 个——和文献报道一致。你看，其实就是两个参数的问题。