R语言的流行病学应用 课件 第18章 孟德尔随机化

第18章孟德尔随机化18.1方法概论18.2基于汇总数据的MR方法18.3示例实践18.4练习案例18.5小结18.1方法概述18.1方法概述18.1.1研究背景1．观察性流行病学推断因果关系的局限性观察性流行病学设计相对简单且易于实施，被广泛用于疾病病因的初步探索，如吸烟与肺癌的关系，但许多观察到的与疾病风险相关的暴露在随机对照试验（randomizedcontrolledtrials，RCT）中未被证实与疾病存在因果关联。观察性流行病学研究在进行因果推断时容易受到选择偏倚、混杂因素和反向因果关系的影响。例如，在补充维生素E与全因死亡的观察性研究中，补充维生素E会受到社会经济地位的影响，而社会经济地位与全因死亡相关。这些因素影响着观察性流行病学在推断因果关系时的可靠性。18.1方法概述18.1.1研究背景2．RCT推断因果关系的局限性RCT通常被认为是进行因果推论的黄金标准。然而在现实中，对一些可能是病因的暴露（如吸烟、饮酒）进行人为干预的随机分配是不道德的，伦理问题使得用RCT研究这些暴露因素几乎不可能。此外，对于需要研究长期暴露和结局关联的RCT，如饮食和疾病的因果推断，需要投入巨大的研究时间和成本，而受试者的依从性也会影响RCT的结果。18.1方法概述18.1.1研究背景3．遗传流行病学推断因果关系：孟德尔随机化孟德尔随机化（mendelianrandomization，MR）研究使用与暴露因素具有强相关关系的遗传变异作为工具变量，推断暴露因素与结局之间的因果关系。遗传变异继承自父母，出生后保持不变，其与结局之间的关联具有时间顺序合理性。此外，遗传变异不会受到环境、社会经济地位、行为因素等传统混杂因素的影响。因此，MR能够降低混杂和反向因果对因果推断的影响。含有不同等位基因的个体可以认为被随机分配了不同的暴露水平，这一过程与RCT中的随机分组类似，形成天然的暴露组和对照组。因此，MR法也被称为“天然的RCT”。但与传统RCT相比，MR研究可以在具有代表性的人群样本中进行，不需要排除标准，也不需要随机分配和干预。18.1方法概论18.1.2MR基本假设MR的本质是利用遗传变异作为工具变量进行因果推断。MR基本假设如图18.1所示，在MR分析中选择的工具变量G需要满足3个基本假设。

（1）关联性：工具变量G与风险因素（暴露）X存在关联；

（2）独立性：工具变量G独立于风险因素X和结局Y之间的混杂因素U；

（3）排他性：工具变量G只通过风险因素X途径来影响结局Y的风险。

图18.1

MR基本假设18.1方法概述18.1.3MR影响因素工具变量满足上述3个基本假设是MR研究进行因果推断的前提。然而，弱工具变量（weakinstrument）、人群分层（populationstratification）、遗传变异的水平多效性（horizontalpleiotropy）和连锁不平衡（linkagedisequilibrium，LD）等因素可能导致工具变量违背MR基本假设，从而导致MR因果推断的结果不稳健。18.1方法概述

18.1方法概述18.1.3MR影响因素2．人群分层人群分层是指遗传变异频率和发病率（或对性状的贡献）在不同遗传背景的人群间存在差异，导致遗传变异与结局之间出现虚假关联。在MR研究中，人群分层可作为混杂因素，导致工具变量与结局的虚假关联，违反MR第2个基本假设，继而导致错误的因果推断。在研究工具变量与暴露、工具变量与结局的关联时，可以通过控制样本来自同一种群的人群并校正潜在的人群结构或者分层分析，控制人群分层的影响。18.1方法概述18.1.3MR影响因素3．水平多效性多效性（pleiotropy）是指一个遗传变异与多个表型相关，可以分为垂直多效性（verticalpleiotropy）和水平多效性（horizontalpleiotropy）。垂直多效性也被称为中介多效性（mediatedpleiotropy）或因果效应，如图18.2所示，虽然工具变量G除了与暴露X1关联，还与暴露X2相关，但暴露X2位于暴露X1与结局的通路上，反映了暴露X1的下游效应，是将暴露X1与结局联系的因果途径。图18.2

垂直多效性示意图18.1方法概述18.1.3MR影响因素3．水平多效性水平多效性也被称为非中介多效性（unmediatedpleiotropy）或生物多效性（biologicalpleiotropy），如图18.3所示。水平多效性可通过“遗传变异-暴露-结局”以外的通路影响结局发生，违背排他性假设，也可能导致独立性假设不成立。假设使用遗传变异G作为工具变量来探索暴露X1与结局之间的因果关联，暴露X1与结局无关联，但遗传变异G可通过暴露X2影响结局，工具变量G违反了排他性假设，可能得到暴露X1与结局存在因果关系的错误结论。对于多效性问题，可以筛选具有生物学意义的工具变量，并可以采用不同的方法进行分析，以期得到稳健的结果。图18.3

水平多效性示意图18.1方法概述18.1.3MR影响因素3．连锁不平衡连锁不平衡是指群体内不同位点等位基因间的非随机关联。假设使用遗传变异G1作为工具变量来探索暴露X1与结局之间的因果关联，如图18.4所示，遗传变异G1与遗传变异G2之间存在连锁不平衡，并且遗传变异G2可通过暴露X2影响结局；G1作为工具变量，违背了MR的排他性假设，可能得到暴露X1与结局存在因果关系的错误结论。对于LD存在的问题，与多效性问题类似，也可以通过筛选具有生物学意义的工具变量，并采用不同的方法进行分析，以期得到稳健的结果。图18.4

连锁不平衡示意图18.1方法概论18.1.4MR分类全基因组关联分析（genome-wideassociationstudy，GWAS）是指在人类全基因组范围内找出与表型存在关联的遗传变异，如单核苷酸多态性（single-nucleotidepolymorphism，SNP）。根据需要使用GWAS个体数据还是汇总数据，MR可分为基于个体数据的MR和基于汇总数据的MR。基于MR研究中暴露与结局的关系，可以分为线性MR和非线性MR。基于MR研究只研究一个暴露还是同时研究多个暴露与结局的关联，可以分为单变量MR和多变量MR。基于MR研究中“基因-暴露”和“基因-结局”关联的样本是否分别来源于两个独立样本，可以分为单样本MR和两样本MR。18.1方法概论18.1.4MR分类GWAS汇总数据易于从公共数据库中获得，单样本MR局限于单个数据来源，检验功效较低，工具变量的选择也比较局限，容易受到潜在混杂因素的影响。而两样本MR的GWAS与暴露、GWAS与结局的关联数据来源于两个独立样本，样本量更大，相对于单样本MR可以获得更大的把握度，因此基于汇总数据的两样本MR被广泛使用。18.2基于汇总数据的MR方法18.2基于汇总数据的MR方法MR研究的步骤包括：（1）选择工具变量，得到工具变量对暴露的效应，只有存在与暴露相关的遗传变异的情况下，才可以进行MR研究，可以从已有的GWAS数据中获取；（2）估计工具变量对结局的作用，工具变量对结局的作用也可以从已有的GWAS数据获取，可以是从所有的研究中估计出来的整体效应，这样可以拒绝单个研究的偏倚；（3）估计暴露和结局的因果效应，可以是单个工具变量的效应，也可以是合并多个工具变量的效应；

（4）进行敏感性分析，检验MR结果的稳健性。18.2基于汇总数据的MR方法基于汇总数据的MR研究中，用于估计因果效应的方法包括Wald方法、逆方差加权法（inversevarianceweighted，IVW）、MR-Egger方法、加权中位数法（weightedmedian）、简单众数法（simplemode）、加权众数法（weightedmode）、稳健校正轮廓得分法（robustadjustedprofilescore，RAPS）、基于汇总数据的广义孟德尔随机化法（generalizedsummary-data-basedMendelianrandomization，GSMR）等。18.2基于汇总数据的MR方法

18.2基于汇总数据的MR方法

18.2基于汇总数据的MR方法

18.2基于汇总数据的MR方法18.2.2逆方差加权法IVW方法的特点是用每个工具变量因果效应估计的方差的倒数作为权重来合并多个工具变量估计的因果效应量，从而得到整合多个工具变量的效应量，只考虑了工具变量与结局之间的估计误差，不考虑工具变量与暴露之间的效应量估计误差。IVW方法的应用前提是所有的工具变量均为有效的工具变量（服从MR的3个假设）且相互之间完全独立，此时IVW方法的统计功效最大。但在实际应用中，工具变量可能存在多效性（违背第3个基本假设），因此该方法在风险因素的GWAS数据样本量较小和工具变量存在多效性时会得到有偏的估计。

18.2基于汇总数据的MR方法

18.2基于汇总数据的MR方法18.2.3MR-Egger方法MR-Egger方法允许所有遗传变异均具有多效性，但是需满足遗传变异多效性对结局的效应独立于遗传变异对暴露因素的效应（instrumentstrengthindependentofdirecteffects，InSIDE）假设。当InSIDE假设成立时，MR-Egger回归可以准确估算因果效应，但是结果的精确度和统计学检验效能将有所降低。但当多效性效应通过暴露和结果的混杂因素起作用时，就会违反InSIDE假设，导致有偏的因果效应估计。

18.2基于汇总数据的MR方法18.2.4加权中位数法当所有的工具变量都是有效的工具变量时（即满足MR的3个基本假设），IVW是有效的MR估计方法，但是即使只有一个无效工具变量，都会导致IVW的估计存在偏倚。因此，已有学者提出基于工具变量的中位效应估计，得到暴露与结局的因果效应估计。基于中位数的方法允许一些（但不是全部）工具变量是无效工具。首先通过Wald方法计算每个工具变量的因果估计值，在未加权中位数（简单中位数）方法中，每个工具变量在分析中获得相等的权重，所有工具变量因果估计值的中位值即简单中位数的因果效应估计。该方法最突出的优势是在高达50%无效工具变量存在的情况下，仍可得到稳健的因果估计。但当单个估计值的精度差异很大时，简单中位数的估计是不恰当的。

18.2基于汇总数据的MR方法

18.2基于汇总数据的MR方法

18.2基于汇总数据的MR方法18.2.6小结：MR方法比较上述几种方法各有优劣，表18.1从是否违反MR基本假设和优缺点等方面对MR方法进行比较。在进行MR研究时，可以采用多种MR方法进行分析，若多种方法因果推断的结论一致，则可认为MR研究的结果稳健而可靠。

表18.1

MR方法比较MR方法放宽MR基本假设优劣势Wald方法无只适用于单个工具变量效应估计，易受到弱工具变量的影响逆方差加权法无统计功效高，易受到工具变量违反MR假设的影响MR-Egger方法允许工具变量存在水平多效性工具变量存在水平多效性时结果较稳健；易受到离群值的影响，易受到违背InSIDE假设的影响，统计功效较低加权中位数法允许少部分工具变量违反MR假设（如多效性）工具变量存在离群值或水平多效性时，MR结果仍然稳健；对增加或减少工具变量敏感，统计功效较低众数法允许大部分工具变量违反MR假设工具变量存在离群值或水平多效性时，MR结果仍然稳健；对带宽参数敏感，相对保守18.3示例实践18.3示例实践首先介绍TwoSampleMR程序包，它拥有三大优势，使之几乎可以满足所有的MR研究需求。（1）提供了大量的计算方法，除了基本IVW方法、MR-Egger方法、加权中位数法和众数法，还提供

了RAPS等算法。（2）提供大量已发表的GWAS数据，使我们可以从公开数据库中快速获取GWAS研究的结果，省时省力。

(3)灵活度非常高，不仅支持用户输入自己的暴露或者结局的数据，还提供了各种方法将TwoSampleMR程序包的分析结果转换为其他R语言程序包可以直接使用的形式。最重要的是，它提供了一套较为完整的分析流程，以帮助用户更快速地掌握MR的研究方法，主要包括读取暴露文件、去除连锁不平衡、提取工具变量在结局文件中的信息、将暴露和结局的效应等位基因统一、进行MR分析和敏感性分析6步。

首先安装TwoSampleMR程序包

（演示）。

18.3示例实践18.3.1暴露数据首先需要提供暴露汇总数据，MR分析所需的默认列名如表18-2所示，其中，前4行是必须的列：表18-2

MR分析所需要的默认列表名列名含义SNP工具变量名称，通常为RSIDbeta效应大小。如果表型是二元的，那么应该使用log(OR)se效应大小的标准误差effect_allele效应等位基因other_allele非效应等位基因eaf效应等位基因频率Phenotype表型chrSNP对应的染色体positionSNP对应的物理位置samplesize样本量ncase病例样本数ncontrol对照样本数pvalP值geneSNP的基因或其他注释18.3示例实践18.3.1暴露数据1．读取数据可以直接从文件读取数据，也可以从已定义的数据，或已有数据库中读取数据，并将其转换为暴露数据对应的格式。1）从文件读取数据若文件的列名与表18-2中的相对应，且为文本文件格式，则可以直接使用read_exposure_data()读取数据（演示）。该函数自动检查数据类型是否包含必要的列。如果文件没有默认列名，则仍可用read_exposure_data()读取文件，但需指定对应的列名（演示）。18.3示例实践18.3.1暴露数据1．读取数据2）从已定义的数据读取如果数据已在R语言中的数据库中存在，则可以使用format_data()将已有的数据进行转换（演示）。如果数据列名不是默认列名，则可以添加snp_col等参数指定对应列名，类似read_exposure_data的对应参数。3）从已有数据库中获取工具变量MRInstruments程序包提供了许多整理好、可直接使用的工具变量（演示）。通过data()访问数据（演示）。available_outcomes()返回数据库中所有可用的研究，可以用head(ao)查看。每个研究对应一个唯一的ID，可以通过指定ID提取需要的研究（演示）。18.3示例实践

18.3示例实践18.3.2结局数据结局数据需要的列与暴露数据类似，在确定暴露因素的工具变量后，需要从结局数据中提取相对应的变量，并将其整理为结局数据对应的格式。1．从IEUGWAS数据库中提取结局数据以身体质量指数（bodymassindex，BMI）与冠心病（coronaryheartdisease）因果关联研究为例，我们需要识别影响BMI的工具变量，并从冠心病的GWAS中提取这些SNP（演示）。CARDIoGRAMplusC4D机构冠心病GWASID为ieu-a-7，我们可以从该GWAS中提取BMI的工具变量（演示）。18.3示例实践

18.3示例实践18.3.2结局数据1．从IEUGWAS数据库中提取结局数据extract_outcome_data()返回工具变量与结局的关联，数据格式的主要列名如下：SNP、beta.outcome、se.outcome、samplesize.outcome、ncase.outcome、ncontrol.outcome、pval.outcome、eaf.outcome、effect_allele.outcome、other_allele.outcome、units.outcome、outcome、consortium.outcome、year.outcome、pmid.outcome、id.outcome、originalname.outcome、proxy.outcome、target_snp.outcome、proxy_snp.outcome、target_a1.outcome、target_a2.outcome、proxy_a1.outcome、proxy_a2.outcome、mr_keep.outcome和data_source.outcome。18.3示例实践18.3.2结局数据2．从本地GWAS汇总数据中提取结局数据假设有一个名为gwas_summary.csv的GWAS汇总数据文件，列名为rsid、effect、SE、a1、a2、a1_freq、p-value、Units、Gene、n，则可以用read_outcome_data()读取结局数据，参数与read_exposure_data()的类似（演示）。3．基于已有数据进行转换如果结局数据已在R语言中的数据库中存在，类似暴露数据，则也可以使用format_data()将已有的数据转换为上述示例中的格式，特别地，type="outcome"，还需要通过指定snps=bmi_exp_dat$SNP提取工具变量。18.3示例实践18.3.3数据预处理基于18.3.1节和18.3.2节得到的暴露和结局数据，需要将暴露和结局的效应等位基因整理为相同的等位基因，并进行合并（演示）。harmonise_data()可以自动处理暴露和结局数据效应等位基因相反、正反链和等位基因为A/T或C/G的SNP问题。18.3示例实践18.3.4进行MR分析1．主要的MR分析继续以18.3.3节得到的dat为例，运行mr(dat)，这将返回一个数据框，包含了一系列不同MR方法估计暴露对结局的因果效应。此外，可以通过mr_method_list()查看MR方法，如果要执行特定的MR方法，则可以在mr()中通过method_list参数指定（演示）。注意：MR假设工具变量首先影响暴露，然后通过暴露影响结局。但有时这很难评估暴露和结局的因果方向，因此可以使用Steiger检验进行测试（演示）。18.3示例实践18.3.4进行MR分析1．主要的MR分析它通过估计SNP对暴露和结局中解释的方差，并检验工具变量对结局的方差解释是否小于暴露。但该测试在暴露和结局中的某些测量误差参数下容易给出不准确的因果方向（如果结局的测量精度低得多，则其所解释的方差比例将被低估）。我们可以通过两种方式来评估推断的因果方向容易产生测量误差的程度。（1）提供暴露和结局的测量误差估计，并获得因果方向的调整估计；（2）对于所有可能的测量误差值，确定支持推断因果方向的参数空间比例。18.3示例实践

(18.7)18.3示例实践18.3.4进行MR分析2．敏感性分析2）水平多效性检验MR-Egger回归中的截距项可以有效指示水平多效性是否使MR结果产生了偏倚（演示）。3）每次剔除一个工具变量每次剔除一个工具变量，使用剩下的工具变量重新进行MR分析，识别是否由单个工具变量导致MR研究的因果关联（演示）。在默认情况下，使用的MR方法是IVW方法，但可以通过method参数更改。18.3示例实践18.3.4进行MR分析3．可视化通过图像展示MR分析的结果。1）散点图可以使用散点图来描述SNP对暴露的影响与SNP对结局的影响之间的关系（演示）。2）森林图使用mr_forest_plot()可将使用不同MR方法的MR估计值与单个工具变量结果进行比较（演示）。要想获取使用不同方法的绘图，可以在mr_singlesnp()中指定（演示）。18.3示例实践18.3.4进行MR分析3．可视化3）每次剔除一个工具变量的MR分析首先通过mr_leaveoneout()计算每次剔除一个工具变量的MR结果（可以用method参数指定MR方法），然后用mr_leaveoneout_plot()展示，检测是否由于某个特定的工具变量导致因果关联的结论（演示）。4）漏斗图基于单个工具变量结果生成漏斗图，漏斗图中的对称性可以衡量MR分析的可靠性。横坐标展示每个工具变量的因果效应估计，以及整合所有工具变量的IVW和MR-Egger估计结果；纵坐标表示每个工具变量因果估计误差的倒数（演示）。18.4练习案例18.4练习案例请利用消化道溃疡疾病(pepticulcerdisease，PUD)和慢性肾疾病(chronickidneydisease，CKD)的汇总数据进行MR研究，以推断PUD是否导致CKD。PUD的汇总数据基于UKBiobank456327个欧洲样本(16666/439661，病例/对照)的GWAS；CKD的汇总数据基于CKDGenConsortium480,698个欧洲样本(41,395/439,303病例/对照)的GWAS。