16S rRNA 扩增子分析全流程¶

一、基础概念与实验设计¶

1.1 16S rRNA 基因简介¶

• 原核生物核糖体小亚基（30S）的组成部分
• 全长 ~1542 bp
• 包含 9 个可变区（V1–V9）+ 保守区交替排列
• 保守区：设计通用引物
• 可变区：区分不同物种
• 为什么选 16S？
• 所有细菌/古菌都有（通用性）
• 进化速率适中（既能区分物种，又不至于太快丢失信号）
• 数据库最完善（Silva、Greengenes2、RDP）

1.2 16S vs 宏基因组¶

• 16S 扩增子
• 仅测 16S rRNA 基因的特定可变区
• 分辨率通常到属级，种级不太可靠
• 只能回答"谁在那里"
• 功能只能预测（PICRUSt2），不够准确
• 成本低、分析简单
• 受引物偏好性影响
• 宏基因组（Shotgun）
• 对所有 DNA 随机打断测序
• 分辨率可达种级/菌株级
• 同时回答"谁在那里 + 在做什么"
• 功能分析直接基于基因（HUMAnN3）
• 成本高、计算资源需求大
• 偏好性小

1.3 引物选择¶

• V3-V4 区（最常用）
• 引物：341F / 805R
• 341F：5'-CCTACGGGNGGCWGCAG-3'
• 805R：5'-GACTACHVGGGTATCTAATCC-3'
• 扩增片段 ~460 bp
• 适用平台：Illumina PE250 / PE300
• 优势：覆盖度广、数据库支持好
• 劣势：对部分古菌覆盖不足
• V4 区
• 引物：515F / 806R（Earth Microbiome Project 标准）
• 515F：5'-GTGYCAGCMGCCGCGGTAA-3'
• 806R：5'-GGACTACNVGGGTWTCTAAT-3'
• 扩增片段 ~253 bp
• 适用平台：Illumina PE150 / PE250
• 优势：EMP 标准、跨研究可比性强
• 劣势：分辨率略低于 V3-V4
• V1-V3 区
• 引物：27F / 534R
• 扩增片段 ~500 bp
• 对口腔/皮肤微生物分辨率较好
• 引物特异性验证
• BLAST 验证
  • blastn -query primer.fa -db nt -task blastn-short -evalue 1000 -word_size 7 -outfmt 6 -num_threads 4
• Bowtie2 验证
  • bowtie2 -x ref_genome -f primers.fa -S primers_aligned.sam --very-sensitive -k 10
  • -k 10：报告最多 10 个比对位置（检查脱靶）
• ⚠️ 注意
• 引物选择直接决定能检测到哪些物种
• 不同引物的结果不宜直接比较
• 引物序列中的简并碱基（N/W/H/V 等）是正常的

1.4 测序设计¶

• 测序平台
• Illumina MiSeq PE300（最常用）
• Illumina MiSeq PE250
• Illumina NovaSeq PE250
• V3-V4 区必须用 PE250 或 PE300（片段 ~460 bp，需要双端重叠）
• V4 区可用 PE150（片段 ~253 bp）
• 测序深度
• 常规分析：10,000–50,000 reads/样本
• 深度分析（稀有物种）：50,000–100,000 reads/样本
• 最低不少于 10,000 reads/样本
• 超过 100,000 reads/样本通常边际效益递减
• 样本量
• 每组最低 ≥ 5 个
• 建议 ≥ 10 个（提高统计效力）
• 有协变量时需更多样本
• 样本保存
• -80°C 冻存（金标准）
• DNA/RNA Shield 室温保存（野外采样）
• 避免反复冻融

1.5 文库质控（上机前三步验证）¶

• 第一步：Qubit 荧光定量
• 测总 dsDNA 浓度
• 包含所有 DNA（有效文库 + 废物）
• 第二步：Bioanalyzer / TapeStation
• 检查片段大小分布
• 应出现单一主峰（如 V3-V4 在 ~550–620 bp 含接头）
• 无接头二聚体峰（~120–170 bp 处不应有峰）
• 第三步：qPCR 定量（KAPA qPCR Kit）
• 只测带有完整接头（P5+P7）的分子浓度
• 这是最准确的定量，直接决定上机 loading 浓度
• ⚠️ 注意
• Qubit 像"称体重"（不分肌肉脂肪）
• Bioanalyzer 像"照 X 光"（看结构）
• qPCR 像"量肌肉量"（只算有用的）
• 三步都做才能保证上机质量

1.6 ASV vs OTU 决策¶

• ASV（Amplicon Sequence Variant）
• 工具：DADA2 / Deblur
• 单碱基分辨率
• 可跨研究比较（序列即 ID）
• 当前主流，推荐使用
• 不需要设定相似度阈值
• OTU（Operational Taxonomic Unit）
• 工具：VSEARCH / UPARSE
• 97% 相似度聚类
• 传统方法，逐渐被 ASV 取代
• 不同研究的 OTU 不可直接比较
• 可能合并不同物种或拆分同一物种
• 推荐：使用 ASV（DADA2）

二、数据预处理（QIIME 2 流程）¶

2.0 环境准备¶

• 安装 QIIME 2
• conda activate qiime2-amplicon-2024.5
• 版本号随年份更新
• 所有 QIIME 2 文件格式
• .qza：数据文件（QIIME 2 Artifact）
• .qzv：可视化文件（用 view.qiime2.org 查看）
• 本质上都是 zip 压缩包，可用 qiime tools export 导出

2.1 数据导入¶

• 准备 manifest 文件
• TSV 格式，三列：sample-id、forward-read、reverse-read
• 路径必须是绝对路径
• 示例
  • sample-id forward-read reverse-read
  • Sample1 /path/to/Sample1_R1.fastq.gz /path/to/Sample1_R2.fastq.gz
• 导入命令
• qiime tools import --type 'SampleData[PairedEndSequencesWithQuality]' --input-path manifest.tsv --output-path paired-end-demux.qza --input-format PairedEndFastqManifestPhred33V2
• 关键参数
• --type：数据类型（双端带质量）
• --input-format：Phred33V2（Illumina 标准质量编码）
• 如果已经 demultiplex：用 Manifest 格式
• 如果未 demultiplex：用 EMPPairedEndSequences 格式
• ⚠️ 注意
• manifest 文件中的路径错误是最常见的导入失败原因
• 确保 sample-id 不含特殊字符（只用字母、数字、短横线）
• Phred33 vs Phred64：现代 Illumina 数据都是 Phred33

2.2 序列质量可视化¶

• 命令
• qiime demux summarize --i-data paired-end-demux.qza --o-visualization paired-end-demux.qzv
• 查看方式
• 上传到 https://view.qiime2.org
• 或 qiime tools view paired-end-demux.qzv
• 关注内容
• Interactive Quality Plot：每个位置的质量分数分布
• 确定截断位置：找到质量开始下降的位置（通常 Q < 20–25）
• 正向 reads 通常质量优于反向 reads
• 记录每个样本的 reads 数量
• ⚠️ 注意
• 这一步是确定 DADA2 截断参数的关键依据
• 反向 reads 质量差是正常的（Illumina 特性）
• 截断后双端 reads 必须仍有足够重叠（≥ 20 bp）

2.3 去引物（Cutadapt）¶

• 为什么要去引物？
• 引物序列不是生物学信号，会干扰后续分析
• 引物区域的简并碱基会引入人为变异
• 命令
• qiime cutadapt trim-paired --i-demultiplexed-sequences paired-end-demux.qza --p-front-f CCTACGGGNGGCWGCAG --p-front-r GACTACHVGGGTATCTAATCC --p-discard-untrimmed --p-no-indels --o-trimmed-sequences trimmed-demux.qza --verbose
• 关键参数
• --p-front-f：正向引物序列（5'→3'）
• --p-front-r：反向引物序列（5'→3'）
• --p-discard-untrimmed：丢弃未找到引物的 reads（推荐开启）
• --p-no-indels：引物匹配时不允许插入缺失
• --p-error-rate 0.1：允许 10% 错配（默认值，通常不需修改）
• 验证
• 再次运行 qiime demux summarize 检查去引物后的数据
• reads 长度应减少引物长度（如 V3-V4 减少 ~17+21=38 bp）
• 保留率应 ≥ 80%
• ⚠️ 注意
• 引物序列必须与实际使用的完全一致（含简并碱基）
• 如果测序中心已去除引物，可跳过此步（但要确认！）
• 不去引物会导致 DADA2 产生虚假 ASV

三、降噪（Denoising）¶

3.1 DADA2 降噪（推荐）¶

• 原理
• 基于错误模型学习，区分真实生物变异和测序错误
• 自动完成：质量过滤 + 去噪 + 双端合并 + 嵌合体去除
• 输出单碱基分辨率的 ASV
• 命令
• qiime dada2 denoise-paired --i-demultiplexed-seqs trimmed-demux.qza --p-trunc-len-f 280 --p-trunc-len-r 200 --p-max-ee-f 2.0 --p-max-ee-r 2.0 --p-trunc-q 2 --p-min-overlap 20 --p-chimera-method consensus --p-n-threads 8 --o-table table.qza --o-representative-sequences rep-seqs.qza --o-denoising-stats stats.qza
• 关键参数详解
• --p-trunc-len-f 280：正向 reads 截断长度
  • 根据质量图确定：在质量中位数降到 Q20–25 之前截断
  • V3-V4 + PE300：正向通常 270–290
• --p-trunc-len-r 200：反向 reads 截断长度
  • 反向质量通常更差，截断更短
  • V3-V4 + PE300：反向通常 200–230
  • 关键约束：trunc-len-f + trunc-len-r - 扩增片段长度 ≥ 20（重叠区）
  • 例：280 + 200 - 460 = 20 bp 重叠 ✓
• --p-max-ee-f 2.0：正向最大期望错误数（默认 2.0）
  • 越小越严格，通常不需修改
• --p-max-ee-r 2.0：反向最大期望错误数
  • 反向质量差时可放宽到 5.0
• --p-trunc-q 2：遇到质量 ≤ 2 的碱基时截断（默认 2）
• --p-min-overlap 20：双端合并最小重叠长度（默认 12，建议 ≥ 20）
• --p-chimera-method consensus：嵌合体检测方法
  • consensus（默认）：逐样本检测，取共识
  • pooled：所有样本合并检测（更灵敏但更慢）
  • none：不检测嵌合体（不推荐）
• --p-n-threads 8：线程数（0 = 使用所有可用线程）
• 三大输出
• table.qza：ASV 特征表（样本 × ASV 计数矩阵）
• rep-seqs.qza：代表序列（每个 ASV 的实际 DNA 序列）
• stats.qza：降噪统计（每个样本的 reads 保留情况）
• 查看降噪统计
• qiime metadata tabulate --m-input-file stats.qza --o-visualization stats.qzv
• 关注列：input → filtered → denoised → merged → non-chimeric
• 最终保留率（non-chimeric / input）应 ≥ 50%
• 如果 merged 步骤损失大 → 截断长度设置不当，重叠不足
• 如果 filtered 步骤损失大 → 原始数据质量差或截断太严格
• ⚠️ 注意
• 截断长度是最关键的参数，必须根据质量图调整
• 不要盲目复制别人的截断长度（不同批次质量不同）
• DADA2 已包含质量过滤，不需要提前用 fastp/Trimmomatic
• 运行时间较长（大数据集可能数小时）

3.2 替代方案：Deblur¶

• 命令
• 先做质量过滤
• qiime quality-filter q-score --i-demux trimmed-demux.qza --o-filtered-sequences filtered.qza --o-filter-stats filter-stats.qza
• 再降噪
• qiime deblur denoise-16S --i-demultiplexed-seqs filtered.qza --p-trim-length 250 --p-sample-stats --o-table table.qza --o-representative-sequences rep-seqs.qza --o-stats stats.qza
• 关键参数
• --p-trim-length 250：统一截断到此长度（只处理单端或已合并的 reads）
• 特点
• 基于正向参考数据库的错误模型
• 只能处理等长序列
• 通常用于单端数据或已合并的双端数据
• 比 DADA2 更保守（假阳性更低，但灵敏度也更低）

3.3 替代方案：VSEARCH OTU 聚类（传统方法）¶

• 命令
• qiime vsearch cluster-features-de-novo --i-table table.qza --i-sequences rep-seqs.qza --p-perc-identity 0.97 --o-clustered-table otu-table.qza --o-clustered-sequences otu-seqs.qza
• 关键参数
• --p-perc-identity 0.97：97% 相似度阈值
• ⚠️ 注意
• 传统方法，不推荐新项目使用
• 如果审稿人要求 OTU，可在 DADA2 ASV 基础上再聚类

四、分类注释¶

4.1 数据库选择¶

• Silva（推荐）
• 当前版本：Silva 138.1（QIIME 2 兼容版）
• 覆盖：细菌 + 古菌 + 真核
• 更新频率高、分类体系规范
• 下载预训练分类器：QIIME 2 官网 Data Resources 页面
• Greengenes2（2022.10）
• 基于全基因组系统发育（GTDB 骨架）
• 与 GTDB 分类体系兼容
• 适合需要与宏基因组结果对比的研究
• RDP（Ribosomal Database Project）
• 经典数据库，更新较慢
• 分类器训练简单
• ⚠️ 注意
• 同一项目内所有样本必须用同一数据库
• 不同数据库的分类结果不宜直接比较

4.2 使用预训练分类器¶

• 下载
• 从 QIIME 2 官网下载与引物区域匹配的预训练分类器
• 例：Silva 138.1 V3-V4 区域（341F/805R）预训练分类器
• 命令
• qiime feature-classifier classify-sklearn --i-classifier silva-138.1-ssu-nr99-341f-805r-classifier.qza --i-reads rep-seqs.qza --p-n-jobs 8 --p-confidence 0.7 --o-classification taxonomy.qza
• 关键参数
• --p-confidence 0.7：置信度阈值（默认 0.7）
  • 低于此置信度的分类层级标记为"未分类"
  • 降低到 0.5 可获得更多注释但假阳性增加
• --p-n-jobs 8：并行线程数
• 可视化
• qiime metadata tabulate --m-input-file taxonomy.qza --o-visualization taxonomy.qzv

4.3 自训练分类器（可选，更精确）¶

• 适用场景：没有现成的区域特异性分类器
• 步骤一：导入参考序列和分类
• qiime tools import --type 'FeatureData[Sequence]' --input-path silva-138.1-ssu-nr99-seqs.fasta --output-path silva-seqs.qza
• qiime tools import --type 'FeatureData[Taxonomy]' --input-format HeaderlessTSVTaxonomyFormat --input-path silva-138.1-ssu-nr99-tax.txt --output-path silva-tax.qza
• 步骤二：提取引物区域
• qiime feature-classifier extract-reads --i-sequences silva-seqs.qza --p-f-primer CCTACGGGNGGCWGCAG --p-r-primer GACTACHVGGGTATCTAATCC --p-min-length 300 --p-max-length 600 --o-reads silva-v3v4-seqs.qza
• 步骤三：训练分类器
• qiime feature-classifier fit-classifier-naive-bayes --i-reference-reads silva-v3v4-seqs.qza --i-reference-taxonomy silva-tax.qza --o-classifier silva-v3v4-classifier.qza
• ⚠️ 注意
• 训练分类器需要较长时间（30 分钟–数小时）和大内存（≥ 32 GB）
• 区域特异性分类器比全长分类器更准确
• 训练一次后可反复使用

4.4 分类结果可视化¶

• 柱状图（Taxa Bar Plot）
• qiime taxa barplot --i-table table.qza --i-taxonomy taxonomy.qza --m-metadata-file metadata.tsv --o-visualization taxa-barplot.qzv
• 交互式：可切换分类层级（门/纲/目/科/属）
• Krona 交互式饼图
• 需先导出为 BIOM 格式再转换

五、数据过滤与清洗¶

5.1 去除非目标序列¶

• 去线粒体
• qiime taxa filter-table --i-table table.qza --i-taxonomy taxonomy.qza --p-exclude mitochondria --o-filtered-table table-no-mito.qza
• 去叶绿体
• qiime taxa filter-table --i-table table-no-mito.qza --i-taxonomy taxonomy.qza --p-exclude chloroplast --o-filtered-table table-no-mito-chloro.qza
• 去未分类（可选）
• qiime taxa filter-table --i-table table-no-mito-chloro.qza --i-taxonomy taxonomy.qza --p-include p__ --o-filtered-table table-filtered.qza
• --p-include p__：只保留至少注释到门级的 ASV
• 同步过滤代表序列
• qiime feature-table filter-seqs --i-data rep-seqs.qza --i-table table-filtered.qza --o-filtered-data rep-seqs-filtered.qza

5.2 低丰度/低频率过滤¶

• 去除低频率 ASV
• qiime feature-table filter-features --i-table table-filtered.qza --p-min-frequency 10 --o-filtered-table table-freq-filtered.qza
• --p-min-frequency 10：总 reads 数 < 10 的 ASV 被移除
• 去除低样本数 ASV（可选）
• qiime feature-table filter-features --i-table table-freq-filtered.qza --p-min-samples 2 --o-filtered-table table-final.qza
• --p-min-samples 2：只出现在 1 个样本中的 ASV 被移除
• 去除低 reads 样本
• qiime feature-table filter-samples --i-table table-final.qza --p-min-frequency 1000 --o-filtered-table table-sample-filtered.qza
• --p-min-frequency 1000：总 reads < 1000 的样本被移除
• ⚠️ 注意
• 过滤顺序：先去非目标 → 再去低丰度 → 最后去低 reads 样本
• 过滤阈值没有绝对标准，需根据数据情况调整
• 过滤前后都要检查特征表摘要

5.3 特征表摘要¶

• 命令
• qiime feature-table summarize --i-table table-final.qza --m-sample-metadata-file metadata.tsv --o-visualization table-summary.qzv
• 关注内容
• 总 ASV 数
• 每个样本的 reads 数分布
• 最小/最大/中位数 reads 数
• 用于确定后续稀释深度（rarefaction depth）

六、构建系统发育树¶

6.1 为什么需要建树？¶

• UniFrac 距离（加权/非加权）需要系统发育树
• Faith's PD（系统发育多样性）需要系统发育树
• 如果只用 Bray-Curtis / Shannon 等非系统发育指标，可跳过

6.2 快速建树流程（推荐）¶

• 命令（一步完成：比对 → 掩码 → 建树 → 有根树）
• qiime phylogeny align-to-tree-mafft-fasttree --i-sequences rep-seqs-filtered.qza --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 --p-n-threads 8
• 内部步骤
• MAFFT 多序列比对
• 掩码（去除高度可变的比对列）
• FastTree 建树（近似最大似然法）
• 中点有根化（Midpoint rooting）
• 四个输出
• aligned-rep-seqs.qza：比对后的序列
• masked-aligned-rep-seqs.qza：掩码后的比对
• unrooted-tree.qza：无根树
• rooted-tree.qza：有根树 → 用于后续多样性分析
• ⚠️ 注意
• ASV 数量很大时（> 50,000）建树可能很慢
• FastTree 是近似方法，精度足够用于多样性分析
• 如需更精确的树：用 RAxML 或 IQ-TREE（但更慢）

七、多样性分析¶

7.1 确定稀释深度（Rarefaction Depth）¶

• 为什么要稀释？
• 16S 扩增子中样本间测序深度差异大
• 深度越大，检测到的物种越多 → 不稀释会引入偏差
• 稀释 = 将所有样本随机抽取到相同 reads 数
• 稀释曲线
• qiime diversity alpha-rarefaction --i-table table-final.qza --i-phylogeny rooted-tree.qza --p-max-depth 50000 --p-steps 20 --m-metadata-file metadata.tsv --o-visualization alpha-rarefaction.qzv
• --p-max-depth：最大稀释深度（设为样本中位数 reads 数附近）
• --p-steps 20：计算 20 个深度点
• 如何选择稀释深度？
• 查看稀释曲线：曲线趋于平台期的最小深度
• 兼顾：保留尽可能多的样本 vs 足够的深度
• 通常选择：最小样本 reads 数的 90%，或曲线平台期起点

• 低于此深度的样本将被丢弃

  7.2 核心多样性一步完成（core-metrics-phylogenetic）¶

  • 命令
  • qiime diversity core-metrics-phylogenetic --i-phylogeny rooted-tree.qza --i-table table-final.qza --p-sampling-depth 10000 --m-metadata-file metadata.tsv --p-n-jobs-or-threads 8 --output-dir core-metrics-results/
  • 关键参数
  • --p-sampling-depth 10000：稀释深度（根据 7.1 稀释曲线确定）
  • 低于此深度的样本自动丢弃
  • 设太高 → 丢失样本多
  • 设太低 → 多样性估计不准
  • --i-phylogeny：有根系统发育树（第六章输出）
  • 一步输出 12 个文件
  • Alpha 多样性（4 个）
  • shannon_vector.qza：Shannon 指数
  • observed_features_vector.qza：观测 ASV 数
  • faith_pd_vector.qza：Faith 系统发育多样性
  • evenness_vector.qza：Pielou 均匀度
  • Beta 多样性距离矩阵（4 个）
  • bray_curtis_distance_matrix.qza：Bray-Curtis 距离
  • jaccard_distance_matrix.qza：Jaccard 距离
  • weighted_unifrac_distance_matrix.qza：加权 UniFrac 距离
  • unweighted_unifrac_distance_matrix.qza：非加权 UniFrac 距离
  • PCoA 排序结果（4 个）
  • 对应上述 4 种距离的 PCoA 坐标
  • 自动生成 Emperor 交互式 3D 可视化（.qzv）
  • ⚠️ 注意
  • 这是 16S 分析中最核心的一步命令
  • 稀释深度的选择直接影响所有下游结果
  • 稀释后的特征表（rarefied_table.qza）也在输出目录中

7.3 Alpha 多样性分析¶

• 四大指标解读
• Shannon 指数
  • 综合丰富度 + 均匀度
  • 值越大 = 多样性越高
  • 对稀有物种敏感
  • 典型范围：肠道 2–6，土壤 5–10
• Observed Features（观测 ASV 数）
  • 最直观的丰富度指标
  • 受测序深度影响最大
  • 稀释后比较才有意义
• Faith's PD（系统发育多样性）
  • 基于系统发育树的总枝长
  • 考虑了物种间的进化关系
  • 需要系统发育树
• Pielou's Evenness（均匀度）
  • 范围 0–1
  • 1 = 所有物种丰度完全相等
  • 0 = 极度不均匀（少数物种占绝对优势）
• 分组统计检验
• qiime diversity alpha-group-significance --i-alpha-diversity core-metrics-results/shannon_vector.qza --m-metadata-file metadata.tsv --o-visualization shannon-significance.qzv
• 对每个 alpha 指标分别运行
• 内置检验方法：Kruskal-Wallis（非参数检验）
  • 两组比较 = 等价于 Wilcoxon Rank-Sum
  • 三组及以上 = Kruskal-Wallis + 两两比较
• 显著性：p < 0.05
• 输出 .qzv 中包含箱线图 + p 值表
• 可视化
• 箱线图（Boxplot）：每组一个箱体
• 标注 p 值或显著性星号（* p<0.05, ** p<0.01, *** p<0.001）
• 每个点代表一个样本（jitter 散点叠加）
• ⚠️ 注意
• 16S 分析中稀释（rarefaction）是标准做法，与宏基因组不同
• 多重比较时需注意 FDR 校正
• 不同指标可能给出不同结论（这是正常的，各有侧重）

7.4 Beta 多样性分析¶

• 四种距离度量
• Bray-Curtis
  • 基于丰度（定量）
  • 不考虑系统发育关系
  • 最常用的非系统发育距离
• Jaccard
  • 基于有无（定性，presence/absence）
  • 不考虑丰度差异
  • 适合关注物种组成差异
• Weighted UniFrac
  • 基于丰度 + 系统发育关系
  • 对优势物种的丰度变化敏感
  • 需要系统发育树
• Unweighted UniFrac
  • 基于有无 + 系统发育关系
  • 对稀有物种的有无变化敏感
  • 需要系统发育树
• 选择建议
  • 关注优势物种变化 → Weighted UniFrac 或 Bray-Curtis
  • 关注稀有物种变化 → Unweighted UniFrac 或 Jaccard
  • 建议至少报告两种距离的结果
• PCoA 可视化
• core-metrics 已自动生成 Emperor 交互式 3D 图
• 也可用 R（ggplot2）绘制 2D PCoA 散点图
• 不同组用不同颜色/形状
• 添加 95% 置信椭圆
• 标注前两个轴的解释方差比例
• 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 --p-pairwise --p-permutations 999 --o-visualization bray-curtis-permanova.qzv
• 关键参数
  • --m-metadata-column Group：分组变量名
  • --p-method permanova：使用 PERMANOVA 检验
  • --p-pairwise：进行两两比较（三组及以上时必须加）
  • --p-permutations 999：置换次数（至少 999，发表建议 9999）
• 结果解读
  • p < 0.05 = 组间群落结构显著不同
  • pseudo-F 值越大 = 组间差异越大
  • R² = 分组变量解释的方差比例（> 0.1 有生物学意义）
• PERMDISP 方差齐性检验
• 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 permdisp --o-visualization bray-curtis-permdisp.qzv
• --p-method permdisp：切换为 PERMDISP 检验
• 结果解读
  • p > 0.05 = 方差齐性 ✓ → PERMANOVA 结果可信
  • p < 0.05 = 方差不齐 ⚠️ → PERMANOVA 显著可能是方差差异导致
• PERMANOVA + PERMDISP 必须配套报告
• ⚠️ 注意
• PERMANOVA 对不平衡设计（各组样本量差异大）敏感
• 有协变量时需在 R 中用 adonis2(dist ~ Covariate + Group)
• QIIME 2 内置的 beta-group-significance 不支持协变量校正
• 需要协变量校正时，导出距离矩阵到 R 中分析

八、差异分析¶

8.0 数据导出（QIIME 2 → R/Python）¶

• 导出特征表（qza → BIOM → TSV）
• qiime tools export --input-path table-final.qza --output-path exported/
• biom convert -i exported/feature-table.biom -o exported/feature-table.tsv --to-tsv
• 输出：feature-table.tsv（行=ASV，列=样本，值=计数）
• 导出分类注释
• qiime tools export --input-path taxonomy.qza --output-path exported/
• 输出：taxonomy.tsv（ASV ID + 分类学路径 + 置信度）
• 导出系统发育树
• qiime tools export --input-path rooted-tree.qza --output-path exported/
• 输出：tree.nwk（Newick 格式）
• 导出代表序列
• qiime tools export --input-path rep-seqs-filtered.qza --output-path exported/
• 输出：dna-sequences.fasta
• 在 R 中构建 phyloseq 对象
• library(phyloseq)
• otu <- otu_table(as.matrix(read.csv("feature-table.tsv", sep="\t", row.names=1, skip=1)), taxa_are_rows=TRUE)
• tax <- tax_table(as.matrix(tax_parsed))
• meta <- sample_data(read.csv("metadata.tsv", sep="\t", row.names=1))
• tree <- read_tree("tree.nwk")
• ps <- phyloseq(otu, tax, meta, tree)
• ⚠️ 注意
• BIOM 导出的 TSV 第一行是注释行（以 # 开头），读取时需 skip=1
• taxonomy.tsv 中的分类路径需要拆分为 Kingdom/Phylum/Class/Order/Family/Genus/Species 七列
• phyloseq 对象是 R 中 16S 下游分析的标准数据结构

8.1 ANCOM-BC2（QIIME 2 内置，推荐）¶

• 原理
• 估计并校正样本间的采样偏差（sampling fraction）
• 基于线性回归模型
• 专门为组分数据设计
• QIIME 2 命令
• qiime composition ancombc --i-table table-final.qza --m-metadata-file metadata.tsv --p-formula Group --p-p-adj-method BH --p-prv-cut 0.1 --p-lib-cut 1000 --p-reference-levels "Group::Control" --o-differentials ancombc-results.qza
• 关键参数
• --p-formula Group：模型公式（分组变量）
  • 多协变量：--p-formula "Group+Age+Sex"
• --p-p-adj-method BH：多重校正方法（Benjamini-Hochberg）
• --p-prv-cut 0.1：流行度过滤（至少 10% 样本中出现）
• --p-lib-cut 1000：文库大小过滤（总 reads < 1000 的样本移除）
• --p-reference-levels：参考组（对照组）
• 输出解读
• lfc（log fold change）> 0：在目标组中富集
• lfc < 0：在目标组中减少
• q_val < 0.05：显著差异
• diff_abn = TRUE：差异丰度特征
• ⚠️ 注意
• ANCOM-BC2 是 ANCOM 的改进版，推荐使用
• 支持协变量校正（优于 LEfSe）
• 对零膨胀数据有专门处理

8.2 ALDEx2（R 包）¶

• 优势：CLR 转换 + Monte Carlo 采样，专门处理组分数据
• 输入：原始计数矩阵（绝对不能用相对丰度！）
• R 代码
• library(ALDEx2)
• counts <- read.csv("feature-table.tsv", sep="\t", row.names=1, skip=1)
• conditions <- c(rep("Control",5), rep("Treatment",5))
• x <- aldex.clr(counts, conditions, mc.samples=128, denom="all")
• results <- aldex.ttest(x, paired.test=FALSE)
• results <- cbind(results, aldex.effect(x))
• aldex.plot(results, type="MA", test="welch")
• 关键参数
• mc.samples=128：Monte Carlo 采样次数（可增加到 1000 提高稳定性）
• denom="all"：使用所有特征计算 CLR（默认）
  • denom="iqlr"：使用四分位间距特征（更稳健）
• paired.test=FALSE：非配对检验
• 筛选标准
• we.eBH < 0.05：Welch's t-test BH 校正后 p 值
• wi.eBH < 0.05：Wilcoxon BH 校正后 p 值
• |effect| > 1：效应量绝对值 > 1
• 推荐取 we.eBH 和 wi.eBH 都显著的特征（更保守）
• 可视化
• MA 图：aldex.plot(results, type="MA")
• MW 图：aldex.plot(results, type="MW")
• 红色点 = 显著差异特征
• ⚠️ 注意
• 输入必须是整数计数矩阵（不能有小数）
• 不支持协变量校正（局限性）
• 适合两组比较；多组用 aldex.kw()
• 16S 数据零值多，CLR 转换前会自动加伪计数（0.5）

8.3 MaAsLin2（R 包，支持协变量）¶

• 优势：灵活的统计模型、支持协变量校正、支持纵向数据
• R 代码
• library(Maaslin2)
• Maaslin2(input_data = "feature-table.tsv", input_metadata = "metadata.tsv", output = "maaslin2_output", fixed_effects = c("Group"), random_effects = c("Subject"), normalization = "TSS", transform = "LOG", min_abundance = 0.0001, min_prevalence = 0.1, analysis_method = "LM")
• 关键参数
• normalization = "TSS"：总和缩放归一化
  • 替代选项：CLR、CSS、NONE
• transform = "LOG"：对数转换
  • 替代选项：LOGIT、AST、NONE
• min_abundance = 0.0001：最低丰度阈值
• min_prevalence = 0.1：最低流行度（至少 10% 样本中出现）
• analysis_method = "LM"：线性模型
  • 替代选项：CPLM（复合泊松）、NEGBIN（负二项）、ZINB（零膨胀负二项）
• fixed_effects：固定效应（分组 + 协变量如年龄、性别、BMI）
• random_effects：随机效应（如纵向研究中的个体 ID）
• 输出解读
• coefficient > 0：在目标组中富集
• coefficient < 0：在目标组中减少
• qval < 0.05：显著差异
• 每个显著特征自动生成散点图
• ⚠️ 注意
• 16S 数据零膨胀严重时，考虑用 ZINB 模型
• 不要对输入做归一化（TSS 由 MaAsLin2 内部完成）
• 样本量太小（< 10/组）时结果不稳定

8.4 LEfSe（经典方法）¶

• 原理：Kruskal-Wallis + Wilcoxon + LDA 效应量
• 输入：相对丰度表（LEfSe 特定格式）
• 在线工具：Galaxy / Huttenhower Lab 网站
• 命令行
• lefse_format_input.py input.txt formatted.in -c 1 -u 2 -o 1000000
• lefse_run.py formatted.in results.res -l 2.0 -a 0.05 -w 0.05
• lefse_plot_res.py results.res results.png --dpi 300
• lefse_plot_cladogram.py results.res cladogram.png --dpi 300
• 关键参数
• -l 2.0：LDA score 阈值（默认 2.0，更严格用 3.0）
• -a 0.05：Kruskal-Wallis p 值阈值
• -w 0.05：Wilcoxon p 值阈值
• 输出
• LDA 效应量柱状图（最经典的 16S 差异分析图）
• 进化分支图（Cladogram）：在分类树上标注差异物种
• ⚠️ 注意
• 每组最低 ≥ 5 个样本
• 不支持协变量校正
• 假阳性率较高，建议与其他方法交叉验证
• 已有 R 包替代：MicrobiomeMarker

8.5 差异分析最佳实践¶

• 至少使用 2 种方法交叉验证
• 推荐组合一：ANCOM-BC2 + ALDEx2
• 推荐组合二：MaAsLin2 + ALDEx2
• 推荐组合三：MaAsLin2 + LEfSe（审稿人常要求 LEfSe）
• 取交集结果 → 高置信差异特征
• 16S 特殊注意事项
• 16S 数据是组分数据（compositional data）：丰度之和恒为 1
• 零膨胀严重：大量 ASV 在多数样本中计数为 0
• 不要直接用 t-test / ANOVA / DESeq2（这些方法假设数据独立，不适合组分数据）
• DESeq2 在 16S 中假阳性率偏高（已有文献证实）
• 可视化
• 火山图：横轴 log2FC / coefficient，纵轴 -log10(q-value)
• 热图：差异 ASV/属在各样本中的丰度（Z-score 标准化）
• 箱线图：单个差异物种在各组中的分布
• LEfSe 柱状图 + 进化分支图

九、功能预测（PICRUSt2）¶

9.1 PICRUSt2 简介¶

• 全称：Phylogenetic Investigation of Communities by Reconstruction of Unobserved States 2
• 原理
• 基于 16S 序列推断物种 → 查找该物种的参考基因组 → 预测其基因含量
• 本质是"借用"已测序基因组的功能信息
• 是"预测"而非"测量"
• 适用场景
• 没有宏基因组数据时的功能分析替代方案
• 初步探索功能差异，指导后续实验设计
• 局限性（必须了解）
• 准确性取决于参考基因组覆盖度
• 人体肠道：预测较准（参考基因组多）
• 土壤/海洋/极端环境：预测不准（参考基因组少）
• 无法检测水平基因转移带来的功能
• 无法区分基因是否真正表达
• 发表时必须说明是"预测"功能

9.2 输入准备¶

• 从 QIIME 2 导出
• 导出特征表
  • qiime tools export --input-path table-final.qza --output-path picrust2_input/
  • 输出：feature-table.biom
• 导出代表序列
  • qiime tools export --input-path rep-seqs-filtered.qza --output-path picrust2_input/
  • 输出：dna-sequences.fasta
• ⚠️ 注意
• 输入必须是未稀释的原始计数表（不要用 rarefied table）
• 代表序列必须与特征表中的 ASV ID 一一对应

9.3 运行 PICRUSt2¶

• 一步完成命令
• picrust2_pipeline.py -s picrust2_input/dna-sequences.fasta -i picrust2_input/feature-table.biom -o picrust2_output/ -p 8 --stratified --per_sequence_contrib
• 关键参数
• -s：代表序列文件（FASTA）
• -i：特征表（BIOM 格式）
• -o：输出目录
• -p 8：线程数
• --stratified：输出分层结果（每个 ASV 对每个功能的贡献）
• --per_sequence_contrib：输出每条序列的功能贡献
• --min_align 0.8：最小比对覆盖度（默认 0.8）
• --max_nsti 2.0：最大 NSTI 阈值（默认 2.0）
• 输出文件
• EC_metagenome_out/pred_metagenome_unstrat.tsv.gz：EC 酶丰度表
• KO_metagenome_out/pred_metagenome_unstrat.tsv.gz：KO 丰度表
• pathways_out/path_abun_unstrat.tsv.gz：MetaCyc 通路丰度表
• marker_predicted_and_nsti.tsv.gz：NSTI 值

9.4 NSTI 质量评估¶

• NSTI = Nearest Sequenced Taxon Index
• 每个 ASV 与最近参考基因组的进化距离
• 值越小 = 预测越可靠
• 阈值
• NSTI < 0.03：预测非常可靠
• NSTI 0.03–0.15：预测较可靠
• NSTI 0.15–0.5：预测可信度一般
• NSTI > 0.5：预测不可靠
• NSTI > 2.0：默认被过滤掉
• 检查方法
• 查看 marker_predicted_and_nsti.tsv.gz
• 计算加权平均 NSTI（按 ASV 丰度加权）
• 如果大量高丰度 ASV 的 NSTI > 0.5 → 功能预测结果不可信
• ⚠️ 注意
• 必须在文章中报告 NSTI 分布
• 高 NSTI 的 ASV 可考虑从分析中移除

9.5 添加功能描述¶

• 命令
• add_descriptions.py -i KO_metagenome_out/pred_metagenome_unstrat.tsv.gz -m KO -o KO_with_descriptions.tsv
• add_descriptions.py -i EC_metagenome_out/pred_metagenome_unstrat.tsv.gz -m EC -o EC_with_descriptions.tsv
• add_descriptions.py -i pathways_out/path_abun_unstrat.tsv.gz -m METACYC -o pathways_with_descriptions.tsv
• -m：描述类型（KO / EC / METACYC）

9.6 功能差异分析¶

• 使用与物种差异分析相同的方法
• MaAsLin2 / ALDEx2 / ANCOM-BC2
• 输入：KO / EC / 通路丰度表
• 方法和参数与第八章完全一致
• KEGG 通路富集分析
• 将差异 KO 映射到 KEGG 通路
• 工具：clusterProfiler（R 包）或 KEGG Mapper 在线工具
• ⚠️ 注意
• 功能预测结果的差异分析同样需要多重校正
• 发表时必须注明"基于 PICRUSt2 预测"

十、可视化与结果整理¶

10.1 物种组成可视化¶

• 堆叠柱状图（Stacked Bar Plot）
• 最经典的 16S 结果图
• 门级（Phylum）：展示整体组成
• 属级（Genus）：展示详细组成（通常取 Top 15–20 属）
• 工具：R（ggplot2 + phyloseq）/ Python（matplotlib）
• 低丰度物种合并为"Others"
• 按分组排列样本
• 热图（Heatmap）
• 展示 Top 30–50 属/种在各样本中的丰度
• Z-score 标准化（按行）
• 行/列聚类（层次聚类）
• 工具：R（pheatmap / ComplexHeatmap）
• Krona 交互式饼图
• 多层级分类的交互式可视化
• 适合展示单个样本或组的分类组成

10.2 多样性可视化¶

• Alpha 多样性箱线图
• 每组一个箱体，叠加散点
• 标注统计检验 p 值
• 工具：R（ggplot2 + ggpubr）
• Beta 多样性 PCoA 散点图
• 2D 散点图 + 95% 置信椭圆
• 标注 PERMANOVA p 值和 R²
• 标注各轴解释方差比例
• 工具：R（ggplot2 + vegan）
• 稀释曲线（Rarefaction Curve）
• 展示测序深度是否足够
• 曲线趋于平台 = 深度足够

10.3 差异分析可视化¶

• 火山图（Volcano Plot）
• 横轴：log2FC 或 coefficient
• 纵轴：-log10(q-value)
• 红色 = 上调，蓝色 = 下调，灰色 = 不显著
• 标注 Top 差异物种名称
• 工具：R（ggplot2 / EnhancedVolcano）
• LEfSe 柱状图
• 横轴：LDA score
• 不同颜色代表不同组
• 经典的 16S 差异分析展示方式
• LEfSe 进化分支图（Cladogram）
• 在分类树上标注差异物种
• 从门到属逐层展示
• 差异物种箱线图
• 单个差异物种在各组中的丰度分布
• 标注 q 值
• Venn 图 / UpSet 图
• 展示不同差异分析方法的交集
• Venn 图：2–3 种方法取交集
• UpSet 图：> 3 种方法时更清晰
  • 工具：R（VennDiagram / UpSetR）
  • 交集特征 = 高置信差异物种

10.4 发表级图表规范¶

• 通用要求
• 分辨率 ≥ 300 dpi
• 格式：PDF（矢量图优先）/ TIFF / PNG
• 字体：Arial 或 Helvetica，≥ 8pt
• 配色：色盲友好配色方案（如 RColorBrewer 的 Set2、Paired）
• 图注（Figure Legend）必须独立可读
• 16S 文章必备图表（审稿人常要求）
• Figure 1：物种组成堆叠柱状图（门级 + 属级）
• Figure 2：Alpha 多样性箱线图（Shannon + Observed Features）
• Figure 3：Beta 多样性 PCoA 散点图（标注 PERMANOVA p 值和 R²）
• Figure 4：差异物种分析（火山图 / LEfSe 柱状图 / 热图）
• Supplementary：稀释曲线、NSTI 分布（如有 PICRUSt2）、完整差异物种列表
• 工具推荐
• R：ggplot2 + cowplot / patchwork（组合多图）
• R：ggsave() 导出高分辨率图
• 在线美化：Adobe Illustrator / Inkscape（调整布局和标注）

十一、常见问题与排错¶

11.1 DADA2 保留率过低（< 50%）¶

• 排查思路
• 检查原始数据质量（MultiQC 报告）
  • Q20 < 80% → 测序质量本身差，联系测序公司
• 检查截断长度是否合理
  • 截断太短 → 双端无法拼接（overlap < 12 bp）
  • 截断太长 → 保留了大量低质量碱基
  • 解决：重新查看质量分布图，调整 --p-trunc-len-f 和 --p-trunc-len-r
• 检查引物是否已去除
  • 未去除引物 → 引物序列被当作错误修剪掉 → reads 大量丢失
  • 解决：先用 Cutadapt 去引物，再跑 DADA2
• 检查 --p-max-ee（最大期望错误数）
  • 默认 2.0，过于严格时可放宽到 3.0 或 5.0
• 检查嵌合体去除比例
  • 嵌合体 > 30% → 可能 PCR 循环数过多（实验问题）
• 正常保留率参考
• 质控过滤：70–90%
• 去噪：60–80%
• 拼接：80–95%（取决于 overlap 长度）
• 嵌合体去除：95–99%
• 总保留率：50–80% 为正常范围

11.2 分类注释大量停留在属级/科级¶

• 原因
• 16S 基因本身分辨率有限（V3-V4 区无法区分很多近缘种）
• 这是 16S 方法的固有局限，不是分析错误
• 数据库中该类群的参考序列不足
• 应对
• 属级分析是 16S 的标准做法，不必强求种级
• 如需种级分辨率 → 考虑宏基因组测序
• 可尝试更换数据库（如 GTDB 对某些类群注释更好）
• 降低置信度阈值（如 0.5）可获得更多种级注释，但假阳性增加
• 在文章中如实报告注释深度分布

11.3 PERMANOVA 不显著（p > 0.05）¶

• 可能原因
• 组间差异确实不大（生物学真实结果）
• 样本量不足（每组 < 5 个样本时统计效力很低）
• 组内变异太大（个体差异掩盖了组间差异）
• 选择的距离度量不合适
• 应对
• 尝试不同距离度量（Bray-Curtis / Weighted UniFrac / Unweighted UniFrac）
  • 不同距离对不同类型的差异敏感
• 检查是否有混杂因素（年龄、性别、批次等）
  • 在 R 中用 adonis2(dist ~ Covariate + Group) 校正协变量
• 增加样本量（最有效的方法）
• 如果确实不显著 → 如实报告，不要 p-hacking
• 可以报告效应量（R²）即使 p 不显著

11.4 样本量不足的影响¶

• Alpha 多样性：检验效力低，难以检测到显著差异
• Beta 多样性：PERMANOVA 置换检验需要足够样本量
• 每组 < 3 个样本 → PERMANOVA 无法给出可靠 p 值
• 每组 5–10 个样本 → 最低要求
• 每组 ≥ 15 个样本 → 较理想
• 差异分析：假阴性率高，假阳性率也可能升高
• 建议
• 16S 研究最低每组 5 个生物学重复
• 理想每组 10–20 个
• 临床研究通常需要更多（≥ 30/组）

11.5 批次效应的识别与处理¶

• 识别
• PCoA 图中样本按批次（而非分组）聚类 → 存在批次效应
• PERMANOVA 检验批次变量：adonis2(dist ~ Batch + Group)
  • 如果 Batch 的 R² > Group 的 R² → 批次效应严重
• 处理
• 最佳方案：实验设计时每批次包含所有分组（平衡设计）
• 分析时校正：在统计模型中加入批次作为协变量
  • MaAsLin2：fixed_effects = c("Batch", "Group")
  • PERMANOVA：adonis2(dist ~ Batch + Group)
• 不推荐：ComBat 等批次校正方法在 16S 组分数据上效果不确定
• ⚠️ 注意
• 如果批次与分组完全混淆（如所有对照在批次 1，所有处理在批次 2）→ 无法区分批次效应和真实效应 → 实验设计失败

11.6 空白对照（Negative Control）中出现大量 reads¶

• 原因
• 试剂污染（DNA 提取试剂盒中的微生物 DNA，即"kitome"）
• 交叉污染（样本间串扰）
• 环境污染（实验操作中引入）
• 处理
• 使用 decontam（R 包）识别并去除污染 ASV
  • 频率法：在空白对照中高丰度但在真实样本中低丰度的 ASV
  • 流行度法：在空白对照中频繁出现的 ASV
• 命令示例
  • library(decontam)
  • contamdf <- isContaminant(ps, method="combined", neg="is.neg", conc="quant_reading")
  • ps.clean <- prune_taxa(!contamdf$contaminant, ps)
• 在文章中报告去污染步骤和去除的 ASV 数量
• ⚠️ 注意
• 低生物量样本（如皮肤、空气、组织）受试剂污染影响最大
• 必须设置空白对照（提取空白 + PCR 空白）
• 高生物量样本（如粪便）受影响较小

十二、发表清单与最佳实践¶

12.1 方法部分必须报告的信息¶

• 样本信息
• 样本类型、采集方法、保存条件
• 样本量（每组 n 数）
• 纳入/排除标准
• DNA 提取
• 试剂盒品牌和型号
• 是否有机械裂解（bead-beating）步骤
• PCR 扩增
• 引物序列（正向 + 反向）
• 扩增区域（如 V3-V4）
• PCR 循环数和退火温度
• 测序
• 平台和型号（如 Illumina NovaSeq 6000）
• 测序模式（如 PE250）
• 测序公司（如适用）
• 生信分析
• QIIME 2 版本号
• DADA2 参数（截断长度、最大期望错误等）
• 分类数据库及版本（如 SILVA 138.1、Greengenes2 2024.09）
• 分类器训练的区域（如 V3-V4 区域特异性分类器）
• 稀释深度
• 差异分析方法及参数
• 多重校正方法（如 BH / FDR）
• 显著性阈值

12.2 数据提交¶

• 原始测序数据
• 提交到 NCBI SRA 或 EBI ENA
• 获取 BioProject 编号（如 PRJNA123456）
• 大多数期刊要求在投稿时提供 accession number
• 提交内容
• 原始 FASTQ 文件（去除引物前或后均可，需说明）
• 样本元数据（metadata）
• BioSample 信息
• 分析代码
• 建议上传到 GitHub 或 Zenodo
• 包含完整的分析脚本和参数
• 提高研究的可重复性

12.3 16S 分析最佳实践总结¶

• 实验设计
• 每组 ≥ 5 个生物学重复（理想 ≥ 10）
• 包含空白对照（提取空白 + PCR 空白）
• 平衡批次设计（每批次包含所有分组）
• 质控
• 使用 DADA2 而非 OTU 聚类（ASV 分辨率更高、可重复性更好）
• 去除线粒体、叶绿体、未注释序列
• 使用 decontam 去除污染（低生物量样本必须做）
• 分析
• Alpha 多样性：稀释后比较，报告多个指标
• Beta 多样性：PERMANOVA + PERMDISP 配套报告
• 差异分析：至少两种方法交叉验证，使用组分数据方法（ANCOM-BC2 / ALDEx2）
• 功能预测：注明是 PICRUSt2 预测，报告 NSTI
• 可视化
• 组成柱状图 + Alpha 箱线图 + PCoA + 差异分析图
• 所有图 ≥ 300 dpi，色盲友好配色
• 报告
• 如实报告所有参数和阈值
• 不显著的结果也要报告
• 提交原始数据和分析代码

12.4 推荐阅读与资源¶

• 官方教程
• QIIME 2 官方教程：最权威的扩增子分析教程
• DADA2 Tutorial：DADA2 R 版教程
• Microbiome Helper：宏基因组/扩增子分析实用工具集
• 核心文献
• Callahan et al., 2016 (Nature Methods)：DADA2 原始论文
• Bolyen et al., 2019 (Nature Biotechnology)：QIIME 2 平台论文
• Segata et al., 2011 (Genome Biology)：LEfSe 原始论文
• Anderson, 2001 (Austral Ecology)：PERMANOVA 方法论文
• Gloor et al., 2017 (Frontiers in Microbiology)：组分数据分析综述（理解为什么不能直接用 t-test）
• Nearing et al., 2022 (Nature Communications)：差异丰度方法系统比较（推荐 ANCOM-BC / ALDEx2）
• R 包生态
• phyloseq：16S 数据的标准 R 数据结构
• microbiome：多样性分析扩展
• vegan：群落生态学经典 R 包（PERMANOVA、NMDS 等）
• ggplot2 + ggsci + ggpubr：发表级可视化