666 微生物组与免疫互作¶

# 分析SCFA产生菌与免疫标志物的关联
import pandas as pd  # 数据处理
import numpy as np  # 数值计算
from scipy import stats  # 统计检验

# 读取数据
microbiome = pd.read_csv("species_abundance.csv", index_col=0)  # 物种丰度
metabolome = pd.read_csv("scfa_levels.csv", index_col=0)  # SCFA浓度
immune = pd.read_csv("immune_markers.csv", index_col=0)  # 免疫指标

# SCFA产生菌列表
scfa_producers = {
    "butyrate": ["Faecalibacterium_prausnitzii",  # 普拉梭菌（最重要的丁酸菌）
                  "Roseburia_intestinalis",  # 肠道罗斯氏菌
                  "Eubacterium_rectale",  # 直肠真杆菌
                  "Anaerostipes_hadrus"],  # 哈德厌氧棒状菌
    "propionate": ["Bacteroides_spp",  # 拟杆菌属
                    "Phascolarctobacterium"],
    "acetate": ["Bifidobacterium_spp",  # 双歧杆菌属
                "Blautia_obeum"]
}

# 计算SCFA产生菌总丰度
for scfa_type, producers in scfa_producers.items():
    present = [p for p in producers if p in microbiome.index]
    if present:
        total_abundance = microbiome.loc[present].sum()  # 总丰度
        # 与SCFA浓度的相关性
        if scfa_type in metabolome.columns:
            r, p = stats.spearmanr(total_abundance, metabolome[scfa_type])
            print(f"{scfa_type}产生菌 vs {scfa_type}浓度: rho={r:.3f}, p={p:.4f}")

# 丁酸与Treg的关联
if "butyrate" in metabolome.columns and "Treg_percent" in immune.columns:
    r, p = stats.spearmanr(metabolome["butyrate"], immune["Treg_percent"])
    print(f"\n丁酸 vs Treg%: rho={r:.3f}, p={p:.4f}")
    # 丁酸促进Treg分化是菌群-免疫轴的经典通路

# 菌群-代谢物-免疫三组学整合分析
from sklearn.cross_decomposition import CCA  # 典型相关分析
from sklearn.preprocessing import StandardScaler  # 标准化

# 1. 数据标准化
scaler = StandardScaler()
X_micro = scaler.fit_transform(microbiome.T)  # 微生物组数据
X_metab = scaler.fit_transform(metabolome)   # 代谢组数据
X_immune = scaler.fit_transform(immune)      # 免疫数据

# 2. 典型相关分析：菌群 vs 免疫
cca = CCA(n_components=2)  # 2个典型变量
cca.fit(X_micro, X_immune)  # 拟合
X_c, Y_c = cca.transform(X_micro, X_immune)  # 变换

# 典型相关系数
for i in range(2):
    r, p = stats.pearsonr(X_c[:, i], Y_c[:, i])
    print(f"典型变量{i+1}: r={r:.3f}, p={p:.4f}")

# 3. 中介分析：菌群 → 代谢物 → 免疫
from statsmodels.regression.linear_model import OLS  # 线性回归
import statsmodels.api as sm

def mediation_analysis(exposure, mediator, outcome, labels):
    """
    简单中介分析
    exposure: 菌群丰度
    mediator: 代谢物浓度
    outcome: 免疫指标
    """
    # 路径a：exposure → mediator
    X_a = sm.add_constant(exposure)
    model_a = OLS(mediator, X_a).fit()
    a = model_a.params[1]  # a路径系数
    p_a = model_a.pvalues[1]

    # 路径b：mediator → outcome（控制exposure）
    X_b = sm.add_constant(np.column_stack([exposure, mediator]))
    model_b = OLS(outcome, X_b).fit()
    b = model_b.params[2]  # b路径系数

    # 路径c'：直接效应
    c_prime = model_b.params[1]

    # 间接效应 = a × b
    indirect = a * b
    # Sobel检验
    se_a = model_a.bse[1]
    se_b = model_b.bse[2]
    se_indirect = np.sqrt(b**2 * se_a**2 + a**2 * se_b**2)
    z_sobel = indirect / se_indirect
    p_sobel = 2 * (1 - stats.norm.cdf(abs(z_sobel)))

    print(f"\n中介分析: {labels[0]} → {labels[1]} → {labels[2]}")
    print(f"  间接效应(a×b): {indirect:.4f}, p(Sobel)={p_sobel:.4f}")
    print(f"  直接效应(c'): {c_prime:.4f}")
    print(f"  中介比例: {abs(indirect/(indirect+c_prime))*100:.1f}%")

# 示例：F. prausnitzii → butyrate → Treg
mediation_analysis(
    microbiome.loc["Faecalibacterium_prausnitzii"].values,
    metabolome["butyrate"].values,
    immune["Treg_percent"].values,
    ["F. prausnitzii", "butyrate", "Treg%"]
)

# 分析菌群对免疫检查点抑制剂(ICI)疗效的影响
def ici_response_analysis(microbiome_df, response_labels):
    """
    比较ICI应答者vs非应答者的菌群差异
    response_labels: "R"(应答) / "NR"(非应答)
    """
    from scipy.stats import mannwhitneyu  # 非参数检验

    results = []
    for species in microbiome_df.index:
        responders = microbiome_df.loc[species, response_labels == "R"]
        non_responders = microbiome_df.loc[species, response_labels == "NR"]

        stat, pval = mannwhitneyu(responders, non_responders,
                                   alternative='two-sided')
        fc = responders.mean() / (non_responders.mean() + 1e-10)  # fold change

        results.append({
            "species": species,
            "mean_R": responders.mean(),
            "mean_NR": non_responders.mean(),
            "fold_change": fc,
            "p_value": pval
        })

    results_df = pd.DataFrame(results)

    # FDR校正
    from statsmodels.stats.multitest import multipletests
    _, results_df["q_value"], _, _ = multipletests(
        results_df["p_value"], method="fdr_bh")

    # 显著差异物种
    sig = results_df[results_df["q_value"] < 0.1].sort_values("fold_change",
                                                                ascending=False)
    print("ICI应答相关物种:")
    print(sig[["species", "fold_change", "q_value"]].to_string())

    # 已知与ICI应答相关的菌属
    # Akkermansia muciniphila ↑ → 抗PD-1应答
    # Faecalibacterium prausnitzii ↑ → 抗CTLA-4应答
    # Bacteroides thetaiotaomicron ↑ → 抗CTLA-4应答

    return sig

# R语言：菌群-免疫-脑轴的网络分析
library(igraph)    # 网络分析
library(ggplot2)   # 绑图
library(Hmisc)     # 相关分析

# 读取三层数据
microbiome <- read.csv("microbiome.csv", row.names = 1)  # 菌群
immune <- read.csv("immune.csv", row.names = 1)  # 免疫指标
brain <- read.csv("brain_markers.csv", row.names = 1)  # 脑功能指标

# 合并数据
all_data <- cbind(t(microbiome), immune, brain)

# 计算相关矩阵（Spearman）
cor_result <- rcorr(as.matrix(all_data), type = "spearman")
r_matrix <- cor_result$r  # 相关系数
p_matrix <- cor_result$P  # p值

# 保留显著相关（FDR<0.05，|r|>0.3）
p_adj <- p.adjust(p_matrix[lower.tri(p_matrix)], method = "BH")
significant <- abs(r_matrix) > 0.3 & p_matrix < 0.05

# 构建网络
adj_matrix <- r_matrix * significant  # 显著相关的邻接矩阵
g <- graph_from_adjacency_matrix(adj_matrix, weighted = TRUE,
                                  mode = "undirected", diag = FALSE)

# 节点分类
V(g)$type <- ifelse(V(g)$name %in% rownames(microbiome), "菌群",
              ifelse(V(g)$name %in% colnames(immune), "免疫", "脑功能"))
V(g)$color <- c("菌群" = "#4CAF50", "免疫" = "#FF9800",
                 "脑功能" = "#2196F3")[V(g)$type]

# 可视化
plot(g,
     vertex.size = degree(g) * 2,  # 节点大小按连接度
     vertex.label.cex = 0.7,
     edge.width = abs(E(g)$weight) * 3,
     edge.color = ifelse(E(g)$weight > 0, "#66BB6A", "#EF5350"),  # 绿=正，红=负
     layout = layout_with_fr(g),
     main = "菌群-免疫-脑轴互作网络")

# 菌群-免疫互作分析流程
1. 菌群分析: 16S/宏基因组 → 物种+功能组成
2. 代谢物检测: LC-MS/GC-MS → SCFA/胆汁酸/色氨酸代谢物
3. 免疫指标: 流式细胞术/CyTOF → 免疫细胞亚群
4. 关联分析: Spearman/CCA/MOFA
5. 因果推断: 中介分析/孟德尔随机化/体外验证

# 关键免疫-菌群互作
丁酸 → 组蛋白去乙酰化 → Treg分化 → 抗炎
吲哚 → AhR受体 → IL-22产生 → 肠道屏障
PSA(脆弱拟杆菌) → TLR2 → IL-10 → 免疫耐受
LPS(革兰阴性菌) → TLR4 → NF-κB → 促炎

# 分析工具
PICRUSt2: 从16S预测功能基因（含SCFA通路）
HUMAnN3: 宏基因组功能通路定量
MaAsLin2: 多变量关联分析
MOFA+: 多组学因子分析