16S 微生物多样性分析全流程

01QIIME2 数据导入数据准备▶

QIIME2 使用 .qza（Artifact）格式统一管理数据，保证来源可追溯。.qzv 是可视化文件，在 view.qiime2.org 打开。

# 导入双端 FASTQ（Casava 1.8 格式） qiime tools import \ --type 'SampleData[PairedEndSequencesWithQuality]' \ --input-path raw_reads/ \ --input-format CasavaOneEightSingleLanePerSampleDirFmt \ --output-path paired-end-demux.qza # 查看质量分布 qiime demux summarize \ --i-data paired-end-demux.qza \ --o-visualization demux-summary.qzv

📖 参数逐行讲解

--type 'SampleData[PairedEndSequencesWithQuality]' 声明数据类型：双端带质量值的测序数据（白话：告诉 QIIME2 "我给你的是双端测序的 FASTQ 文件"）

--input-path raw_reads/ 原始 FASTQ 文件所在目录（把所有样本的 R1/R2 文件放在这个文件夹里）

--input-format CasavaOneEightSingleLanePerSampleDirFmt Illumina Casava 1.8 命名格式（文件名必须像 SampleName_S1_L001_R1_001.fastq.gz，QIIME2 靠文件名识别样本）

--output-path paired-end-demux.qza 输出 .qza 格式（QIIME2 专用压缩包，包含数据 + 元信息 + 完整来源记录，保证可复现）

qiime demux summarize 生成质量分布可视化（.qzv），用于确定后续 DADA2 的截断位置（看 Q30 从哪里开始掉）

--i-data 输入参数前缀 --i- = input artifact；--o- = output；--p- = parameter；--m- = metadata

.qza 文件本质

.qza = ZIP 压缩包，内部结构： ├── data/ ← 实际数据（FASTQ/BIOM/FASTA） ├── metadata.yaml ← 数据类型、格式 └── provenance/ ← 完整处理历史可用 qiime tools peek / export 查看/导出

QIIME2 核心设计是"可复现性"——每个 .qza 记录完整数据来源和处理链路。

02DADA2 去噪 → ASV 特征表核心步骤▶

DADA2 通过统计模型区分测序错误和真实变异，输出精确到单碱基的 ASV（扩增子序列变体）。

# DADA2 去噪（截断位置根据质量图确定） qiime dada2 denoise-paired \ --i-demultiplexed-seqs paired-end-demux.qza \ --p-trunc-len-f 240 \ --p-trunc-len-r 200 \ --p-trim-left-f 17 \ --p-trim-left-r 21 \ --p-n-threads 8 \ --o-table table.qza \ --o-representative-sequences rep-seqs.qza \ --o-denoising-stats stats.qza

📖 参数逐行讲解

--i-demultiplexed-seqs 输入上一步导入的 .qza 文件（已拆分好的双端序列）

--p-trunc-len-f 240 正向 reads 在第 240bp 截断（根据质量图，找 Q30 开始下降的位置，把后面质量差的尾巴砍掉）

--p-trunc-len-r 200 反向 reads 在第 200bp 截断（R2 质量通常比 R1 差，所以截得更短；但要保证 R1+R2 截断后仍有足够重叠区拼接）

--p-trim-left-f 17 正向去掉前 17bp = 正向引物长度（如 341F 是 17bp；引物序列不是目标区域，必须去掉）

--p-trim-left-r 21 反向去掉前 21bp = 反向引物长度（如 805R 是 21bp；不同引物对长度不同，要查自己实验用的引物）

--p-n-threads 8 多线程加速（DADA2 的错误模型学习阶段最耗时，8 线程在 100 个样本约需 2-4 小时）

--o-table table.qza 输出特征表（ASV x 样本的计数矩阵，白话：每个样本里每种细菌有多少条序列）

--o-representative-sequences rep-seqs.qza 输出代表序列（每个 ASV 的 DNA 序列，用于后续物种注释和建树）

--o-denoising-stats stats.qza 输出去噪统计（每步过滤了多少 reads，用来判断参数设得好不好）

去噪统计 — 真实输出

sample-id input filtered denoised merged non-chimeric SRR001 85432 78965 76234 71892 68453 SRR002 92145 85678 83456 79234 75890 每列：原始 → 质量过滤 → 去噪 → 双端合并 → 去嵌合体最终保留率 60-80% 正常，<50% 要警惕

ASV vs OTU — 面试必考

维度	OTU（97% 聚类）	ASV（精确去噪）
原理	按 97% 相似度聚类	统计模型区分错误和真实变异
分辨率	种水平勉强	株/亚种水平
可重现	不同运行结果不同	完全可重复
跨研究	不能直接比较	同序列全球统一
趋势	逐渐淘汰	当前主流

"OTU 把 97% 相似的混在一起，像把所有柯基都叫'狗'；ASV 精确到单碱基，能区分不同品种的柯基。"

03物种分类注释（SILVA / GTDB）分类▶

# Naive Bayes 分类器 qiime feature-classifier classify-sklearn \ --i-classifier silva-138-99-nb-classifier.qza \ --i-reads rep-seqs.qza \ --o-classification taxonomy.qza # 物种柱状图 qiime taxa barplot \ --i-table table.qza \ --i-taxonomy taxonomy.qza \ --m-metadata-file metadata.tsv \ --o-visualization taxa-bar-plots.qzv

📖 参数逐行讲解

qiime feature-classifier classify-sklearn 用 scikit-learn 的 Naive Bayes 分类器给 ASV 做物种注释（白话：拿未知序列去数据库里"查户口"）

--i-classifier silva-138-99-nb-classifier.qza 预训练分类器文件（基于 SILVA 138 数据库，99% 相似度聚类；需从 QIIME2 官网下载对应版本，~2GB）

--i-reads rep-seqs.qza 输入 DADA2 产出的代表序列（就是上一步的 ASV 序列）

--o-classification taxonomy.qza 输出分类结果（每个 ASV 对应 "域;门;纲;目;科;属;种" 七级分类 + 置信度分数）

qiime taxa barplot 生成交互式物种组成柱状图（.qzv），可在 view.qiime2.org 按不同分类层级查看

--i-table table.qza 输入特征表（提供每个 ASV 在每个样本中的丰度）

--i-taxonomy taxonomy.qza 输入上面的分类注释结果（提供每个 ASV 的物种名）

--m-metadata-file metadata.tsv 样本元数据（分组信息，让柱状图可以按组排列和着色）

注释结果格式

Feature ID Taxon Confidence 4b5eeb30... d__Bacteria;p__Firmicutes;...;g__Roseburia;s__intestinalis 0.99 a2f7dc88... d__Bacteria;p__Bacteroidota;...;g__Bacteroides;s__vulgatus 0.97 Confidence < 0.7 → 截断到更高层级

数据库	版本	特点
SILVA	138.2 (2024)	16S 首选，引用最高，CC-BY 4.0
GTDB	R232 (2026)	基因组系统发育，901K 基因组，最新分类
Greengenes2	2024.09	全长 16S + 基因组，UniFrac 分析用

GTDB 重命名：Firmicutes → Bacillota，Proteobacteria → Pseudomonadota
16S V3-V4 区种水平准确率~70%，属水平~90%

04Alpha 多样性 — 样本内部统计▶

# 核心多样性分析（含稀释） qiime diversity core-metrics-phylogenetic \ --i-phylogeny rooted-tree.qza \ --i-table table.qza \ --p-sampling-depth 10000 \ --m-metadata-file metadata.tsv \ --output-dir core-metrics-results/ # Alpha 分组显著性检验 qiime diversity alpha-group-significance \ --i-alpha-diversity core-metrics-results/shannon_vector.qza \ --m-metadata-file metadata.tsv \ --o-visualization shannon-significance.qzv

📖 参数逐行讲解

qiime diversity core-metrics-phylogenetic 一站式核心多样性分析（同时计算 Alpha + Beta 多样性，先抽平再算）

--i-phylogeny rooted-tree.qza 输入有根系统发育树（用于计算 Faith PD 和 UniFrac 距离；树由 qiime phylogeny align-to-tree-mafft-fasttree 生成）

--i-table table.qza 输入 ASV 特征表（DADA2 产出的计数矩阵）

--p-sampling-depth 10000 抽平深度（每个样本随机抽 10000 条 reads；低于此数的样本被丢弃。选值看 qiime feature-table summarize 输出，取最低样本量的 90% 左右）

--m-metadata-file metadata.tsv 样本元数据（分组、临床信息，用于生成分组 PCoA 图）

--output-dir core-metrics-results/ 输出目录（包含 shannon_vector.qza、observed_features_vector.qza、faith_pd_vector.qza、bray_curtis_distance_matrix.qza 等十几个文件）

qiime diversity alpha-group-significance 对 Alpha 多样性做分组间统计检验（Kruskal-Wallis 非参数检验）

--i-alpha-diversity shannon_vector.qza 输入某个 Alpha 指标的向量（可换成 faith_pd_vector.qza、observed_features_vector.qza 等）

指标	含义	白话
Shannon	H = -Σ(p_i × ln(p_i))	物种多+分布均匀→值高。最常用
Simpson	D = 1 - Σ(p_i²)	随机抓两个是不同物种的概率
Chao1	S + F1²/(2×F2)	估计你没检测到的稀有种
Faith PD	系统发育树总枝长	考虑进化关系的唯一 Alpha 指标

QIIME2 输出示例

Kruskal-Wallis (all groups) — Shannon H = 8.234, p-value = 0.0041 Healthy: 4.21 ± 0.38 | T2D: 3.68 ± 0.52 → T2D 组多样性显著降低

Alpha 看"样本内"：Shannon 越高=菌群越丰富均匀=通常更健康。检验用 Kruskal-Wallis（多组）或 Wilcoxon（两组）。

稀释争议：Schloss 2024 支持 rarefaction，McMurdie 2014 反对。共识：Alpha 分析时稀释合理

05Beta 多样性 — 样本间差异统计▶

# PERMANOVA 检验 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 \ --o-visualization bray-curtis-significance.qzv

📖 参数逐行讲解

qiime diversity beta-group-significance 对 Beta 多样性距离矩阵做分组间统计检验（白话：检验两组菌群整体结构是否有显著差异）

--i-distance-matrix bray_curtis_distance_matrix.qza 输入距离矩阵（来自 core-metrics-results/，可换成 weighted_unifrac、unweighted_unifrac 等）

--m-metadata-file metadata.tsv 样本元数据文件（提供分组信息）

--m-metadata-column Group 指定元数据中的分组列名（如 Group、Treatment、Disease 等，必须是分类变量）

--p-method permanova 统计方法选 PERMANOVA（置换多元方差分析，非参数，不要求正态分布；也可选 anosim、permdisp）

--o-visualization 输出 .qzv 可视化文件（包含 p 值、检验统计量、组间箱线图）

距离	原理	进化关系	适用
Bray-Curtis	丰度差异	否	最通用首选
Weighted UniFrac	丰度+进化距离	是	需系统发育树
Unweighted UniFrac	有无+进化距离	是	关注稀有种
Jaccard	有无（0/1）	否	简单比较

PERMANOVA 真实输出

adonis2(dist_bray ~ Group, data=metadata, permutations=999) Df SumOfSqs R2 F Pr(>F) Group 1 0.342 0.058 2.84 0.012 * Residual 162 5.486 0.942

必须配合 betadisper：验证方差齐性，否则 PERMANOVA 显著性可能是假的。

报告格式："PERMANOVA on Bray-Curtis, R²=0.058, F=2.84, p=0.012, 999 permutations。betadisper p=0.18（方差齐性成立）。"

06差异物种分析（LEfSe / ANCOM / DESeq2）差异▶

方法	原理	优点	注意
LEfSe	KW + LDA 效应量	出图直观，引用最高	不控制混杂因素
ANCOM-BC	对数线性 + 偏差校正	专为组成数据设计	零值多时不稳
DESeq2	负二项 + 经验贝叶斯	处理 count，自带 FDR	输入必须原始 count
MaAsLin2	多变量关联模型	可加协变量	计算慢

LEfSe 输出

Feature Group LDA p_value g__Faecalibacterium Healthy 4.23 0.001 g__Roseburia Healthy 3.87 0.003 g__Escherichia T2D 3.92 0.002 LDA > 2.0 且 p < 0.05 → 有生物学意义

DESeq2 输出

baseMean log2FC lfcSE stat pvalue padj ASV_001 181.56 2.20 0.43 5.13 2.9e-07 0.0001 padj < 0.05 = BH 校正后显著 log2FC > 0 = 疾病组富集, < 0 = 健康组富集

三选一："发文章用 LEfSe + ANCOM-BC 双验证。有混杂因素用 MaAsLin2。要 count 级用 DESeq2。"

维度	16S 扩增子	宏基因组 Shotgun
原理	只测 16S rRNA 片段	随机打断测所有 DNA
分辨率	通常到属	到种/株
功能信息	无	有（基因、通路）
成本	~200元/样本	~2000元/样本
数据量	~10万 reads	~5000万 reads
分析	QIIME2 一站式	多工具组合
PCR偏好	有	无

一图看懂 16S 分析流程

ASV vs OTU — 面试必考

16S vs 宏基因组测序

面试高频追问