微生物组代谢预测（PICRUSt2 / Tax4Fun / MIMOSA）¶

# QIIME2 DADA2去噪
qiime dada2 denoise-paired \
  --i-demultiplexed-seqs reads.qza \
  --p-trunc-len-f 250 --p-trunc-len-r 200 \
  --p-n-threads 16 \
  --o-table table.qza \
  --o-representative-sequences rep-seqs.qza \
  --o-denoising-stats stats.qza

# 去除线粒体和叶绿体
qiime taxa filter-table \
  --i-table table.qza \
  --i-taxonomy taxonomy.qza \
  --p-exclude mitochondria,chloroplast \
  --o-filtered-table table_filtered.qza

# 导出为BIOM格式（PICRUSt2输入）
qiime tools export --input-path table_filtered.qza --output-path exported/
qiime tools export --input-path rep-seqs.qza --output-path exported/
# 得到 exported/feature-table.biom 和 exported/dna-sequences.fasta

# PICRUSt2内置的序列放置（通常不需要单独运行）
place_seqs.py \
  -s exported/dna-sequences.fasta \
  -o placed_seqs.tre \
  -p 16 \
  --intermediate placement_working

# QIIME2方式（使用q2-picrust2插件）
qiime fragment-insertion sepp \
  --i-representative-sequences rep-seqs.qza \
  --i-reference-database sepp-refs-silva-128.qza \
  --p-threads 16 \
  --o-tree insertion-tree.qza \
  --o-placements insertion-placements.qza

# 基因家族预测（通常作为picrust2_pipeline的一部分自动运行）
hsp.py \
  -i 16S \
  -t placed_seqs.tre \
  -o KO_predicted.tsv.gz \
  -p 16 \
  -n \
  --observed_trait_table /path/to/picrust2/default_files/prokaryotic/KO.tsv.gz

# EC number预测
hsp.py \
  -i 16S \
  -t placed_seqs.tre \
  -o EC_predicted.tsv.gz \
  -p 16 \
  -n \
  --observed_trait_table /path/to/picrust2/default_files/prokaryotic/EC.tsv.gz

# 宏基因组预测
metagenome_pipeline.py \
  -i exported/feature-table.biom \
  -m KO_predicted.tsv.gz \
  -o metagenome_out \
  --strat_out \
  --max_nsti 2

# 输出文件：
# metagenome_out/pred_metagenome_unstrat.tsv.gz  (非分层，样本×基因)
# metagenome_out/pred_metagenome_strat.tsv.gz    (分层，样本×基因×ASV贡献)
# metagenome_out/seqtab_norm.tsv.gz              (标准化后的ASV表)
# metagenome_out/weighted_nsti.tsv.gz            (加权NSTI)

# 通路丰度推断
pathway_pipeline.py \
  -i metagenome_out/pred_metagenome_strat.tsv.gz \
  -o pathways_out \
  --intermediate pathway_working \
  -p 16

# 输出文件：
# pathways_out/path_abun_unstrat.tsv.gz    (通路丰度，非分层)
# pathways_out/path_abun_strat.tsv.gz      (通路丰度，分层)
# pathways_out/path_cov_unstrat.tsv.gz     (通路覆盖度)

# Tax4Fun2 R包使用
library(Tax4Fun2)

# 1. 构建参考数据库（只需一次）
buildReferenceData(
  path_to_tax4fun2_data = "Tax4Fun2_ReferenceData_v2",
  path_to_user_data = "user_reference_genomes/",  # 可选：自定义参考基因组
  path_to_SILVA = "SILVA_138_SSURef_NR99_tax_silva.fasta",
  num_threads = 16
)

# 2. 运行功能预测
# 输入：ASV的FASTA文件和丰度表
result <- runRefBlast(
  path_to_otus = "rep_seqs.fasta",
  path_to_reference_data = "Tax4Fun2_ReferenceData_v2",
  path_to_temp_folder = "tax4fun2_temp",
  database_mode = "Ref99NR",
  num_threads = 16
)

# 预测功能丰度
predictions <- makeFunctionalPrediction(
  path_to_otu_table = "asv_table.txt",  # 样本×ASV丰度表
  path_to_reference_data = "Tax4Fun2_ReferenceData_v2",
  path_to_temp_folder = "tax4fun2_temp",
  min_identity_to_reference = 0.97,
  normalize_by_copy_number = TRUE
)

# 提取KO丰度和通路丰度
ko_table <- predictions$KO_table
pathway_table <- predictions$pathway_table

write.table(ko_table, "tax4fun2_KO_predictions.tsv", sep="\t", quote=FALSE)
write.table(pathway_table, "tax4fun2_pathway_predictions.tsv", sep="\t", quote=FALSE)

library(ALDEx2)

# 读取功能预测表
ko_table <- read.table("metagenome_out/pred_metagenome_unstrat.tsv.gz",
                       header=TRUE, row.names=1, sep="\t")

# 样本分组信息
groups <- c(rep("Control", 10), rep("Treatment", 10))

# ALDEx2差异分析
aldex_result <- aldex(
  reads = round(ko_table),  # ALDEx2需要整数
  conditions = groups,
  mc.samples = 128,
  test = "t",
  effect = TRUE,
  include.sample.summary = FALSE
)

# 筛选显著差异功能
sig_kos <- aldex_result[aldex_result$wi.eBH < 0.05 & abs(aldex_result$effect) > 1, ]
sig_kos <- sig_kos[order(sig_kos$wi.eBH), ]

cat("显著差异KO数:", nrow(sig_kos), "\n")

# MIMOSA2分析
library(mimosa2)

# 输入数据准备
# 1. 微生物丰度表（属/种水平）
species_table <- read.table("species_abundance.tsv", header=TRUE, row.names=1)

# 2. 代谢物浓度表
metabolite_table <- read.table("metabolite_concentrations.tsv", header=TRUE, row.names=1)

# 3. 运行MIMOSA
mimosa_result <- run_mimosa(
  species_data = species_table,
  metabolite_data = metabolite_table,
  kegg_data = "kegg_reaction_data.rds",  # KEGG反应数据
  num_permutations = 1000
)

# 4. 查看显著关联
significant_pairs <- mimosa_result$results %>%
  filter(p_adjusted < 0.05) %>%
  arrange(p_adjusted)

# 5. 特定代谢物的贡献分解
butyrate_contributors <- get_contributors(
  mimosa_result,
  metabolite = "C00246"  # KEGG ID for butyrate
)

library(pathview)

# 准备数据：KO编号为行名，效应量为值
ko_data <- sig_kos$effect
names(ko_data) <- rownames(sig_kos)

# 可视化特定通路（如丁酸代谢 map00650）
pathview(
  gene.data = ko_data,
  pathway.id = "00650",
  species = "ko",
  out.suffix = "butyrate_metabolism",
  kegg.native = TRUE,
  limit = list(gene=c(-2, 2)),
  low = list(gene="blue"),
  high = list(gene="red")
)

# 批量可视化多个通路
pathways_of_interest <- c("00010", "00020", "00030", "00650", "00680")
for (pw in pathways_of_interest) {
  pathview(gene.data=ko_data, pathway.id=pw, species="ko",
           out.suffix=paste0("pathway_", pw))
}

# 验证：比较PICRUSt2预测与宏基因组实测
predicted <- read.table("pred_metagenome_unstrat.tsv.gz", header=TRUE, row.names=1)
observed <- read.table("actual_metagenome_ko.tsv", header=TRUE, row.names=1)

# 对共同样本和共同KO计算相关
common_samples <- intersect(colnames(predicted), colnames(observed))
common_kos <- intersect(rownames(predicted), rownames(observed))

correlations <- sapply(common_samples, function(s) {
  cor(predicted[common_kos, s], observed[common_kos, s], method="spearman")
})

cat("预测-实测相关 (Spearman):\n")
cat("  中位数:", median(correlations), "\n")
cat("  范围:", range(correlations), "\n")

# NSTI分布
nsti <- read.table("metagenome_out/weighted_nsti.tsv.gz", header=TRUE)
cat("\n加权NSTI:\n")
cat("  中位数:", median(nsti$weighted_NSTI), "\n")
cat("  >2的样本:", sum(nsti$weighted_NSTI > 2), "/", nrow(nsti), "\n")

#!/bin/bash
# ====================================================
# 微生物组代谢功能预测完整流程 (PICRUSt2)
# ====================================================

set -euo pipefail
THREADS=16
OUTDIR="functional_prediction"
mkdir -p ${OUTDIR}

# ---- PICRUSt2一站式流程 ----
echo "=== PICRUSt2完整流程 ==="
picrust2_pipeline.py \
  -s exported/dna-sequences.fasta \
  -i exported/feature-table.biom \
  -o ${OUTDIR}/picrust2 \
  -p ${THREADS} \
  --stratified \
  --per_sequence_contrib \
  --max_nsti 2 \
  --verbose

echo "PICRUSt2输出文件:"
ls -la ${OUTDIR}/picrust2/

# 输出说明：
# KO_metagenome_out/     - KEGG Orthology预测
# EC_metagenome_out/     - EC number预测
# pathways_out/          - MetaCyc通路预测

# ---- 添加功能描述 ----
echo "=== 添加功能描述 ==="
add_descriptions.py \
  -i ${OUTDIR}/picrust2/KO_metagenome_out/pred_metagenome_unstrat.tsv.gz \
  -m KO \
  -o ${OUTDIR}/picrust2/KO_with_descriptions.tsv.gz

add_descriptions.py \
  -i ${OUTDIR}/picrust2/pathways_out/path_abun_unstrat.tsv.gz \
  -m METACYC \
  -o ${OUTDIR}/picrust2/pathways_with_descriptions.tsv.gz

echo "=== 完成 ==="

#!/usr/bin/env Rscript
# ====================================================
# 功能预测结果分析与可视化
# ====================================================

library(tidyverse)
library(ALDEx2)
library(ggplot2)
library(pheatmap)

# ---- 1. 读取PICRUSt2结果 ----
ko_table <- read.table("functional_prediction/picrust2/KO_with_descriptions.tsv.gz",
                       header=TRUE, sep="\t", row.names=1, comment.char="")
pathway_table <- read.table("functional_prediction/picrust2/pathways_with_descriptions.tsv.gz",
                           header=TRUE, sep="\t", row.names=1, comment.char="")

# 分离描述列
ko_desc <- ko_table$description
ko_table$description <- NULL
pathway_desc <- pathway_table$description
pathway_table$description <- NULL

# ---- 2. 差异功能分析 (ALDEx2) ----
metadata <- read.table("metadata.tsv", header=TRUE, sep="\t")
groups <- metadata$group[match(colnames(ko_table), metadata$sample_id)]

# KO水平差异分析
aldex_ko <- aldex(
  reads = round(as.data.frame(ko_table)),
  conditions = groups,
  mc.samples = 128,
  test = "t",
  effect = TRUE
)

sig_kos <- aldex_ko %>%
  rownames_to_column("KO") %>%
  filter(wi.eBH < 0.05, abs(effect) > 0.5) %>%
  arrange(wi.eBH)

cat(sprintf("显著差异KO: %d\n", nrow(sig_kos)))

# 通路水平差异分析
aldex_pw <- aldex(
  reads = round(as.data.frame(pathway_table)),
  conditions = groups,
  mc.samples = 128,
  test = "t",
  effect = TRUE
)

sig_pathways <- aldex_pw %>%
  rownames_to_column("pathway") %>%
  filter(wi.eBH < 0.05, abs(effect) > 0.5) %>%
  arrange(wi.eBH)

cat(sprintf("显著差异通路: %d\n", nrow(sig_pathways)))

# ---- 3. 可视化 ----
# 差异通路柱状图
top_pathways <- sig_pathways %>% head(20)
top_pathways$description <- pathway_desc[match(top_pathways$pathway, rownames(pathway_table))]

ggplot(top_pathways, aes(x=reorder(description, effect), y=effect, fill=effect>0)) +
  geom_col() +
  coord_flip() +
  scale_fill_manual(values=c("TRUE"="coral", "FALSE"="steelblue"),
                    labels=c("TRUE"="Treatment ↑", "FALSE"="Control ↑")) +
  labs(x="", y="Effect Size (ALDEx2)", title="显著差异代谢通路 (Top 20)") +
  theme_minimal(base_size=12) +
  theme(legend.title=element_blank())

ggsave("functional_prediction/diff_pathways.pdf", width=10, height=8)

# NSTI分布图
nsti <- read.table("functional_prediction/picrust2/KO_metagenome_out/weighted_nsti.tsv.gz",
                   header=TRUE)
ggplot(nsti, aes(x=weighted_NSTI)) +
  geom_histogram(bins=30, fill="steelblue", color="white") +
  geom_vline(xintercept=2, color="red", linetype="dashed") +
  labs(x="Weighted NSTI", y="样本数", title="预测可靠性评估 (NSTI分布)") +
  annotate("text", x=2.1, y=Inf, label="NSTI=2 阈值", vjust=2, color="red") +
  theme_minimal()

ggsave("functional_prediction/nsti_distribution.pdf", width=7, height=5)