宏基因组与微生物组

引言

宏基因组（metagenomics）研究的是「一份环境样本里的全部遗传物质」——肠道、土壤、海水、皮肤，任何一个微生物群落。它与常规基因组学的根本区别在于：样本不是一个个体，而是一个混合群体，且组成比例本身就是研究的核心变量。测序结果不是「某个物种的序列」，而是「哪些物种、各占多少」的相对组成。

这带来一个贯穿全流程的统计难题：组成性（compositionality）。测序只给出相对丰度（所有物种加起来是 100%），某个物种「变多」可能只是其他物种变少造成的假象。基于比例数据的统计检验（t 检验、Pearson 相关）在数学上不成立，必须用专门的组成性数据分析方法（CLR 变换、ALDEx2、ANCOM-BC）。这是微生物组分析中最容易出错、也最容易被审稿人抓住的地方。

工程上的难点还有三个：一是数据量巨大且混杂，一份粪便样本测序后包含几百到几千个物种、数万个基因，且夹杂大量宿主 DNA；二是参考数据库不完整，环境样本里常有大量「未培养微生物」无法分类；三是低生物量样本（如血液、肺部灌洗液、洁净环境拭子）的信噪比极低，试剂污染可能比真实信号还强。

本文按「问题定义 → 16S 与 WMS 取舍 → 实验设计 → 质控 → 分类 → 丰度 → 功能 → 多样性 → 差异分析」的顺序展开。命令基于 QIIME2 2024.10、Kraken2 2.1.3、MetaPhlAn 4.1 与 HUMAnN 3.9。读完你应该能独立设计并完成一个微生物组项目，并避开最常见的统计陷阱。

目录

1. 宏基因组的问题定义与数据特点
2. 16S 扩增子与全宏基因组的取舍
3. 实验设计、采样与低生物量污染控制
4. 质控与宿主序列去除
5. 物种分类：Kraken2 与 Bracken
6. 丰度估计：MetaPhlAn 与标记基因法
7. 功能注释：HUMAnN 与通路推断
8. 多样性分析：alpha、beta 与统计检验
9. 差异丰度分析与组成性陷阱

1. 宏基因组的问题定义与数据特点

一个微生物组项目的典型问题形式是：「处理组与对照组的肠道菌群组成有差异吗？差异体现在哪些物种或功能通路上？」注意问题本身是关于群落结构的，而非关于某个个体的基因组。这决定了分析的两个基本特征：以「物种/功能单元」为行、「样本」为列的丰度表（OTU/ASV 表或物种表）是全流程的中心产物；所有统计都在这个丰度表上进行。

数据特征决定了后续所有方法选择：

• 相对丰度而非绝对丰度：测序深度是人为设定的，无法直接推断「某物种的绝对数量」。除非额外做 qPCR 或加入已知量的内标（spike-in），否则只能谈比例。
• 高度偏态：少数物种占据 90% 以上的丰度（如肠道里的拟杆菌与厚壁菌），长尾里是大量低丰度物种。这个分布形态会让很多统计方法失效。
• 稀疏：一份样本里测到 500 个物种，跨样本合并后可能是 5000 个，多数物种在多数样本里是 0。稀疏让相关性分析极不稳定。
• 参考不完整：环境样本中可分类到物种的比例常只有 50%-70%，其余是「未分类细菌」。

理解这四点，就能理解为什么微生物组分析会衍生出一整套自己的统计方法学——它不能直接套用转录组或基因组分析的现成工具。

OTU 与 ASV：丰度表怎么来

16S 分析的起点是把读段聚成「分类单元」。历史上有两种做法：

• OTU（Operational Taxonomic Unit）：按 97% 相似度聚类，同一簇的序列算作一个单元。优点是容忍测序错误，缺点是「97% 阈值」是人为设定，不同研究、不同聚类参数给出的 OTU 边界不一致，跨研究不可比。
• ASV（Amplicon Sequence Variant）：用 DADA2 / Deblur 做错误模型校正后，每个独特的真实序列就是一个单元，精确到单碱基分辨率。优点是结果可跨研究复现（同一个序列在任何数据集里都是同一个 ASV），当前已是主流。

# QIIME2 + DADA2 生成 ASV 表（16S V3-V4 示例）
qiime dada2 denoise-paired \
  --i-demultiplexed-seqs demux.qza \
  --p-trim-left-f 17 --p-trim-left-r 21 \
  --p-trunc-len-f 250 --p-trunc-len-r 230 \
  --p-n-threads 32 --o-table table.qza \
  --o-representative-sequences rep-seqs.qza --o-denoising-stats stats.qza

--p-trunc-len 的选择有讲究：截断过短会丢信息，过长会把质量下降的尾部保留进来。判据是看 demux.qzv 里的质量分布图，在质量中位数掉到 30 左右的位置截断。这一步没做好，DADA2 的错误模型会把大量低质量读段误判成「新 ASV」，导致 ASV 数量虚高。

2. 16S 扩增子与全宏基因组的取舍

两种主流技术路线，选择决定了后续所有分析的可能性边界：

16S rRNA 扩增子测序：用保守区引物（如 V3-V4 区的 341F/805R）扩增 16S 基因的可变区，只测这一个标记基因。成本极低（单样本几十元到几百元），通量高，适合大样本队列。局限是只能给出「属」级别的分类精度（V3-V4 通常到属，部分到种）、无法做功能分析、引物偏好会系统性扭曲某些类群的丰度。

全宏基因组测序（WMS / shotgun）：把样本里所有 DNA 打碎后全部测序。能给出种/株级分类、能直接分析功能基因与代谢通路、能组装出 MAG（宏基因组组装基因组）。代价是成本高 10-50 倍、数据量大、宿主 DNA 污染会浪费大量测序量。

选型口诀：「只想知道谁在那儿」用 16S，「想知道它们在干什么」用 WMS。预算允许时，先做一批 16S 摸清群落结构，再挑代表性样本做 WMS 深入机制，是性价比最高的组合策略。

引物与测序区域的选择

16S 有 9 个可变区（V1-V9），常用组合是 V3-V4、V4、V4-V5。选择不是随意的：

• V4 单区（515F/806R）：最通用，Earth Microbiome Project 的标准，跨研究可比性最好；
• V3-V4：分辨率略高，国内实验室常用，但与其他研究的 V4 数据不可直接合并；
• V1-V3：对某些类群（如双歧杆菌）分辨更好，但读长要求高。

关键约束是扩增子长度必须能被 2×250 或 2×300 的双端测序覆盖并重叠，否则无法拼接。更要紧的是「跨研究合并」：不同引物扩增的同一物种，其 ASV 序列完全不同，绝不能把不同引物区域的数据合到一个丰度表里。这是 meta 分析中最常见的致命错误之一。

WMS 没有引物偏好问题，但文库制备的片段化与 PCR 循环数仍会引入偏差（GC 偏好、低丰度物种丢失）。建库时用低循环数（< 12）并做技术重复，能显著改善可复现性。

3. 实验设计、采样与低生物量污染控制

微生物组研究里，实验设计的重要性远高于分析——采样一旦出错，再好的分析也救不回来。

样本量与统计功效：群落数据的高变异性决定了需要的样本量比常规实验更大。经验值：检测中等效应量（如 β 多样性差异）通常需要每组 20-30 个样本；队列研究建议 100 例以上。样本量不足是微生物组论文被质疑最多的问题。

采样与保存：粪便样本的室温放置时间会显著改变群落组成（兼性厌氧菌增殖）。规范做法是采样后立即 −80°C 冻存，或用保存液（如 DNA/RNA Shield）室温稳定。采样部位（粪便不同区段、口腔不同部位）也会带来差异，必须全程一致。

低生物量样本的污染控制是这类研究的头号陷阱。当样本里的微生物 DNA 总量极低（如血液、脑脊液、支气管灌洗液、洁净环境拭子）时，试剂与耗材自带的背景 DNA（「kitome」）会成为主要信号：

低生物量样本的污染控制要求
  1. 每批次必须做阴性对照：采样空白（未使用的拭子）+ 提取空白（无样本的提取管）+ 无模板 PCR 对照
  2. 对照与真实样本必须同批处理、同批测序（不能只测一次对照了事）
  3. 数据分析时先看对照里有什么：对照中占主导的属（常见如 Ralstonia、
     Bradyrhizobium、Cutibacterium、Halomonas）必须从所有样本中扣除或标记
  4. 若某物种在对照中的丰度与真实样本相当，该物种的结论一律不可信

工具方面，decontam（R 包）用「对照样本」与「DNA 浓度」两种信息做统计去污染，是当前的标准做法。它的原理是利用「污染物在低浓度样本中占比更高」这一模式，通过频率法或浓度回归识别污染序列。但工具不能替代对照——没有阴性对照的数据，任何去污染都是猜测。

# decontam 的两种模式（R）
library(decontam); library(phyloseq)

# 频率法：靠「是否为阴性对照」这个标签
contamdf.freq <- isContaminant(ps, method = "frequency", conc = NULL,
                               neg = "is_neg_control", threshold = 0.1)

# 浓度法：靠「样本 DNA 浓度」——污染物丰度与浓度负相关
contamdf.conc <- isContaminant(ps, method = "prevalence",
                               conc = "DNA_conc", threshold = 0.5)

ps_clean <- prune_taxa(!contamdf.freq$contaminant, ps)

两种方法各有前提：频率法要求对照与真实样本同批处理，浓度法要求每个样本都测了 DNA 浓度（Qubit 或 qPCR）。判据的直觉是——污染物的丰度不应与真实生物量正相关，而真实物种的丰度通常随生物量上升。若某物种在浓度最高的样本里也最丰富，它更可能是真实的。

元数据的设计

元数据表（metadata）的质量决定分析的上限。至少要记录：样本 ID（与测序数据严格对应）、分组变量、批次（提取批次、测序 lane、建库日期）、技术参数（DNA 浓度、提取试剂盒批次）、以及所有可能的混杂变量（年龄、性别、BMI、用药、饮食）。批次信息尤其关键——如果处理组恰好都来自同一个提取批次，那么「处理效应」与「批次效应」在数学上无法区分，这个项目从设计阶段就失败了。

4. 质控与宿主序列去除

WMS 数据的质控流程与常规测序相似（接头去除、低质量剪裁、长度过滤），可参考 FASTQ 质控与预处理 。微生物组特有的两个步骤：

宿主序列去除。宿主 DNA 在某些样本里占比极高（组织活检样本中宿主可占 99%），不去除会浪费大量计算并污染下游分类：

# 用宿主参考基因组比对并剔除（Bowtie2 快速模式，人宿主示例）
bowtie2 -x host_index -1 R1.fq.gz -2 R2.fq.gz -p 32 --very-sensitive \
  --un-conc-gz host_removed_%.fq.gz -S /dev/null

# 或用 KneadData 把质控与去宿主串成一条流水线
kneaddata --input1 R1.fq.gz --input2 R2.fq.gz \
  --reference-db host_bowtie2_index --output kneaddata_out \
  --trimmomatic /path/to/trimmomatic --run-trimmomatic \
  --processes 32

注意 --un-conc-gz 输出的是「未能比对到宿主」的读段，这才是你要的；很多人误取了比对上的文件。KneadData 把 Trimmomatic 与 Bowtie2 串起来，是常用的一站式方案。

去重与低复杂度过滤。宏基因组数据里常有完全相同的读段（PCR 重复或高丰度物种），工具如 fastp 的去重功能可用；低复杂度序列（如 poly-A、微卫星）会干扰组装与分类，可用 bbduk 的 entropy 过滤剔除。

去污染（decontamination）与去宿主是两件事，不要混淆：前者针对试剂污染（用阴性对照），后者针对样本自身的宿主 DNA（用宿主参考比对）。两者都要做，顺序上一般先去宿主（减少数据量），再去污染（在丰度表层面）。

5. 物种分类：Kraken2 与 Bracken

Kraken2 是当前最主流的 k-mer 分类器：把参考数据库里所有序列切成 k-mer（默认 k=35）建哈希表，然后把每条读段的所有 k-mer 查表，取「最长公共祖先」（LCA）作为分类结果。它的优势是极快（百万读段/分钟级）且对未分类序列友好（能给出「只在属级确定」的结果）。

# 建库（标准库约 50-100 GB，含细菌/古菌/病毒/人类）
kraken2-build --standard --db $DBNAME --threads 32

# 分类
kraken2 --db $DBNAME --threads 32 --paired R1.fq.gz R2.fq.gz \
  --report sample.report --output sample.kraken2 --confidence 0.1

# Bracken：把属级结果重新分配到种级，并估计丰度
bracken -d $DBNAME -i sample.report -o sample.bracken -r 150 -l S

--confidence 参数是最关键的调节旋钮。它要求「一条读段上至少有一定比例的 k-mer 支持同一分类」，默认 0，调到 0.1-0.5 能显著降低假阳性（代价是未分类比例上升）。低生物量或高宿主污染样本建议用 0.2-0.5。

Kraken2 的原始输出是「读段计数」，不是「物种丰度」——一个基因组大的物种天然会被更多读段命中。Bracken（Bayesian Reestimation of Abundance with KrakEN）做的事就是「按基因组大小与读长做重估计」，把读段计数转成更接近真实的相对丰度。不做 Bracken 直接拿 Kraken2 的计数当丰度，是极常见的错误，会让基因组大的物种（如大肠杆菌）系统性高估。

一个补充认知：Kraken2 用的是 LCA 策略，因此「属级确定、种级不定」的读段会被归到属。这实际上是有价值的保守信息，但在做种级差异分析时会被算作「属」的丰度，需要意识到这一点。

6. 丰度估计：MetaPhlAn 与标记基因法

MetaPhlAn 走的是完全不同的路线：它不用全基因组比对，而是用物种特异的标记基因（marker genes）（约 100 万个，覆盖 2.6 万个物种）来估计相对丰度。这带来几个优势：数据库小（几 GB 而非上百 GB）、速度快、丰度估计更接近真实（因为标记基因在物种间是单拷贝的，天然做了「基因组大小归一化」）。

# MetaPhlAn 4：直接输出物种级相对丰度表
metaphlan sample_R1.fq.gz,sample_R2.fq.gz \
  --input_type fastq --nproc 32 \
  --bowtie2out sample.bowtie2.bz2 \
  -o sample_profile.txt

# 合并多样本（输出可读入下游统计工具的表）
merge_metaphlan_tables.py *_profile.txt > merged_abundance.txt

两种方法的定位差异：

实践建议：大样本队列用 MetaPhlAn（快、可比性好），探索性/环境样本用 Kraken2（能发现未知）。两者结果常不一致，这是方法学差异而非 bug——严谨研究会用两者交叉验证，报告一致的部分。若要做菌株级分辨（strain-level），两者都不够，需要专门工具（如 StrainPhlAn、inStrain）。

7. 功能注释：HUMAnN 与通路推断

物种组成只是「谁在」，功能才是「在干什么」。功能分析有两条路：

直接测定（WMS）：HUMAnN 3 把读段比对到 UniRef90 蛋白库，得到基因家族丰度，再映射到 MetaCyc 通路并做「通路覆盖度」校正：

# HUMAnN 3：输入是去宿主后的读段与 MetaPhlAn 的物种剖面
humann --input sample.fastq.gz \
  --output humann_out --threads 32 \
  --nucleotide-database chocophlan --protein-database uniref90 \
  --taxonomic-profile sample_profile.txt

# 标准化：把基因家族丰度转成「每百万」（CPM）后再比较
humann_renorm_table --input sample_genefamilies.tsv \
  --output sample_genefamilies_cpm.tsv --units cpm

HUMAnN 的关键设计是分层（stratified）输出：它不只给「某通路的丰度」，还给出「该通路由哪些物种贡献」。这对机制解释极有价值——同样是「丁酸合成通路升高」，由谁贡献意味着完全不同的生物学含义。

预测（16S）：PICRUSt2 用 16S 序列推断功能潜力。它的前提是「16S 系统发育与基因组功能相关」，精度有限，且对参考中没有的类群无效。当前共识是：PICRUSt2 只能用于假设生成，不能作为功能证据；有 WMS 数据时不要用预测。

一个常见误区是「把基因家族丰度直接相加得到通路丰度」。通路的活性取决于「限速步骤」，简单的丰度求和会掩盖关键酶的缺失。HUMAnN 的「通路覆盖度」（pathway coverage）就是对这一点的修正——只有当通路上的大部分反应都被检出时，才算这条通路「完整」。

通路数据库的选择同样影响结论：MetaCyc（HUMAnN 默认）偏重代谢反应网络的精细刻画，KEGG 覆盖更广、更常用于跨研究比较，GO 偏功能术语而非通路。同一个基因集映射到不同数据库会得到不同粒度的结果，报告时必须写明用的是哪个数据库与哪个版本。跨研究比较功能丰度时，若数据库版本不同，差异可能纯属注释差异。

还有一个实践细节：功能丰度也需要做组成性处理。HUMAnN 输出的 CPM 仍是相对值，做差异分析时同样不能用 t 检验直接比较，要沿用第 9 节的组成性方法。这一点常被忽略——很多人对物种表做了 ANCOM-BC，却对功能表直接跑 t 检验，导致结论不一致。

8. 多样性分析：alpha、beta 与统计检验

多样性分析是微生物组最经典的产出，也最容易误用：

Alpha 多样性描述「单个样本内部的丰富度与均匀度」：

关键前提：alpha 多样性必须先做「测序深度抽平」（rarefaction）或使用深度校正方法。观测物种数随测序深度单调上升，不抽平就直接比较「样本 A 比 B 更丰富」是没有意义的——差异可能纯粹来自测序量。抽平到最小深度是传统做法，但有争议（丢弃数据），替代方案是使用对深度不敏感的指标或深度作为协变量。

Beta 多样性描述「样本之间的差异」，需要先算距离矩阵再降维：

# QIIME2 中的 beta 多样性流程
qiime diversity beta --i-table table.qza --p-metric braycurtis \
  --o-distance-matrix bray.qza
qiime diversity pcoa --i-distance-matrix bray.qza --o-pcoa pcoa.qza
qiime diversity adonis --i-distance-matrix bray.qza \
  --m-metadata-file metadata.tsv --p-formula "Treatment + Batch" \
  --o-visualization adonis.qzv

距离指标的选择影响结论：Bray-Curtis（考虑丰度，对优势物种敏感）、UniFrac（考虑系统发育，分加权/非加权）、Jaccard（只看有无）。UniFrac 需要系统发育树，对 16S 分析是常规选择；WMS 常直接用 Bray-Curtis。

统计检验用 PERMANOVA（adonis2） 检验「分组是否解释了群落差异」，用 PERMDISP 检验「组内离散度是否不同」。这两个必须一起报告——PERMANOVA 显著可能只是因为某一组的样本更分散（离散度差异），而非真正的组成差异。这是 beta 多样性分析最常被忽略的一步。

9. 差异丰度分析与组成性陷阱

差异丰度分析（找出「处理组显著富集的物种」）是微生物组论文的核心结果，也是陷阱最密集的地方。

陷阱一：组成性。丰度表是相对数据，所有物种的比例之和恒为 1，因此物种之间天然存在负相关。用 t 检验或 DESeq2（为计数数据设计）直接分析，会得到大量假阳性。正确方法是专门为组成性数据设计的方法：

陷阱二：稀疏与零膨胀。大量物种在多数样本里是 0，导致检验功效极低。过滤策略（保留在至少 10%-20% 样本中检出的物种）是常规做法，但过滤阈值会显著影响结果，必须在方法里如实报告。

陷阱三：多重检验与假发现。一次分析要检验几千个物种，必须做 FDR 校正（Benjamini-Hochberg）。报告时给出校正后的 q 值而非原始 p 值。

陷阱四：混杂变量。批次、测序深度、DNA 提取方法、饮食、用药（尤其抗生素）都会强烈影响群落。分析时必须把这些作为协变量纳入模型，或在设计阶段就做好配对/分层。一个只报告「处理 vs 对照显著」而不控制混杂的差异丰度结果，几乎不可信。

一个务实的工作流是：先用 beta 多样性确认「整体上有差异」，再用差异丰度方法找「哪些物种驱动了这个差异」，最后用功能分析解释「这些物种的变化意味着什么」。三步递进，且每一步都要报告效应量而非只报 p 值。

权衡取舍

常见坑清单

1. 不做 Bracken 直接用 Kraken2 计数：基因组大的物种被系统性高估；必须做重估计。
2. 用 DESeq2 分析相对丰度表：组成性导致大量假阳性；改用 ANCOM-BC 或 ALDEx2。
3. 不抽平就比 alpha 多样性：测序深度差异被误读为丰富度差异；先抽平或用深度校正模型。
4. 低生物量样本无阴性对照：试剂污染被当成真实信号；每批必须带采样、提取、PCR 三种对照。
5. PERMANOVA 显著就下结论：可能只是组内离散度不同；必须补 PERMDISP。
6. 去宿主取错输出文件：Bowtie2 的 --un-conc-gz 才是非宿主读段；取反了等于丢掉全部微生物数据。
7. 忽略引物偏好（16S）：某些门类（如 TM7、疣微菌）扩增效率低被系统性低估；跨研究比较需注意。
8. 样本量与功效不足：群落数据方差大，每组 5 例几乎检不出效应；按功效分析定样本量。
9. 忽略抗生素/饮食等混杂：真实效应被混杂掩盖或放大；纳入协变量或设计阶段配对。
10. 未分类序列直接丢弃：环境样本中 30%-50% 未分类，丢弃会系统性低估多样性；如实报告未分类比例。

小结

微生物组分析的核心特殊性在于「数据是相对的、群落是整体的」。组成性、稀疏性、参考不完整这三个特征，决定了它不能照搬转录组或基因组分析的现成方法学。从实验设计（对照、样本量、采样一致性）到统计分析（CLR、PERMANOVA、ANCOM-BC），每一步都需要针对群落数据做专门处理。

工程上最该建立的三个习惯：一是「先看对照」，任何低生物量研究在分析真实样本前，必须先搞清对照里有什么；二是「先看整体再看局部」，用 beta 多样性确认有差异，再找具体物种，避免在噪声里挑 p 值；三是「报告效应量与置信区间」，微生物组数据方差大，只报 p 值极易产生不可重复的结果。

下一步可以看 单细胞测序数据分析 对比另一种「高维稀疏数据」的分析思路，或 RNA-seq 转录组分析 了解差异表达方法的设计逻辑——ANCOM-BC 与 DESeq2 的差异正来自「相对 vs 绝对」这个根本区别。如果你的差异丰度结果「每个物种都显著」，先检查是否忘了组成性校正。

常见问题

Q：16S 和 WMS 的结果能直接比较吗？
不能。两者的分类精度、引物偏好、丰度估计方式都不同，同一份样本给出不一致的属级丰度是正常的。跨研究比较时至少要统一分类数据库与丰度估计方法。

Q：为什么我的样本里有一半读段无法分类？
这通常不是错误。环境样本中大量微生物没有参考基因组（未培养微生物），标准数据库覆盖不全。可以尝试更全的数据库（如 GTDB）、或用组装 + 基因预测的方式获取「未分类但可用」的序列信息。参见 序列比对与算法 了解比对敏感度对分类率的影响。