微生物多样性分析（Microbial Diversity Analysis）¶
模块编号：02 | 彭文强 面试知识库 | 2026届 生信分析工程师
1. 一句话概述¶
微生物多样性分析通过扩增子测序（16S/ITS）或宏基因组测序数据，评估样本中微生物群落的组成（谁在那里）、结构（各有多少）和差异（组间有什么不同），是微生物组研究的核心分析内容。
2. 核心知识点¶
2.1 扩增子分析基础¶
2.1.1 常用标记基因¶
| 标记基因 | 目标生物 | 常用区段 | 典型引物 | 测序平台 |
|---|---|---|---|---|
| 16S rRNA | 细菌/古菌（Bacteria/Archaea） | V3-V4 区（最常用） | 341F/806R | Illumina PE250/PE300 |
| 16S rRNA | 细菌/古菌 | V4 区 | 515F/806R（Earth Microbiome Project） | Illumina PE150 |
| ITS | 真菌（Fungi） | ITS1 / ITS2 | ITS1F/ITS2R | Illumina PE250/PE300 |
| 18S rRNA | 真核微生物（Eukaryota） | V4 / V9 区 | TAReuk454FWD1/TAReukREV3 | Illumina PE250 |
为什么选 V3-V4 区？ - V3-V4 区长度约 460 bp，包含足够的变异信息用于物种分类 - Illumina PE250/PE300 可以完全覆盖，双端拼接后质量好 - 数据库中 V3-V4 区的参考序列最全，注释准确率高
2.1.2 16S 扩增子 vs 宏基因组¶
| 对比项 | 16S 扩增子（Amplicon） | 宏基因组（Metagenomics） |
|---|---|---|
| 原理 | PCR 扩增特定标记基因 | 直接打断全部 DNA 测序 |
| 分辨率 | 属级（genus）为主，种级（species）有限 | 可达种/株级（species/strain） |
| 功能信息 | 无直接功能信息（可用 PICRUSt2 预测） | 直接获得功能基因信息 |
| 引物偏好 | 有（primer bias），部分物种扩增不到 | 无引物偏好 |
| 成本 | 低（约 50-100 元/样本） | 高（约 500-2000 元/样本） |
| 数据量 | 5-10 万条 reads/样本 | 5-20 Gb/样本（通常 10 Gb） |
| 宿主污染 | 低（引物特异性高） | 需去宿主（如人、小鼠） |
| 分析复杂度 | 低，流程标准化 | 高，涉及组装、binning 等 |
| 适用场景 | 大规模群落调查、临床筛查 | 深入功能研究、新物种发现 |
| 数据库依赖 | SILVA / Greengenes2 / UNITE | NCBI nr / GTDB / UniRef |
| 定量方式 | 相对丰度（relative abundance） | 相对丰度 / 可结合 spike-in 做绝对定量 |
💡 实习经验：在百迈客，16S 项目量远大于宏基因组项目，因为客户预算有限且很多课题只需要群落组成信息。宏基因组通常在需要功能分析或高分辨率物种鉴定时才会选用。
2.2 OTU vs ASV ★¶
| 对比项 | OTU（Operational Taxonomic Unit） | ASV（Amplicon Sequence Variant） |
|---|---|---|
| 定义 | 按序列相似度（通常 97%）聚类的操作分类单元 | 经去噪（denoising）得到的精确序列变体 |
| 核心方法 | 聚类（clustering）：UPARSE、CD-HIT、VSEARCH | 去噪（denoising）：DADA2、Deblur、UNOISE3 |
| 序列分辨率 | 97% 相似度，约属级 | 100% 精确，单碱基差异可区分 |
| 代表序列 | 聚类中心序列（centroid） | 真实的生物序列（biological sequence） |
| 跨研究可比性 | 差，不同研究聚类结果不可直接比较 | 好，同一序列在不同研究中一致 |
| 对噪音处理 | 聚类会将测序错误吸收进 OTU | 统计模型主动识别并去除测序错误 |
| 特征表大小 | 较小（聚类后合并） | 较大（保留更多真实变体） |
| 主流趋势 | 传统方法，逐渐被 ASV 取代 | 当前推荐方法 |
| 代表工具 | USEARCH/VSEARCH + QIIME 1 | DADA2 + QIIME 2 |
DADA2 去噪原理简述： 1. 学习测序错误模型（error model）：基于每个 run 的数据，拟合错误率与质量值的关系 2. 样本推断（sample inference）：利用错误模型，区分真实生物序列和测序错误 3. 合并双端序列（merge paired reads）：要求至少 12 bp 重叠，无错配 4. 去除嵌合体（remove chimeras）：用 consensus 方法识别嵌合序列
📌 面试要记住：OTU 是 97% 聚类，ASV 是 100% 去噪；现在主流推荐 ASV（DADA2），因为分辨率更高、可跨研究比较。
2.3 Alpha 多样性（样本内多样性）★¶
Alpha 多样性（alpha diversity）衡量单个样本内微生物群落的多样性，回答"这个样本里的微生物种类有多丰富、分布有多均匀"。
2.3.1 常用指标详解¶
| 指标 | 英文名 | 衡量维度 | 值域范围 | 直觉理解 |
|---|---|---|---|---|
| Observed species | Observed features | 丰富度（richness） | >= 0，整数 | 直接数有多少种 |
| Chao1 | Chao1 index | 丰富度估计 | >= Observed | 加上对稀有种的补偿估计 |
| ACE | ACE index | 丰富度估计 | >= Observed | 类似 Chao1，用丰度分布估计 |
| Shannon | Shannon-Wiener index | 丰富度 + 均匀度 | 0 ~ ln(S) | 越大 = 种类多且分布均匀 |
| Simpson | Simpson's index (D) | 优势度 | 0 ~ 1 | 越大 = 优势种越突出（低多样性） |
| Inverse Simpson | 1/D | 多样性 | 1 ~ S | 越大 = 多样性越高 |
| Pielou's evenness | Pielou's J | 均匀度 | 0 ~ 1 | 越接近 1 = 越均匀 |
| Good's coverage | Good's coverage | 测序覆盖度 | 0 ~ 1 | 越接近 1 = 覆盖越充分 |
| Faith's PD | Faith's phylogenetic diversity | 系统发育多样性 | >= 0 | 进化枝越多 = 多样性越高 |
2.3.2 各指标详细说明¶
Observed species（观察到的物种数） - 最直观的丰富度指标，直接统计样本中出现的 ASV/OTU 数量 - 缺点：严重依赖测序深度，测序浅的样本值偏低 - 需要配合稀释曲线（rarefaction curve）评估
Chao1 指数 - 公式直觉：Chao1 = 已观察到的种数 + 对未观察到的稀有种的估计 - 核心思想：如果样本中有很多"只出现 1 次"（singleton）或"只出现 2 次"（doubleton）的物种，说明还有很多稀有种没被测到 - Chao1 = S_obs + F1^2 / (2 * F2)，其中 F1 = singleton 数，F2 = doubleton 数 - 值增大：说明样本有更多物种（包括估计的稀有种） - 值减小/接近 Observed：说明大部分物种已被充分采样
ACE 指数（Abundance-based Coverage Estimator） - 与 Chao1 类似，但利用丰度 <= 10 的物种分布来估计 - 对稀有物种的估计更稳健 - 通常与 Chao1 趋势一致
Shannon 指数（Shannon-Wiener Index） ★ - H = -sum(pi * ln(pi))，其中 pi 是第 i 个物种的相对丰度 - 同时考虑丰富度（种类数）和均匀度（分布是否均匀） - 值范围：0（只有 1 个物种）到 ln(S)（所有物种等比例） - 典型值：肠道微生物 2-5，土壤微生物 5-10 - H 增大 = 物种多且分布均匀 = 高多样性 - H 减小 = 物种少或某些种极度优势 = 低多样性
Simpson 指数 ★ - D = sum(pi^2)，衡量"随机取两个个体属于同一物种"的概率 - ⚠️ 注意方向：D 越大 = 优势种越突出 = 多样性越低 - 这与 Shannon 指数方向相反！Shannon 越大多样性越高，Simpson 的 D 越大多样性越低 - 因此常用 1-D（Gini-Simpson）或 1/D（Inverse Simpson）来与 Shannon 方向一致 - 1-D 值范围：0（完全单一）到 1-1/S（非常多样）
⚠️ 易混淆：Shannon 越大 = 多样性越高；Simpson 的 D 越大 = 多样性越低。面试时一定注意这个方向差异！软件输出的"Simpson"可能是 D、1-D 或 1/D，要看清楚文档。
Good's coverage（Good 覆盖度） - 估计测序对群落的覆盖程度 - 公式：C = 1 - (F1 / N)，F1 = singleton 数，N = 总序列数 - 值 > 0.97 通常认为测序深度足够 - 值 < 0.90 建议增加测序深度
Faith's PD（Faith 系统发育多样性） - 计算样本中所有物种在系统发育树上覆盖的总枝长 - 既考虑物种数，又考虑物种间的进化距离 - 优点：进化距离越远的物种贡献越大 - 需要构建系统发育树（QIIME 2 中用 fragment-insertion 或 align-to-tree-mafft-fasttree）
2.3.3 Alpha 多样性统计检验¶
| 比较场景 | 统计方法 | R/QIIME2 函数 |
|---|---|---|
| 两组比较 | Wilcoxon 秩和检验（非参数） | wilcox.test() / qiime diversity alpha-group-significance |
| 多组比较 | Kruskal-Wallis 检验 | kruskal.test() / qiime diversity alpha-group-significance |
| 多组两两比较 | Dunn 检验（事后检验） | dunn.test::dunn.test() |
| 配对样本 | Wilcoxon 符号秩检验 | wilcox.test(paired=TRUE) |
2.4 Beta 多样性（样本间差异）★¶
Beta 多样性（beta diversity）衡量不同样本之间微生物群落组成的差异程度，回答"这些样本的微生物群落有多不一样"。
2.4.1 常用距离/相异度¶
| 距离指标 | 英文 | 考虑因素 | 值域 | 适用场景 |
|---|---|---|---|---|
| Bray-Curtis | Bray-Curtis dissimilarity | 物种丰度差异 | 0~1 | 最常用，关注丰度变化 |
| Jaccard | Jaccard distance | 物种有无（presence/absence） | 0~1 | 只关心组成，不关心丰度 |
| Weighted UniFrac | Weighted UniFrac | 丰度 + 进化关系 | 0~1 | 丰度差异 + 系统发育 |
| Unweighted UniFrac | Unweighted UniFrac | 有无 + 进化关系 | 0~1 | 稀有种差异 + 系统发育 |
| Euclidean | Euclidean distance | 丰度绝对差 | >= 0 | 较少用于微生物组 |
| Aitchison | Aitchison distance | 组分数据（CLR 变换后） | >= 0 | 组分数据分析推荐 |
2.4.2 核心距离详解¶
Bray-Curtis 距离 ★ - 公式直觉：两个样本共有物种的丰度差异 / 总丰度 - BC = 1 - 2 * sum(min(xi, yi)) / (sum(xi) + sum(yi)) - 0 = 两样本完全相同 - 1 = 两样本没有任何共有物种 - 最常用的 beta 多样性距离，适合大多数微生物组分析 - 对丰度敏感，优势种影响大
UniFrac 距离 ★ - 独特之处：考虑物种间的进化关系（phylogenetic relationship） - 需要系统发育树作为输入
| 类型 | Weighted UniFrac | Unweighted UniFrac |
|---|---|---|
| 考虑因素 | 物种进化距离 + 丰度 | 仅物种进化距离（有/无） |
| 敏感性 | 对优势种（高丰度物种）变化敏感 | 对稀有种（低丰度/独有物种）变化敏感 |
| 适用场景 | 群落结构整体差异 | 物种组成差异（不论丰度） |
| 计算 | 进化枝长 * 丰度权重 | 仅计算独有枝长比例 |
📌 选择建议： - 关注整体群落变化（优势种主导）→ Weighted UniFrac 或 Bray-Curtis - 关注稀有种或物种组成差异 → Unweighted UniFrac 或 Jaccard - 没有系统发育树 → Bray-Curtis（不需要树） - 有树且想考虑进化 → UniFrac
Jaccard 距离 - 只考虑物种是否存在（有/无），不考虑丰度 - J = 1 - |A ∩ B| / |A ∪ B| - 适合关注"共有了哪些物种"而非"各有多少"
2.4.3 可视化方法¶
| 方法 | 英文 | 原理 | 优缺点 |
|---|---|---|---|
| PCoA | Principal Coordinates Analysis | 基于距离矩阵的降维 | 最常用，保留距离关系 |
| NMDS | Non-metric Multidimensional Scaling | 非参数排序，保留秩次 | stress < 0.2 可接受 |
| PCA | Principal Component Analysis | 基于原始数据的降维 | 适合 CLR 变换后的数据 |
| CCA/RDA | Constrained Analysis | 约束排序，结合环境因子 | 解释环境因子的影响 |
| UPGMA 树 | UPGMA clustering | 层次聚类 | 展示样本间层次关系 |
PCoA vs NMDS： - PCoA：线性方法，保留实际距离，轴有百分比解释度 - NMDS：非线性方法，保留距离的秩次关系，用 stress 值评估拟合 - 一般优先 PCoA，NMDS 在非线性关系明显时更好
2.4.4 统计检验 ★¶
| 方法 | 英文 | 用途 | 工具/函数 |
|---|---|---|---|
| PERMANOVA | Permutational MANOVA | 检验组间群落结构是否有显著差异 | vegan::adonis2() / qiime diversity beta-group-significance |
| ANOSIM | Analysis of Similarities | 类似 PERMANOVA，非参数 | vegan::anosim() |
| PERMDISP | Permutational Dispersion | 检验组间方差是否均一（PERMANOVA 前提） | vegan::betadisper() |
| Mantel test | Mantel test | 检验距离矩阵间的相关性 | vegan::mantel() |
PERMANOVA（adonis） ★ - 最常用的 beta 多样性统计检验 - 原理：将总变异分解为组间和组内，用置换检验（permutation test）评估显著性 - 前提假设：各组方差齐性（需先做 PERMDISP 检验） - p < 0.05 说明组间差异显著 - R^2 值反映分组因素解释了多少变异
2.5 差异分析 ★¶
差异分析用于找出哪些特定微生物在不同组间有显著差异，即寻找"标志物"（biomarker）。
2.5.1 常用方法对比¶
| 方法 | 全称 | 原理 | 输入 | 优势 | 局限 |
|---|---|---|---|---|---|
| LEfSe | LDA Effect Size | Kruskal-Wallis + LDA | 相对丰度 | 直观、图形好看、广泛引用 | 假阳性偏高、不校正 |
| DESeq2 | Differential Expression (adapted) | 负二项分布模型 | raw counts | 统计严谨、自动 FDR 校正 | 对零膨胀敏感 |
| ANCOM-BC | ANCOM with Bias Correction | 组分数据分析 | raw counts | 处理组分性偏差 | 计算复杂 |
| Wilcoxon | Wilcoxon rank-sum test | 非参数秩检验 | 相对丰度 | 简单稳健 | 需人工多重校正 |
| MaAsLin2 | Microbiome Multivariable Association | 多变量线性模型 | 相对丰度/counts | 可纳入协变量 | 参数选择影响大 |
| ALDEx2 | ANOVA-Like Differential Expression 2 | CLR + 贝叶斯 | raw counts | 处理组分数据好 | 保守（假阴性偏高） |
2.5.2 LEfSe 详解 ★¶
LEfSe（Linear Discriminant Analysis Effect Size） 是微生物组差异分析中最广泛使用的工具之一。
分析步骤： 1. Kruskal-Wallis 检验：在所有组间检验每个分类单元是否有显著差异（p < 0.05） 2. Wilcoxon 秩和检验：对通过步骤 1 的分类单元，做两两组间比较，确认差异一致 3. LDA 打分：用线性判别分析计算效应值（effect size），默认 LDA > 2 为显著
结果解读： - LDA 柱状图：每个柱子代表一个显著差异的分类单元，长度 = 效应值大小，颜色 = 富集在哪个组 - 进化分支图（cladogram）：在分类树上标注差异物种及其归属组 - LDA > 2：弱差异；LDA > 3：中等差异；LDA > 4：强差异
⚠️ 注意：LEfSe 不做多重检验校正，在物种数多时可能有假阳性。高分文章建议配合 DESeq2 或 ANCOM-BC 交叉验证。
2.5.3 多重检验校正¶
| 方法 | 英文 | 严格程度 | 适用场景 |
|---|---|---|---|
| Bonferroni | Bonferroni correction | 最严格 | 少量比较 |
| FDR (BH) | Benjamini-Hochberg | 中等 | 最常用，控制假发现率 |
| FDR (BY) | Benjamini-Yekutieli | 较严格 | 相关检验较多时 |
⚠️ 易错：做多组比较时，必须进行多重检验校正。不做校正的 p 值会大幅膨胀假阳性。用 p.adjust(p, method="BH") 校正。
2.6 分析流程总览¶
原始数据 (FASTQ)
|
v
[1] 质量控制 (QC) ─── fastp / Trimmomatic / Cutadapt（去引物）
|
v
[2] 去噪/聚类 ─── DADA2（推荐）/ VSEARCH / Deblur
|
v
[3] 特征表 (Feature Table) + 代表序列 (Rep Seqs)
|
v
[4] 物种注释 ─── classify-sklearn (QIIME2) → SILVA / Greengenes2 / UNITE
|
v
[5] 构建进化树 ─── MAFFT 比对 → FastTree 建树
|
v
[6] Alpha 多样性 ─── Shannon, Chao1, Observed, Faith's PD...
|
v
[7] Beta 多样性 ─── Bray-Curtis, UniFrac → PCoA/NMDS → PERMANOVA
|
v
[8] 差异分析 ─── LEfSe / DESeq2 / ANCOM-BC
|
v
[9] 功能预测 ─── PICRUSt2（16S→功能预测）/ Tax4Fun2
|
v
[10] 可视化 & 报告 ─── 柱状图、热图、网络图、Venn图
3. 实战命令/代码¶
3.1 QIIME 2 扩增子分析流程¶
# ============================================================
# QIIME 2 16S 扩增子分析流程
# 环境：conda activate qiime2-amplicon-2024.5
# ============================================================
# ---------- 0. 数据导入 ----------
# 将 FASTQ 文件导入 QIIME 2 格式（.qza）
# manifest 文件包含样本名和 FASTQ 路径的对应关系
qiime tools import \
--type 'SampleData[PairedEndSequencesWithQuality]' \
--input-path manifest.tsv \ # 样本清单文件
--output-path paired-end-demux.qza \ # 输出的 qza 文件
--input-format PairedEndFastqManifestPhred33V2 # 输入格式
# 查看导入数据的质量分布（生成可视化文件）
qiime demux summarize \
--i-data paired-end-demux.qza \
--o-visualization paired-end-demux.qzv # 用 view.qiime2.org 查看
# ---------- 1. DADA2 去噪 ----------
# 根据质量分布图确定截断位置
# 📌 trunc-len 参数需要根据质量图调整，确保双端拼接有 >= 12 bp 重叠
qiime dada2 denoise-paired \
--i-demultiplexed-seqs paired-end-demux.qza \
--p-trim-left-f 0 \ # 正向读段左端裁剪长度（去引物后设 0）
--p-trim-left-r 0 \ # 反向读段左端裁剪长度
--p-trunc-len-f 230 \ # 正向读段截断到 230 bp（根据质量图）
--p-trunc-len-r 220 \ # 反向读段截断到 220 bp
--p-n-threads 16 \ # 使用 16 个线程
--o-table table.qza \ # 输出特征表（ASV 丰度矩阵）
--o-representative-sequences rep-seqs.qza \ # 输出代表序列
--o-denoising-stats dada2-stats.qza # 输出去噪统计
# 查看去噪统计（检查每步保留了多少序列）
qiime metadata tabulate \
--m-input-file dada2-stats.qza \
--o-visualization dada2-stats.qzv
# 查看特征表摘要（样本数、ASV 数、测序深度）
qiime feature-table summarize \
--i-table table.qza \
--o-visualization table.qzv \
--m-sample-metadata-file metadata.tsv
# ---------- 2. 物种注释 ----------
# 使用预训练的 Naive Bayes 分类器
# 📌 分类器要与引物区段和 SILVA 版本匹配
qiime feature-classifier classify-sklearn \
--i-classifier silva-138-99-nb-classifier-V3V4.qza \ # SILVA 138 V3-V4 分类器
--i-reads rep-seqs.qza \
--p-n-jobs 16 \
--o-classification taxonomy.qza
# 查看注释结果
qiime metadata tabulate \
--m-input-file taxonomy.qza \
--o-visualization taxonomy.qzv
# 生成物种组成柱状图（不同分类水平）
qiime taxa barplot \
--i-table table.qza \
--i-taxonomy taxonomy.qza \
--m-metadata-file metadata.tsv \
--o-visualization taxa-bar-plots.qzv
# ---------- 3. 构建进化树 ----------
# 用于计算 UniFrac 距离
qiime phylogeny align-to-tree-mafft-fasttree \
--i-sequences rep-seqs.qza \
--o-alignment aligned-rep-seqs.qza \ # MAFFT 多序列比对
--o-masked-alignment masked-aligned.qza \ # 屏蔽高变异位点
--o-tree unrooted-tree.qza \ # 无根树
--o-rooted-tree rooted-tree.qza \ # 有根树（用于 UniFrac）
--p-n-threads 16
# ---------- 4. Alpha + Beta 多样性核心分析 ----------
# 📌 sampling-depth 根据 table.qzv 中最小测序深度确定
# 一般取最小深度附近，排除测序太浅的样本
qiime diversity core-metrics-phylogenetic \
--i-phylogeny rooted-tree.qza \
--i-table table.qza \
--p-sampling-depth 10000 \ # 稀释深度，需根据数据调整
--m-metadata-file metadata.tsv \
--p-n-jobs-or-threads 16 \
--output-dir core-metrics-results # 输出目录，包含多个结果文件
# 输出目录包含：
# - shannon_vector.qza Shannon 指数
# - observed_features_vector.qza Observed species
# - faith_pd_vector.qza Faith's PD
# - evenness_vector.qza Pielou 均匀度
# - bray_curtis_distance_matrix.qza Bray-Curtis 距离矩阵
# - weighted_unifrac_distance_matrix.qza Weighted UniFrac
# - unweighted_unifrac_distance_matrix.qza Unweighted UniFrac
# - bray_curtis_emperor.qzv PCoA 可视化（Emperor）
# ... 等
# ---------- 5. Alpha 多样性组间比较 ----------
# Kruskal-Wallis 检验各组 Shannon 指数差异
qiime diversity alpha-group-significance \
--i-alpha-diversity core-metrics-results/shannon_vector.qza \
--m-metadata-file metadata.tsv \
--o-visualization shannon-group-significance.qzv
# Chao1 指数需要单独计算（core-metrics 不包含 Chao1）
qiime diversity alpha \
--i-table table.qza \
--p-metric chao1 \
--o-alpha-diversity chao1_vector.qza
qiime diversity alpha-group-significance \
--i-alpha-diversity chao1_vector.qza \
--m-metadata-file metadata.tsv \
--o-visualization chao1-group-significance.qzv
# ---------- 6. Beta 多样性统计检验 ----------
# PERMANOVA 检验（基于 Bray-Curtis）
qiime diversity beta-group-significance \
--i-distance-matrix core-metrics-results/bray_curtis_distance_matrix.qza \
--m-metadata-file metadata.tsv \
--m-metadata-column Group \ # 分组列名
--p-method permanova \ # 使用 PERMANOVA
--p-pairwise \ # 两两比较
--o-visualization bray-curtis-permanova.qzv
# PERMANOVA 检验（基于 Weighted UniFrac）
qiime diversity beta-group-significance \
--i-distance-matrix core-metrics-results/weighted_unifrac_distance_matrix.qza \
--m-metadata-file metadata.tsv \
--m-metadata-column Group \
--p-method permanova \
--p-pairwise \
--o-visualization weighted-unifrac-permanova.qzv
# ---------- 7. 稀释曲线 ----------
qiime diversity alpha-rarefaction \
--i-table table.qza \
--i-phylogeny rooted-tree.qza \
--p-max-depth 50000 \ # 最大稀释深度
--p-steps 20 \ # 采样步数
--m-metadata-file metadata.tsv \
--o-visualization alpha-rarefaction.qzv
3.2 R 语言分析代码¶
# ============================================================
# R 语言微生物多样性分析
# 需要安装：vegan, ggplot2, phyloseq, ape
# ============================================================
# ---------- 加载包 ----------
library(vegan) # 群落生态学核心包
library(ggplot2) # 画图
library(ape) # 进化树操作
# ---------- 读取数据 ----------
# 特征表：行=样本，列=物种/ASV
otu_table <- read.csv("otu_table.csv", row.names = 1, check.names = FALSE)
# 分组信息
metadata <- read.csv("metadata.csv", row.names = 1)
# ---------- Alpha 多样性计算 ----------
# Shannon 指数
shannon <- diversity(otu_table, index = "shannon")
# Simpson 指数（注意：vegan 返回的是 1-D，即 Gini-Simpson）
simpson <- diversity(otu_table, index = "simpson")
# Observed species（观测到的物种数）
observed <- specnumber(otu_table)
# Chao1 指数（使用 estimateR 函数，需要整数 count 数据）
# estimateR 返回多个估计值，包含 S.chao1
chao1_result <- t(estimateR(round(otu_table))) # 需要整数输入
chao1 <- chao1_result[, "S.chao1"]
# ACE 指数
ace <- chao1_result[, "S.ACE"]
# Pielou 均匀度 = Shannon / ln(Observed)
pielou <- shannon / log(observed)
# 汇总为数据框
alpha_div <- data.frame(
Sample = rownames(otu_table),
Group = metadata[rownames(otu_table), "Group"], # 从 metadata 获取分组
Shannon = round(shannon, 3),
Simpson = round(simpson, 3),
Observed = observed,
Chao1 = round(chao1, 1),
ACE = round(ace, 1),
Pielou = round(pielou, 3)
)
# ---------- Alpha 多样性箱线图 ----------
# Shannon 指数箱线图
p_shannon <- ggplot(alpha_div, aes(x = Group, y = Shannon, fill = Group)) +
geom_boxplot(outlier.shape = NA, alpha = 0.7) + # 箱线图，隐藏离群点
geom_jitter(width = 0.2, size = 1.5, alpha = 0.6) + # 添加散点
theme_bw() + # 白色背景主题
labs(
title = "Shannon Diversity Index",
x = "Group",
y = "Shannon Index"
) +
theme(legend.position = "none") + # 隐藏图例（颜色已区分）
stat_compare_means(method = "wilcox.test") # 添加 p 值（需 ggpubr 包）
ggsave("alpha_shannon_boxplot.pdf", p_shannon, width = 6, height = 5)
# ---------- Alpha 多样性统计检验 ----------
# 两组比较：Wilcoxon 秩和检验
wilcox_shannon <- wilcox.test(
Shannon ~ Group,
data = alpha_div
)
cat("Shannon Wilcoxon p-value:", wilcox_shannon$p.value, "\n")
# 多组比较：Kruskal-Wallis 检验
kw_shannon <- kruskal.test(Shannon ~ Group, data = alpha_div)
cat("Shannon Kruskal-Wallis p-value:", kw_shannon$p.value, "\n")
# 多组事后两两比较：Dunn 检验（需 dunn.test 包）
# library(dunn.test)
# dunn.test(alpha_div$Shannon, alpha_div$Group, method = "bh")
# ---------- 稀释曲线（Rarefaction Curve） ----------
# 使用 vegan 的 rarecurve 函数
pdf("rarefaction_curve.pdf", width = 8, height = 6)
rarecurve(
round(otu_table), # 需要整数 count
step = 200, # 每隔 200 条序列采样一次
col = as.numeric(factor(metadata$Group)), # 按组着色
lwd = 1.5,
xlab = "Number of Sequences", # X 轴：测序深度
ylab = "Observed Species", # Y 轴：观测物种数
main = "Rarefaction Curves"
)
# 曲线趋于平台 = 测序深度足够
# 曲线仍在上升 = 测序深度不足
dev.off()
# ---------- Beta 多样性分析 ----------
# 1. 计算 Bray-Curtis 距离矩阵
bc_dist <- vegdist(otu_table, method = "bray")
# 2. 计算 Jaccard 距离矩阵（基于有无）
jaccard_dist <- vegdist(otu_table, method = "jaccard", binary = TRUE)
# 3. PCoA 分析（基于 Bray-Curtis）
pcoa_result <- cmdscale(bc_dist, k = 2, eig = TRUE) # 取前 2 个轴
# 计算各轴解释度
eig_percent <- round(pcoa_result$eig / sum(pcoa_result$eig) * 100, 1)
# 整理 PCoA 数据
pcoa_df <- data.frame(
PC1 = pcoa_result$points[, 1],
PC2 = pcoa_result$points[, 2],
Group = metadata[rownames(pcoa_result$points), "Group"]
)
# 4. PCoA 画图
p_pcoa <- ggplot(pcoa_df, aes(x = PC1, y = PC2, color = Group)) +
geom_point(size = 3, alpha = 0.8) + # 样本散点
stat_ellipse(level = 0.95, linewidth = 0.8) + # 95% 置信椭圆
theme_bw() +
labs(
title = "PCoA based on Bray-Curtis Distance",
x = paste0("PC1 (", eig_percent[1], "%)"), # 轴标签含解释度
y = paste0("PC2 (", eig_percent[2], "%)")
) +
theme(
plot.title = element_text(hjust = 0.5), # 标题居中
legend.position = "right"
)
ggsave("pcoa_bray_curtis.pdf", p_pcoa, width = 7, height = 5)
# ---------- NMDS 分析 ----------
nmds_result <- metaMDS(
otu_table,
distance = "bray",
k = 2, # 2 维
trymax = 100 # 最大迭代次数
)
cat("NMDS stress:", nmds_result$stress, "\n")
# stress < 0.05：非常好
# stress < 0.1：好
# stress < 0.2：可接受
# stress > 0.2：不可靠
# ---------- PERMANOVA 统计检验 ----------
# 📌 最重要的 beta 多样性统计检验
permanova_result <- adonis2(
otu_table ~ Group, # 公式：特征表 ~ 分组因子
data = metadata,
method = "bray", # 使用 Bray-Curtis 距离
permutations = 999 # 置换次数
)
print(permanova_result)
# 关注：
# - Pr(>F) 即 p 值，< 0.05 说明组间差异显著
# - R2 值：分组因素解释了多少比例的变异
# - F 值：越大说明组间差异相对组内越大
# PERMDISP：检验方差齐性（PERMANOVA 的前提）
betadisp_result <- betadisper(bc_dist, metadata$Group)
permutest(betadisp_result)
# p > 0.05：方差齐性假设成立，PERMANOVA 结果可靠
# p < 0.05：方差不齐，PERMANOVA 显著可能是方差差异导致
# ---------- ANOSIM 检验 ----------
anosim_result <- anosim(
otu_table,
metadata$Group,
distance = "bray",
permutations = 999
)
print(anosim_result)
# R 值范围 -1 ~ 1
# R > 0.75：组间差异大
# R > 0.5：组间差异中等
# R > 0.25：组间差异小
# R ~ 0：无差异
3.3 LEfSe 分析流程¶
# ============================================================
# LEfSe 分析流程
# 环境：conda activate lefse
# 输入：lefse_input.txt（第一行为分组信息，后续行为物种丰度）
# ============================================================
# 格式说明：lefse_input.txt 格式
# 第 1 行：class（分组信息）
# 第 2 行（可选）：subclass（亚组信息）
# 第 3 行（可选）：subject（个体ID）
# 后续行：每行一个物种的相对丰度
# 步骤 1：格式化输入文件
# 将文本格式转换为 LEfSe 内部格式
lefse_format_input.py \
lefse_input.txt \ # 输入文件
lefse_formatted.in \ # 格式化后的文件
-c 1 \ # 第 1 行是 class（主分组）
-s 2 \ # 第 2 行是 subclass（可选，没有则去掉 -s）
-u 3 \ # 第 3 行是 subject（可选）
-o 1000000 # 归一化到百万
# 步骤 2：运行 LEfSe 分析
# Kruskal-Wallis + Wilcoxon + LDA
lefse_run.py \
lefse_formatted.in \ # 输入
lefse_results.res \ # 输出结果
-l 2.0 \ # LDA 阈值（默认 2.0，可调高到 3.0/4.0 更严格）
-a 0.05 \ # Kruskal-Wallis p 值阈值
-w 0.05 # Wilcoxon p 值阈值
# 步骤 3：绘制 LDA 柱状图
lefse_plot_res.py \
lefse_results.res \
lefse_lda_bar.pdf \ # 输出柱状图
--format pdf \
--dpi 300
# 步骤 4：绘制进化分支图（cladogram）
lefse_plot_cladogram.py \
lefse_results.res \
lefse_cladogram.pdf \
--format pdf \
--dpi 300
# 查看显著差异物种列表
# 筛选 LDA > 2 且有分组归属的结果
grep -v "^$" lefse_results.res | awk -F'\t' '$4 != ""' | sort -k4 -nr
# 输出列含义：
# 第 1 列：物种分类
# 第 2 列：log 最高均值对应的组
# 第 3 列：该物种的 LDA 得分
# 第 4 列：该物种在哪个组中富集（空 = 不显著）
3.4 功能预测（PICRUSt2）¶
# ============================================================
# PICRUSt2：基于 16S 预测功能基因组成
# 环境：conda activate picrust2
# 注意：PICRUSt2 结果是"预测"，不能替代宏基因组
# ============================================================
# 一步法运行（推荐）
picrust2_pipeline.py \
-s rep-seqs.fasta \ # ASV 代表序列（FASTA 格式）
-i feature_table.biom \ # 特征表（BIOM 格式）
-o picrust2_output \ # 输出目录
-p 16 \ # 线程数
--stratified # 输出每个 ASV 对功能的贡献
# 输出文件：
# picrust2_output/
# ├── EC_metagenome_out/ EC 酶注释结果
# ├── KO_metagenome_out/ KEGG KO 注释结果
# ├── pathways_out/ MetaCyc 通路丰度
# └── ...
# 将 KO 丰度表转换为 KEGG 通路丰度
# 可用于后续 STAMP 或 R 画图
4. 面试常问点¶
★ 什么是 alpha 多样性和 beta 多样性？区别是什么？¶
参考答案：
Alpha 多样性衡量的是单个样本内微生物群落的多样性，主要关注物种丰富度（有多少种）和均匀度（分布是否均匀）。常用指标包括 Shannon 指数、Simpson 指数、Chao1 指数等。
Beta 多样性衡量的是不同样本之间微生物群落的差异程度，回答"两个样本的微生物组成有多大不同"。常用 Bray-Curtis 距离、UniFrac 距离等来量化，用 PCoA 或 NMDS 来可视化，用 PERMANOVA 来统计检验。
简单来说：alpha 看"一个样本有多丰富"，beta 看"两个样本有多不同"。
💡 加分回答：在百迈客做 16S 项目时，我们标准分析流程中 alpha 和 beta 多样性都是必做的。alpha 多样性用箱线图展示组间差异，beta 多样性用 PCoA 图展示样本聚类情况。客户最关心的往往是 PCoA 图上两组是否明显分开，以及 PERMANOVA 的 p 值是否显著。
★ Shannon 和 Simpson 指数的区别？¶
参考答案：
两者都是 alpha 多样性指标，但侧重不同：
Shannon 指数： - 同时考虑物种丰富度和均匀度 - 对稀有物种更敏感（因为用 ln 变换，稀有种的贡献不会被压缩太多） - 值越大，多样性越高
Simpson 指数（D）： - 更侧重优势种的影响 - 表示"随机取两个个体属于同一物种的概率" - D 越大，优势种越突出，多样性越低（方向与 Shannon 相反） - 实际使用时常用 1-D 或 1/D，使方向与 Shannon 一致
最重要的区别是对稀有种的敏感度：Shannon 更敏感，Simpson 更关注优势种。如果样本间稀有种差异大，Shannon 更能体现；如果优势种变化是关注点，Simpson 更合适。
★ Bray-Curtis 和 UniFrac 的区别？什么时候用哪个？¶
参考答案：
核心区别：是否考虑进化关系。
- Bray-Curtis：只基于物种丰度计算，不考虑物种间的进化距离。两个物种无论亲缘关系远近，丰度差异的权重相同。
- UniFrac：利用系统发育树，考虑物种间的进化距离。亲缘关系远的物种差异贡献更大。
选择建议：
| 场景 | 推荐距离 | 原因 |
|---|---|---|
| 有系统发育树，想考虑进化关系 | Weighted/Unweighted UniFrac | 进化信息是额外维度 |
| 没有系统发育树（如功能基因分析） | Bray-Curtis | 不需要树 |
| 关注优势种的丰度变化 | Weighted UniFrac / Bray-Curtis | 两者都对丰度敏感 |
| 关注稀有种/独有种的差异 | Unweighted UniFrac | 只看有无，不看丰度 |
| 通用/不确定 | Bray-Curtis | 最常用，发表最多 |
实际项目中通常同时做多个距离，放在补充材料或正文中展示一致性。
💡 加分回答：在百迈客的报告中，我们一般同时给 Bray-Curtis 和 Weighted UniFrac 的 PCoA 图，让客户看两种距离下的聚类是否一致。如果两种距离结果一致，说明结论更可靠。
★ OTU 和 ASV 的区别？现在推荐用哪个？¶
参考答案：
| 对比 | OTU | ASV |
|---|---|---|
| 方法 | 97% 相似度聚类 | DADA2 等去噪，保留精确序列 |
| 分辨率 | 低（约属级） | 高（单碱基可区分） |
| 可重复性 | 低，不同研究聚类结果不同 | 高，同一序列跨研究可比 |
现在推荐 ASV（DADA2），原因： 1. 单碱基分辨率，能区分更多真实变体 2. 可跨研究比较（同一 ASV 序列全球通用） 3. 不会将不同物种错误合并（97% 聚类可能把近缘种合并） 4. QIIME 2 默认使用 DADA2，已成行业标准
但 OTU 在老文献中仍然很常见，读文献时需要理解。
★ 如何判断两组样本的微生物群落有显著差异？¶
参考答案：
判断方法分两个层面：
整体群落层面（beta 多样性）： 1. 计算距离矩阵（如 Bray-Curtis） 2. 做 PCoA 可视化，看两组是否在图上分开 3. 用 PERMANOVA（adonis2 函数）做统计检验 - p < 0.05：组间差异显著 - R^2：分组因素解释了多少变异 4. 先做 PERMDISP 检验方差齐性（PERMANOVA 的前提）
单个物种层面（差异分析）： 1. 用 LEfSe 找哪些物种在不同组间有显著差异 2. 用 DESeq2 或 ANCOM-BC 进行更严谨的差异丰度分析 3. 注意做多重检验校正（FDR）
💡 实习经验回答：在百迈客做项目时，标准分析报告中我们先用 PCoA 展示整体差异趋势，再用 PERMANOVA 给出 p 值，最后用 LEfSe 找具体的差异物种。客户看报告时最在意的三个数字：PERMANOVA 的 p 值、R^2 值，以及 LEfSe 中 LDA 最高的物种是什么。
★ LEfSe 分析的原理和结果怎么解读？¶
参考答案：
原理（三步法）： 1. Kruskal-Wallis 检验：在所有组间筛选差异显著的物种（p < 0.05） 2. Wilcoxon 秩和检验：对通过第一步的物种做两两组间比较，确认差异方向一致 3. LDA 打分：用线性判别分析计算效应值，量化差异的生物学意义
结果解读： - LDA 柱状图：每个柱子 = 一个显著差异物种，长度 = 效应值大小（LDA score），颜色 = 富集在哪个组 - LDA > 2：差异显著（默认阈值） - LDA > 3 或 > 4：差异更强 - 进化分支图（cladogram）：在分类树上用不同颜色标注各组的标志物种，从门到属逐层展示
局限性： - 不做多重检验校正，物种数多时假阳性偏高 - 建议配合 DESeq2 或 ANCOM-BC 交叉验证
★ 稀释曲线的作用是什么？怎么判断测序深度够不够？¶
参考答案：
作用： 评估测序深度是否足够捕获样本中的物种多样性。
原理： 随机从测序数据中抽取不同数量的序列，计算对应的物种数（observed species），绘制"序列数-物种数"曲线。
判断标准： - 曲线趋于平台：说明测序深度足够，继续增加测序量也不会发现新物种 → 测序深度 OK - 曲线仍在上升：说明还有很多物种没被测到，测序深度不足 → 需要增加测序量 - 定量指标：Good's coverage > 0.97 通常认为足够
面试加分： - 稀释曲线不仅评估测序深度，还能初步比较样本间的物种丰富度 - 如果不同组的曲线平台高度不同，说明这些组的物种丰富度确实有差异 - 在做 alpha 多样性分析前，通常需要将所有样本稀释（rarefy）到相同测序深度，以消除测序量差异的影响
💡 实习经验：在百迈客，客户有时会问"测序够不够深"，我们就看稀释曲线。16S 项目通常 5-10 万条 clean reads 就足够了，Good's coverage 一般都在 0.99 以上。
★ 你在百迈客做 16S 项目时，分析流程是怎样的？¶
参考答案：
标准 16S 分析流程：
- 数据质控：用 fastp 或 Cutadapt 去引物、去低质量碱基、去接头
- DADA2 去噪：在 QIIME 2 中运行 DADA2，得到 ASV 特征表和代表序列
- 物种注释：用 classify-sklearn 分类器比对 SILVA 138 数据库，得到各 ASV 的分类信息
- Alpha 多样性：计算 Shannon、Simpson、Chao1、Observed species 等指标，做组间比较（Wilcoxon/Kruskal-Wallis）
- Beta 多样性：计算 Bray-Curtis 和 UniFrac 距离，做 PCoA 可视化，PERMANOVA 检验
- 物种组成：绘制各分类水平（门、纲、目、科、属）的物种组成柱状图
- 差异分析：用 LEfSe 找标志物种，必要时用 DESeq2 补充
- 功能预测（客户需要时）：用 PICRUSt2 预测 KEGG 通路
- 出报告：整理结果图表，撰写分析报告
💡 关键细节： - 引物序列必须去干净，否则影响 DADA2 去噪和物种注释 - DADA2 的 trunc-len 参数要根据质量图调整，保证拼接重叠 >= 12 bp - 稀释深度要根据最小测序量确定，太低会丢样本，太高会丢信息 - 报告中要给 Good's coverage，证明测序深度充分
5. 易错/易混淆点¶
⚠️ Shannon vs Simpson 方向相反¶
- Shannon 越大 = 多样性越高
- Simpson D 越大 = 多样性越低（优势种越突出）
- 软件输出的 "Simpson" 可能是 D、1-D 或 1/D，需查看文档
- R 的 vegan 包
diversity(x, "simpson")返回的是 1-D（方向与 Shannon 一致） - QIIME 2 的 Simpson 也是 1-D
⚠️ Weighted vs Unweighted UniFrac 的选择¶
- Weighted UniFrac：考虑丰度，对优势种敏感
- 适合：关注整体群落结构变化
- 可能被少数高丰度物种主导
- Unweighted UniFrac：不考虑丰度，只看有/无
- 适合：关注稀有种或独有种的差异
- 可能被测序噪音（低丰度假阳性）影响
- 建议：两者都做，结果一致 = 结论稳健
⚠️ 稀释（Rarefaction）vs 不稀释的争议¶
传统做法（稀释）： - 将所有样本随机抽取到相同测序深度 - 优点：消除测序量差异的影响，简单公平 - 缺点：丢弃大量数据（深测序样本的数据被浪费）
反对方（不稀释）： - McMurdie & Holmes (2014) 提出稀释是"不可接受的"，建议用统计模型（如 DESeq2）处理 - 理由：稀释丢弃有效数据，降低统计功效
当前共识： - Alpha 多样性：稀释仍然是主流做法，因为丰富度指标对测序深度极其敏感 - Beta 多样性：可以不稀释，用相对丰度或 CSS/TMM 标准化 - 差异分析：用 DESeq2、ANCOM-BC 等模型方法，不需要稀释 - 实际操作：QIIME 2 的 core-metrics 默认稀释
⚠️ 相对丰度 vs 绝对丰度¶
- 扩增子测序得到的是相对丰度（relative abundance），不是绝对丰度
- 相对丰度是"组分数据"（compositional data），总和恒为 100%
- 一个物种看起来"增加"可能只是因为其他物种"减少"了
- 要获得绝对丰度，需要：
- 加入内标（spike-in，如已知浓度的合成 DNA）
- 结合 qPCR 定量总细菌量
- 使用流式细胞术计数
📌 面试中被问"16S 结果能代表绝对丰度吗"，答"不能，16S 只能给相对丰度，要绝对丰度需要 spike-in 或 qPCR 辅助。"
⚠️ p 值校正（多重检验问题）¶
- 同时检验多个物种的差异时，假阳性概率会膨胀
- 例如：检验 100 个物种，p < 0.05，期望有 5 个假阳性
- 必须做多重检验校正：
- FDR（BH 法）：控制假发现率，最常用
- Bonferroni：最严格，适合少量比较
- R 代码：
p.adjust(p_values, method = "BH") - 校正后的 p 值称为 q 值（q-value）或 adjusted p-value
⚠️ 其他常见错误¶
| 错误 | 正确做法 |
|---|---|
| 不去引物序列直接去噪 | Cutadapt 去引物后再 DADA2 |
| DADA2 的 trunc-len 设太短导致无法拼接 | 确保 F + R - 扩增长度 >= 12（重叠区） |
| 用错分类器（V4 区分类器注释 V3-V4 数据） | 分类器必须与引物区段匹配 |
| PERMANOVA 显著就说差异大 | 还要看 R^2，R^2 小说明分组解释力弱 |
| 只报 p 值不报效应量 | p 值 + R^2（PERMANOVA）/ LDA score（LEfSe） |
| 用 Euclidean 距离做群落分析 | 微生物组数据用 Bray-Curtis 或 UniFrac |
| 不检查方差齐性就做 PERMANOVA | 先 betadisper + permutest 检验 |
6. 术语速查表（中英对照）¶
| 中文 | 英文 | 缩写 |
|---|---|---|
| 扩增子测序 | Amplicon sequencing | - |
| 操作分类单元 | Operational Taxonomic Unit | OTU |
| 扩增子序列变体 | Amplicon Sequence Variant | ASV |
| 阿尔法多样性 | Alpha diversity | - |
| 贝塔多样性 | Beta diversity | - |
| 香农指数 | Shannon-Wiener index | H' |
| 辛普森指数 | Simpson's index | D |
| 物种丰富度 | Species richness | - |
| 均匀度 | Evenness | J |
| 稀释曲线 | Rarefaction curve | - |
| 主坐标分析 | Principal Coordinates Analysis | PCoA |
| 非度量多维标度 | Non-metric Multidimensional Scaling | NMDS |
| 置换多元方差分析 | Permutational MANOVA | PERMANOVA |
| 线性判别分析效应量 | LDA Effect Size | LEfSe |
| 差异丰度分析 | Differential abundance analysis | - |
| 物种注释 | Taxonomic classification | - |
| 系统发育树 | Phylogenetic tree | - |
| 相对丰度 | Relative abundance | - |
| 绝对丰度 | Absolute abundance | - |
| 测序深度 | Sequencing depth | - |
| 去噪 | Denoising | - |
| 嵌合体 | Chimera | - |
| 多重检验校正 | Multiple testing correction | FDR |
| 假发现率 | False Discovery Rate | FDR |
| 组分数据 | Compositional data | - |
| 中心对数比变换 | Centered Log-Ratio transform | CLR |
7. 推荐阅读¶
| 资源 | 说明 |
|---|---|
| QIIME 2 官方教程 | 最权威的扩增子分析教程 |
| DADA2 Tutorial | DADA2 R 版教程 |
| Microbiome Helper | 宏基因组/扩增子分析实用工具集 |
| Callahan et al., 2016 (Nature Methods) | DADA2 原始论文 |
| Segata et al., 2011 (Genome Biology) | LEfSe 原始论文 |
| Anderson, 2001 (Austral Ecology) | PERMANOVA 方法论文 |
| Gloor et al., 2017 (Frontiers in Microbiology) | 组分数据分析综述 |
最后更新：2026-04-27 | 面试前重点复习 ★ 标记的问答和 ⚠️ 易混淆点