孟德尔随机化：脑卒中研究新钥匙
孟德尔随机化，作为一种基于遗传变异的因果推断方法，正逐渐成为医学研究领域的有力工具，在脑卒中研究中发挥着关键作用 。它的基本原理源于孟德尔遗传定律，即个体在受精过程中，等位基因会从父母那里随机分配给后代，这种随机分配类似于随机对照试验中的随机分组 。通过巧妙利用这一特性，孟德尔随机化将遗传变异作为 “自然发生” 的工具变量，以此来探索暴露因素与疾病之间的因果关系。
在传统的观察性研究中，混杂因素和反向因果关系是难以避免的难题 。比如，在研究生活习惯与脑卒中的关系时，年龄、性别、遗传背景、生活环境等众多混杂因素可能会干扰研究结果，使得我们难以准确判断生活习惯是否真正导致了脑卒中 。而且，疾病的发生往往是一个复杂的过程，很难确定是生活习惯导致了脑卒中，还是脑卒中影响了生活习惯，即存在反向因果关系 。而孟德尔随机化分析则能有效规避这些问题 。由于遗传变异在受孕时就已确定，不受后天环境和生活方式的影响，因此可以作为独立于其他因素的工具变量 。这就意味着，孟德尔随机化能够在一定程度上克服传统研究方法的局限性，为我们提供更为可靠的因果推断 。
以研究吸烟与脑卒中的因果关系为例，传统观察性研究可能会受到个体生活环境、职业压力、饮食结构等多种混杂因素的干扰 。而孟德尔随机化研究则可以选取与吸烟行为相关的遗传变异作为工具变量 。这些遗传变异在个体出生时就已确定，不会受到后天环境因素的影响 。通过分析这些遗传变异与脑卒中发生风险之间的关系，我们可以更准确地判断吸烟是否真的会增加脑卒中的发病风险 。再比如，在研究血脂水平与脑卒中的关系时，孟德尔随机化分析能够排除其他心血管危险因素的干扰，直接探究血脂水平对脑卒中的因果效应 。这种独特的优势使得孟德尔随机化在脑卒中病因学研究中具有重要价值，为我们深入了解脑卒中的发病机制提供了新的视角 。
搭建 Python 分析环境
在进行脑卒中孟德尔随机化分析之前，我们需要搭建一个合适的 Python 环境，确保代码能够顺利运行。以下是详细的环境搭建步骤：
安装 Python：Python 的最新版安装包可以在其官方网站下载，它可用于多个平台，如 Windows、Mac OS X 和 Linux 。根据自己电脑的操作系统进行选择，下载时不推荐选择最新版，可能会存在兼容问题，这里以 Windows 系统下 Python3.8.6 64bit 版本的安装包为例 。下载好安装包之后，双击打开，启动安装程序。在安装界面中，需要勾选 “添加 Python 到系统的环境变量中”，这一步至关重要，如果没有勾选，将不能在 cmd 运行窗口使用 python 和 pip 命令 。接着选择自定义安装，可以设置安装路径，安装路径不要包含汉字并且尽量不要包含空格 。进入自定义安装界面后，默认勾选的选项不用动，点击 “next” 进入下一步 。在后续界面中，勾选 “为所有用户安装”，设置好 Python 的安装路径，点击 “Install” 开始安装 。等待几分钟，安装成功后点击 “Close” 关闭页面即可 。
安装必要的库：
NumPy：这是一个用于科学计算的基础库，提供了高性能的多维数组对象和一系列用于处理数组的函数 。安装方法如下：打开命令提示符（Win+R 键，输入 cmd 回车），在命令提示符中输入 “pip install numpy”，然后回车，等待安装完成 。安装完成后，可以输入 “pip list” 查看 numpy 是否安装成功 。也可以在命令提示符中输入 python，进入 Python 交互环境，再输入 “import numpy”，如果没有报错，说明 numpy 库可以被调用，安装成功 。
pandas：这是一个强大的数据分析和处理库，提供了快速、灵活、明确的数据结构，旨在简单、直观地处理关系型、标记型数据 。安装命令为在命令提示符中输入 “pip install pandas” 。安装完成后，可在 Python 交互环境中输入 “import pandas as pd”，如果没有报错，再输入 “print (pd.version)”，若能显示出版本号，则表示 pandas 已成功安装 。
SciPy：这是一个基于 NumPy 的科学计算库，包含了优化、线性代数、积分、插值、特殊函数等模块 。在命令提示符中输入 “pip install scipy” 即可进行安装 。安装完成后，在 Python 交互环境中输入 “from scipy import some_module”（some_module 为 SciPy 中的某个模块，如 optimize、stats 等），若不报错，则安装成功 。
除了上述核心库之外，还可以根据具体需求安装其他辅助库，例如绘图库 matplotlib 用于结果可视化，seaborn 能提供更美观的绘图风格 。这些库将为我们进行脑卒中孟德尔随机化分析提供有力的支持，确保我们能够高效地处理和分析数据，挖掘出有价值的信息 。
代码实现全解析
数据准备
在开始孟德尔随机化分析之前，首先需要准备好相关数据，这些数据通常来自大规模的基因组关联研究（GWAS）。数据准备阶段的主要任务是读取和预处理数据，确保其适合后续的分析。
import pandas as pd

读取暴露因素数据，这里假设数据文件为exposure_data.csv，包含SNP、beta_exposure、se_exposure等列

exposure_data = pd.read_csv(‘exposure_data.csv’)

读取结局数据，假设数据文件为outcome_data.csv，包含SNP、beta_outcome、se_outcome等列

outcome_data = pd.read_csv(‘outcome_data.csv’)

数据清洗，例如去除缺失值

exposure_data = exposure_data.dropna()
outcome_data = outcome_data.dropna()

格式转换，将数据类型转换为合适的类型，如将SNP列转换为字符串类型

exposure_data[‘SNP’] = exposure_data[‘SNP’].astype(str)
outcome_data[‘SNP’] = outcome_data[‘SNP’].astype(str)

在上述代码中，我们使用pandas库的read_csv函数读取暴露因素数据和结局数据 。接着，通过dropna方法去除数据中的缺失值，以保证数据的完整性 。最后，将SNP列的数据类型转换为字符串类型，以便后续处理 。这些步骤都是为了确保数据的质量，避免因数据问题导致分析结果出现偏差 。
工具变量选择
工具变量的选择是孟德尔随机化分析的关键步骤，直接影响到分析结果的准确性和可靠性 。在脑卒中研究中，我们通常选择与血脂、载脂蛋白等相关的单核苷酸多态性（SNP）作为工具变量 。这些 SNP 与暴露因素（如血脂水平）密切相关，且满足孟德尔随机化的假设条件，即与混杂因素无关，且仅通过暴露因素影响结局 。

选择与血脂相关的SNP作为工具变量，假设这些SNP存储在一个列表中

lipid_snps = [‘rs12345’, ‘rs67890’, ‘rs54321’]

从暴露因素数据和结局数据中筛选出工具变量对应的行

instrumental_variables_exposure = exposure_data[exposure_data[‘SNP’].isin(lipid_snps)]
instrumental_variables_outcome = outcome_data[outcome_data[‘SNP’].isin(lipid_snps)]

在这段代码中，我们首先定义了一个包含与血脂相关的 SNP 的列表lipid_snps 。然后，使用isin方法从暴露因素数据和结局数据中筛选出这些 SNP 对应的行，得到工具变量数据集instrumental_variables_exposure和instrumental_variables_outcome 。这样，我们就完成了工具变量的选择过程 。在实际应用中，选择工具变量时需要综合考虑多个因素，如 SNP 与暴露因素的关联强度、SNP 在人群中的频率等 。通常，我们会选择与暴露因素在传统的全基因组显著性水平上相关（P < 5×10-8）的 SNP 作为工具变量 。同时，还需要通过生物学知识和统计方法验证这些 SNP 是否满足孟德尔随机化的假设条件，以确保工具变量的有效性 。
分析方法实现
加权逆方差法（IVW）是孟德尔随机化分析中常用的方法之一，它通过对每个工具变量的效应估计值进行加权，来综合评估暴露因素与结局之间的因果关系 。下面是使用 Python 实现加权逆方差法的核心代码：
import numpy as np

计算每个工具变量的权重，权重为1/标准误的平方

weights = 1 / instrumental_variables_exposure[‘se_exposure’]**2

计算加权后的暴露效应估计值

weighted_beta_exposure = instrumental_variables_exposure[‘beta_exposure’] * weights

计算加权后的结局效应估计值

weighted_beta_outcome = instrumental_variables_outcome[‘beta_outcome’] * weights

计算因果效应估计值，即加权后的结局效应估计值之和除以加权后的暴露效应估计值之和

causal_effect = np.sum(weighted_beta_outcome) / np.sum(weighted_beta_exposure)

计算因果效应的标准误

se_causal_effect = np.sqrt(1 / np.sum(weights))

在上述代码中，我们首先计算每个工具变量的权重，权重为 1 / 标准误的平方，这样可以使效应估计值更精确的工具变量在分析中具有更大的权重 。接着，分别计算加权后的暴露效应估计值和结局效应估计值 。然后，通过将加权后的结局效应估计值之和除以加权后的暴露效应估计值之和，得到因果效应估计值 。最后，计算因果效应的标准误，用于评估因果效应估计值的可靠性 。这些步骤完整地实现了加权逆方差法的计算过程，通过合理地运用权重，能够更准确地推断暴露因素与结局之间的因果关系 。
结果解读
代码运行完成后，我们需要对结果进行解读，以获取有价值的信息 。在孟德尔随机化分析中，关键的结果指标包括因果关系强度和统计学显著性 。

设定显著性水平

alpha = 0.05

判断因果关系是否显著

if causal_effect > 0 and se_causal_effect < causal_effect / (1.96 * np.sqrt(len(weights))):
print(f"暴露因素与脑卒中之间存在正向因果关系，因果效应估计值为{causal_effect}，在{alpha}水平上显著。“)
elif causal_effect < 0 and se_causal_effect < -causal_effect / (1.96 * np.sqrt(len(weights))):
print(f"暴露因素与脑卒中之间存在负向因果关系，因果效应估计值为{causal_effect}，在{alpha}水平上显著。”)
else:
print(f"未发现暴露因素与脑卒中之间存在显著的因果关系，因果效应估计值为{causal_effect}。")

在这段代码中，我们首先设定了显著性水平alpha为 0.05 。然后，根据因果效应估计值和标准误的关系来判断因果关系是否显著 。如果因果效应估计值大于 0，且标准误小于因果效应估计值除以 1.96 倍的权重平方根（对应 95% 置信区间），则表明暴露因素与脑卒中之间存在正向因果关系且在设定的显著性水平上显著 。反之，如果因果效应估计值小于 0，且标准误满足相应条件，则表明存在负向因果关系且显著 。如果不满足上述条件，则说明未发现显著的因果关系 。通过这样的判断，我们可以从分析结果中准确地解读出暴露因素与脑卒中之间的因果关系情况，为进一步的研究和决策提供有力的依据 。
案例实操：加深理解
为了让大家更深入地理解孟德尔随机化分析在脑卒中研究中的应用，我们来看一个实际案例。假设我们要研究血脂水平与脑卒中之间的因果关系，以血脂相关的 SNP 作为工具变量 。
数据准备
import pandas as pd

读取暴露因素数据，这里假设数据文件为exposure_data.csv，包含SNP、beta_exposure、se_exposure等列

exposure_data = pd.read_csv(‘exposure_data.csv’)

读取结局数据，假设数据文件为outcome_data.csv，包含SNP、beta_outcome、se_outcome等列

outcome_data = pd.read_csv(‘outcome_data.csv’)

数据清洗，例如去除缺失值

exposure_data = exposure_data.dropna()
outcome_data = outcome_data.dropna()

格式转换，将数据类型转换为合适的类型，如将SNP列转换为字符串类型

exposure_data[‘SNP’] = exposure_data[‘SNP’].astype(str)
outcome_data[‘SNP’] = outcome_data[‘SNP’].astype(str)

在这个案例中，我们从exposure_data.csv文件中读取了暴露因素数据，其中beta_exposure表示每个 SNP 与血脂水平关联的效应估计值，se_exposure表示对应的标准误 。从outcome_data.csv文件中读取了结局数据，beta_outcome表示每个 SNP 与脑卒中关联的效应估计值，se_outcome表示对应的标准误 。通过数据清洗和格式转换，确保数据质量和格式符合后续分析要求 。
工具变量选择

选择与血脂相关的SNP作为工具变量，假设这些SNP存储在一个列表中

lipid_snps = [‘rs12345’, ‘rs67890’, ‘rs54321’]

从暴露因素数据和结局数据中筛选出工具变量对应的行

instrumental_variables_exposure = exposure_data[exposure_data[‘SNP’].isin(lipid_snps)]
instrumental_variables_outcome = outcome_data[outcome_data[‘SNP’].isin(lipid_snps)]

我们事先确定了与血脂相关的 SNP 列表lipid_snps ，然后从暴露因素数据和结局数据中筛选出这些 SNP 对应的行，得到工具变量数据集instrumental_variables_exposure和instrumental_variables_outcome 。这一步确保了我们后续分析所使用的工具变量是与血脂密切相关的，符合孟德尔随机化分析的要求 。
分析方法实现
import numpy as np

计算每个工具变量的权重，权重为1/标准误的平方

weights = 1 / instrumental_variables_exposure[‘se_exposure’]**2

计算加权后的暴露效应估计值

weighted_beta_exposure = instrumental_variables_exposure[‘beta_exposure’] * weights

计算加权后的结局效应估计值

weighted_beta_outcome = instrumental_variables_outcome[‘beta_outcome’] * weights

计算因果效应估计值，即加权后的结局效应估计值之和除以加权后的暴露效应估计值之和

causal_effect = np.sum(weighted_beta_outcome) / np.sum(weighted_beta_exposure)

计算因果效应的标准误

se_causal_effect = np.sqrt(1 / np.sum(weights))

在这部分代码中，我们按照加权逆方差法的原理进行计算 。首先计算每个工具变量的权重，权重与标准误的平方成反比，这意味着标准误越小，权重越大，即效应估计值更精确的工具变量在分析中具有更大的权重 。接着分别计算加权后的暴露效应估计值和结局效应估计值 。通过加权处理，综合考虑了每个工具变量的信息 。最后，计算出因果效应估计值和标准误 。
结果解读

设定显著性水平

alpha = 0.05

判断因果关系是否显著

if causal_effect > 0 and se_causal_effect < causal_effect / (1.96 * np.sqrt(len(weights))):
print(f"暴露因素与脑卒中之间存在正向因果关系，因果效应估计值为{causal_effect}，在{alpha}水平上显著。“)
elif causal_effect < 0 and se_causal_effect < -causal_effect / (1.96 * np.sqrt(len(weights))):
print(f"暴露因素与脑卒中之间存在负向因果关系，因果效应估计值为{causal_effect}，在{alpha}水平上显著。”)
else:
print(f"未发现暴露因素与脑卒中之间存在显著的因果关系，因果效应估计值为{causal_effect}。")

运行上述代码后，我们设定显著性水平alpha为 0.05 。根据因果效应估计值和标准误的关系来判断因果关系是否显著 。如果因果效应估计值大于 0，且标准误满足一定条件，表明血脂水平与脑卒中之间存在正向因果关系，即血脂水平升高会增加脑卒中的发病风险 。反之，如果因果效应估计值小于 0 且标准误满足相应条件，则存在负向因果关系 。若不满足上述条件，则说明在当前数据和分析方法下，未发现血脂水平与脑卒中之间存在显著的因果关系 。通过这个案例实操，我们完整地展示了孟德尔随机化分析的过程，希望能帮助大家更好地掌握这一分析方法在脑卒中研究中的应用 。