当孟德尔随机化遇上中介分析：用Python+Statsmodels拆解疾病因果通路（避坑指南）-平芜编程栈

当孟德尔随机化遇上中介分析：用Python+Statsmodels拆解疾病因果通路（避坑指南）

在生物医学研究中，确定因果关系远比发现相关性更具挑战性。想象一下，你发现某种蛋白质水平与疾病风险显著相关——但这究竟是因为蛋白质导致了疾病，还是疾病状态影响了蛋白质表达？又或者两者都被某个隐藏因素所驱动？这正是孟德尔随机化（Mendelian Randomization, MR）结合中介分析大显身手的场景。

对于已经掌握基础MR技术的研究者来说，将这种方法扩展到多步因果链分析（如"基因→蛋白质→代谢物→疾病"）时，往往会遇到工具变量重叠、效应量传递偏差、共线性干扰等实际问题。本文将以Python技术栈为核心，结合pandas、statsmodels和TwoSampleMR等工具，带你构建一套可诊断、可调试的中介效应分析工作流。我们将重点解决以下痛点：

如何验证工具变量在中介分析中的跨层级有效性？
当直接效应和间接效应符号相反时，如何正确解释**"遮掩效应"**？
使用多变量回归控制混杂时，怎样避免方差膨胀导致的假阳性？

1. 工具变量的跨层级验证

中介分析的核心是建立"暴露→中介→结局"的因果链，而每个箭头都需要独立的工具变量支持。但在实际操作中，研究者常犯的错误是假设同一组SNP能完美服务于不同层级的分析。

1.1 工具变量的层级特异性检验

理想的工具变量应在每一环节都满足三大假设：

相关性：SNP与暴露/中介的强关联（F>10）
独立性：SNP与混杂因素无关联（Hansen's J检验p>0.05）
排他性：SNP仅通过目标变量影响下游（MR-Egger截距检验）

用Python实现自动化验证：

from TwoSampleMR import harmonise_data, mr import pandas as pd def validate_iv_strength(snps, exposure_df, mediator_df, outcome_df): # 第一步：验证暴露-中介环节(X→M) xm_data = harmonise_data(exposure_df[exposure_df['snp'].isin(snps)], mediator_df) xm_results = mr(xm_data, method='ivw') # 第二步：验证中介-结局环节(M→Y) my_data = harmonise_data(mediator_df[mediator_df['snp'].isin(snps)], outcome_df) my_results = mr(my_data, method='ivw') # 返回各环节F统计量和p值 return pd.DataFrame({ '环节': ['X→M', 'M→Y'], 'F值': [xm_results.f_statistic[0], my_results.f_statistic[0]], 'p值': [xm_results.p_value[0], my_results.p_value[0]] })

注意：当同一SNP在不同环节的效应方向不一致时（如X→M为正效应而M→Y为负效应），需检查等位基因对齐情况，这可能是链翻转(flip strand)问题导致的假象。

1.2 工具变量重叠的处理策略

当暴露和中介共享部分工具变量时，会导致效应量估计偏差。以下是三种常见场景的解决方案：

场景	问题表现	解决方案
完全独立IV	SNP_X∩SNP_M=∅	直接进行两阶段分析
部分重叠	SNP_X∩SNP_M≠∅	采用MVMR（多变量MR）控制交叉影响
完全重叠	SNP_X=SNP_M	需要引入第三方工具变量

使用MVMR控制交叉影响的示例代码：

import statsmodels.api as sm def run_mvmr(exposure_effects, mediator_effects, outcome_effects): # 准备设计矩阵 X = pd.DataFrame({ 'exposure_beta': exposure_effects, 'mediator_beta': mediator_effects }) X = sm.add_constant(X) # 添加截距项 y = outcome_effects # 检查方差膨胀因子(VIF) vif = pd.DataFrame() vif["变量"] = X.columns vif["VIF"] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])] # 拟合模型 model = sm.OLS(y, X).fit() return model.summary(), vif

2. 效应量传递与中介计算

中介效应的量化看似简单（β_indirect = β_XM × β_MY），但在实际应用中存在多个技术陷阱。

2.1 效应量标准化的一致性

不同数据库的效应量单位可能不同，需要进行标准化处理：

连续变量：转换为每标准差变化(SD)的效应

# 将OR值转换为对数尺度 df['beta'] = np.log(df['or']) # 计算每SD变化对应的beta df['beta_sd'] = df['beta'] / df['unit_sd']

二分类变量：使用log(OR)作为效应量

# 确保OR值大于0 assert (df['or'] > 0).all() df['beta'] = np.log(df['or'])

2.2 中介比例的非常规情况

当中介比例出现以下特殊值时，需要特别注意：

>100%：通常意味着存在遮掩效应(suppression effect)，即直接效应和间接效应方向相反
<0：提示中介因子可能起保护作用，抵消了暴露的部分风险
不显著但β_XM和β_MY均显著：可能是样本重叠导致的假阳性

计算中介效应及其置信区间的完整流程：

from scipy.stats import norm def bootstrap_mediation(xm_beta, xm_se, my_beta, my_se, n_bootstrap=1000): # 模拟抽样分布 xm_samples = norm.rvs(loc=xm_beta, scale=xm_se, size=n_bootstrap) my_samples = norm.rvs(loc=my_beta, scale=my_se, size=n_bootstrap) indirect_samples = xm_samples * my_samples # 计算95% CI ci_lower = np.percentile(indirect_samples, 2.5) ci_upper = np.percentile(indirect_samples, 97.5) return { '间接效应': xm_beta * my_beta, '95%CI下限': ci_lower, '95%CI上限': ci_upper }

3. 共线性诊断与解决方案

在中介分析中，暴露和中介变量常存在共线性，导致回归系数不稳定。以下是关键诊断指标：

3.1 方差膨胀因子(VIF)计算

from statsmodels.stats.outliers_influence import variance_inflation_factor def check_vif(exposure, mediator): X = pd.DataFrame({'exposure': exposure, 'mediator': mediator}) X = sm.add_constant(X) vif = pd.DataFrame() vif["变量"] = X.columns vif["VIF"] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])] return vif

经验阈值：VIF>5表示存在中度共线性，>10则需要采取矫正措施

3.2 共线性处理方案对比

方法	原理	适用场景	Python实现
主成分回归	将相关变量转换为正交成分	高度线性相关	`sklearn.decomposition.PCA`
岭回归	增加L2正则化约束	中等共线性	`sklearn.linear_model.Ridge`
弹性网络	L1+L2正则化组合	共线性+变量选择	`sklearn.linear_model.ElasticNet`
工具变量法	利用外生变量估计	存在有效工具变量	`linearmodels.iv`

以岭回归为例的代码实现：

from sklearn.linear_model import Ridge from sklearn.preprocessing import StandardScaler def ridge_adjustment(X, y, alpha=1.0): # 标准化数据 scaler = StandardScaler() X_scaled = scaler.fit_transform(X) # 拟合模型 model = Ridge(alpha=alpha) model.fit(X_scaled, y) # 返回标准化系数 return pd.DataFrame({ '变量': ['截距'] + X.columns.tolist(), '系数': [model.intercept_] + model.coef_.tolist() })

4. 结果解释与可视化

正确的统计结果需要结合生物学背景才能产生科学价值。以下是常见误区和解决方案：

4.1 中介效应方向解读

当中介分析结果出现以下模式时：

同号效应（β_XM×β_MY与β_XY同号）：典型的中介通路
异号效应：可能存在遮掩效应或竞争通路
中介比例>100%：提示存在未被测量的抑制因子

使用森林图展示各环节效应量：

import matplotlib.pyplot as plt def plot_mediation_results(xm_effect, my_effect, xy_effect): fig, ax = plt.subplots(figsize=(10, 6)) # 绘制效应量 effects = [xm_effect, my_effect, xy_effect] labels = ['X→M效应', 'M→Y效应', 'X→Y总效应'] ax.errorbar( x=effects, y=labels, xerr=[effect*0.2 for effect in effects], # 假设SE为效应量的20% fmt='o', capsize=5 ) # 标注中介比例 mediation_prop = (xm_effect * my_effect) / xy_effect * 100 ax.annotate( f'中介比例: {mediation_prop:.1f}%', xy=(xy_effect, 2), xytext=(5, 5), textcoords='offset points' ) ax.axvline(x=0, color='grey', linestyle='--') ax.set_title('中介效应分解结果') return fig

4.2 敏感性分析报告

完整的敏感性分析应包含以下要素：

异质性检验（Cochran's Q）

from statsmodels.stats.anova import anova_lm def cochran_q_test(residuals, predictors): model = sm.OLS(residuals**2, predictors).fit() return anova_lm(model)['F'][0], anova_lm(model)['Pr(>F)'][0]

水平多效性检验（MR-Egger截距）

def mregger_test(beta_exposure, beta_outcome, se_outcome): X = sm.add_constant(beta_exposure) model = sm.WLS(beta_outcome, X, weights=1/se_outcome**2).fit() return model.params[0], model.pvalues[0] # 返回截距和p值

留一法分析（Leave-one-out）

def loo_analysis(snps, beta, se): results = [] for i in range(len(snps)): mask = [True]*len(snps) mask[i] = False loo_beta = np.average(beta[mask], weights=1/se[mask]**2) results.append(loo_beta) return results

在实际分析中，我们常遇到工具变量数量不足的问题。这时可以考虑使用基因聚合评分（Gene Aggregate Score）方法，将同一基因区域的多个SNP合并：

def calculate_gas(snps, beta, genotypes): """ 计算基因聚合评分 :param snps: SNP ID列表 :param beta: 各SNP的效应量 :param genotypes: 样本基因型矩阵 (n_samples × n_snps) :return: 各样本的加权风险评分 """ assert len(snps) == len(beta) assert genotypes.shape[1] == len(snps) # 标准化基因型 (0,1,2 → -1,0,1) std_geno = (genotypes - 1) / 1.0 # 计算加权评分 gas = np.dot(std_geno, beta) return gas

最后要强调的是，统计上的中介效应只是假设生成工具，真正的因果机制需要实验验证。一个完整的研究闭环应该包含：