16S 扩增子测序分析——从原理到面试¶

5' ─[C1]─[V1]─[C2]─[V2]─[C3]─[V3]─[C4]─[V4]─[C5]─[V5]─[C6]─[V6]─[C7]─[V7]─[C8]─[V8]─[C9]─[V9]─ 3'
     保守  可变  保守  可变  保守  可变  保守  可变  保守  可变  保守  可变  保守  可变  保守  可变  保守  可变

你研究的是什么生物？
├── 细菌/古菌 → 16S rRNA（最成熟、数据库最全）
├── 真菌 → ITS（真菌鉴定的金标准，UNITE 数据库）
├── 真核微生物（原生生物、藻类） → 18S rRNA
└── 什么都想看 → 考虑鸟枪法宏基因组（shotgun metagenomics）

三问法则：

Q1: 你研究什么环境的菌？
    → 人体（肠道/口腔/皮肤）→ V3-V4 或 V4
    → 环境（土壤/水体）→ V4 或 V4-V5
    → 口腔 → V1-V2

Q2: 你需要兼顾古菌吗？
    → 需要 → V4（515F/806R，EMP 方案）
    → 只看细菌 → V3-V4 也行

Q3: 你的测序平台读长是多少？
    → PE150 → 选 V4（~253 bp，拼得上）
    → PE250 → V3-V4（~460 bp，刚好能拼）
    → PE300 → 都行

白话：引物的选择就像选考试范围——V3-V4 考的范围大、题目多、
区分度高，是最通用的方案；V4 单区段短小精悍、兼容性好；
V1-V2 对某些特殊菌群有独到优势。没有"最好的"引物，
只有"最适合你实验目的"的引物。

核心思想：
把序列按相似度 ≥97% 聚成一堆，每堆叫一个 OTU。

经典算法：
1. UCLUST / VSEARCH —— 选一条代表序列（centroid），其他序列和它比对
   如果相似度 ≥ 97%，就归入同一个 OTU
2. 为每个 OTU 选一条代表序列（rep-seq）去做分类注释

步骤：
原始序列 → 质控 → 去嵌合体 → 按 97% 相似度聚类 → OTU 表 → 分类注释

问题：
- 97% 阈值是人为规定的，没有生物学依据（不同菌种之间 16S 差异可能 <3%）
- 不同实验产生的 OTU 编号不通用（OTU_001 在 A 实验和 B 实验是不同的菌）
- 测序错误会产生大量假 OTU（一个碱基错误就可能产生新 OTU）

白话：OTU 聚类就像"长得差不多的人编成一组"——你和你表弟可能被编成一组，
但你们其实是不同的人。而且每个村的编组方式不一样，跨村比较没意义。

核心思想：
不是"按相似度聚类"，而是用统计模型学习测序错误的模式，
然后把错误纠正回去，推断出真实的生物序列（ASV）。

四步流程：

1. 学习错误模型（learnErrors）
   - 统计每种碱基替换的频率（A→G 多常见？T→C 多常见？）
   - 建立一个"错误率 vs 质量值"的关系模型
   白话：先学习"这台测序仪通常犯什么样的错误"

2. 去噪（dada）
   - 对每条序列，计算它是"真实生物序列"还是"测序错误产生的假序列"的概率
   - 如果一条序列和某条丰度更高的序列只差一两个碱基，
     且差异碱基的质量值低 → 大概率是测序错误 → 纠正回去
   白话：对每条序列做"有罪推定"——你和高丰度序列只差一个碱基？
         那这个碱基质量还很低？八成你就是人家的测序错误版本。

3. 合并双端（mergePairs）
   - PE 测序的正向和反向读段有重叠区
   - 重叠区的序列必须一致才合并
   白话：两张从不同方向拍的半张照片，重叠部分对得上才拼成完整照片

4. 去嵌合体（removeBimeraDenovo）
   - 检测并移除 PCR 过程中产生的嵌合体序列
   白话：PCR 有时候会把两条不同的 DNA "缝"在一起产生假序列，
         这一步就是把这些"缝合怪"找出来删掉

最终产出：ASV 表（每个 ASV 是一条精确到单碱基的去噪序列）

主流方法：

1. Naive Bayes 分类器（QIIME2 默认 = sklearn 分类器）
   - 原理：用已知物种的 16S 序列训练一个贝叶斯分类器
   - 对每条 ASV 序列，计算它属于每个物种的后验概率
   - 取概率最高的作为分类结果
   - 白话：就像垃圾邮件过滤器——你给它看过一万封垃圾邮件的特征，
           新来一封邮件它就能算出"这封 80% 概率是垃圾邮件"

2. BLAST/VSEARCH 比对
   - 把 ASV 序列直接比对到参考数据库
   - 按最佳匹配结果赋予分类
   - 白话：拿"身份证号"去户籍数据库查，找最像的那个

3. 置信度阈值
   - 通常设置 0.7（70%），低于此阈值的分类结果标为"unassigned"
   - 白话：查出来像但不确定的，宁可说"不知道"也不瞎猜

分类层级（从粗到细）：
界（Kingdom）→ 门（Phylum）→ 纲（Class）→ 目（Order）
→ 科（Family）→ 属（Genus）→ 种（Species）

16S 的通常极限：能可靠注释到属（Genus）级别，种级别常常不确定。

你的情况（T2D 肠道菌群）→ SILVA 138.2 或 Greengenes2

日常分析决策树：
├── 只研究细菌 → SILVA 138.2（最稳妥）或 Greengenes2（更新更先进）
├── 细菌+真核微生物 → SILVA（同时覆盖 16S 和 18S）
├── 只研究真菌 → UNITE
└── 需要和宏基因组数据对比 → Greengenes2（同一棵系统发育树）

白话：数据库就像户籍系统的"花名册"。SILVA 是最全的花名册，
Greengenes2 是最新用"基因组族谱"重新编排的花名册，
RDP 是老牌花名册，UNITE 是只登记真菌的专用花名册。

白话原理（三步走）：

第 1 步：查族谱、找近亲
   输入你的 ASV 序列，在已知基因组的参考树上找到最近的亲戚
   白话：你不知道张三会什么技能，但你知道他表哥张大会修车、
         他堂弟张二会做饭——按照"一家人技能差不多"来猜

第 2 步：预测基因家族
   根据近亲的已知基因组，用系统发育的方法推断你的菌可能有哪些基因
   参考数据库：
   - KEGG KO（KEGG Orthology）—— 基因功能分类
   - EC 酶编号 —— 酶活性分类
   - MetaCyc 代谢通路

第 3 步：推断通路丰度
   用 MinPath 算法，从预测的基因集合中推断哪些代谢通路是完整的
   白话：有了原料清单（基因），推断能做出什么菜（代谢通路）

典型实验设计（发一篇好文章的标配）：

1. 先用 16S 做"大筛"（样本量大、成本低）
   → 找到组间差异最显著的菌群模式（如 T2D 组 Firmicutes/Bacteroidetes 比值升高）

2. 再用宏基因组做"精查"（选重要样本深入分析）
   → 看差异菌到底是功能基因不同（如产丁酸的基因少了）
   → 该 T2D 项目用的就是这一步

3. 结合两者讲故事：
   16S 告诉你"谁变了" → 宏基因组告诉你"它变了之后能干嘛"

白话：16S 是"海选"，宏基因组是"复赛"。
海选便宜、筛得快但只看脸（物种组成）；
复赛贵、慢但看综合实力（功能能力）。

关联点 1：该项目为什么选宏基因组而不是 16S？
"我们需要菌株级别的分辨率来训练机器学习模型，
 16S 只能到属级别，信息量不够，所以选了鸟枪法。"

关联点 2：如果经费有限你会怎么做？
"如果经费有限，我会先用 16S 做大量样本的菌群结构调查，
 用 QIIME2 的 DADA2 流程生成 ASV 表做多样性分析，
 再挑出差异显著的样本子集做宏基因组验证功能。"

关联点 3：你的 Kraken2 结果能和 16S 数据比较吗？
"直接比较要注意方法学差异——16S 有引物偏好性、只能看细菌，
 鸟枪法能看所有微生物。但如果用 Greengenes2 数据库，
 它基于全基因组系统发育树统一了 16S 和鸟枪法的分类体系，
 使得两种方法的结果可以在同一棵树上对照。"

关联点 4：T2D 研究中 16S 发现了什么经典结论？
"大量 16S 研究发现 T2D 患者肠道中产丁酸菌（如 Roseburia、
 Faecalibacterium）丰度降低，这与我们宏基因组中看到的
 丁酸合成通路基因下调是一致的。16S 提供了流行病学证据，
 宏基因组提供了机制证据。"

标准答案：
16S rRNA 基因是几乎所有细菌和古菌都有的一个约 1500bp 的基因。
它由保守区和可变区交替组成——保守区用来设计通用引物，一次 PCR
就能捕获样本中所有细菌的这个基因；可变区（通常选 V3-V4 或 V4）
序列因物种而异，通过高通量测序读取这些序列，再和参考数据库比对，
就能鉴定样本中有哪些菌、各占多少比例。

加分点：
"现在主流方法是用 DADA2 降噪得到 ASV，而不是传统的 97% OTU 聚类。
 ASV 精确到单碱基分辨率，而且不同实验间可以直接比较。"

标准答案：
OTU 是传统方法，把相似度 ≥97% 的序列聚成一组，每组取一条代表序列。
问题是：(1) 97% 阈值是人为设定的；(2) 测序错误会产生大量假 OTU；
(3) 不同实验的 OTU 编号不通用。

ASV 是现在的主流方法（DADA2/Deblur），通过统计模型学习测序错误模式，
把错误纠正回去，得到精确到单碱基的真实生物序列。
ASV 的优势是：精度更高、可重现、跨实验可比、假阳性更低。

白话版：OTU 是"长得像的编一组"，ASV 是"每个人都精确识别"。

标准答案：
16S 测序本身不能直接做功能分析，因为它只测了一个标记基因，
不包含功能基因的信息。但可以用 PICRUSt2 做间接的功能预测——
它根据 ASV 序列在参考系统发育树上的位置，利用近缘物种的已知基因组
来推断样本中可能存在的功能基因和代谢通路。

但必须指出其局限性：PICRUSt2 假设近缘菌功能相似，对于参考数据库中
没有近亲的菌、或经历了水平基因转移的菌，预测准确度会下降。
如果需要准确的功能信息，应该做鸟枪法宏基因组测序。

标准答案：
V3-V4 区的综合表现最均衡：
(1) 变异度够高——能区分大多数细菌到属甚至种级别；
(2) 扩增长度约 460bp，和 Illumina PE250 测序匹配，双端能拼接；
(3) 对人体相关菌群（特别是肠道菌群）的覆盖度最好；
(4) 文献积累最多，便于和已发表数据比较。

V4 单区段（EMP 方案）也是很好的选择，特别是需要兼顾古菌、
或者测序预算只能做 PE150 的时候。
选择哪个 V 区取决于研究对象、测序平台和预算。

标准答案：
两者不是替代关系，而是互补关系：
- 16S 扩增子：成本低、通量高，适合大样本量的菌群结构调查，
  但只能看"有哪些菌"，分辨率到属级别
- 鸟枪法宏基因组：能看物种+功能+甚至菌株变异，但成本高、
  数据量大、分析复杂

实际应用中，经典策略是"16S 先筛、宏基因组精查"：
先用 16S 在大量样本中找到关键的菌群差异模式，
再对重要样本做宏基因组探究功能机制。

结合自身经历：
"该 T2D 项目直接用的鸟枪法，因为我们关注功能基因层面的差异。
 但如果预算有限或样本量很大，我会建议先做 16S 筛选。"

原始数据（FASTQ）
    │
    ▼
① 质控（FastQC/MultiQC 查看质量 → Cutadapt 去引物 → Trimmomatic/fastp 去低质量）
    │
    ▼
② 降噪/聚类（DADA2 生成 ASV 表 或 VSEARCH 聚类生成 OTU 表）
    │
    ▼
③ 分类注释（Naive Bayes 分类器 + SILVA/Greengenes2 数据库 → 物种组成表）
    │
    ▼
④ 多样性分析
   ├── Alpha 多样性：Shannon, Simpson, Chao1, Observed ASVs（单样本丰富度/均匀度）
   └── Beta 多样性：Bray-Curtis, UniFrac（样本间差异）→ PCoA/NMDS 可视化
    │
    ▼
⑤ 差异分析（LEfSe, ANCOM-BC, DESeq2 → 找组间差异菌）
    │
    ▼
⑥ 功能预测（PICRUSt2 → KEGG/MetaCyc 通路）
    │
    ▼
⑦ 可视化与报告（R/ggplot2, phyloseq, MicrobiomeAnalyst）