文章前言 | INTRODICTION
孟德尔随机化作为一种前沿的因果推断工具，为生物信息学和遗传学研究提供了强有力的方法论支持。在生物信息学和遗传学研究中，孟德尔随机化（Mendelian Randomization，MR）作为一种创新且有效的因果推断方法，正在帮助研究人员解开复杂的表型与疾病关系之谜。通过利用基因变异作为工具变量，孟德尔随机化提供了一种避免混杂因素和反向因果关系的新途径。
那么，什么是孟德尔随机化？它究竟有什么特别之处令其迅速走红？本篇公众号将带您深入了解这一方法，并手把手教您如何使用R包TwosampleMR进行实际分析。让我们一起走进孟德尔随机化的世界吧！
Part.01
孟德尔随机化解读
孟德尔随机化的概念源于19世纪遗传学的奠基人——格里高尔·孟德尔（Gregor Mendel）。孟德尔通过对豌豆植物的杂交实验，发现了基因遗传的基本规律，这为现代遗传学奠定了基础。现代的孟德尔随机化方法借鉴了他的遗传学理论，将基因变异视为一种“随机分配”的工具变量，用以推断因果关系。
随着基因组学和统计学的发展，孟德尔随机化在20世纪末和21世纪初得到了进一步的发展。著名统计学家Fisher发现，通过利用基因变异这一天然的“实验”，可以有效减少混杂因素和反向因果关系的干扰，从而获得更加可靠的因果推断结果。
CONCEPT
孟德尔随机化的核心是孟德尔第二定律：分离定律，即当DNA在配子形成时，基因对的不同等位基因会独立地分配到不同的配子中。这点类似于随机对照实验（Randomized controlled trials，RCT）中的随机分组。
例如，在血清钙水平对冠心病风险影响的研究中（图1），RCT受试者被随机分为两组，一组为对照组，另一组则是使用了钙补充剂的实验组。而根据孟德尔第二定律，MR研究人群根据是否继承血清钙水平相关的遗传变异，也被“随机化”分为了两组。
在排除混杂因素的条件下，如果RCT中的随机组平均血清钙水平较高的同时也具有较高的冠心病风险，则表明钙水平与冠心病风险存在因果关系。因此，按照同样的推理，如果MR发现影响血清钙水平的遗传变异与冠心病风险的差异有关，则它提供了血清钙对冠心病有因果影响的证据。
图1. 随机对照试验RCT和孟德尔随机化MR研究设计的比较（以血清钙水平对冠心病的因果关系为例）
3 HYPOTHESE
孟德尔随机化实际上是通过工具变量（Instrumental variables）来揭示事物之间的因果关系。在生物学研究中，这个工具变量是指遗传变异，即利用全基因组关联研究（Genome-wide association study，GWAS）数据SNP（Single-Nucleotide Polymorphism）作为工具变量来分析表型和疾病之间的因果关系。
为了确保结果的有效性，MR必须满足3个核心假设（图2）：
1)关联性：遗传变异必须与暴露因素X强相关；
2)独立性：遗传变异不能与任何可能的混淆因素U相关；
3)排他性：遗传变异不能与结局变量Y直接相关。
图2.孟德尔随机化3大假设
上述例子中，假如我们发现血清钙水平较高的患者更容易得冠心病，这只能表明高血清钙和冠心病同时出现的概率较大（相关性较强），不能证明血清钙水平高是冠心病的风险因素（相关性不代表因果关系）。因为可能存在其他因素导致冠心病的发生，或者也许冠心病才是导致血清钙水平高的原因。
对此，为了探究高血清钙（暴露因素X）是否会导致冠心病（结局变量Y）的发生，我们引入遗传变异SNP作为工具变量，探究高血清钙和冠心病的因果关系：
1) 关联性：弱工具变量的使用会导致估计出现偏移，所以遗传变异必须与高血清钙存在稳定的强相关关系。这在全基因组水平上表现为：该SNP在高血清钙上的p值小于5 × 10−8。
2) 独立性：遗传变异必须与高血压高血脂等其他因素相互独立，没有相关性。这在全基因组水平上表现为：除了高血清钙以外，该SNP在其他因素上的p值大于5 × 10−8。
3) 排他性：遗传变异不能与冠心病直接相关。如果该SNP能直接导致冠心病，那么相当于跳过了高血清钙这一暴露因素直接导致了结局，我们无法证明高血清钙和冠心病有因果关系。
Part.02
TwoSampleMR分析孟德尔随机化
在MR中，两样本孟德尔随机化（Two sample Mendelian randomisation，2SMR）是仅使用GWAS概括性数据来估计暴露对结局的因果效应的方法。下图为TwoSampleMR官方使用说明中的示意图，介绍了MR分析的工作流程和可能使用到的命令参数，简单概括为：
1) 选取与暴露因素显著相关的SNP，并保证它们的独立性（不与其他因素相关）；
2) 从结局变量的GWAS中提取上述SNP；
3) 对同一位点的暴露变量和结局变量的SNP数据进行预处理并整合，确保它们在同一方向上；
4) 进行孟德尔随机化分析、敏感性分析、结果可视化。
图3.TwoSampleMR分析流程
INSTALL
我们可以在R上安装TwoSampleMR包，需要从GitHub安装下载：
install.packages("remotes")remotes::install_github("MRCIEU/TwoSampleMR")
Part.03
GWAS数据获取
使用TwoSampleMR包进行孟德尔随机化分析时，通常需要用到两个GWAS数据集，一个是暴露因素的SNP数据，另一个是结局变量的SNP数据。GWAS数据可以从数据库或文献中获取，常用的数据库有：IEU OpenGWAS（图4）、GWAS-Catalog、FINNGEN、UK-Biobank，其中IEU数据库可以直接通过TwoSampleMR包进行下载。
图4. IEU OpenGWAS数据库
EXPOSURE SNP
例如，我们想分析上述例子中血清钙和冠心病之间的因果关系，为了便于演示，血清钙的SNP数据将从本地整理读取，冠心病SNP直接通过TwoSampleMR下载。
TwoSampleMR需要读取暴露因素SNP的data frame文件，每行为一个SNP，文件信息如下表所示（表1）：
表1. 暴露SNP文件信息
整理好文件后，使用read_exposure_data函数读取文件：
exposure <- read_exposure_data(filename = 'exposure.csv',sep = ",",snp_col = "SNP",beta_col = "beta.exposure",se_col = "se.exposure",effect_allele_col = "effect_allele.exposure",other_allele_col = "other_allele.exposure",eaf_col = "eaf.exposure",pval_col = "pval.exposure",chr_col = "chr.exposure",pos_col = "pos.exposure",samplesize_col = "samplesize.exposure")
读取暴露SNP后，并不是所有血清钙的SNP都能纳入后续分析，根据之前提到的孟德尔随机化三大假设，我们将对血清钙SNP进行筛选：
1) 根据MR的关联性假设，工具变量必须与暴露因素强相关，所以去除p值大于5 × 10−8的SNP：
#过滤掉pval大于5×10^−8的SNPexposure <- exposure[exposure$pval.exposure < 5e-08,]
2) 连锁不平衡（Linkage Disequilibrium，LD）是指遗传变异的非独立分离。同一染色体上邻近的遗传变异可以一起遗传，如果等位基因频率相似，则会导致它们之间存在相关性。为了确保SNP的独立性，我们需要去除连锁不平衡：
#去除连锁不平衡（去掉在10000kb范围内与r2大于0.001的SNP，默认参数）exposure <- clump_data(exposure, clump_kb = 10000, clump_r2 = 0.001, pop = "EUR")
3) 利用PhenoScanner网站去除混杂因素的影响，所选SNP的p值不得在其他因素上小于5 × 10−8。
4) 此外，还可通过筛选MAF >1%、F检验值>10等过滤条件，进一步去除弱工具变量。
OUTCOME SNP
过滤完暴露SNP后，我们现在要从结局GWAS数据中提取与筛选后的暴露数据相匹配的SNP，这里直接使用IEU数据库中冠心病的数据集：
outcome_data <- extract_outcome_data(snps = exposure$SNP, # 通过暴露SNP的rs id去提取结局的SNPoutcomes = 'ieu-a-7', # 冠心病在IEU数据库的编号proxies = FALSE, # 不使用代理SNPmaf_threshold = 0.01, # 最小等位基因频率（maf）阈值)
Part.04
孟德尔随机化分析
在获取了暴露因素和结局变量的SNP后，我们需要对两份数据进行预处理，以确保效应等位的一致性。这一步骤对于后续分析至关重要，因为它能确保我们在分析过程中比较的是同一种等位基因的效应。
HARMONIZATION
为了简化这一复杂的预处理过程，我们可以使用R包中的harmonise_data函数。该函数可以自动检查和调整暴露因素和结局数据中的效应等位，使其保持一致，从而避免了手动处理可能带来的错误和不一致。
data <- harmonise_data(exposure_dat = exposure, #暴露数据outcome_dat = outcome_data #结局数据)
具体来说，harmonise_data函数会对暴露数据和结局数据进行比对，并根据需要调整等位基因的方向和符号。这样，我们就能确保在后续分析中，两组数据之间的比较是有效且一致的（图5）。
图5. harmonise_data函数整合暴露和结局SNP
MR ANALYSIS
整合好上述数据后，就可以使用mr函数进行孟德尔随机化分析了：
result <- mr(data)result# id.exposure id.outcome outcome exposure method nsnp b se pval#1 mpEOi1 ieu-a-7 mpEOi1 || id:ieu-a-7 exposure MR Egger 9 0.9844 0.92156 0.320868#2 mpEOi1 ieu-a-7 mpEOi1 || id:ieu-a-7 exposure Weighted median 9 0.1575 0.12749 0.216598#3 mpEOi1 ieu-a-7 mpEOi1 || id:ieu-a-7 exposure Inverse variance weighted 9 0.7836 0.31012 0.011516#4 mpEOi1 ieu-a-7 mpEOi1 || id:ieu-a-7 exposure Simple mode 9 0.4765 0.19407 0.039591#5 mpEOi1 ieu-a-7 mpEOi1 || id:ieu-a-7 exposure Weighted mode 9 0.1523 0.08771 0.120797
mr函数默认使用5种分析方法：MR Egger、Weighted median、Inverse variance weighted、Simple mode和Weighted mode。
上述结果中，Inverse variance weighted和Simple mode的p值小于0.05，表明在这两种MR方法中，血清钙与冠心病存在因果关系。
可以使用mr_method_list()查看当前支持的统计方法。如果需要指定其他方法，可使用method_list参数：
mr_method_list()# obj name PubmedID Description use_by_default heterogeneity_test#1 mr_wald_ratio Wald ratio TRUE FALSE#2 mr_two_sample_ml Maximum likelihood FALSE TRUE#3 mr_egger_regression MR Egger 26050253 TRUE TRUE#4 mr_egger_regression_bootstrap MR Egger (bootstrap) 26050253 FALSE FALSE#5 mr_simple_median Simple median FALSE FALSE#6 mr_weighted_median Weighted median TRUE FALSE#7 mr_penalised_weighted_median Penalised weighted median FALSE FALSE#8 mr_ivw Inverse variance weighted TRUE TRUE#9 mr_ivw_radial IVW radial FALSE TRUE#10 mr_ivw_mre Inverse variance weighted (multiplicative random effects) FALSE FALSE#11 mr_ivw_fe Inverse variance weighted (fixed effects) FALSE FALSE#12 mr_simple_mode Simple mode TRUE FALSE#13 mr_weighted_mode Weighted mode TRUE FALSE#14 mr_weighted_mode_nome Weighted mode (NOME) FALSE FALSE#15 mr_simple_mode_nome Simple mode (NOME) FALSE FALSE#16 mr_raps Robust adjusted profile score (RAPS) FALSE FALSE#17 mr_sign Sign concordance test Tests for concordance of signs between exposure and outcome FALSE FALSE#18 mr_uwr Unweighted regression Doesn't use any weights FALSE TRUEresult_othermethod <- mr(data, method_list = c("mr_weighted_median", "mr_simple_mode", "mr_weighted_mode"))
Part.05
敏感性分析
敏感性分析是一种用于评估系统或模型对输入变量变化的响应的技术。在MR中，敏感性分析的目的是检验结果是否对某些假设或数据处理方法的变化敏感，从而验证因果推断是否稳固。具体而言，敏感性分析可以评估工具变量的有效性，检查是否有潜在的偏倚或混杂因素影响结果，评估不同模型的稳健性，以及确认是否存在某一个工具变量严重影响结局变量。
敏感性分析主要有异质性（heterogeneity）检验、水平多效性（horizontal pleiotropy）检验和留一法（leave-one-out）检验。
HETEROGENEITY
在MR中，理想的工具变量应该只通过暴露因素影响结局变量。如果工具变量不仅通过暴露因素，而且还通过其他途径影响结局变量，这种现象称为“水平多效性”或“水平多重效应”。水平多效性会导致工具变量违背孟德尔随机化3大假设之一的独立性，这会导致因果效应的估计出现偏倚。
# 异质性检验mr_heterogeneity(data)# id.exposure id.outcome outcome exposure method Q Q_df Q_pval#1 mpEOi1 ieu-a-7 mpEOi1 || id:ieu-a-7 exposure MR Egger 199.1175 7 1.76E-39#2 mpEOi1 ieu-a-7 mpEOi1 || id:ieu-a-7 exposure Inverse variance weighted 200.6679 8 4.62E-39# 也可以指定方法mr_heterogeneity(data, method_list = c("mr_ivw", "mr_simple_mode"))
如果Q_pval小于0.05说明存在异质性，这时候可以剔除结局数据中某些P值非常小的SNP，或者直接使用随机效应模型（mr_ivm_mre方法）来估计MR效应量。
PLEIOTROPY
在MR中，异质性检验旨在评估不同工具变量（遗传变异）对因果效应估计的影响是否一致。具体来说，它检测不同遗传变异在其对暴露因素和结局之间的因果关系中的作用是否存在显著差异。如果不同的工具变量对因果效应的估计存在显著差异，这可能表明工具变量在不同的生物学背景或环境中发挥着不同的作用，从而影响研究结论的解释。
# 水平多效性检验mr_pleiotropy_test(data)# id.exposure id.outcome outcome exposure egger_intercept se pval#1 mpEOi1 ieu-a-7 mpEOi1 || id:ieu-a-7 exposure -0.01420588 0.06084932 0.8220844
上述结果中，pval远大于显著性水平0.05，说明在我们的数据中没有找到显著的水平多效性。
LEAVE ONE OUT
留一法检验通过将每一个工具变量从分析中排除，然后使用剩余的工具变量重新计算因果效应估计，观察剔除每个工具变量后结果是否发生变化，来检查结果对单个工具变量的敏感性。它可以帮助识别和理解工具变量对因果效应估计的贡献，以及潜在的异常值或影响较大的工具变量。这种方法有助于提高因果推断的可靠性，确保研究结论不被单一工具变量过度影响。
# 留一法# 如果剔除了某一个SNP后，结果改变很大，说明存在某一个SNP对结果影响很大，这是我们不愿意看到的情况。# 理想的情况应该是剔除每个SNP后，总体的误差线变化不大result_loo <- mr_leaveoneout(data)plot <- mr_leaveoneout_plot(result_loo)plot
图6. 留一图
Part.06
可视化
最后，我们可以将MR结果通过可视化的方式进行展示。
DOT PLOT
使用mr_scatter_plot函数对上述MR分析结果绘制散点图：
# 散点图p1 <- mr_scatter_plot(result, data)p1# 选择两个具有显著性的结果result_add_method <- mr(data, method_list = c("mr_ivw", "mr_simple_mode"))p1_1 <- mr_scatter_plot(result_add_method, data)p1_1
图7. 散点图
FOREST PLOT
使用mr_singlesnp函数获得单个SNP的结果后，在使用mr_forest_plot函数绘制森林图：
# 森林图用来展示每个 SNP 的因果效应估计值及其置信区间，主要看总效应情况。# 同上，可用method_list指定方法result_single <- mr_singlesnp(data)p2 <- mr_forest_plot(result_single)p2
图8. 森林图
FUNNEL PLOT
使用mr_singlesnp函数获得单个SNP的结果后，在使用mr_funnel_plot函数绘制漏斗图：
# 漏斗图，用于展示异质性result_single <- mr_singlesnp(data)p3 <- mr_funnel_plot(result_single)p3
图9. 漏斗图
文章结语 | SUMMARY
孟德尔随机化作为一种前沿的因果推断工具，为生物信息学和遗传学研究提供了强有力的方法论支持。通过利用基因变异作为工具变量，MR不仅能够有效减少混杂因素和反向因果关系的干扰，还能在无需昂贵和时间耗费巨大的随机对照试验（RCT）的情况下，揭示暴露因素与结局之间的真实因果关系。这种方法的快速性、成本效益以及对罕见疾病研究的适用性，使其在现代科学研究中逐渐成为不可或缺的工具。
通过本文的介绍，无论您是生物信息学的研究人员，还是对生物研究感兴趣的学者，我们希望大家对孟德尔随机化有了一个清晰的认识，从基本概念到实际操作，都有了系统的了解。我们希望通过本文，您能顺利应用这些方法，深入挖掘数据中的潜在因果关系，为您的研究成果增添更多价值。如果您有任何问题或想要进一步探讨的内容，请随时联系我们。
关注纽科，我们将持续为您带来更多前沿的科研资讯和实用的分析技巧，一起探索更多可能！
参考文献：
Larsson, S. C., Butterworth, A. S., & Burgess, S. (2023). Mendelian randomization for cardiovascular diseases: principles and applications. European heart journal, 44(47), 4913-4924.
关于我们
纽科生物提供专业的生物信息学数据分析和高通量测序服务。目前，公司已经和四川大学、复旦大学、上海交通大学、中山医院、华中科技大学等多所医院、高校的研究团队建立了长期良好的合作关系，提供高品质的数据分析和测序服务，帮助客户在European Heart Journal、Circulation Research、Nature Communications等多个著名杂志期刊上发表高水平科研文章，欢迎各位老师前来咨询。
