181_微生物组功能通路差异¶

# 方案1：HUMAnN3（宏基因组数据）
# 安装
conda install -c bioconda humann

# 下载数据库
humann_databases --download chocophlan full /path/to/databases
humann_databases --download uniref uniref90_diamond /path/to/databases

# 运行HUMAnN3
humann --input sample.fastq.gz \
       --output humann_output \
       --threads 16 \
       --nucleotide-database /path/to/databases/chocophlan \
       --protein-database /path/to/databases/uniref

# 输出文件
# sample_genefamilies.tsv     # 基因家族丰度
# sample_pathabundance.tsv    # 通路丰度
# sample_pathcoverage.tsv     # 通路覆盖度

# 合并多个样本
humann_join_tables -i humann_output/ -o all_pathabundance.tsv --file_name pathabundance

# 标准化（推荐CPM）
humann_renorm_table -i all_pathabundance.tsv -o all_pathabundance_cpm.tsv -u cpm

# 分离物种特异和社区整体通路
humann_split_stratified_table -i all_pathabundance_cpm.tsv -o stratified/
# stratified/all_pathabundance_cpm_stratified.tsv   # 按物种分
# stratified/all_pathabundance_cpm_unstratified.tsv # 总体通路

# 方案2：PICRUSt2（16S数据）
# 安装
conda install -c bioconda -c conda-forge picrust2

# 从BIOM文件运行
picrust2_pipeline.py \
    -s study_seqs.fna \
    -i study_seqs.biom \
    -o picrust2_output \
    -p 16 \
    --stratified \
    --per_sequence_contrib

# 输出
# picrust2_output/pathways_out/path_abun_unstrat.tsv  # MetaCyc通路丰度
# picrust2_output/KO_metagenome_out/pred_metagenome_unstrat.tsv  # KO丰度
# picrust2_output/EC_metagenome_out/pred_metagenome_unstrat.tsv  # EC丰度

library(Maaslin2)

# 读取通路丰度数据（样本为行，通路为列）
pathway_data <- read.table(
  "stratified/all_pathabundance_cpm_unstratified.tsv",
  header = TRUE, sep = "\t", row.names = 1,
  check.names = FALSE, comment.char = ""
)
pathway_data <- t(pathway_data)  # 转置为样本×通路

# 元数据
metadata <- data.frame(
  Sample = rownames(pathway_data),
  Group = c(rep("Control", 10), rep("Treatment", 10)),
  Age = rnorm(20, 50, 10),
  Sex = sample(c("M", "F"), 20, replace = TRUE),
  row.names = rownames(pathway_data)
)

# 运行MaAsLin2
results <- Maaslin2(
  input_data = as.data.frame(pathway_data),
  input_metadata = metadata,
  output = "maaslin2_output",
  fixed_effects = c("Group", "Age", "Sex"),
  reference = c("Group,Control"),
  normalization = "TSS",
  transform = "LOG",
  analysis_method = "LM",
  max_significance = 0.05,
  min_abundance = 0.0001,
  min_prevalence = 0.1,
  cores = 4
)

# 查看显著通路
sig_pathways <- results$results[results$results$qval < 0.05, ]
sig_pathways <- sig_pathways[order(sig_pathways$qval), ]
print(sig_pathways[, c("feature", "metadata", "coef", "pval", "qval")])

library(ALDEx2)

# 准备数据（需要count或整数数据）
# 如果是CPM数据，需要转回近似counts
pathway_counts <- round(pathway_data * 1000)  # 近似转换

# 分组
groups <- metadata$Group

# 运行ALDEx2
aldex_result <- aldex(
  pathway_counts,
  conditions = groups,
  mc.samples = 128,      # Monte Carlo采样次数
  test = "t",             # 使用Welch's t检验
  effect = TRUE,          # 计算effect size
  denom = "all"           # 使用全部特征作为参考
)

# 查看显著结果
sig <- aldex_result[aldex_result$we.eBH < 0.05, ]
sig <- sig[order(sig$effect, decreasing = TRUE), ]

# 效应量图
aldex.plot(aldex_result, type = "MA", test = "welch")
aldex.plot(aldex_result, type = "MW", test = "welch")

library(ggplot2)
library(dplyr)
library(pheatmap)

# 1. 条形图（Top差异通路）
top_pathways <- sig_pathways %>%
  head(20) %>%
  mutate(direction = ifelse(coef > 0, "Enriched in Treatment", "Enriched in Control"))

ggplot(top_pathways, aes(x = reorder(feature, coef), y = coef, fill = direction)) +
  geom_bar(stat = "identity") +
  coord_flip() +
  theme_minimal() +
  labs(x = "Pathway", y = "Coefficient", title = "Differential Pathways") +
  scale_fill_manual(values = c("blue", "red"))

# 2. 热图
sig_pathway_names <- head(sig_pathways$feature, 30)
heatmap_data <- pathway_data[, sig_pathway_names]

# 标准化用于热图
heatmap_scaled <- scale(heatmap_data)

annotation_row <- data.frame(
  Group = metadata$Group,
  row.names = rownames(metadata)
)

pheatmap(t(heatmap_scaled),
         annotation_col = annotation_row,
         clustering_distance_rows = "correlation",
         clustering_distance_cols = "correlation",
         show_colnames = FALSE,
         fontsize_row = 8,
         main = "Differential Pathways Heatmap")

# 3. 火山图
library(EnhancedVolcano)
volcano_data <- results$results[results$results$metadata == "Group", ]
EnhancedVolcano(volcano_data,
                lab = volcano_data$feature,
                x = 'coef',
                y = 'qval',
                pCutoff = 0.05,
                FCcutoff = 0.5,
                title = 'Pathway Differential Analysis')

# 将MetaCyc通路按上层分类汇总
# HUMAnN3提供通路重新分组功能

# 命令行：将通路映射到更高层次
# humann_regroup_table -i pathabundance.tsv -g uniref90_go -o go_grouped.tsv
# humann_regroup_table -i genefamilies.tsv -g uniref90_ko -o ko_grouped.tsv

# R中手动分类
# MetaCyc通路有层次结构，可以按超级通路分组
pathway_hierarchy <- data.frame(
  pathway = c("PWY-5100", "PWY-6471", "GLYCOLYSIS"),
  superclass = c("Amino Acid Biosynthesis", "Nucleotide Biosynthesis", "Energy Metabolism")
)

# 按超级分类汇总差异通路
superclass_summary <- sig_pathways %>%
  left_join(pathway_hierarchy, by = c("feature" = "pathway")) %>%
  group_by(superclass) %>%
  summarise(n_pathways = n(), mean_effect = mean(abs(coef)))

# HUMAnN3完整流程
for sample in *.fastq.gz; do
    humann --input $sample --output humann_out/ --threads 8
done
humann_join_tables -i humann_out/ -o joined_pathways.tsv --file_name pathabundance
humann_renorm_table -i joined_pathways.tsv -o pathways_cpm.tsv -u cpm
humann_split_stratified_table -i pathways_cpm.tsv -o split/

# PICRUSt2快速流程
picrust2_pipeline.py -s seqs.fna -i table.biom -o picrust2_out -p 16

# KEGG通路映射
humann_regroup_table -i genefamilies.tsv -g uniref90_ko -o ko.tsv