602 肠道菌群与代谢组联合分析¶

早期整合（Early Integration）：先把两套数据拼成一张大表，再做分析
├── 优点：简单直接
├── 缺点：特征数爆炸，噪声大
└── 适用：样本量大、数据维度相近

中期整合（Intermediate Integration）：各自降维后在潜在空间中关联  ← 最推荐
├── 优点：保留各自数据结构
├── 缺点：方法选择复杂
├── 代表：CCA、MOFA、MintTea
└── 适用：大多数研究

晚期整合（Late Integration）：各自分析完再比较结果
├── 优点：各自流程独立
├── 缺点：可能丢失跨组学信号
└── 适用：快速探索

# ========== 16S扩增子数据 ==========
# 用QIIME2处理得到特征表（ASV表）
qiime dada2 denoise-paired \
    --i-demultiplexed-seqs paired-end-demux.qza \  # 输入：去多路复用的双端测序数据
    --p-trunc-len-f 240 \                          # 正向读段截断长度
    --p-trunc-len-r 200 \                          # 反向读段截断长度
    --o-table table.qza \                          # 输出：特征表（ASV丰度矩阵）
    --o-representative-sequences rep-seqs.qza      # 输出：代表序列

# 导出为BIOM格式，方便后续R处理
qiime tools export \
    --input-path table.qza \      # 输入：QIIME2特征表
    --output-path exported-table   # 输出目录

# 转换为TSV文本格式
biom convert \
    -i exported-table/feature-table.biom \  # 输入：BIOM格式
    -o asv_table.tsv \                      # 输出：TSV文本格式
    --to-tsv                                # 指定输出格式为TSV

# ========== 宏基因组数据 ==========
# 用HUMAnN3做功能注释，得到通路丰度表
humann \
    --input merged_reads.fastq.gz \           # 输入：质控后的合并读段
    --output humann_output/ \                 # 输出目录
    --threads 16 \                            # 使用16个线程
    --nucleotide-database chocophlan_db \     # 核酸比对数据库
    --protein-database uniref90_db            # 蛋白比对数据库

# 合并多样本通路丰度表
humann_join_tables \
    --input humann_output/ \      # 输入：各样本的通路丰度文件
    --output merged_pathabundance.tsv \  # 输出：合并后的通路丰度表
    --file_name pathabundance           # 匹配文件名中包含pathabundance的文件

# ========== R语言：代谢组数据预处理 ==========

# 加载必要的包
library(metabolomics)  # 代谢组分析包
library(MetaboAnalystR)  # MetaboAnalyst的R版本

# 读取代谢组数据（通常是峰面积矩阵）
metab_data <- read.csv("metabolomics_data.csv",  # 读取CSV格式的代谢组数据
                        row.names = 1)             # 第一列作为行名（样本名）

# 缺失值处理：用最小值的1/5填充（常用策略）
metab_data[metab_data == 0] <- NA  # 先把0替换为NA
for(col in colnames(metab_data)) {  # 遍历每一列（每个代谢物）
  min_val <- min(metab_data[[col]], na.rm = TRUE)  # 计算该代谢物的最小值
  metab_data[[col]][is.na(metab_data[[col]])] <- min_val / 5  # 用最小值/5填充NA
}

# 数据标准化：log2转换 + Pareto缩放
metab_log <- log2(metab_data)  # log2转换，减小数值差异

# Pareto缩放（除以标准差的平方根）
metab_scaled <- scale(metab_log,           # 输入：log转换后的数据
                      center = TRUE,        # 中心化（减均值）
                      scale = apply(metab_log, 2, function(x) sqrt(sd(x))))  # 除以SD的平方根

# 检查数据维度
cat("微生物组样本数:", nrow(asv_table), "\n")  # 打印微生物组样本数
cat("代谢组样本数:", nrow(metab_scaled), "\n")  # 打印代谢组样本数
# 两者的样本数和样本名必须一致！

# ========== 确保两套数据样本完全匹配 ==========
# 这一步非常关键，是新手最容易忽略的

# 找到两个数据集的共同样本
common_samples <- intersect(rownames(asv_table),   # 微生物组的样本名
                            rownames(metab_scaled))  # 代谢组的样本名

cat("共同样本数:", length(common_samples), "\n")  # 打印共同样本数

# 只保留共同样本，并按相同顺序排列
asv_matched <- asv_table[common_samples, ]       # 微生物组取共同样本
metab_matched <- metab_scaled[common_samples, ]  # 代谢组取共同样本

# 过滤低丰度特征（微生物组）
# 至少在20%的样本中出现的ASV才保留
prevalence_threshold <- 0.2  # 普遍性阈值
keep_asv <- colSums(asv_matched > 0) / nrow(asv_matched) >= prevalence_threshold
asv_filtered <- asv_matched[, keep_asv]  # 只保留高普遍性的ASV

cat("过滤后ASV数:", ncol(asv_filtered), "\n")  # 打印过滤后的ASV数量

# ========== Procrustes分析：两套数据的整体相关性 ==========
library(vegan)  # 群落生态学包

# 分别做PCoA降维
# 微生物组用Bray-Curtis距离
asv_dist <- vegdist(asv_filtered, method = "bray")  # 计算Bray-Curtis距离矩阵
asv_pcoa <- cmdscale(asv_dist, k = 2)               # PCoA降维到2维

# 代谢组用Euclidean距离
metab_dist <- vegdist(metab_matched, method = "euclidean")  # 计算欧氏距离矩阵
metab_pcoa <- cmdscale(metab_dist, k = 2)                   # PCoA降维到2维

# Procrustes分析
pro <- protest(asv_pcoa, metab_pcoa,  # 输入两个PCoA结果
               permutations = 999)      # 置换检验999次

# 查看结果
cat("M2 统计量:", pro$ss, "\n")           # M2越小，两套数据越相关
cat("相关系数:", sqrt(1 - pro$ss), "\n")  # 转换为相关系数
cat("P值:", pro$signif, "\n")             # P<0.05表示显著相关

# 可视化Procrustes结果
plot(pro,                          # 绘制Procrustes图
     main = "Microbiome vs Metabolome",  # 标题
     ar.col = "blue",              # 箭头颜色
     xlab = "Dimension 1",        # X轴标签
     ylab = "Dimension 2")        # Y轴标签

# ========== Mantel检验：两个距离矩阵的全局相关性 ==========

# 计算距离矩阵
micro_dist <- vegdist(asv_filtered, method = "bray")       # 微生物组Bray-Curtis距离
metab_dist <- vegdist(metab_matched, method = "euclidean")  # 代谢组欧氏距离

# Mantel检验
mantel_result <- mantel(micro_dist, metab_dist,  # 输入两个距离矩阵
                        method = "spearman",       # 使用Spearman相关
                        permutations = 9999)       # 置换检验9999次

# 输出结果
cat("Mantel统计量 r:", mantel_result$statistic, "\n")  # r值：相关性强度
cat("P值:", mantel_result$signif, "\n")                # P<0.05为显著
# r > 0.25 且 P < 0.05 通常认为有意义的关联

# ========== 微生物-代谢物相关网络 ==========
library(Hmisc)     # 用于计算带P值的相关矩阵
library(corrplot)  # 相关性可视化

# 合并数据做相关分析
combined <- cbind(asv_filtered, metab_matched)  # 列合并微生物和代谢物

# 计算Spearman相关矩阵（带P值校正）
cor_result <- rcorr(as.matrix(asv_filtered),  # 微生物数据
                    as.matrix(metab_matched),  # 代谢物数据
                    type = "spearman")          # Spearman秩相关

# 提取微生物-代谢物之间的相关系数
n_micro <- ncol(asv_filtered)   # 微生物特征数
n_metab <- ncol(metab_matched)  # 代谢物特征数
cor_matrix <- cor_result$r[1:n_micro, (n_micro+1):(n_micro+n_metab)]  # 取交叉部分
pval_matrix <- cor_result$P[1:n_micro, (n_micro+1):(n_micro+n_metab)]  # 取P值交叉部分

# BH校正P值
pval_adj <- matrix(p.adjust(pval_matrix, method = "BH"),  # BH多重检验校正
                   nrow = nrow(pval_matrix))

# 筛选显著相关（|r| > 0.5 且 调整后P < 0.05）
sig_pairs <- which(abs(cor_matrix) > 0.5 & pval_adj < 0.05,  # 筛选条件
                   arr.ind = TRUE)                              # 返回行列索引

cat("显著相关对数:", nrow(sig_pairs), "\n")  # 打印显著相关的微生物-代谢物对数

# 热图可视化（只展示显著的）
library(pheatmap)
# 筛选在至少1个显著对中出现的微生物和代谢物
sig_micro <- unique(sig_pairs[, 1])  # 有显著关联的微生物
sig_metab <- unique(sig_pairs[, 2])  # 有显著关联的代谢物

pheatmap(cor_matrix[sig_micro, sig_metab],  # 只画显著部分
         color = colorRampPalette(c("blue", "white", "red"))(100),  # 颜色梯度
         main = "Microbe-Metabolite Correlations",  # 标题
         fontsize_row = 8,    # 行标签字号
         fontsize_col = 8)    # 列标签字号

# ========== MMVEC：用神经网络学习微生物-代谢物共现概率 ==========
# MMVEC比简单相关更适合组成型数据，因为它基于条件概率而非相关系数

# 安装MMVEC
pip install mmvec  # 安装MMVEC包

# 运行MMVEC
mmvec paired-omics \
    --microbe-file asv_table.biom \       # 输入：微生物丰度表（BIOM格式）
    --metabolite-file metab_table.biom \  # 输入：代谢物丰度表（BIOM格式）
    --summary-dir mmvec_output/ \         # 输出目录
    --epochs 1000 \                       # 训练轮数
    --batch-size 32 \                     # 批次大小
    --latent-dim 3                        # 潜在维度数（降维到3维）

# 输出解读：
# ranks.tsv  — 微生物对代谢物的条件概率排序（核心结果）
# biplot.qza — 可以导入QIIME2做可视化的双标图

# ========== MOFA2：多组学因子分析 ==========
# MOFA把多组学数据分解为少数几个"因子"，每个因子代表一种跨组学的变异模式

library(MOFA2)  # 加载MOFA2包

# 准备输入数据（列表格式）
data_list <- list(
  microbiome = t(asv_filtered),    # 微生物组（转置：特征在行，样本在列）
  metabolome = t(metab_matched)    # 代谢组（转置：特征在行，样本在列）
)

# 创建MOFA对象
mofa <- create_mofa(data_list)  # 从列表创建MOFA对象

# 设置参数
data_opts <- get_default_data_options(mofa)    # 获取默认数据选项
model_opts <- get_default_model_options(mofa)  # 获取默认模型选项
model_opts$num_factors <- 10                    # 设置提取10个因子

train_opts <- get_default_training_options(mofa)  # 获取默认训练选项
train_opts$maxiter <- 1000                         # 最大迭代次数
train_opts$seed <- 42                              # 随机种子

# 准备并训练模型
mofa <- prepare_mofa(mofa,                   # 准备MOFA对象
                     data_options = data_opts,
                     model_options = model_opts,
                     training_options = train_opts)
mofa <- run_mofa(mofa, outfile = "mofa_model.hdf5")  # 运行MOFA并保存模型

# 查看各因子解释的方差
plot_variance_explained(mofa)  # 可视化：每个因子在各组学中解释了多少方差

# 查看因子与分组的关联
plot_factor(mofa,                     # 可视化因子值
            factor = 1,                # 查看第1个因子
            color_by = "group")        # 按分组着色（如疾病vs对照）

# 提取因子权重（各特征对因子的贡献）
weights <- get_weights(mofa, views = "all", factors = "all")  # 获取所有权重

# ========== MintTea：疾病相关跨组学模块发现 ==========
# MintTea是2024年发表的新工具，专门找跨组学的"疾病模块"

# 安装
pip install minttea  # 安装MintTea

# 准备输入文件
# microbiome.tsv: 行=样本，列=微生物特征
# metabolome.tsv: 行=样本，列=代谢物特征
# metadata.tsv: 行=样本，列=分组信息（必须包含case/control列）

# ========== Python：运行MintTea ==========
from minttea import MintTea  # 导入MintTea

# 初始化
mt = MintTea(
    omic_files=["microbiome.tsv", "metabolome.tsv"],  # 输入：各组学数据文件
    metadata_file="metadata.tsv",                       # 输入：样本元数据
    label_column="group",                               # 分组列名
    case_label="T2D",                                   # 病例组标签
    control_label="Control"                             # 对照组标签
)

# 运行分析
results = mt.run(
    n_components=5,      # 提取5个成分
    n_permutations=100,  # 置换检验100次
    n_consensus=50       # 共识分析50次
)

# 查看结果：跨组学模块
for module in results.modules:            # 遍历每个模块
    print(f"模块: {module.name}")          # 模块名称
    print(f"微生物特征: {module.microbiome_features}")  # 模块中的微生物
    print(f"代谢物特征: {module.metabolome_features}")  # 模块中的代谢物
    print(f"疾病关联P值: {module.pvalue}")              # 与疾病的关联P值

# ========== MetOrigin：区分代谢物来源 ==========
# MetOrigin可以判断代谢物是宿主产的、微生物产的、还是共代谢产的
# 在线工具：http://metorigin.met-bioinformatics.cn/

# 也可以用KEGG数据库手动判断
# 思路：
# 1. 将代谢物映射到KEGG通路
# 2. 检查该通路在微生物和宿主中的分布
# 3. 仅在微生物中有的通路 → 微生物来源
# 4. 两者都有的通路 → 共代谢

# R语言示例：用KEGGREST查询
library(KEGGREST)  # KEGG API接口

# 查询丁酸（butyrate）相关通路
butyrate_pathways <- keggFind("compound", "butyrate")  # 搜索丁酸
print(butyrate_pathways)  # 打印相关通路

# 查看通路详情
pathway_info <- keggGet("map00650")  # 获取丁酸代谢通路详情
cat("通路名称:", pathway_info[[1]]$NAME, "\n")  # 打印通路名

========== 方法选择决策树 ==========
问题：两套数据整体相关吗？
  ├── Mantel检验（P值 + r值）
  └── Procrustes分析（可视化 + P值）

问题：哪些微生物和哪些代谢物相关？
  ├── 组成型数据 → MMVEC（推荐）
  ├── 非组成型数据 → Spearman相关 + BH校正
  └── 复杂关系 → 稀疏CCA（sCCA）

问题：有没有跨组学的疾病模块？
  ├── MintTea（有case/control分组时）
  ├── MOFA2（探索性分析）
  └── iCluster（有多个组学层时）

问题：代谢物是谁产的？
  ├── MetOrigin（在线工具）
  └── KEGG通路映射（手动判断）

========== 关键R包 ==========
整合分析：mixOmics, MOFA2, mmvec
相关分析：Hmisc, psych, correlation
可视化  ：pheatmap, ggplot2, igraph
距离计算：vegan (vegdist)

========== 关键Python包 ==========
mmvec：微生物-代谢物条件概率
minttea：疾病相关模块发现
scikit-learn：PCA/PLS等降维

========== 常用数据库 ==========
HMDB：人类代谢组数据库（代谢物鉴定）
KEGG：代谢通路数据库
MetaCyc：代谢通路百科
MassBank：质谱标准品数据库
SILVA：16S参考数据库