代谢组与微生物组整合分析¶

# === 微生物组-代谢组整合分析 ===
import pandas as pd
import numpy as np
from scipy.stats import spearmanr
import seaborn as sns
import matplotlib.pyplot as plt

# 1. 读取数据
# 微生物属水平相对丰度（CLR变换后）
microbiome = pd.read_csv("genus_clr.csv", index_col=0)  # 行=样本，列=菌属
# 代谢物浓度（log变换后）
metabolome = pd.read_csv("metabolites_log.csv", index_col=0)  # 行=样本，列=代谢物

# 确保样本顺序一致
common_samples = microbiome.index.intersection(metabolome.index)
microbiome = microbiome.loc[common_samples]
metabolome = metabolome.loc[common_samples]

# 2. 计算Spearman相关矩阵
n_microbes = microbiome.shape[1]
n_metabolites = metabolome.shape[1]
corr_matrix = np.zeros((n_microbes, n_metabolites))  # 相关系数矩阵
pval_matrix = np.zeros((n_microbes, n_metabolites))   # p值矩阵

for i, mic in enumerate(microbiome.columns):
    for j, met in enumerate(metabolome.columns):
        rho, pval = spearmanr(microbiome[mic], metabolome[met])
        corr_matrix[i, j] = rho
        pval_matrix[i, j] = pval

# 3. FDR校正
from statsmodels.stats.multitest import multipletests
pvals_flat = pval_matrix.flatten()
_, padj_flat, _, _ = multipletests(pvals_flat, method="fdr_bh")
padj_matrix = padj_flat.reshape(pval_matrix.shape)

# 4. 筛选显著相关（|rho|>0.4, FDR<0.05）
significant = (np.abs(corr_matrix) > 0.4) & (padj_matrix < 0.05)
n_sig = significant.sum()
print(f"显著相关对数: {n_sig}")

# 5. 绘制相关性热图（只显示显著的）
corr_df = pd.DataFrame(corr_matrix,
                        index=microbiome.columns,
                        columns=metabolome.columns)
# 只保留至少有一个显著相关的行和列
sig_rows = significant.any(axis=1)
sig_cols = significant.any(axis=0)
corr_sig = corr_df.loc[sig_rows, sig_cols]

plt.figure(figsize=(14, 10))
sns.heatmap(corr_sig, cmap="RdBu_r", center=0, vmin=-0.8, vmax=0.8,
            xticklabels=True, yticklabels=True)
plt.title("Microbiome-Metabolome Spearman Correlation")
plt.tight_layout()
plt.savefig("microbe_metabolite_heatmap.png", dpi=150)

# === 使用mixOmics的DIABLO做多组学整合 ===
library(mixOmics)

# 1. 准备数据列表
X <- list(
  microbiome = as.matrix(microbiome_data),  # 微生物丰度矩阵
  metabolome = as.matrix(metabolome_data)   # 代谢物浓度矩阵
)
Y <- factor(group_labels)  # 分组标签（如Disease/Control）

# 2. 参数调优
design <- matrix(0.1, ncol=2, nrow=2)  # 设置组学间关联权重
diag(design) <- 0

# 3. 运行DIABLO
result <- block.splsda(X, Y, ncomp=2, design=design)

# 4. 可视化
plotDiablo(result, ncomp=1)  # 组学间相关得分图
circosPlot(result, cutoff=0.7)  # 环形图展示跨组学关联