第3节：16S扩增子分析基础

这节我们来聊聊16S扩增子分析最底层的那些事儿。说白了，所有下游的物种组成、多样性、差异分析，质量好不好，都取决于这第一步的质控和聚类做得扎不扎实。我个人习惯把这一阶段叫“地基工程”——地基没打好，上面的房子早晚要裂。

3.1 原始数据质控

从测序公司拿到手的数据，一般就是双端FASTQ文件。我刚开始做项目时，拿到数据就急着往下跑流程，结果跑出来一堆奇怪的物种——后来才发现是接头没切干净，引入了很多假阳性。所以现在我的第一件事永远是：先看看数据长什么样。

3.1.1 质量评估：FastQC

FastQC是我个人的首选工具。它不负责修数据，只负责告诉你数据哪有问题。你想想看，数据只有几千条序列的时候还好排查，但一旦上了几百万条，肉眼基本没法判断。FastQC能帮你快速定位到低质量区间、接头残留、GC含量异常等问题。

基本用法很简单：

fastqc sample1_R1.fastq sample1_R2.fastq -o qc_results/

跑完之后会生成HTML报告。我一般重点看这么几个模块：

模块	关注点
Per base sequence quality	碱基位置上Q30比例，尾部质量下降是常态
Per sequence GC content	正常16S应该接近理论值，异常双峰提示污染
Adapter Content	接头残留——有的话必须切掉
Overrepresented sequences	可能来自引物二聚体或非目标扩增

嗯，这里要注意：FastQC只是看，不修。发现问题之后，你得用专门的工具来处理。

3.1.2 质量修剪与过滤：Trimmomatic

Trimmomatic支持多种修剪策略，我常用的是这一步：先切除接头，再滑动窗口裁剪低质量碱基，最后丢掉过短的序列。配置写成一个参数文件会更方便复现：

java -jar trimmomatic-0.39.jar PE \
  sample1_R1.fastq sample1_R2.fastq \
  sample1_R1_trimmed.fastq sample1_R1_unpaired.fastq \
  sample1_R2_trimmed.fastq sample1_R2_unpaired.fastq \
  ILLUMINACLIP:adapters.fa:2:30:10 \
  SLIDINGWINDOW:4:20 \
  MINLEN:150

这几个参数是什么意思？我解释一下：

ILLUMINACLIP：根据接头文件（adapters.fa）识别并切除接头，容许2个错配
SLIDINGWINDOW:4:20：每4个碱基算一次平均质量，低于Q20就从该位置剪掉
MINLEN:150：过滤掉修剪后短于150bp的序列——太短会导致后续比对困难

我曾经踩过的坑：双端测序的两条read长度不一致时，Trimmomatic会把短的那条对应的另一端标记为未配对。如果你后续要用DADA2合并双端，一定要保留未配对的输出文件，虽然它在合并步骤用不上，但可以作为备份。

跑完Trimmomatic后，我会再跑一次FastQC看看修剪效果。质量曲线尾部应该明显抬高，接头峰消失。

3.2 OTU/ASV 聚类

数据干净了，接下来是关键一步：如何从海量序列中归纳出代表物种的特征。这里有两个流派：传统OTU和现代ASV。我想重点聊聊后者。

3.2.1 DADA2：从错误推断真实序列

DADA2的核心思路是：测序是有错误的，但错误模型是可推断的。它不把所有序列按97%相似度硬性聚类，而是先估计每条序列来自真实生物序列的概率，再把高度相似的序列合并为ASV（Amplicon Sequence Variant）。

你在R里执行DADA2的大致流程是这样的：

library(dada2)
path <- "path/to/trimmed_fastqs"
fns <- list.files(path, pattern="_R1_trimmed.fastq")

# 从FASTQ中推断错误率
err <- learnErrors(fns, nbases=1e8, multithread=TRUE)

# 去噪：核心步骤
dadaFs <- dada(fns, err=err, pool="pseudo")

# 合并双端
mergers <- mergePairs(dadaFs, fns, reverse=sapply(fns, function(x) sub("_R1_", "_R2_", x)))

# 生成ASV表格和代表序列
seqtab <- makeSequenceTable(mergers)
seqtab.nochim <- removeBimeraDenovo(seqtab, method="consensus", multithread=TRUE)

这里我想强调一点：DADA2的pool="pseudo"参数，我推荐用这个而非默认的"pool=FALSE"。伪池化能让低丰度但更可靠的序列被保留下来，对检测稀有物种很关键。

3.2.2 Deblur：快速处理的另一个选择

Deblur是QIIME 2内置的ASV方法。它的逻辑跟DADA2不一样：先用一个固定阈值去噪，再根据序列间的碱基差异做进一步清理。速度比DADA2快，适合大规模数据集。但说实话，我对比过几次，DADA2在低丰度序列的准确性上稍微占优。所以我现在做正式项目首选DADA2，只有需要快速出初稿时才用Deblur。

你可能会问：OTU还有必要用吗？我个人建议：如果数据库支持，尽量用ASV。ASV在跨研究比较中的可重复性远高于OTU，毕竟传统的97%阈值是人定的，不同批次间的OTU无法直接比较。

3.3 物种注释

拿到ASV或OTU的代表序列后，得给它安个名字——也就是注释到哪个分类水平（门纲目科属种）。这一步依赖参考数据库的质量。我遇到过最糟心的情况：注释结果里一半是“Unclassified”，后来发现是数据库版本不对。

3.3.1 SILVA 数据库

SILVA是目前最主流的16S参考数据库。它覆盖了细菌、古菌和真核生物的SSU rRNA序列，经过人工校对，质量很高。我用的版本是SILVA 138.1，专门提取了V3-V4高变区的参考序列用于训练分类器。

在R中使用DADA2的assignTaxonomy函数时，需要提前下载训练好的贝叶斯分类器模型：

# 下载 SILVA 138.1 训练好的分类器
download.file("https://zenodo.org/record/4587955/files/silva_nr99_v138.1_train_set.fa.gz",
              destfile="silva_nr99_v138.1_train_set.fa.gz")

# 物种注释
taxa <- assignTaxonomy(seqtab.nochim, 
                       refFasta="silva_nr99_v138.1_train_set.fa.gz",
                       minBoot=80, multithread=TRUE)

minBoot=80 表示自展支持度低于80%时视为未分类。我试过降到50，分类率高了但错误率明显上升——后来说服自己，宁可少而准，不要多而杂。

3.3.2 Greengenes 数据库

Greengenes是另一个常用数据库，尤其在16S功能预测（比如PICRUSt）中必不可少。但它最后一次更新是2013年的13_8版本，已经停止维护了。有些新的分类单元在里面找不到。

数据库	更新状态	适用场景
SILVA	持续更新，最新为138.1	通用注释，精度高
Greengenes 13_8	停更于2013年	功能预测（PICRUSt）、历史项目对比
GTDB	快速发展	基因组分类，适合深度分析

我的个人经验是：大部分项目用SILVA就够。除非你要和2013年之前的数据做比较，或者需要PICRUSt功能预测，那时候才考虑Greengenes。

3.4 本章小结与流程总览

把上面的步骤串起来，一个典型的16S扩增子上游分析管线就是这样的：

核心原则回顾：

质控是产出的基石：不做FastQC直接往下跑，等于盲人开车
ASV优于OTU：分辨率更高，跨研究可重现
SILVA + DADA2是目前通行的黄金组合

做完整套流程后，你手里会有这三个核心产出：ASV特征表（行是样本，列是ASV）、代表序列FASTA、以及物种注释表格。接下来的α/β多样性分析、差异丰度分析，都离不开它们。

我自己在团队里推动过一个经验：一定要把每一步的参数记录下来。哪一版数据库、Trimmomatic用的什么滑动窗口阈值、DADA2的pool模式是什么——这些信息将来写论文时都是必须的。我见过不少项目因为参数记录缺失，被迫重跑整个流程。别让自己也掉进这个坑里。