微生物多样性分析（Microbial Diversity Analysis）¶

原始数据 (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图

# ============================================================
# QIIME 2 16S 扩增子分析流程
# 环境：conda activate qiime2-amplicon-2026.1
# （注意：QIIME 2 版本号格式为 年份.月份，最新稳定版为 2026.1）
# （2025.4 起文档站点迁移到 https://amplicon-docs.qiime2.org）
# ============================================================

# ---------- 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.2-99-nb-classifier-V3V4.qza \  # SILVA 138.2 V3-V4 分类器（138.2 为当前最新维护版）
  --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

# ============================================================
# 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：无差异

# ============================================================
# 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 列：该物种在哪个组中富集（空 = 不显著）

# ============================================================
# 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 画图