一、引言：多组学时代的因果推断困境与破局之道

随着高通量测序技术的飞速发展，基因组、转录组、蛋白质组、代谢组等多维度数据的整合分析已成为生物信息学研究的核心范式。然而，多组学数据的复杂性也带来了严峻挑战：传统关联分析（如相关性检验、差异表达分析）往往只能揭示变量间的统计关联，无法区分因果关系与混杂因素导致的虚假关联，这使得许多基于关联分析的分子机制研究陷入 “重复验证困难”“临床转化受阻” 的困境。

例如，在肿瘤研究中，通过转录组数据分析发现基因 A 的表达量与患者预后显著相关，但后续功能实验却证实这种关联是由肿瘤微环境中的炎症因子（混杂因素）介导的，基因 A 本身并非驱动预后的因果因子。类似地，在复杂疾病的 biomarker 筛选中，约 70% 的关联分析发现的候选分子无法通过临床验证，核心原因在于缺乏对因果关系的精准推断。

孟德尔随机化（Mendelian Randomization, MR）作为一种基于 “遗传变异为工具变量” 的因果推断方法，为解决这一问题提供了新思路。其核心逻辑是：遗传变异（如 SNP）在减数分裂时随机分配，不受后天环境和表型的影响，且仅通过特定暴露因素（如基因表达、代谢物水平）影响结局（如疾病状态），因此可模拟 “自然随机对照试验”，有效排除混杂因素，推断暴露与结局的因果关系。

但传统 MR 分析仍存在局限：仅能分析单一层级的暴露因素（如单一基因、单一代谢物），无法整合多组学数据揭示复杂的分子调控网络；且对非线性因果关系、多中介效应的建模能力不足。而人工智能（AI）技术（如深度学习、图神经网络）具备强大的高维数据处理和复杂关系建模能力，能够有效弥补传统 MR 的缺陷。

本文将系统介绍 “孟德尔随机化 + AI 融合建模” 的核心逻辑、技术流程，并提供 R/Python 双语言实战代码，帮助读者掌握多组学数据中的因果机制挖掘方法，为科研项目和临床转化研究提供有力支撑。

二、核心理论基础：孟德尔随机化与 AI 融合的底层逻辑

2.1 孟德尔随机化（MR）的核心原理与关键假设

MR 的核心是利用遗传变异作为 “工具变量（Instrumental Variable, IV）”，推断暴露（Exposure）与结局（Outcome）的因果关系。其成立需满足三大关键假设：

1. 相关性假设：工具变量（SNP）与暴露因素显著相关（即 SNP 对暴露有强预测力）；
2. 独立性假设：工具变量与结局的混杂因素（如年龄、性别、生活习惯）无关联；
3. 排他性假设：工具变量仅通过暴露因素影响结局，不直接影响结局或通过其他路径影响结局。

根据暴露和结局的类型，MR 可分为不同亚型，适用于多组学场景的主要包括：

• 两样本 MR：暴露数据与结局数据来自不同人群（但需满足人群遗传背景一致性），可利用公共数据库（如 GTEx、UK Biobank、FinnGen）快速开展分析，是目前应用最广泛的 MR 类型；
• 多暴露 MR：同时分析多个暴露因素（如多个基因、多个代谢物）与结局的因果关系，适用于多组学数据整合；
• 中介 MR：分析暴露通过中介变量（如基因→蛋白质→代谢物）影响结局的因果路径，可揭示分子调控的级联效应。

2.2 AI 技术在 MR 融合建模中的核心作用

AI 技术并非替代 MR，而是通过以下方式增强 MR 的分析能力：

1. 高维工具变量筛选：利用机器学习（如 LASSO、随机森林）从全基因组 SNP 中筛选出对多组学暴露（如基因表达谱、代谢组谱）有强预测力的工具变量，解决传统 MR 中 “工具变量筛选依赖先验知识” 的问题；
2. 复杂因果关系建模：通过深度学习（如神经网络、注意力机制）捕捉暴露与结局的非线性因果关系、多因素交互作用，突破传统 MR 仅能建模线性关系的局限；
3. 多组学网络构建：利用图神经网络（GNN）整合基因组、转录组、蛋白质组等多维度数据，构建 “遗传变异→多组学暴露→结局” 的因果网络，揭示分子机制的系统性调控关系；
4. 因果效应量化与解释：通过可解释 AI（XAI）技术（如 SHAP 值、注意力权重）解释 AI 模型中各变量的因果贡献度，解决传统 AI “黑箱” 问题，增强结果的生物学可解释性。

2.3 孟德尔随机化 + AI 融合建模的核心流程

融合建模的核心思路是：以 MR 的因果推断逻辑为框架，利用 AI 技术解决 MR 在高维数据处理、复杂关系建模中的不足，最终实现 “因果推断 + 机制解析” 的双重目标。具体流程如下：

1. 数据预处理：整合多组学暴露数据（基因组、转录组、蛋白质组等）、结局数据（疾病状态、表型数据）和遗传变异数据，进行质量控制和标准化；
2. 工具变量筛选：利用 AI 算法（如 LASSO、随机森林）从全基因组 SNP 中筛选出与多组学暴露显著相关的工具变量，并验证 MR 三大假设；
3. 因果关系初步推断：基于筛选的工具变量，通过传统 MR 方法（如 Inverse Variance Weighted, IVW）初步验证暴露与结局的因果关系；
4. AI 融合建模：构建 “工具变量→多组学暴露→结局” 的 AI 模型（如深度学习模型、GNN 模型），捕捉复杂因果关系和调控网络；
5. 结果验证与解释：通过交叉验证、敏感性分析验证模型稳定性，利用可解释 AI 技术解析关键因果路径和核心分子；
6. 生物学验证：结合功能实验（如 CRISPR-Cas9 基因编辑、细胞实验）验证核心因果关系，形成 “数据分析→机制提出→实验验证” 的闭环。

三、实战准备：数据来源与软件环境搭建

3.1 核心数据来源

融合建模需要三类核心数据，推荐优先使用公共数据库资源，降低数据采集成本：

1. 遗传变异数据：全基因组关联研究（GWAS）汇总数据（如 NHGRI-EBI GWAS Catalog、FinnGen、UK Biobank），包含 SNP 位点、等位基因、效应值（beta）、P 值等信息；
2. 多组学暴露数据：
   • 转录组：GTEx（基因型 - 组织表达数据库）、GEUVADIS（淋巴母细胞系表达数据）；
   • 蛋白质组：Somalogic、Olink Explore（血浆蛋白质组数据）；
   • 代谢组：Metabolomics GWAS Server、UK Biobank 代谢组数据；
3. 结局数据：疾病状态（如肿瘤、心血管疾病）、表型数据（如身高、血糖），可来自 GWAS Catalog、FinnGen、UK Biobank 等数据库。

数据格式要求：

• 遗传变异数据：需包含 SNP ID（rs 号）、染色体位置、等位基因（效应等位基因 / 非效应等位基因）、效应值（beta）、标准误（se）、P 值、样本量；
• 多组学暴露数据：矩阵格式，行为样本，列为分子特征（基因、蛋白质、代谢物），需包含样本的遗传信息（如 SNP 基因型）；
• 结局数据：二分类变量（如疾病与否）或连续变量（如表型数值），需包含样本的协变量信息（如年龄、性别）。

3.2 软件环境搭建

3.2.1 R 语言环境

推荐使用 R 4.2.0 及以上版本，核心包安装：

# 安装核心包
install.packages(c("TwoSampleMR", "MRInstruments", "MRmediation", "glmnet", "randomForest", "xgboost"))
install.packages(c("torch", "rTorch", "ggraph", "igraph", "shapr", "DALEX"))
# 加载包
library(TwoSampleMR)  # 两样本MR分析
library(MRInstruments)  # 工具变量筛选
library(MRmediation)  # 中介MR分析
library(glmnet)  # LASSO回归（工具变量筛选）
library(randomForest)  # 随机森林（工具变量筛选）
library(xgboost)  # XGBoost（预测模型）
library(torch)  # 深度学习框架
library(rTorch)  # R与PyTorch接口
library(ggraph)  # 因果网络可视化
library(shapr)  # SHAP值解释
library(DALEX)  # 模型解释

3.2.2 Python 环境

推荐使用 Python 3.8 及以上版本，核心库安装：

pip install pandas numpy scipy matplotlib seaborn scikit-learn torch torchvision torchaudio dgl dglgo shap networkx

核心库说明：

• pandas/numpy/scipy：数据处理与统计分析；
• scikit-learn：机器学习算法（LASSO、随机森林）；
• torch：深度学习框架；
• dgl/dglgo：图神经网络构建；
• shap：SHAP 值解释；
• networkx/matplotlib/seaborn：可视化。

四、实战流程：孟德尔随机化 + AI 融合建模（R/Python 双语言）

本节将以 “转录组 + 蛋白质组多组学暴露→肺癌风险” 为例，详细演示融合建模的完整流程，包含数据预处理、工具变量筛选、MR 因果推断、AI 融合建模、结果解释与验证等关键步骤。

4.1 步骤 1：数据预处理与质量控制

4.1.1 数据加载与整合（R 语言）

# 1. 加载数据（示例数据，实际需替换为真实数据或公共数据库数据）
# 遗传变异数据（GWAS汇总数据，结局：肺癌）
gwas_lung_cancer <- read.csv("gwas_lung_cancer_summary.csv", stringsAsFactors = F)
# 转录组数据（GTEx肺组织）
transcriptome <- read.csv("gtex_lung_transcriptome.csv", row.names = 1)
# 蛋白质组数据（Olink肺相关蛋白质组）
proteome <- read.csv("olink_lung_proteome.csv", row.names = 1)
# 样本信息（包含SNP基因型、协变量）
sample_info <- read.csv("sample_info.csv", stringsAsFactors = F)

# 2. 数据质量控制
# 遗传变异数据QC：过滤低质量SNP（P值<5e-8，样本量>5000）
gwas_qc <- gwas_lung_cancer[gwas_lung_cancer$p_value < 5e-8 & gwas_lung_cancer$sample_size > 5000, ]
# 过滤等位基因信息缺失的SNP
gwas_qc <- gwas_qc[!is.na(gwas_qc$effect_allele) & !is.na(gwas_qc$other_allele), ]

# 多组学数据QC：过滤缺失值>20%的分子特征
transcriptome_qc <- transcriptome[, apply(transcriptome, 2, function(x) sum(is.na(x))/nrow(transcriptome) < 0.2)]
proteome_qc <- proteome[, apply(proteome, 2, function(x) sum(is.na(x))/nrow(proteome) < 0.2)]

# 3. 数据整合：确保样本匹配
common_samples <- intersect(intersect(rownames(transcriptome_qc), rownames(proteome_qc)), sample_info$sample_id)
transcriptome_final <- transcriptome_qc[common_samples, ]
proteome_final <- proteome_qc[common_samples, ]
sample_info_final <- sample_info[match(common_samples, sample_info$sample_id), ]

# 4. 标准化：多组学数据z-score标准化
transcriptome_norm <- t(scale(t(transcriptome_final)))  # 按基因标准化
proteome_norm <- t(scale(t(proteome_final)))  # 按蛋白质标准化

# 5. 合并多组学暴露数据
exposure_data <- cbind(transcriptome_norm, proteome_norm)

4.1.2 数据预处理（Python 语言）

import pandas as pd
import numpy as np
from sklearn.preprocessing import StandardScaler

# 1. 加载数据
gwas_lung_cancer = pd.read_csv("gwas_lung_cancer_summary.csv")
transcriptome = pd.read_csv("gtex_lung_transcriptome.csv", index_col=0)
proteome = pd.read_csv("olink_lung_proteome.csv", index_col=0)
sample_info = pd.read_csv("sample_info.csv")

# 2. 遗传变异数据QC
gwas_qc = gwas_lung_cancer[(gwas_lung_cancer["p_value"] < 5e-8) & (gwas_lung_cancer["sample_size"] > 5000)]
gwas_qc = gwas_qc.dropna(subset=["effect_allele", "other_allele"])

# 3. 多组学数据QC：过滤缺失值>20%的特征
transcriptome_qc = transcriptome.dropna(thresh=0.8*len(transcriptome), axis=1)
proteome_qc = proteome.dropna(thresh=0.8*len(proteome), axis=1)

# 4. 样本匹配
common_samples = list(set(transcriptome_qc.index) & set(proteome_qc.index) & set(sample_info["sample_id"]))
transcriptome_final = transcriptome_qc.loc[common_samples]
proteome_final = proteome_qc.loc[common_samples]
sample_info_final = sample_info[sample_info["sample_id"].isin(common_samples)].set_index("sample_id")

# 5. 标准化
scaler = StandardScaler()
transcriptome_norm = pd.DataFrame(
    scaler.fit_transform(transcriptome_final),
    index=transcriptome_final.index,
    columns=transcriptome_final.columns
)
proteome_norm = pd.DataFrame(
    scaler.fit_transform(proteome_final),
    index=proteome_final.index,
    columns=proteome_final.columns
)

# 6. 合并多组学暴露数据
exposure_data = pd.concat([transcriptome_norm, proteome_norm], axis=1)

4.2 步骤 2：工具变量筛选（AI 增强 MR）

传统 MR 工具变量筛选依赖单 SNP - 暴露关联分析，效率低且易遗漏弱关联但联合预测力强的 SNP。此处采用 LASSO 回归（AI 算法）筛选对多组学暴露有强预测力的工具变量，同时满足 MR 假设。

4.2.1 工具变量筛选（R 语言）

# 1. 提取样本的SNP基因型数据（假设sample_info_final中包含SNP基因型，0/1/2编码）
snp_genotype <- sample_info_final[, grep("rs", colnames(sample_info_final))]  # 筛选SNP列
snp_genotype <- na.omit(snp_genotype)  # 过滤基因型缺失的样本

# 2. 匹配暴露数据与SNP数据的样本
common_samples_snp_exposure <- intersect(rownames(snp_genotype), rownames(exposure_data))
snp_genotype_final <- snp_genotype[common_samples_snp_exposure, ]
exposure_data_final <- exposure_data[common_samples_snp_exposure, ]

# 3. LASSO回归筛选工具变量：以多组学暴露为因变量，SNP为自变量
# 由于暴露是多维度的，采用多响应LASSO回归
x <- as.matrix(snp_genotype_final)
y <- as.matrix(exposure_data_final)

# 交叉验证确定LASSO惩罚系数
cv_lasso <- glmnet::cv.glmnet(
  x = x,
  y = y,
  family = "mgaussian",  # 多响应变量
  alpha = 1,  # LASSO回归（alpha=1）
  nfolds = 10  # 10折交叉验证
)

# 提取最优惩罚系数对应的模型
best_lasso <- glmnet::glmnet(
  x = x,
  y = y,
  family = "mgaussian",
  alpha = 1,
  lambda = cv_lasso$lambda.min
)

# 4. 筛选非零系数的SNP（工具变量）
iv_coef <- coef(best_lasso)
iv_snp <- c()
for (i in 1:length(iv_coef)) {
  coef_mat <- as.matrix(iv_coef[[i]])
  non_zero_snp <- rownames(coef_mat)[coef_mat[, 1] != 0]
  iv_snp <- c(iv_snp, non_zero_snp)
}
iv_snp <- unique(iv_snp)[iv_snp != "(Intercept)"]  # 去重并排除截距项

# 5. 验证MR假设：工具变量与暴露的相关性（相关性假设）
iv_exposure_corr <- cor(x[, iv_snp], y)
iv_exposure_corr_mean <- apply(iv_exposure_corr, 1, function(x) mean(abs(x)))  # 计算每个SNP与所有暴露的平均相关性
iv_snp_final <- names(iv_exposure_corr_mean[iv_exposure_corr_mean > 0.1])  # 筛选平均相关性>0.1的SNP

# 6. 输出工具变量列表
write.csv(data.frame(snp_id = iv_snp_final), "instrumental_variables.csv", row.names = F)

4.2.2 工具变量筛选（Python 语言）

import pandas as pd
import numpy as np
from sklearn.linear_model import MultiTaskLassoCV
from sklearn.preprocessing import StandardScaler

# 1. 提取SNP基因型数据（0/1/2编码）
snp_genotype = sample_info_final.filter(regex="^rs", axis=1)  # 筛选SNP列
snp_genotype = snp_genotype.dropna()  # 过滤基因型缺失的样本

# 2. 样本匹配
common_samples_snp_exposure = list(set(snp_genotype.index) & set(exposure_data.index))
snp_genotype_final = snp_genotype.loc[common_samples_snp_exposure]
exposure_data_final = exposure_data.loc[common_samples_snp_exposure]

# 3. 多任务LASSO回归筛选工具变量
x = snp_genotype_final.values
y = exposure_data_final.values

# 标准化自变量
scaler_x = StandardScaler()
x_scaled = scaler_x.fit_transform(x)

# 多任务LASSO交叉验证（适合多响应变量）
mt_lasso = MultiTaskLassoCV(
    cv=10,  # 10折交叉验证
    random_state=42,
    n_jobs=-1
)
mt_lasso.fit(x_scaled, y)

# 4. 筛选非零系数的SNP（工具变量）
coef_matrix = mt_lasso.coef_  # 系数矩阵：(n_exposures, n_snps)
non_zero_coef_mask = np.any(coef_matrix != 0, axis=0)  # 至少对一个暴露有非零系数的SNP
iv_snp = snp_genotype_final.columns[non_zero_coef_mask]

# 5. 验证相关性假设：计算SNP与暴露的平均相关性
iv_snp_indices = [list(snp_genotype_final.columns).index(snp) for snp in iv_snp]
x_iv = x_scaled[:, iv_snp_indices]
iv_exposure_corr = np.corrcoef(x_iv.T, y.T)[:len(iv_snp), len(iv_snp):]
iv_exposure_corr_mean = np.mean(np.abs(iv_exposure_corr), axis=1)  # 每个SNP与所有暴露的平均相关性
iv_snp_final = iv_snp[iv_exposure_corr_mean > 0.1]  # 筛选平均相关性>0.1的SNP

# 6. 输出工具变量列表
pd.DataFrame({"snp_id": iv_snp_final}).to_csv("instrumental_variables.csv", index=False)

4.3 步骤 3：传统 MR 分析（因果关系初步验证）

利用筛选的工具变量，通过两样本 MR 分析初步验证多组学暴露与肺癌风险的因果关系。此处以 “转录组暴露” 为例，蛋白质组及联合暴露的分析流程类似。

4.3.1 两样本 MR 分析（R 语言）

library(TwoSampleMR)

# 1. 构建暴露数据（转录组基因表达）
# 提取转录组暴露数据
transcriptome_exposure <- exposure_data_final[, colnames(exposure_data_final) %in% colnames(transcriptome_norm)]
# 构建暴露数据集（格式：SNP ID、效应等位基因、非效应等位基因、效应值、标准误、P值）
exposure_mr <- data.frame()
for (gene in colnames(transcriptome_exposure)) {
  # 计算每个SNP与该基因的回归系数（效应值）和标准误
  for (snp in iv_snp_final) {
    model <- lm(transcriptome_exposure[, gene] ~ snp_genotype_final[, snp])
    beta <- coef(model)[2]
    se <- summary(model)$coefficients[2, 2]
    p_val <- summary(model)$coefficients[2, 4]
    exposure_mr <- rbind(exposure_mr, data.frame(
      snp = snp,
      exposure = gene,
      effect_allele = gwas_qc$effect_allele[match(snp, gwas_qc$snp_id)],
      other_allele = gwas_qc$other_allele[match(snp, gwas_qc$snp_id)],
      beta = beta,
      se = se,
      pval = p_val
    ))
  }
}

# 2. 构建结局数据（肺癌GWAS汇总数据）
outcome_mr <- gwas_qc[, c("snp_id", "effect_allele", "other_allele", "beta", "se", "p_value")]
colnames(outcome_mr) <- c("snp", "effect_allele", "other_allele", "beta", "se", "pval")
outcome_mr$outcome <- "lung_cancer"

# 3. 两样本MR分析：匹配暴露与结局的SNP
mr_data <- harmonise_data(
  exposure_dat = exposure_mr,
  outcome_dat = outcome_mr,
  action = 2  # 自动匹配等位基因
)

# 4. MR分析：采用IVW（主要方法）、MR-Egger、Weighted Median等多种方法验证
mr_results <- mr(
  mr_data,
  method_list = c("ivw", "mr_egger_regression", "weighted_median", "simple_mode")
)

# 5. 敏感性分析：验证结果稳健性
# 异质性检验（I²统计量）
heterogeneity <- mr_heterogeneity(mr_data, method_list = c("ivw", "mr_egger_regression"))
# 水平多效性检验（MR-Egger截距项）
pleiotropy <- mr_pleiotropy_test(mr_data)

# 6. 输出MR结果
write.csv(mr_results, "mr_results.csv", row.names = F)
write.csv(heterogeneity, "mr_heterogeneity.csv", row.names = F)
write.csv(pleiotropy, "mr_pleiotropy.csv", row.names = F)

# 7. 可视化：森林图（展示各SNP的因果效应）
p_forest <- mr_forest_plot(mr_results, mr_data)
ggsave("mr_forest_plot.pdf", p_forest, width = 12, height = 8)

4.3.2 两样本 MR 分析（Python 语言）

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
import seaborn as sns

# 1. 构建暴露数据（转录组基因表达）
transcriptome_exposure = exposure_data_final.filter(regex="^ENSG", axis=1)  # 假设基因ID以ENSG开头
exposure_mr = []

for gene in transcriptome_exposure.columns:
    for snp in iv_snp_final:
        # 线性回归计算SNP对基因表达的效应值和标准误
        x_snp = snp_genotype_final[snp].values.reshape(-1, 1)
        y_gene = transcriptome_exposure[gene].values
        model = stats.linregress(x_snp.flatten(), y_gene)
        beta = model.slope
        se = model.stderr
        p_val = model.pvalue
        
        # 获取SNP的等位基因信息
        snp_info = gwas_qc[gwas_qc["snp_id"] == snp].iloc[0] if snp in gwas_qc["snp_id"].values else None
        effect_allele = snp_info["effect_allele"] if snp_info is not None else "NA"
        other_allele = snp_info["other_allele"] if snp_info is not None else "NA"
        
        exposure_mr.append({
            "snp": snp,
            "exposure": gene,
            "effect_allele": effect_allele,
            "other_allele": other_allele,
            "beta": beta,
            "se": se,
            "pval": p_val
        })

exposure_mr = pd.DataFrame(exposure_mr)

# 2. 构建结局数据
outcome_mr = gwas_qc[["snp_id", "effect_allele", "other_allele", "beta", "se", "p_value"]].copy()
outcome_mr.columns = ["snp", "effect_allele", "other_allele", "beta", "se", "pval"]
outcome_mr["outcome"] = "lung_cancer"

# 3. 匹配暴露与结局的SNP（等位基因对齐）
def harmonize_alleles(exposure_df, outcome_df):
    # 合并暴露和结局数据
    merged = pd.merge(
        exposure_df,
        outcome_df,
        on="snp",
        suffixes=("_exposure", "_outcome")
    )
    
    # 等位基因对齐：如果效应等位基因不同，翻转暴露的beta值
    merged["beta_exposure_harmonized"] = merged.apply(
        lambda x: -x["beta_exposure"] if x["effect_allele_exposure"] != x["effect_allele_outcome"] else x["beta_exposure"],
        axis=1
    )
    
    return merged

mr_data = harmonize_alleles(exposure_mr, outcome_mr)

# 4. MR分析：IVW方法（主要方法）
def ivw_mr(mr_data):
    # 计算权重（1/se²）
    mr_data["weight"] = 1 / (mr_data["se_outcome"] ** 2)
    # IVW效应值 = sum(beta_exposure * beta_outcome * weight) / sum(beta_exposure² * weight)
    numerator = np.sum(mr_data["beta_exposure_harmonized"] * mr_data["beta_outcome"] * mr_data["weight"])
    denominator = np.sum((mr_data["beta_exposure_harmonized"] ** 2) * mr_data["weight"])
    ivw_beta = numerator / denominator
    # 标准误
    ivw_se = np.sqrt(1 / denominator)
    # P值
    ivw_pval = 2 * (1 - stats.norm.cdf(abs(ivw_beta / ivw_se)))
    return pd.DataFrame({
        "method": "ivw",
        "beta": [ivw_beta],
        "se": [ivw_se],
        "pval": [ivw_pval]
    })

mr_results = ivw_mr(mr_data)

# 5. 敏感性分析：异质性检验（I²统计量）
def heterogeneity_test(mr_data):
    # 计算每个SNP的效应值
    mr_data["snp_effect"] = mr_data["beta_outcome"] / mr_data["beta_exposure_harmonized"]
    # 加权方差
    mr_data["weight"] = 1 / (mr_data["se_outcome"] ** 2)
    weighted_mean = np.sum(mr_data["snp_effect"] * mr_data["weight"]) / np.sum(mr_data["weight"])
    # 总异质性
    q_total = np.sum(mr_data["weight"] * (mr_data["snp_effect"] - weighted_mean) ** 2)
    # 自由度
    df = len(mr_data) - 1
    # I²统计量
    i2 = (q_total - df) / q_total * 100 if q_total > df else 0
    return pd.DataFrame({"i2": [i2], "q_statistic": [q_total], "df": [df]})

heterogeneity = heterogeneity_test(mr_data)

# 6. 可视化：森林图
plt.figure(figsize=(12, 8))
sns.barplot(x="snp_effect", y="snp", data=mr_data, palette="viridis")
plt.axvline(x=mr_results["beta"].iloc[0], color="red", linestyle="--", label=f"IVW effect (beta={mr_results['beta'].iloc[0]:.3f})")
plt.xlabel("Causal Effect (beta)")
plt.ylabel("SNP")
plt.title("MR Forest Plot: Transcriptome Exposure vs Lung Cancer Risk")
plt.legend()
plt.tight_layout()
plt.savefig("mr_forest_plot.pdf")

# 7. 输出结果
mr_results.to_csv("mr_results.csv", index=False)
heterogeneity.to_csv("mr_heterogeneity.csv", index=False)

4.4 步骤 4：AI 融合建模（多组学因果网络构建）

通过图神经网络（GNN）构建 “工具变量→多组学暴露→结局” 的因果网络，捕捉复杂的调控关系和非线性因果效应。GNN 适合处理网络结构数据，能够有效整合多组学节点和边（因果关系）信息。

4.4.1 图神经网络因果建模（R 语言）

library(torch)
library(dgl)
library(igraph)
library(ggraph)

# 1. 构建图数据结构
# 节点：工具变量（SNP）、多组学暴露（基因、蛋白质）、结局（肺癌）
nodes_snp <- data.frame(id = paste0("snp_", iv_snp_final), type = "snp")
nodes_exposure <- data.frame(id = colnames(exposure_data_final), type = ifelse(colnames(exposure_data_final) %in% colnames(transcriptome_norm), "gene", "protein"))
nodes_outcome <- data.frame(id = "lung_cancer", type = "outcome")
nodes <- rbind(nodes_snp, nodes_exposure, nodes_outcome)

# 边：SNP→暴露（基于LASSO回归系数）、暴露→结局（基于MR因果效应）
# SNP→暴露边
edges_snp_exposure <- data.frame()
for (snp in iv_snp_final) {
  for (exposure in colnames(exposure_data_final)) {
    coef <- coef(best_lasso)[[which(colnames(exposure_data_final) == exposure)]][snp, 1]
    if (coef != 0) {
      edges_snp_exposure <- rbind(edges_snp_exposure, data.frame(
        source = paste0("snp_", snp),
        target = exposure,
        weight = abs(coef)
      ))
    }
  }
}

# 暴露→结局边（基于MR结果）
mr_exposure_outcome <- mr_results[mr_results$pval < 0.05, ]  # 筛选显著因果效应的暴露
edges_exposure_outcome <- data.frame(
  source = mr_exposure_outcome$exposure,
  target = "lung_cancer",
  weight = abs(mr_exposure_outcome$beta)
)

# 合并边
edges <- rbind(edges_snp_exposure, edges_exposure_outcome)

# 2. 构建DGL图
# 节点ID映射
node_ids <- setNames(0:(nrow(nodes)-1), nodes$id)
edges$source_id <- node_ids[edges$source]
edges$target_id <- node_ids[edges$target]

# 创建图
g <- dgl_graph()
g <- dgl_add_nodes(g, nrow(nodes))
g <- dgl_add_edges(g, edges$source_id, edges$target_id)
dgl_edges_attr_set(g, "weight", as.matrix(edges$weight))

# 3. 节点特征：SNP基因型、暴露表达量、结局状态
# SNP节点特征：基因型均值
snp_features <- apply(snp_genotype_final[, iv_snp_final], 2, mean)
# 暴露节点特征：表达量均值
exposure_features <- apply(exposure_data_final, 2, mean)
# 结局节点特征：1（肺癌）
outcome_features <- 1

# 合并节点特征
node_features <- c(snp_features, exposure_features, outcome_features)
node_features <- matrix(node_features, nrow = nrow(nodes), ncol = 1)  # 1维特征

# 4. 构建GNN模型（GCN）
model <- nn_module(
  "GCN",
  initialize = function(in_feats, hidden_feats, out_feats) {
    self$conv1 <- dglnn_gcn_conv(in_feats, hidden_feats)
    self$conv2 <- dglnn_gcn_conv(hidden_feats, out_feats)
  },
  forward = function(g, features) {
    x <- torch_relu(self$conv1(g, features))
    x <- self$conv2(g, features)
    return(x)
  }
)

# 初始化模型
in_feats <- ncol(node_features)
hidden_feats <- 64
out_feats <- 1  # 结局预测（肺癌风险）
net <- model(in_feats, hidden_feats, out_feats)

# 5. 训练模型
# 数据准备
features <- torch_tensor(node_features, dtype = torch_float32)
labels <- torch_tensor(rep(1, nrow(nodes)), dtype = torch_float32)  # 简化标签，实际需根据结局数据调整

# 优化器
optimizer <- torch_optim_adam(net$parameters(), lr = 0.01)

# 训练循环
for (epoch in 1:100) {
  optimizer$zero_grad()
  logits <- net(g, features)
  loss <- nnf_mse_loss(logits, labels)
  loss$backward()
  optimizer$step()
  if (epoch %% 10 == 0) {
    cat(sprintf("Epoch %d, Loss: %.4f\n", epoch, loss$item()))
  }
}

# 6. 提取节点重要性（因果贡献度）
node_importance <- as.array(logits$detach())
nodes$importance <- node_importance

# 7. 可视化因果网络
p_causal_network <- ggraph(igraph::graph_from_data_frame(edges, vertices = nodes), layout = "kk") +
  geom_edge_link(aes(edge_width = weight), alpha = 0.6) +
  geom_node_point(aes(size = importance, color = type)) +
  geom_node_text(aes(label = id), size = 3, repel = TRUE) +
  scale_color_manual(values = c("snp" = "blue", "gene" = "green", "protein" = "orange", "outcome" = "red")) +
  theme_minimal() +
  labs(title = "Causal Network: MR + GNN Fusion Model", edge_width = "Edge Weight", size = "Node Importance")

ggsave("causal_network.pdf", p_causal_network, width = 16, height = 12)

4.4.2 图神经网络因果建模（Python 语言）

import torch
import torch.nn as nn
import torch.optim as optim
import dgl
import dgl.nn as dglnn
import networkx as nx
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import seaborn as sns

# 1. 构建图数据结构
# 节点：SNP、基因、蛋白质、结局
nodes_snp = pd.DataFrame({"id": [f"snp_{snp}" for snp in iv_snp_final], "type": "snp"})
nodes_gene = pd.DataFrame({"id": colnames(transcriptome_norm), "type": "gene"})
nodes_protein = pd.DataFrame({"id": colnames(proteome_norm), "type": "protein"})
nodes_outcome = pd.DataFrame({"id": ["lung_cancer"], "type": "outcome"})
nodes = pd.concat([nodes_snp, nodes_gene, nodes_protein, nodes_outcome], ignore_index=True)

# 边：SNP→暴露（基于LASSO系数）、暴露→结局（基于MR结果）
# SNP→暴露边
edges_snp_exposure = []
coef_matrix = mt_lasso.coef_  # (n_exposures, n_snps)
exposure_names = exposure_data_final.columns
snp_names = snp_genotype_final.columns

for i, exposure in enumerate(exposure_names):
    for j, snp in enumerate(snp_names):
        if snp in iv_snp_final and coef_matrix[i, j] != 0:
            edges_snp_exposure.append({
                "source": f"snp_{snp}",
                "target": exposure,
                "weight": abs(coef_matrix[i, j])
            })

# 暴露→结局边（基于MR显著结果）
mr_significant = mr_results[mr_results["pval"] < 0.05]
edges_exposure_outcome = [
    {"source": exp, "target": "lung_cancer", "weight": abs(beta)}
    for exp, beta in zip(mr_significant["exposure"], mr_significant["beta"])
]

# 合并边
edges = pd.DataFrame(edges_snp_exposure + edges_exposure_outcome)

# 2. 构建DGL图
# 节点ID映射
node_id_map = {node: idx for idx, node in enumerate(nodes["id"])}
edges["source_id"] = edges["source"].map(node_id_map)
edges["target_id"] = edges["target"].map(node_id_map)

# 创建DGL图
g = dgl.graph((edges["source_id"].values, edges["target_id"].values))
g.edata["weight"] = torch.tensor(edges["weight"].values, dtype=torch.float32)

# 3. 节点特征
# SNP特征：基因型均值
snp_features = snp_genotype_final[iv_snp_final].mean(axis=0).values
# 暴露特征：表达量均值
exposure_features = exposure_data_final.mean(axis=0).values
# 结局特征：1
outcome_features = np.array([1])

# 合并特征
node_features = np.concatenate([snp_features, exposure_features, outcome_features])
node_features = torch.tensor(node_features.reshape(-1, 1), dtype=torch.float32)  # (n_nodes, 1)

# 4. 构建GCN模型
class GCN(nn.Module):
    def __init__(self, in_feats, hidden_feats, out_feats):
        super(GCN, self).__init__()
        self.conv1 = dglnn.GraphConv(in_feats, hidden_feats, activation=nn.ReLU())
        self.conv2 = dglnn.GraphConv(hidden_feats, out_feats)
    
    def forward(self, g, features):
        x = self.conv1(g, features)
        x = self.conv2(g, x)
        return x

# 初始化模型
in_feats = node_features.shape[1]
hidden_feats = 64
out_feats = 1
model = GCN(in_feats, hidden_feats, out_feats)

# 5. 训练模型
labels = torch.ones(node_features.shape[0], dtype=torch.float32)  # 简化标签
optimizer = optim.Adam(model.parameters(), lr=0.01)
criterion = nn.MSELoss()

model.train()
for epoch in range(100):
    optimizer.zero_grad()
    logits = model(g, node_features)
    loss = criterion(logits, labels.unsqueeze(1))
    loss.backward()
    optimizer.step()
    if (epoch + 1) % 10 == 0:
        print(f"Epoch {epoch+1}, Loss: {loss.item():.4f}")

# 6. 提取节点重要性
node_importance = logits.detach().numpy().flatten()
nodes["importance"] = node_importance

# 7. 可视化因果网络
plt.figure(figsize=(16, 12))
nx_graph = nx.from_pandas_edgelist(edges, "source", "target", ["weight"])
# 添加节点属性
node_attrs = nodes.set_index("id")["type"].to_dict()
nx.set_node_attributes(nx_graph, node_attrs, "type")
node_importance_attr = nodes.set_index("id")["importance"].to_dict()
nx.set_node_attributes(nx_graph, node_importance_attr, "importance")

# 绘制图
pos = nx.kamada_kawai_layout(nx_graph)  # 布局
# 节点颜色
node_colors = {"snp": "blue", "gene": "green", "protein": "orange", "outcome": "red"}
colors = [node_colors[data["type"]] for _, data in nx_graph.nodes(data=True)]
# 节点大小
sizes = [data["importance"] * 1000 + 100 for _, data in nx_graph.nodes(data=True)]
# 边宽度
edge_weights = [data["weight"] * 2 for _, _, data in nx_graph.edges(data=True)]

nx.draw_networkx_nodes(nx_graph, pos, node_color=colors, node_size=sizes, alpha=0.8)
nx.draw_networkx_edges(nx_graph, pos, width=edge_weights, alpha=0.6)
nx.draw_networkx_labels(nx_graph, pos, font_size=8)
plt.title("Causal Network: MR + GNN Fusion Model", fontsize=16)
plt.axis("off")
plt.tight_layout()
plt.savefig("causal_network.pdf")

4.5 步骤 5：结果解释与生物学验证

4.5.1 可解释 AI（XAI）解析核心因果路径

利用 SHAP 值解释 AI 模型中各分子的因果贡献度，筛选核心调控节点：

import shap
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# 1. 准备模型输入数据（SNP+多组学暴露）
x_input = pd.concat([snp_genotype_final, exposure_data_final], axis=1)
y_output = sample_info_final["lung_cancer_status"]  # 假设样本信息中包含肺癌状态（0/1）

# 2. 训练可解释模型（以XGBoost为例）
import xgboost as xgb
model_xgb = xgb.XGBClassifier(n_estimators=100, random_state=42)
model_xgb.fit(x_input, y_output)

# 3. 计算SHAP值
explainer = shap.TreeExplainer(model_xgb)
shap_values = explainer.shap_values(x_input)

# 4. 可视化SHAP值
# 特征重要性图
plt.figure(figsize=(12, 8))
shap.summary_plot(shap_values, x_input, plot_type="bar")
plt.title("SHAP Feature Importance: Core Causal Factors", fontsize=14)
plt.tight_layout()
plt.savefig("shap_feature_importance.pdf")

# 依赖图：展示核心特征与结局的关系
core_feature = "ENSG00000123456"  # 假设排名第一的核心基因
plt.figure(figsize=(10, 6))
shap.dependence_plot(core_feature, shap_values, x_input)
plt.title(f"SHAP Dependence Plot: {core_feature} vs Lung Cancer Risk", fontsize=14)
plt.tight_layout()
plt.savefig("shap_dependence_plot.pdf")

# 5. 筛选核心因果分子（SHAP值绝对值前20）
shap_importance = pd.DataFrame({
    "feature": x_input.columns,
    "shap_importance": np.abs(shap_values).mean(axis=0)
}).sort_values("shap_importance", ascending=False).head(20)

print("Core Causal Factors (Top 20):")
print(shap_importance)

4.5.2 生物学验证建议

1. 功能实验验证核心因果关系：
   • 基因编辑：利用 CRISPR-Cas9 敲除 / 过表达核心基因，验证其对肺癌细胞增殖、迁移、侵袭的影响；
   • 蛋白质相互作用：通过 Co-IP、质谱分析验证核心蛋白质与下游分子的相互作用；
   • 动物模型：构建裸鼠移植瘤模型，验证核心分子对肿瘤生长的调控作用。
2. 临床样本验证：
   • 收集肺癌患者与健康人样本，检测核心分子的表达水平，验证其与疾病诊断、预后的相关性；
   • 回顾性队列研究：分析核心分子表达量与肺癌发病风险的因果关联。

五、最佳实践与常见问题解决方案

5.1 最佳实践

1. 数据质量控制：
   • 遗传变异数据：严格过滤低质量 SNP（P 值 <5e-8、样本量> 5000），确保工具变量的有效性；
   • 多组学数据：过滤缺失值 > 20% 的特征，进行标准化处理，减少批次效应影响；
2. 工具变量筛选：
   • 结合多种 AI 算法（LASSO、随机森林、XGBoost）筛选工具变量，提高稳定性；
   • 严格验证 MR 三大假设，尤其是排他性假设（可通过 MR-Egger 截距项检验、敏感性分析验证）；
3. AI 模型选择：
   • 低维数据：选择 LASSO、随机森林等传统机器学习算法，兼顾效率与可解释性；
   • 高维数据 / 网络数据：选择图神经网络（GNN）、深度学习模型，捕捉复杂调控关系；
4. 结果验证：
   • 采用多种 MR 方法（IVW、MR-Egger、Weighted Median）交叉验证因果关系；
   • 结合可解释 AI 技术（SHAP、注意力机制）解析模型，增强生物学可解释性；
   • 必须通过功能实验验证核心因果关系，形成 “数据分析→机制提出→实验验证” 的闭环。

5.2 常见问题解决方案

1. 工具变量与暴露相关性弱：
   • 扩大 SNP 筛选范围（如降低 GWAS P 值阈值至 1e-5）；
   • 采用多 SNP 联合预测（如 LASSO、随机森林），提高工具变量对暴露的预测力；
2. MR 结果存在异质性 / 水平多效性：
   • 异质性：采用 MR-Egger、Weighted Median 等稳健方法，或排除异质性强的 SNP；
   • 水平多效性：使用 MR-PRESSO 方法检测并校正异常值，或采用多工具变量 MR 分析；
3. AI 模型过拟合：
   • 数据增强：通过 Bootstrap 抽样扩充样本量；
   • 正则化：在 AI 模型中加入 L1/L2 正则化、Dropout 层；
   • 交叉验证：采用 10 折交叉验证评估模型稳定性；
4. 结果生物学可解释性差：
   • 结合生物学数据库（如 GO、KEGG、STRING）注释核心分子的功能；
   • 利用可解释 AI 技术（SHAP、CAM）解析模型，明确核心调控路径；
   • 优先选择在疾病相关通路中已被报道的分子进行验证。

六、总结与展望

孟德尔随机化（MR）作为因果推断的核心工具，能够有效排除混杂因素，揭示多组学暴露与疾病结局的因果关系；而 AI 技术则具备强大的高维数据处理和复杂关系建模能力，两者融合可实现 “因果推断 + 机制解析” 的双重目标。本文通过系统介绍 “MR + AI 融合建模” 的核心理论、实战流程，并提供 R/Python 双语言代码，帮助读者掌握多组学数据中的因果机制挖掘方法。

随着多组学数据的积累和 AI 技术的发展，未来融合建模将向以下方向发展：

1. 多模态 AI 模型：整合基因组、转录组、蛋白质组、影像组学等多模态数据，构建更全面的因果网络；
2. 动态因果建模：结合时间序列多组学数据，构建动态 MR-AI 模型，揭示疾病进展过程中的因果关系变化；
3. 临床转化应用：基于融合模型筛选疾病诊断、预后的 biomarker，开发靶向治疗药物，推动个性化医疗发展。