QIIME 2 扩增子全流程分析（16S / ITS）

适用场景

当用户提出"用 QIIME 2 跑我的扩增子数据""做 16S 多样性分析""ITS 真菌菌群分析"时，完成从原始 FASTQ 到统计结论的完整流程：导入 → 去噪 → 物种注释 → 建树 → 稀释曲线 → 多样性与统计检验 → 差异丰度 → 溯源与导出。

QIIME 2 相对 R（DADA2/phyloseq）流程的核心优势是每个 .qza 产物内嵌完整 provenance，可在分析结束后反向生成可执行脚本。请把溯源当作流程的一部分主动使用，而不是可选项——见下方「数据溯源与复现」章节。

版本兼容性与命名变更（重要）

QIIME 2 每季度发布一次（2026.1 / 2026.4 / 2026.7 / 2026.10…），发行版名称与环境文件名近两年发生过实质变更，不要照抄旧教程：

| 变更 | 生效版本 | 说明 | |---|---|---| | 框架 Q2F 更名为 rachis | 2026.1 起过渡 | Python 包层面的重命名 | | amplicon 发行版更名为 qiime2 | 2026.4 起 | 旧名 amplicon 仅适用于历史版本 | | 环境文件前缀 qiime2- → rachis- | 2026.4 起 | 仅适用于 >=2026.4 及开发版 | | macOS 不再构建 moshpit/pathogenome | 2026.1 之后 | 这两个发行版请用 HPC 或 Docker | | q2-dada2 n_reads_learn 参数弃用 | 2026.10 计划移除 | 见 2026.7 公告 |

执行前必做：不要硬编码版本号或 .yml URL。先运行以下命令确认真实环境，再据此调整命令与参数：

qiime info                      # 确认发行版、版本号、已安装插件
qiime <插件名> --help           # 确认该插件的可用 action
qiime <插件> <action> --help    # 确认参数名（--p-* 参数改名较频繁）

安装命令请从官方 quickstart 获取当前版本对应的写法，按平台（Linux/WSL、macOS Intel、macOS Apple Silicon）选择： https://library.qiime2.org/quickstart/amplicon

若命令报 Error: no such option 或 Plugin error，先用 --help 核对实际参数名并调整，不要反复重试同一条命令。

核心概念

Artifact（.qza）：数据 + 元数据 + 完整 provenance 的压缩包。所有中间产物都是 .qza。
Visualization（.qzv）：交互式可视化，用 qiime tools view 本地打开，或上传至 https://view.qiime2.org（该站点仅在浏览器本地解析，不上传到服务器）。
Metadata（.tsv）：样本分组信息，几乎每个统计步骤都要用。
半随机性：DADA2、稀释抽样等步骤含随机过程，如需严格复现请固定 --p-random-seed（若该 action 支持）。

流程总览

FASTQ ──导入──> demux.qza ──去噪──> table.qza + rep-seqs.qza
                                        │
                    ┌───────────────────┼───────────────────┐
                    ↓                   ↓                   ↓
              taxonomy.qza        rooted-tree.qza      稀释曲线（选深度）
                    │                   │                   │
                    └───────────────────┴───────────────────┘
                                        ↓
                            core-metrics-phylogenetic
                                        ↓
                      alpha/beta 显著性检验 + ANCOM-BC 差异丰度
                                        ↓
                              provenance 溯源 + 导出

步骤 0：准备输入文件

manifest.tsv（制表符分隔，路径必须为绝对路径）：

sample-id	forward-absolute-filepath	reverse-absolute-filepath
sample1	/abs/path/sample1_R1.fastq.gz	/abs/path/sample1_R2.fastq.gz
sample2	/abs/path/sample2_R1.fastq.gz	/abs/path/sample2_R2.fastq.gz

metadata.tsv（首列必须是 sample-id；可选第二行为类型声明行）：

sample-id	Group	Timepoint	Age
#q2:types	categorical	categorical	numeric
sample1	Treatment	Day0	34
sample2	Control	Day0	41

导入前先校验元数据，避免后续统计步骤才报错：

qiime tools inspect-metadata metadata.tsv
qiime metadata tabulate --m-input-file metadata.tsv --o-visualization metadata.qzv

元数据列类型陷阱：形如 1、2、3 的分组编号会被识别为 numeric，导致 beta-group-significance 报错。用 #q2:types 行显式声明为 categorical。

步骤 1：导入数据

# 双端，manifest 方式（最通用）
qiime tools import \
    --type 'SampleData[PairedEndSequencesWithQuality]' \
    --input-path manifest.tsv \
    --output-path demux.qza \
    --input-format PairedEndFastqManifestPhred33V2

# 单端
qiime tools import \
    --type 'SampleData[SequencesWithQuality]' \
    --input-path manifest_se.tsv \
    --output-path demux.qza \
    --input-format SingleEndFastqManifestPhred33V2

# Casava 1.8 目录格式（文件名形如 L2S357_15_L001_R1_001.fastq.gz）
qiime tools import \
    --type 'SampleData[PairedEndSequencesWithQuality]' \
    --input-path casava_dir/ \
    --output-path demux.qza \
    --input-format CasavaOneEightSingleLanePerSampleDirFmt

# 质量总览 —— 用它决定截断位置，必看
qiime demux summarize --i-data demux.qza --o-visualization demux.qzv

步骤 2：去除引物（16S 与 ITS 都建议做）

若测序数据仍含引物序列，必须先去除，否则会被误当作生物学变异保留：

qiime cutadapt trim-paired \
    --i-demultiplexed-sequences demux.qza \
    --p-front-f CCTACGGGNGGCWGCAG \
    --p-front-r GACTACHVGGGTATCTAATCC \
    --p-discard-untrimmed \
    --o-trimmed-sequences demux-trimmed.qza

上例为 16S V3-V4（341F/805R）。请按用户实际引物替换，不要照抄。

步骤 3A：DADA2 去噪（16S 主流路径）

截断长度是全流程最关键的人工决策。打开 demux.qzv 的交互式质量图，选择质量中位数开始明显下降前的位置。硬性约束：

trunc_len_f + trunc_len_r - overlap >= 扩增子长度 + 12 正反读长截断后必须保留 至少 12 bp 重叠，否则合并率会骤降至接近 0。

qiime dada2 denoise-paired \
    --i-demultiplexed-seqs demux-trimmed.qza \
    --p-trunc-len-f 240 \
    --p-trunc-len-r 200 \
    --p-trim-left-f 0 \
    --p-trim-left-r 0 \
    --p-max-ee-f 2 \
    --p-max-ee-r 2 \
    --p-pooling-method independent \
    --p-chimera-method consensus \
    --p-n-threads 0 \
    --o-table table.qza \
    --o-representative-sequences rep-seqs.qza \
    --o-denoising-stats denoising-stats.qza

# 必看：检查每一步的留存率
qiime metadata tabulate \
    --m-input-file denoising-stats.qza \
    --o-visualization denoising-stats.qzv

去噪统计的诊断标准——出问题时先看这里，不要直接往下跑：

| 现象 | 原因 | 处理 | |---|---|---| | merged 比 denoised 骤降（<50%） | 重叠不足 | 增大 trunc-len，或改用单端 | | filtered 比 input 骤降 | 质量差或截断过松 | 缩短 trunc-len 或放宽 max-ee | | non-chimeric 骤降（<60%） | 引物未去除 | 回到步骤 2 | | 各样本留存率差异极大 | 批次/建库问题 | 检查是否需剔除低质量样本 |

步骤 3B：ITS 真菌路径（与 16S 的关键差异）

ITS 区域长度高度可变（可从 <200 bp 到 >500 bp），因此：

不能按固定位置截断 —— 必须设 --p-trunc-len-* 0，否则会系统性丢弃长 ITS 变体，造成物种组成偏倚。
先用 ITSxpress 提取 ITS 区，去除两侧保守的 18S/5.8S/28S 侧翼。
分类器必须用 UNITE，不能用 SILVA/Greengenes。

# 1) 提取 ITS2 区域（--p-taxa F 表示真菌；ITS1 请改 --p-region ITS1）
qiime itsxpress trim-pair-output-unmerged \
    --i-per-sample-sequences demux.qza \
    --p-region ITS2 \
    --p-taxa F \
    --p-threads 4 \
    --o-trimmed itsxpress-trimmed.qza

# 2) 去噪 —— 注意 trunc-len 必须为 0
qiime dada2 denoise-paired \
    --i-demultiplexed-seqs itsxpress-trimmed.qza \
    --p-trunc-len-f 0 \
    --p-trunc-len-r 0 \
    --p-max-ee-f 2 \
    --p-max-ee-r 2 \
    --p-n-threads 0 \
    --o-table table.qza \
    --o-representative-sequences rep-seqs.qza \
    --o-denoising-stats denoising-stats.qza

ITS 的另一个差异：ITS 序列无法可靠地做多序列比对，因此不要构建系统发育树，也不要使用 UniFrac 等系统发育距离。改用 qiime diversity core-metrics （非 phylogenetic 版本），指标限于 Bray-Curtis、Jaccard、Shannon、Observed features。

步骤 3C：Deblur 替代路径（仅 16S）

qiime quality-filter q-score \
    --i-demux demux.qza \
    --o-filtered-sequences demux-filtered.qza \
    --o-filter-stats filter-stats.qza

qiime deblur denoise-16S \
    --i-demultiplexed-seqs demux-filtered.qza \
    --p-trim-length 250 \
    --p-sample-stats \
    --o-representative-sequences rep-seqs.qza \
    --o-table table.qza \
    --o-stats deblur-stats.qza

Deblur 仅对 16S 有专用参考（denoise-16S），且需先合并双端。默认优先用 DADA2。

步骤 4：物种注释

分类器必须与扩增子区域匹配。全长分类器用于 V4 数据会降低属水平准确度，建议对特定引物区做 extract-reads 后重训练。

# 从官方 Data Resources 下载当前版本的预训练分类器（勿硬编码版本）：
#   https://library.qiime2.org/data-resources
#   16S -> SILVA 或 Greengenes2；ITS -> UNITE

qiime feature-classifier classify-sklearn \
    --i-classifier silva-nb-classifier.qza \
    --i-reads rep-seqs.qza \
    --p-n-jobs -1 \
    --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-barplot.qzv

针对自己引物区训练分类器（显著提升准确度）：

qiime feature-classifier extract-reads \
    --i-sequences silva-seqs.qza \
    --p-f-primer GTGCCAGCMGCCGCGGTAA \
    --p-r-primer GGACTACHVGGGTWTCTAAT \
    --p-trunc-len 0 \
    --o-reads ref-seqs-v4.qza

qiime feature-classifier fit-classifier-naive-bayes \
    --i-reference-reads ref-seqs-v4.qza \
    --i-reference-taxonomy silva-tax.qza \
    --o-classifier classifier-v4.qza

注意：sklearn 分类器与训练时的 scikit-learn 版本强绑定。跨 QIIME 2 版本使用会报版本不匹配错误，必须下载对应版本的分类器或重新训练。

步骤 5：过滤（按需）

# 移除线粒体与叶绿体污染（植物/宿主样本常见）
qiime taxa filter-table \
    --i-table table.qza \
    --i-taxonomy taxonomy.qza \
    --p-exclude mitochondria,chloroplast \
    --o-filtered-table table-filtered.qza

# 移除极低频特征（例如总频次 < 10 的 ASV）
qiime feature-table filter-features \
    --i-table table-filtered.qza \
    --p-min-frequency 10 \
    --o-filtered-table table-final.qza

# 查看每样本测序深度 —— 用于下一步选采样深度
qiime feature-table summarize \
    --i-table table-final.qza \
    --m-sample-metadata-file metadata.tsv \
    --o-visualization table-final.qzv

步骤 6：系统发育树（仅 16S；ITS 跳过）

qiime phylogeny align-to-tree-mafft-fasttree \
    --i-sequences rep-seqs.qza \
    --p-n-threads auto \
    --o-alignment aligned-rep-seqs.qza \
    --o-masked-alignment masked-aligned-rep-seqs.qza \
    --o-tree unrooted-tree.qza \
    --o-rooted-tree rooted-tree.qza

步骤 7：稀释曲线 —— 选择采样深度（不可跳过）

--p-sampling-depth 是一个权衡：depth 越高每样本信息越充分，但低于该深度的样本会被直接丢弃。必须先看曲线再定值，不要沿用 10000 这类默认数字。

qiime diversity alpha-rarefaction \
    --i-table table-final.qza \
    --i-phylogeny rooted-tree.qza \
    --p-max-depth 20000 \
    --p-steps 20 \
    --m-metadata-file metadata.tsv \
    --o-visualization alpha-rarefaction.qzv

判读方法：

上半部分曲线在何处趋于平台（说明该深度已捕获大部分多样性）。
下半部分显示各分组在每个深度还剩多少样本。
取「曲线已平台」且「保留样本数可接受」的最小深度。
结合 table-final.qzv 的每样本深度表，明确告知用户哪些样本会被丢弃。

步骤 8：多样性分析与统计检验

# 16S（含系统发育指标：Faith PD、UniFrac）
qiime diversity core-metrics-phylogenetic \
    --i-phylogeny rooted-tree.qza \
    --i-table table-final.qza \
    --p-sampling-depth 10000 \
    --m-metadata-file metadata.tsv \
    --output-dir core-metrics-results

# ITS 或无树时（仅非系统发育指标）
qiime diversity core-metrics \
    --i-table table-final.qza \
    --p-sampling-depth 10000 \
    --m-metadata-file metadata.tsv \
    --output-dir core-metrics-results

Alpha 多样性组间检验（Kruskal-Wallis）

for METRIC in shannon observed_features faith_pd evenness; do
    [ -f core-metrics-results/${METRIC}_vector.qza ] || continue
    qiime diversity alpha-group-significance \
        --i-alpha-diversity core-metrics-results/${METRIC}_vector.qza \
        --m-metadata-file metadata.tsv \
        --o-visualization ${METRIC}-significance.qzv
done

# 与连续变量的相关性（Spearman）
qiime diversity alpha-correlation \
    --i-alpha-diversity core-metrics-results/shannon_vector.qza \
    --m-metadata-file metadata.tsv \
    --p-method spearman \
    --o-visualization shannon-correlation.qzv

Beta 多样性组间检验（PERMANOVA）

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 \
    --p-permutations 999 \
    --o-visualization weighted-unifrac-permanova.qzv

PERMANOVA 显著必须配套检验组内离散度——否则无法区分「组间中心不同」与「组内变异度不同」，这是扩增子分析中最常见的解读错误：

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 permdisp \
    --o-visualization weighted-unifrac-permdisp.qzv

判读：PERMANOVA 显著 + PERMDISP 不显著 → 可解释为组间群落组成差异；两者均显著 → 差异可能部分来自离散度不同，需在结论中明确说明。

多因素与重复测量

# 控制协变量的 PERMANOVA（adonis 公式语法）
qiime diversity adonis \
    --i-distance-matrix core-metrics-results/weighted_unifrac_distance_matrix.qza \
    --m-metadata-file metadata.tsv \
    --p-formula "Group+Age+Timepoint" \
    --p-permutations 999 \
    --o-visualization adonis.qzv

# 纵向数据：线性混合效应模型
qiime longitudinal linear-mixed-effects \
    --m-metadata-file metadata.tsv \
    --m-metadata-file core-metrics-results/shannon_vector.qza \
    --p-metric shannon_entropy \
    --p-group-columns Group \
    --p-state-column Timepoint \
    --p-individual-id-column SubjectID \
    --o-visualization lme-shannon.qzv

adonis 公式中项的顺序影响结果（序贯平方和），把要控制的协变量放在前面。

主坐标分析可视化

qiime emperor plot \
    --i-pcoa core-metrics-results/weighted_unifrac_pcoa_results.qza \
    --m-metadata-file metadata.tsv \
    --o-visualization weighted-unifrac-emperor.qzv

步骤 9：差异丰度（ANCOM-BC）

旧的 qiime composition ancom 与 add-pseudocount 已弃用，使用 ANCOM-BC （自带偏倚校正与零值处理，无需手工加伪计数）：

qiime taxa collapse \
    --i-table table-final.qza \
    --i-taxonomy taxonomy.qza \
    --p-level 6 \
    --o-collapsed-table table-l6.qza

qiime composition ancombc \
    --i-table table-l6.qza \
    --m-metadata-file metadata.tsv \
    --p-formula Group \
    --o-differentials ancombc-l6.qza

qiime composition da-barplot \
    --i-data ancombc-l6.qza \
    --p-significance-threshold 0.05 \
    --o-visualization ancombc-l6.qzv

--p-level：2=门 3=纲 4=目 5=科 6=属 7=种。属水平（6）最常用； ASV 水平（不 collapse）零值过多，统计效力低。

差异丰度输入不要用稀释后的表，直接用过滤后的原始计数表 —— ANCOM-BC 内部自行处理测序深度差异。

数据溯源与复现（QIIME 2 的核心价值）

每个 .qza/.qzv 都内嵌了生成它的完整命令链、参数、软件版本与引用。分析结束后应主动产出溯源材料，这是 QIIME 2 相对 R 流程的最大优势。

# 快速查看类型、UUID、数据格式
qiime tools peek core-metrics-results/weighted_unifrac_distance_matrix.qza

# 校验完整性（MD5）—— 文件被改动过会在此报出
qiime tools validate table-final.qza

# 【最重要】从任一结果反推出可执行的 bash 脚本
qiime tools replay-provenance \
    --in-fp ancombc-l6.qzv \
    --out-fp reproduce-analysis.sh

# 生成 Python API 版本
qiime tools replay-provenance \
    --in-fp ancombc-l6.qzv \
    --p-usage-driver python3 \
    --out-fp reproduce-analysis.py

# 自动汇总本次分析所有方法的文献引用（写论文方法学部分直接可用）
qiime tools replay-citations \
    --in-fp ancombc-l6.qzv \
    --out-fp citations.bib

# 一键打包完整复现材料（脚本 + 引用 + 环境清单），投稿补充材料用
qiime tools replay-supplement \
    --in-fp . \
    --p-recurse \
    --o-out-fp reproducibility-supplement.zip

给用户交付结果时，默认同时产出 replay-provenance 脚本与 replay-citations，并说明：.qza 文件本身即是完整的方法学记录，无需另行维护流程文档。

溯源工具由 provenance-lib 提供，自 QIIME 2 2023.5 起已包含在发行版中，无需单独安装。

导出到 R / Python

qiime tools export --input-path table-final.qza --output-path exported/
biom convert -i exported/feature-table.biom -o feature-table.tsv --to-tsv

qiime tools export --input-path taxonomy.qza --output-path exported/
qiime tools export --input-path rooted-tree.qza --output-path exported/

导出到 phyloseq 时注意：需将 taxonomy.tsv 表头改为 #OTUID taxonomy confidence 后才能用 biom add-metadata 合并。

导出即脱离 provenance 追踪 —— 导出后的下游分析请自行记录代码版本。

常见错误速查

| 报错/现象 | 原因 | 处理 | |---|---|---| | Plugin error from dada2 无细节 | 需看完整日志 | 加 --verbose，读 /tmp/qiime2-q2cli-err-*.log | | 合并率极低 | 重叠不足 | 增大 trunc-len，或改单端分析 | | All features were filtered | 采样深度高于所有样本 | 看 table.qzv 重选 depth | | not a categorical column | 分组列被识别为数值 | 元数据加 #q2:types 行 | | 分类器报版本不匹配 | scikit-learn 版本绑定 | 下载对应版本分类器或重训练 | | ITS 物种组成异常偏倚 | 误用了固定截断 | trunc-len 设 0，见步骤 3B | | --p-* 参数不存在 | 跨版本参数改名 | qiime <插件> <action> --help 核对 |

执行守则

先 qiime info，确认版本与插件，再据实际版本调整命令，不要照抄本文档的版本假设。
两个决策点必须让用户看图后再定：截断长度（看 demux.qzv）、采样深度（看 alpha-rarefaction.qzv）。不要擅自使用默认值就往下跑。
每一步产出 .qzv 并提示用户查看，尤其 denoising-stats.qzv —— 去噪失败继续往下跑毫无意义。
区分 16S 与 ITS 路径：ITS 不截断、不建树、用 UNITE、不用 UniFrac。
PERMANOVA 显著时必须补 PERMDISP，并在结论中如实说明。
收尾必做溯源：replay-provenance + replay-citations。
长时间步骤（DADA2、建树、classify-sklearn）建议后台运行并明确告知预期耗时。

东方天意QIIME2 微生物组扩增子分析全流程