【方法学特辑】孟德尔随机化（Mendelian Randomization）保姆级教程！！！
孟德尔随机化（Mendelian Randomization, 以下简称MR）保姆级教程
写在前面
这是一篇超详细超啰嗦的孟德尔随机化（Mendelian Randomization）教程。
因为字太多，所以可能看一半就累坏了，又因为字太多，囊括了小编组所有成员在实践中踩过的坑和走过的弯路，希望能给第一次接触孟德尔随机化的你有所帮助！
全文分为十一个章节，这里列出来当目录用，如果哪天突然忘记了一些相关细节需要回顾，照着目录查也会方便一些~
一、 MR的概念
二、 MR的原理
三、 核心假设违背的影响与检测方法
四、 工具变量的选择与验证
五、 因果关系估计方法
六、 敏感性分析
七、 结果可视化
八、 GWAS数据介绍
九、 文章导读
十、 R语言示例代码（以本地TSMR为例）
十一、 MR的类别与应用
楔子：今天干嘛要学孟德尔随机化？
所谓学以致用学以致用，能用上的玩意才有学的意义。这个方法名字倒是相当有趣，但我们总得知道大费周章地学了干嘛？
话不多说，先画三个大饼：孟德尔随机化简直好极了！
（一）超级热门的科研选题
（这是在PubMed以“孟德尔随机化”为主题词语，搜索到的发文数量随年份的变化）
从PubMed上的检索结果可以很明显地看出来，近几年，孟德尔随机化的发文数量呈现井喷式增长，成为了医学界的“新晋顶流”；而孟德尔随机化与其他研究结合的思路更是高分文章频出；2024年，孟德尔随机化的发文数量达到了近7000篇。可见，在孟德尔随机化依然是一支充满无限可能的“潜力股”，等着大家去不断发掘~换句话说，如果你现在就开始学习，赶上这个浪潮正合时宜！
（二）小白也可以快速上手的科研思路
做孟德尔随机化的出发点，简单的来说，就是探究不同事物之间的因果关系。它所依托的严密逻辑和假设让研究结果既可靠又成熟，通过最基础的假设，却能够衍生出一系列丰富的研究思路。对于初学者而言，它不仅是一个快速入门的绝佳工具，也是锻炼和完善科研思维的大好机会。
同时，孟德尔随机化对代码的简单易于操作，要求不算苛刻，相比较于传统的队列研究耗时久远，一个突发奇想的idea + 一段R语言代码，就可以快速地开始一次小型MR研究，听起来就很让人心动~
（三）本科生友好的科研选择
不同于那些需要在实验室苦心培养细胞，或者长年累月的广泛调查的实验，我们可以从许多免费公开的网站获取到多样且海量的数据。试想一下，当你的同学还在为收不到问卷数据头秃时，你已在国际数据库挖到宝藏，实现了科研的“弯道超车”，这难道不香吗？不用做实验、易于上手操作……上述的种种优势都使孟德尔随机化成为了刚刚接触科研、动手能力仍需提高的本科生们入门的不二之选，可以说，孟德尔随机化是一个十分“本科生友好”的科研方法。
一、孟德尔随机化的概念
在介绍MR是什么之前，我们需要从最基础的部分开始，稍微看看一些别的概念~
1.1 病因
简单来说，病因是指能使人群发病概率增加的因素。当这些因素不存在时，人群疾病发生频率会下降。比如，吸烟是肺癌的病因，因为吸烟人群的肺癌发病率远高于不吸烟人群。由此可见病因和疾病密切相关，研究病因对研究疾病来说具有很大的现实意义。
1.2 病因研究方法的分类与传统流行病学方法的弊端
1.2.1传统流行病学方法（简单了解一下就可以啦）
（图源：预防医学-第七版）
a. RCTs系统综述：多个RCT研究进行综合分析的方法
b. RCT（Randomized Controlled Trial）：随机对照试验，医学研究中的“金标准”
c. 队列研究：跟踪一组人群随时间的变化来研究特定因素与结果之间的关系，前瞻性地观察暴露与结局的关系。（前瞻性：在研究开始时确定研究对象，并在随后的时间中跟踪（或称前瞻）这些对象）
d. 病例对照研究：回顾性地比较某种疾病的病例与未患该疾病的受试者，确定可能导致该疾病的暴露因素
e. 生态学研究：用医疗保健生态学模型描述社区人群健康状况、就医行为与医疗卫生资源的关联
从上到下，这些方法的观察性依次增强，实验性依次减弱，证据等级也逐渐减弱。
但是在研究中，大家发现这些方法都有难以避免的弊端：
实验性强的研究虽然置信度高，但是要求高，控制严，难度大，还可能存在伦理问题
观察性强的研究存在反向因果（搞不清究竟是a导致了b还是b导致了a）、混杂因素（一些我们不想研究的，却能够影响结果的因素）等等问题，置信度低
那么为了找到一种可以有效规避这些问题的研究方法，大家找到了我们今天的主角~
1.2.2 孟德尔随机化（简称MR）（证据等级位于RCT和队列研究之间）
1.3.1 孟德尔随机化的定义
为什么MR这个方法的名称跟那位在修道院种豌豆的伟大生物学家有关？
大家在高中生物课都学过孟德尔定律，也就是亲代的等位基因会在受精时随机分配给子代。这是一个带大家回忆一下孟德尔定律的图：
而这个过程特别像随机对照试验中的随机分组：这张图展示了随机对照实验（RCT）和孟德尔随机化（MR）的相似之处。
我们可以假设，个体从父母随机继承的基因变异是固定的，不受后天混杂因素影响，并且这些基因变异在人群中随机分布，有这个基因和没有这个基因的人在其他与这个基因无关的方面都差不多，就相当于一个天然的控制变量，也形成了一个模拟对照试验。孟德尔随机化建立在这个假设上，依靠模拟随机对照试验减少混杂因素影响，所以孟德尔定律正是孟德尔随机化的理论依据。
孟德尔随机化的基本思想是：利用与暴露因素强相关的遗传变异作为工具变量，从而推断暴露因素和结局变量之间的因果关系
相信大家看得有点懵，所以要解释一下这句话：
举个例子，如果某个特定的遗传变异（姑且称之为a）与高胆固醇水平显著相关，那么携带变异a的人可能天生就有较高的胆固醇水平，也就是说，这个遗传变异a和高胆固醇的性状是强相关的，可以被用作研究高胆固醇（作为暴露因素）和其他表型之间的因果关系的工具变量。
接下来的研究中，如果我们发现携带遗传变异a的个体在结局变量（比如二型糖尿病）上表现出显著的差异（例如，携带a变异的人明显有更高的概率患有二型糖尿病），那么，这可能表明高胆固醇与二型糖尿病之间存在因果关系。
这里，我们特别介绍一下一个经常会用到的两个名词——SNP和GWAS：
SNP（单核苷酸多态性）是指基因组中，出现在单个核苷酸位置上的变异。这种变异可以是单个碱基的转换、颠换，或由碱基的插入和缺失引起，比如某个碱基对从A-T突变成了C-G，造成了基因突变。（看看下图，大概就是这样）
GWAS是一种数据类型。研究员通过对大量人的整个基因组进行扫描，检测他们基因组中数百万个遗传变异位点，分析这些位点与特定性状或疾病之间的关联，得出的结果就是GWAS数据。GWAS的作用是发现遗传变异（尤其是SNP）与特定疾病或者性状的联系，关于GWAS数据，我们接下来会有一栏专门进行介绍哒。
SNP是基因组中最常见的遗传变异形式之一，也是GWAS中最常用的研究对象。
1.3.2 孟德尔随机化的优势
1. 减少混杂因素的干扰：由于遗传变异在受精时随机分配，特别类似于随机对照试验的随机分组，也就可以有效减少混杂因素的影响。
2. 提供因果推断的有力证据：与传统流行病学方法相比，mr更能提供因果关系的证据。显然，个体的遗传变异是天生的，在发生的顺序上早于疾病，因此一般来说是遗传变异导致了疾病，而不是疾病导致了遗传变异。这样就可以排除由于逆向因果问题所带来的偏倚。
1.3.3 孟德尔随机化的应用
它在现实应用中主要用评估某些因素（比如吸烟喝酒，或者药物靶点）对健康结果的潜在因果影响。当暴露因素对人体有害时，在现实世界中进行随机对照试验不符合医学伦理，此时应用孟德尔随机化来进行研究，相较于实验性研究更加可行的同时，也比观察性研究置信度更高哦。
二、孟德尔随机化的原理
2.1 核心假设
先来看一些最基本的理论知识，这将是建设对于MR理解的第一步。对于MR，存在三个核心假设，这些假设是因果推断有效性的关键，也是MR存在的基础。
关联性假设（Relevance）
工具变量（遗传变异）必须与暴露因素强相关。例如，如果研究BMI与糖尿病的因果关系，选择与BMI强相关的遗传变异（SNP）作为工具变量。如果要量化“强”这个程度，就要进行一些统计学指标的计算，这点我们接下来再讲。
独立性假设（Independence）
工具变量必须独立于混杂因素。工具变量无法直接通过影响混杂因素进而影响结局，从而保证了暴露必然是结局的因素之一。
排他性假设（Exclusion Restriction）
工具变量只能通过暴露因素影响结局，而不能通过其他途径影响结局。如果我们要研究吸烟与肺癌的关系，工具变量不能直接影响肺癌的发生，也不能通过非暴露因素间接影响结局（ps：这点与独立性假设类似），只能通过影响吸烟行为来间接影响肺癌。
下图更直观地体现出这三个假设（IV为工具变量instrumental variable的缩写）：
a.IV与暴露（exposure）之间存在强关联；
b.IV与任何混杂因素（confounder）之间没有关联；
c.IV对结果（outcome）的影响只能通过暴露（exposure）实现，不能通过其他途径直接影响结果（outcome）
这里出现了一个新的名词叫confounder——混杂因素。混杂因素更常被描述为存在一个开放的“后门路径”，就像一只老鼠偷偷钻进了我们因果关系的机房，悄咪咪咬断了一根电线……我们却没有意识到它的存在，这就会导致意料之外的错误。它同时影响暴露和结局，但却不在暴露与结局因果路径上，所以当混杂因素存在但是却被忽略的时候，就会导致结果出现偏差。
如果任何一个核心假设被违背，因果效应的估计可能会出现偏倚。因此，在进行MR分析时，需要仔细验证这些假设是否成立。
基于这三个核心假设，在选择MR工具变量的时候会有一系列筛选标准，相关内容请看下文工具变量的选择与验证部分~
2.2 进行MR的前提条件
1. 选择同质人群。同质人群是指由具有相似特征或共同背景的个体组成的群体，比如研究用到的两个数据都来自欧洲人，那这两个数据对应的人群互为同质人群。选择同质人群进行研究可以有效避免人群分层的现象（详见下文核心假设的违背一栏），避免不必要的偏倚。
2. 纳入统计的人群要符合孟德尔群体遗传的Hardy-Weinberg定律（这个定律高中学过，就是在接近理想条件下，基因频率和基因型频率在每一代间保持恒定）
3. 纳入研究的基因型要在人群中有一定的变异率，不能特别特别罕见（也就是MAF（最小等位基因频率）至少大于1%），样本量不能太小。不过一般来说，低质量的基因型在GWAS数据中已经被筛选掉了，不用怎么担心~
4. 明确基因功能、基因——基因和基因——环境交互作用，这样才能尽量排除和控制基因多效性和连锁不平衡等对效应估计的影响（后头咱再来解释什么是基因多效性和连锁不平衡吧）
三、核心假设违背的影响与检测方法
这里的原理简单了解就行，重点关注怎么利用这些方法，解决什么问题
3.1 水平多效性（Horizontal Pleiotropy）
3.1.1. 定义：水平多效性是指工具变量在通过暴露因素（exposure）影响结局（outcome）的同时，也会通过其他途径影响结局，也就是，作为工具变量的SNPs背着暴露因素，偷偷也跟结局变量扯上关系了。
需要提到的是，多效性分为两种，水平多效性和垂直多效性。但只有水平多效性是需要避免的。垂直多效性指的是一个基因变异如SNP通过影响一个可遗传的表型，进而间接影响另一个相关表型的现象。垂直多效性是mr的本质，在这里不需要考虑。
在这篇文章中，我们将水平多效性统一简称为多效性。
3.1.2 影响：违反排他性假设
3.1.3 检测方法：
a. MR-Egger回归：用于检测和校正水平多效性
原理：假设基因变异（工具变量）对暴露和结局的效应分别为bx和by，若不存在多效性，在线性回归的结果图中，所有SNP的by与bx的关系线应该是通过原点；若存在多效性，by与bx的关系会出现截距。MR-Egger通过估计截距来判断是否存在多效性及方向和大小。简单来说，可以通过看生成的图像的截距来判断多效性。
b. MR-PRESSO：检测和校正异常值，这个非常非常重要。
具体包括这三个部分（稍微了解就可以啦）
全局测试：用于检测所有工具变量（SNPs）中是否存在水平多效性。若全局测试的P值小于0.05，表明存在水平多效性。
异常值测试：用于识别和剔除那些可能由于水平多效性而影响MR分析结果的异常SNP。会逐一检查每个SNP，若某个SNP的P值小于0.05，则认为该SNP是异常值，可能需要从分析中剔除。
扭曲测试：用于评估在剔除异常值前后结果是否存在差异。若扭曲测试的P值小于0.05，表明剔除异常值后，因果效应的估计发生了显著变化，提示结果可能受到了异常值的显著影响。
3.2连锁不平衡（LD）
3.2.1 定义：连锁不平衡是指遗传位点之间存在高度相关性。某些遗传变异可能位于基因组的同一区域内，相互关联，倾向于共同遗传，也就是它们老是黏在一起，不怎么进行基因重组，导致不符合孟德尔定律。
3.2.2 影响：违背了排他性假设。
3.2.3 解决方法：
a. LD clumping：进行“修剪”，根据算法在两个高度相关的snp中保留一个更有代表性的，删掉另一个，确保工具变量之间没有高度相关性。
通过设置r²阈值（一般为0.001）以及窗口（在GWAS数据中一般为10000kb）筛选独立的工具变量。
那么r²和窗口是什么呢？
r²的取值介于0～1之间，r² = 1表示两个SNP是完全的连锁不平衡关系，意味着这两个SNP在遗传过程中总是一起遗传，不会因为基因重组而分开；r² = 0则表示两个SNP是完全连锁平衡的关系，即它们是完全随机分配的，彼此之间没有连锁关系。后者则是我们做孟德尔随机化时所需要的，一般来讲0.001是一个常规阈值。
窗口是指在基因组上用于检测连锁不平衡的区域范围。窗口的大小通常以碱基对（kb）为单位，在窗口内相关的SNP都会被“修剪“掉。如果设置clump_kb = 10000，意味着以每个snp为中心，对其上下游分别10000kb的区域进行修剪操作。
注意啦，任意满足r²或者kb其中之一修剪条件，就会被修剪。
3.3人群分层
3.3.1 定义：人群分层是指研究人群由于存在不同的亚群结构，而导致遗传和环境因素在这些亚群中的分布不均衡，进而可能影响研究结果的现象。
举个例子，西欧人和东亚人都是人，虽然因为没有生殖隔离所以不能算两个物种，但由于饮食习惯，气候原因等因素的共同影响，这两个人类亚群的遗传变异可能存在很大区别，在同一个研究因果的MR中，可能会导致不同的结果。
3.3.2 影响：
a. 产生虚假关联：可能使遗传标记SNP与疾病之间出现假阳性关联，也就是我们以为结果是阳性，实则本应是阴性的。
b. 降低统计效能：如果把人群不做区分，一起来看，这样会增加数据变异性，使检测真实遗传效应的难度加大，需要更大样本量才能达到相同的统计效力，导致研究效率降低。
3.3.3 检测方法：
一般来说，我们直接研究时会看GWAS数据源是全欧洲人群还是亚洲人群还是混杂的，这个时候尽量选取单一人群进行研究，就在一定程度上可以降低人群分层的影响。（当然也有其他方法，大家感兴趣的话可以自行了解）
3.4 collider（了解即可）
3.4.1 定义：Collider 是一个变量，它受到两个其他变量的共同影响。当如暴露因素和结局能够独立影响第三个变量Collider时，调整这个第三个变量可Collider能会引入偏倚。
举个例子，假设你正在研究吸烟（暴露）和肺癌（结局）之间的关系。同时，还有一个变量，寿命。按照常识我们知道，吸烟的人寿命可能更短，患肺癌的人寿命也可能更短。这时，寿命就是一个 Collider，因为它同时受到吸烟和肺癌的影响。
如果我们在研究中只关注长寿人群（即可以理解为调整了寿命这个变量），就会出现偏倚：在长寿人群中，吸烟者可能更少，因为吸烟者往往很难进入长寿人群。同时，患肺癌的人也可能更少，因为肺癌患者寿命也较短。
因此，如果你只研究长寿人群，可能会发现吸烟和肺癌之间的关系减弱了，这就是 Collider 偏倚。
3.4.2 影响：违背了排他性假设。
3.4.3 检测方法：
a）敏感性分析：通过控制可能的collider，观察暴露与结局之间的变化进行调整
b）因果图分析：此方法可直接观测出是否有可能的collider，方便采取策略
仔细想想会发现collider和confounder是有区别的。collider是受到暴露和结局共同影响的因素，而confounder则是同时影响暴露和结局的因素。所以对待它们的处理方式也就是不一样的啦。
3.5 弱工具变量
3.5.1 定义：弱工具变量是指与暴露因素的关联程度较弱的工具变量。一般而言，如果工具变量与暴露因素之间的相关系数较小，例如F统计量小于10，或者第一阶段回归的R²值较低（了解即可），通常认为该工具变量可能是弱工具变量。
需要注意的是，p值小并不等价于F值大，也就是说即便一个SNP对应的p值可能很小，它也可能是弱工具变量，仍然需要用F统计量检验进行判断。
3.5.2 影响：违背了相关性假设。
3.5.3 检测方法
F统计量检验（一般来说都用这个方法）：F统计量是一种用于比较不同组之间方差差异的统计量，在基于线性回归的工具变量分析中，F统计量用于衡量工具变量与暴露变量的关联强度。和高中我们学的R²的概念有点类似，F统计量越大，关联越强，越符合条件
计算公式是：F=（R²/k）/【（1-R²）/（n-k-1）】
（不过大多统计软件会直接计算出F的值，所以其实不用背的……）
计算工具变量与暴露变量回归模型中的F统计量，若F统计量小于10，通常提示可能存在弱工具变量问题。
四、工具变量的选择与验证
现在我们已经基本了解了MR的原理，不过在正式进行MR分析之前，还有一些准备工作需要完成……
为了提高孟德尔随机化结果的稳健性并减少效应估计的偏差，在每一组MR开始前，都需要进行一系列对工具变量的筛选，让这些工具变量的各项数据都符合三大假设的要求，才能保证因果推断的有效性和准确性。
4.1 筛选标准：基于三大假设
a. 强相关性：工具变量必须与暴露因素强相关。这一步的实现通常包含两次筛选。
第一是针对该SNP对暴露的P值进行限制。普遍而言，选取p＜5e-08的工具变量是更加稳健的（也就是更保守），但这可能会导致能用的工具变量很少，故有时候也以5e-06、1e-05等作为阈值。
第二，部分工具变量可能具有较弱的对表型变异的解释效果，也就是弱工具变量问题。此时我们需要对F统计量进行阈值设定，一般选取F＞10的工具变量以去除弱工具变量（参考上文弱工具变量一栏相关内容）。
b. 独立性：工具变量必须独立于混杂因素。混杂因素当然是相对的，对于每个研究都不同。在研究开始之前，我们必须先了解对于我们当前的研究而言，什么是所谓的“混杂因素”。
通常而言，我们可以通过临床经验、流病调查、文献综述等方式，得知除了当前敲定的暴露以外，还有什么表型会导致结局的发生，这就是我们寻找的“混杂因素”。其次，通过NCBI等网站（或通过接入相关网站的网页端口），我们能够得到与当前留下的各个工具变量相关的多种表型，若某个工具变量涉及已确定的混杂因素，则需将该工具变量去除。
c. 排他性：工具变量只能通过暴露因素影响结局。这通常也被描述为“没有”水平多效性，表明我们希望工具变量和结局“没”强相关性。为实现这一目标，我们应当对工具变量对结局的值进行过滤，通常设置为p＞5e-05。
4.2 其他筛选：
Steiger Filtering：除了需要满足三大假设，我们还应当注意，由工具变量呈现的因果关系也可能出现因果倒置的问题，也就是我们不知道事实上究竟是“暴露导致了结局”还是“结局导致了暴露”。Steiger Filtering 的核心思想是通过比较SNP对暴露和结果的解释变异（r²）来判断因果方向。一个有效的工具变量应该对暴露的解释变异大于对结果的解释变异。如果某个SNP对结果的解释变异大于对暴露的解释变异，则该SNP不适合用作工具变量。为此，我们还需要在Harmonization后，孟德尔随机化之前，添加Steiger Filtering，确保不会出现反向因果。
五、因果关系估计方法
5.1 为了方便大家阅读数据和结果，这里有一些重要的参数的解释
a. beta（β）：表示某个变异（例如SNP）与表型之间的效应大小。β值反映了每个等位基因的效应。在数学上，它通常是线性回归中的斜率slope（具体是什么跟回归模型的种类有关）。
b. se：表示β值估计的精确度，标准误差越小，表示β值估计的精确度越高。标准误差是β值估计的不确定性度量。
c. pval：表示β值是否显著不同于零的概率。
d. OR：表示为比值比，即表达暴露对结局的影响。如果OR大于1，说明暴露增加了某个结局发生的概率；如果小于1，则说明暴露降低了结局发生的概率。（β值和OR值可以很简单地互相换算，β=ln(OR)）
e. samplesize：也就是N，意思是样本量，在Steiger Filtering时用得到，以及在共定位时也需要使用得到
i. case：病例数
ii. control：阴性数（健康人的数量）
这么多的概念……这么多的别名……没关系的没关系的，看多了就熟练了（
5.2 方法学汇总
还坚持看到这里的朋友们，恭喜，这才是MR的本体部分！
MR的本质是计算暴露对结局的影响，也就是OR值。
MR的方法随着研究的深入有着日新月异的变化，也不断地在开发出新的算法，解决了先前存在的问题和缺陷。但请注意，每一种方法都有优缺性，他们或许具有很强的可靠性但过于严格，有可能具有简便性但结果却不稳健。因此，在真正实践MR时，我们要根据研究目的、研究流程、研究具体情况进行分析，选择合适恰当的方法。
（虽然90%的情况下我们只会看IVW的结果）
a. Wald Ratio：最基础的默认方法，适用于单个工具变量。
b. IVW（Inverse variance weighted）：最重要的方法，检验效能高，但是受多效性影响比较严重
c. MR-RAPS：适用于工具变量较弱且存在多效性的分析。这个方法允许纳入些弱工具变量进行稳健统计估计。该方法可以作为增加工作量的补充，结果与IVW结果相近。
d. MR-PRESSO：（参见上文水平多效性一栏对于检测和校正异常值的介绍）本质而言，MR PRESSO不算一个因果效应的估计方法，而应该作为一种检验水平多样性的方法。
e. Weighted Median：常用于存在多效性或离群值的场景，要求至少有一半是工具变量。优点是对于离群值很稳健。缺点是对于遗传变异的增加或者减少很敏感。
f. MR-Egger：适用于存在多效性或工具变量可能无效的情况。优点是可以在较弱的假设下（InSIDE假设）提供一致的因果效应估计，缺点是统计效能较低。
当然还有很多其他的MR方法，比如：
weighted mode
MR-robust
cML-MA
debiased IVW
对同一个数据应用不同的MR方法结果可能各不相同，它们各自都有对应的优缺点，因此还是要在研究中因地制宜地进行选择。
六、 敏感性分析
现在我们已经可以初步通过mr的结果判断因果关系了，但是结果的出现并不意味着研究结束，我们还需要“旁侧敲击”地确认一下结果是不是可靠的。
6.1 目的
在MR分析之后还需要判断结果是否可靠，并探究潜在的偏倚，比如是否存在基因的多效性或数据的异质性，同时评估某个工具变量有没有显著影响结局变量从而导致偏倚。
6.2 处理方法
a. 基因多效性检验：水平多效性是指当暴露因素为0时，工具变量对结局仍有效应，也就是就算不存在暴露因素了，工具变量也会通过混杂因素对结局起作用。直观来看，就是看散点图中的y轴截距为0（后面会有对散点图的介绍）。如果水平多效性检验显著，说明工具变量不止与暴露相关，还与其他混杂因素相关，因此mr结果很大程度不可靠。
我们可以通过MR-PRESSO判断是否存在水平多效性。与进行因果推断不同，我们所希望presso的结果是P>0.05，只有这样才能说明没有多效性
b. 异质性检验：由于两样本随机化的暴露和结局分别来自于不同的样本，可能存在人群不同等异质性，因此需要通过异质性检验说明。与多效性检验一样，我们所希望的是P>0.05，若P值＜0.05说明存在异质性。
c. 留一法（Leave-One-Out Analysis）：逐一排除工具变量再进行计算，观察结果是否稳定。理想的情况是每一次剔除SNP后结果变化都不大（在留一法森林图中，看到单独删除各个SNP后以及原来结果的beta值以及置信区间都在0的同一侧）
如果有一次剔除某个SNP（称为a）后效应值变化很大，说明这个SNP（a）可能通过影响别的SNP（称为b）对整个因果关系进行了影响，说明SNP（a）可能其实违反了三大假设，需要剔除SNP（a）。
需要注意的是，在进行敏感性分析时，由于我们的零假设是不存在水平多效性以及异质性，所以我们希望p>0.05（所谓p值指的是零假设其实是对的，但我们做出的结论是零假设不成立这种错误发生的概率，常以0.05为一个关键值辅助我们进行判断），越大越好从而证明我们的假设成立。
这点与进行因果推断时不同。因为在进行因果推断时零假设是暴露与结局之间是独立的，但我们希望证明因果关系是存在的，p<0.05说明我们因认为暴露与结局之间有因果关系而犯错的概率小于0.05（这已经非常小了），可以认为有因果关系。
七、结果可视化
得到了结果之后，为了更直观地展现给别人看，我们需要把结果画成图，也就是结果的可视化。
MR分析可视化一般有四种图像，分别为散点图、森林图、漏斗图、留一法图
7.1 散点图
散点图向我们展示了作为工具变量的SNP（也就是IV） 与暴露 (Exposure) 、结局(Outcome)之间的关系。
其中每一个点都表示一个SNP，而x轴表示与暴露相关的遗传变异的效应大小(Effect Size)，y轴表示与结局相关的遗传变异的效应大小。每个点所发出的水平线与垂直线代表了其在x轴与y轴上的置信区间。
图上还有根据不同因果效应估计方法绘制的拟合回归直线，其斜率代表了暴露对结局的因果效应的估计，对应效应Beta值（beta>0代表暴露是结局的风险因素，beta<0代表暴露是结局的保护因素）
7.2 森林图
森林图的x轴是暴露对结局的影响（beta），y轴则为每一个rsid对应的SNP，每条水平实线表示用我们选择的MR方法（比如Wald ratio）基于单个SNP估算的结果: 完全在0左边的表示该SNP推算出暴露增加会降低结局风险；完全在0右边的表示暴露增加会提高结局风险；跨过0的表示结果不显著。
单个SNP的结果可能与总体结果相悖，但别紧张，对我们而言综合所有SNP得到的最底下的红线，才是我们希望得到的最终结果，也是最为关键的。
7.3 漏斗图
漏斗图可以帮助我们来评估MR研究结果是否存在偏倚, 我们重点需要观察漏斗图整体上是否是对称的。如果存在严重离群的SNP（高度不对称）则说明可能存在离群值，可以去除后重新进行MR分析。
每个点是一个SNP, 其中x轴表示该遗传变异对结局变量的效应大小, y轴表示相应的标准误差或权重(1/标准误)。
图中的实线为为使用不同方法方法获得的Beta效应值结果。
7.4 留一法图
留一法图形式与森林图相似，但不同之处在于留一法图用于展示每个独立遗传变异对整体MR分析结果的影响。顾名思义，这个方法是逐个剔除n个SNP中的单个SNP，然后用剩下的（n-1）个SNP重新做MR分析。
Y轴对应每个剔除的rsid号，和全部变量都不剔除的ivw 方法，X轴对应 ivw结果的具体数值，黑点和红点是beta效应值，线段代表beta值的置信区间
一般而言最高标准是黑点红点和所有的线段都在0的同侧，最低标准是黑点和红点在 0这个线的一侧，否则说明某个SNP对MR结果存在严重的偏倚，影响了结果的稳健性，可以去除该SNP重新分析。
下图简单展示了一般森林图和留一法森林图的区别，大家可以按照我们上面的文案找找不同喔~（图源CSDN，侵删）
那么至此为止MR的基本流程我们就介绍完了，但是显然还有许多令人疑惑的地方……比如说，我们的数据从哪儿来的？
八、GWAS数据
8.1什么是GWAS数据
8.1.1 GWAS数据是MR分析中最常用的工具变量来源。
众所周知，不同人的基因并不完全相同，相应的，每个人在疾病或者某些性状上的发生风险也有不同，由此，GWAS数据就是分析这之间关联的数据。
简单来说，GWAS数据就是把很多人的所有基因进行扫描，分析其中SNP位点的不同，通过一些算法，凭借海量样本来计算总结一些基因与疾病之间的关系。
GWAS是生信研究中最为关键的数据种类之一，同时也是我们做MR分析中最常用的工具变量来源。所以，了解一下GWAS数据的组成和结构是非常必要的！
8.1.2 数据类型
a.基因型数据：是GWAS数据的核心部分，记录了每个样本在众多SNP位点上的基因型信息，比如有些人在某个位点是A，有些人是T等等
b.表型数据：与个体的外在表现或疾病状态相关，包括身高、体重等生理性状，以及是否患某种疾病、疾病的严重程度等疾病相关表型。
c.样本信息：包含样本的来源、性别、年龄等基本信息。
例：BMI数据（基因数据）
（ieu open GWAS 网站：IEU OpenGWAS project）
这里我们可以看到这个数据库有很多信息（不是每个数据都这么详尽的，因此需要我们来进行筛选数据来源），简单介绍几个的比较重要的：
trait，也就是性状，表型，这里就是BMI，
Dataset：ukb-b-19953，这里是这个数据库的编号，相当于数据库的“身份证”，
Year是日期的意思，我们尽量选取日期较新的（在保证数据符合要求的前提下），
Population是人群的意思，这里就是欧洲人群，
sample size样本数量，
Number of SNPs，所包含的SNP的数量，
build，这里代表的是一种数据格式，我们一般选取后面标注GRCh37中含“37”的（还有一种是GRCh38，这种格式处理起来会比较困难一点，因此尽量选择37的）
8.2 GWAS数据那些乱七八糟的列名到底都是什么意思
在数据处理完导入数据后，打开（使用View函数）大概率可以看到这样的界面
（需要注意的是，每个下载下来的数据源的命名方式可能不同，甚至可能缺失关键的列，所以会有清洗数据这一步，也就是自己选择出需要的数据列，并改成我们想要的列名）
在进行数据清洗之前，要记得“colnames（）”或者“head（）”一下，先了解一下数据哦。
我们可以看到有很多列名，其中只有部分是我们需要的。在本地MR中，我们一般需要7个数据（理论上最少4个：SNPs，effect-allele，se，beta）:
SNPs：图中的variant id，可以理解成SNP的编号，有可能叫SNPID等等各种名字，这时候就需要具体看看数据内容，一般都是“rs＋数字串”的形式
effect allele：效应等位基因简称ea，效应等位基因是一种遗传变异，它会增加某些情况（癌症、一般疾病）的可能性，或者更有可能产生特定特征。
other allele：其他等位基因简称oa，其他位点上与之对应的另外一个等位基因，
effect allele frequency：效应等位基因频率，简称eaf，一种衡量遗传物质变异指数的指标，效应等位基因在群体中特定基因的所有等位基因中的比例
standard error：标准误简称se，可以理解成样本中不同组平均值的标准差，
beta: β，等位基因1的效应大小估计，也称为等位基因的效应值。
P-value: 有时候简称P，是目标SNP与表型是否显著关联。小的P-value被视为零假设（目标SNP与表型无关系）可能不成立，即备择假设成立，也就是目标SNP与表型有关系的证据。（零假设与备择假设部分参考医学统计学）
其他非必须的列名简单介绍：
phenotype：SNP在所影响的表型
chr和pos：Chromosome (Chr) and position (Pos) 染色体和编码位置，利用染色体编号和位置编号确定SNP；但参考基因版本不同，同一位点chr:pos会不同。如果数据的SNPs缺失，可以通过代码由chr和pos找到我们需要的SNPs噢，不用担心~
sample size: 样本量大小
gene:基因或其他SNP的注释
有些数据之间是可以相互计算得到的（比如beta, se, p-value），网上有很多这方面相关的帖子，大家可以自行搜索。
到这里，我们对于孟德尔随机化的基本知识讲解就结束了。现在让我们一起看看一篇早期的TSMR（双变量孟德尔随机化）高分文章，了解一下如何用孟德尔随机化产出paper吧！
九、文章导读
（原文链接在此：https://pubmed.ncbi.nlm.nih.gov/30673084/
）
这篇文章研究的是甲状腺功能遗传决定因素与房颤关系。
9.1 一些背景知识：
甲状腺是什么咱们就不介绍了。文章中提到并且用上了多个甲状腺功能的遗传决定因素作为研究对象。之所以选用这些遗传决定因素进行研究，是因为过往研究可能表明这些因素在甲状腺功能评价中比较重要。
“甲状腺功能的调节复杂，涉及甲状腺、垂体、下丘脑及反馈机制。FT₄水平通过反馈调节垂体抑制促甲状腺激素的产生。在甲状腺和外周组织中，FT₄可转化为活性游离三碘甲状腺原氨酸（FT₃），可通过循环中 FT₃:FT₄比值评估这一过程。因此，FT₄水平只是甲状腺功能的一个衡量指标，其与房颤的关联即使存在因果关系，也可能并非直接联系。”
↑↑不是我说的，是作者在文章中自己写的，所以一个遗传决定因素不够研究，两者之间的因果关系需要通过多个遗传决定因素进行研究。
作者选定的相应的遗传决定因素和来源我们在下面列出来：
FT₄：来自 Teumer 等（2018）的研究。
TSH：来自 Teumer 等（2018）的研究。
FT₃:FT₄比值：来自 Panicker 等（2008）的研究。
甲状腺功能减退：来自 Eriksson 等（2012）的研究。
甲状腺过氧化物酶抗体（TPOAb）：来自 Medici 等（2014）的研究。
甲状腺功能亢进：来自 Gudmundsson 等（2012）的研究。
房颤，全称心房颤动，是一种常见的快速心律失常，随着年龄增长发生率成倍增加，其中无器质性心脏病患者占3%～11%。房颤的发生与年龄和基础疾病类型有关。不过这不是咱们关注的重点，就当补充一点医学知识好了
那问题来了，为什么作者会“突发奇想”研究一下甲状腺功能遗传因素和房颤的关系呢？他在文章中给出的解释是这样的：
简单来说，就是他们在以前的研究者做的观察性研究中发现，参考范围内游离甲状腺素（FT₄）增加与房颤风险增加有关，而明显和亚临床甲状腺功能减退与房颤风险降低有关……吧啦吧啦，总之就是二者有关。
但游离甲状腺素FT₄究竟是甲亢 - 房颤关联的因果因素还是生物标志物尚不明确，相关的各种决定因素和房颤的潜在关系都还没有人搞清楚过，也就是我们之前提到的，观察性研究有很明显的弊端，研究不够透彻，所以作者们在这里挖掘到了潜在的研究点。
9.2 首先要关注的是研究数据来源：
这图截自paper开头部分的abstract，如果需要更详细的数据信息可以看正文Participants和Summary-Level Analysis开头部分的介绍，因为太长所以这里就不截取了。当然，如果想要找到数据最原始的来源，就应该去文章末尾的引用部分找到作者是从哪里找到的数据，然后再进行溯源。
整理成中文大概内容如下：
研究水平（Study-level）数据来源：11 项队列研究，均为欧洲血统人群。共纳入 56,912 名参与者，其中房颤病例 7,679 例（2,093 例现患，5,586 例新发），对照 49,233 例。
汇总水平（Summary-level）数据来源：AFGen 联盟（房颤遗传学联盟）的 GWAS（全基因组关联研究）汇总数据。涉及 66 项研究，均为欧洲血统人群，共 537,409 人，其中房颤病例 55,114 例（7,672 例新发，47,442 例现患），对照 482,295 例。
这两个数据，分别对应着文章中的两个MR研究方法——研究水平的分析和汇总水平的孟德尔随机化分析。
9.3 我们把重点放在方法学部分
这个研究分为两个部分：研究水平分析和汇总水平的MR
9.3.1 观察性研究：对纳入的 11 项队列研究进行了分析，这些队列研究属于观察性研究范畴。通过对这些队列研究数据的整理和分析，作者得以探究甲状腺功能相关指标（如游离甲状腺素FT4、促甲状腺激素等）与心房颤动（AF）之间的关系。
结果呢？在多变量调整的观察性研究水平数据的荟萃分析中，结果表明，FT4与房颤风险存在一定关联，而促甲状腺激素与房颤风险的关联不明显。
因为荟萃分析并不是今天的重点，所以我们先不讲啦，大家有兴趣自己查查看。
这一段放在这里其实是为了说明方法并不是孤立使用的，在真实的研究中我们往往要结合多种方法，除了孟德尔随机化分析外，还可以结合观察性研究数据、实验研究结果或其他统计方法来验证因果关系。研究先通过荟萃分析发现甲状腺功能指标与房颤的关联，进而引发对因果关系的思考，促使进行 MR 分析。同时，MR 的结果可对荟萃分析结果进行验证和补充。如果这些不同来源的证据能够相互印证，形成一个“证据三角形”，那么研究结论的可靠性和说服力就会显著增强。
9.3.2汇总水平的孟德尔随机化分析（Summary-level MR）
（1）数据来源：GWAS汇总数据，来源就不再讲一次了。
（2）工具变量：选择与甲状腺功能指标（FT4、FT3:FT4比值、促甲状腺激素、甲状腺功能减退症、甲状腺过氧化物酶抗体和甲状腺功能亢进症）相关的GWAS显著SNPs。
这些甲状腺功能指标是怎么选出来的，我们上文已经提到过了，是作者翻过往研究了解相关的甲状腺功能通路找到的相对重要的指标。
那跟这些指标显著相关的SNP又是怎么选出来的呢？这就要回顾一下工具变量选择的知识了：我们通过在与以上功能指标相关的GWAS数据中，筛选P值小于某个阈值（通常是5e-8，正文中没有提到）的SNP，同时文中还特意提到了使用了选择F＞10的SNP作为工具变量，从而避免了弱工具变量造成的偏倚。
总之，通过这样的筛选方式，作者找到了可以用于当mr工具变量的SNP们，并正式开始研究
9.3.3 分析方法
采用TSMR
使用逆方差加权固定效应模型（IVW-FE）和随机效应模型（IVW-RE）计算汇总因果效应。
那么，怎么样的结果才能算是阳性，也就是说明了存在显著关联呢？一般我们认为结果P小于0.05就算是阳性，但是文章中却选取了P＜0.008（也就是0.05/6的结果）作为阈值：
这样的做法被称为Bonferroni校正。
要介绍这个我们需要提到一些统计学知识：
在统计分析中，当进行多次假设检验时，会增加犯第一类错误（假阳性）的概率。为了控制这种错误率，需要进行多重比较校正。常见的校正方法包括：
Bonferroni校正：将显著性水平（α）除以检验的次数，在这个研究中我们检验了六次（看上面那张图作者提到了FT4等等六个比较内容），所以最终的显著性阈值就变成了0.05/6=0.008。不过Bonferroni比较严格，在比较次数多的时候又可能会增加假阴性概率，这时候就需要用别的方法。
Benjamini-Hochberg校正：又称FDR。比Bonferroni校正更宽松，适用于比较次数较多的情况。
Holm校正：一种逐步校正方法，结合了Bonferroni校正的优点。
跳过这个对统计学校正的题外话，回到我们的研究方法中来，然后——
使用MR-Egger回归，用于检测和校正水平多效性。
使用加权中位数分析（Weighted Median Analysis）进行对多效性稳健的估计。
留一法分析（Leave-One-Out Analysis）：评估单个SNP对汇总效应的影响。
漏斗图（Funnel Plot）：评估工具变量的分布是否对称，判断是否存在偏倚。
MR完毕，我们将结果可视化如图：
（由于这个研究研究了很多种甲状腺功能决定因素和房颤的关系，所以结果并不止这一张图，我们这里只放促甲状腺激素-房颤的图意思意思，其他研究对的结果大家自己看原文哈）
如大家所见，森林图是结果可视化中最最重要的图，为了说明的清晰性，作者把原数据和森林图放在了一起，这样阅读更加方便。
我们看到，森林图中结果的绿色方形代表的是优势比OR，以及其代表百分之九十五置信区间（95%CI）的穿过方形的小线条。由于OR和beta具有换算关系（β=lnOR），所以方形和小线条没碰到代表OR=1的判定线，就等价于β不为零，也就是可以认为这张图中的促甲状腺激素含量和房颤确实有统计学上的因果关系。并且，因为方形在判定线左边，所以结论就是——促甲状腺激素每增加 1 个标准差，与房颤呈负相关！
肯定会有朋友问留一法图和漏斗图跑哪儿去了，没看见啊？
答案是由于不太重要所以被作者扔进了补充材料。惨。）
补充材料去哪里找呢？文章中并没有给出链接，不过，搜索这篇文章发表的期刊对应的网站，找到它，就可以找到补充材料：
这篇文章来源于期刊JAMA Cardiology（网址：https://jamanetwork.com/journals/jamacardiology
），输入文章标题或 DOI（10.1001/jamacardio.2018.4635）进行检索，在文章页面中查找补充材料的链接，一般会有 “Supplement”“Supplementary Material” 等字样的明显标识，找就完了。
9.4 关于研究结果
原文冗杂，大家自己慢慢读，结果大致就是下面这些（
9.4.1 研究水平分析（多变量调整的观察性研究水平数据的荟萃分析）
FT4 与房颤关系：FT4 对新发房颤的合并风险比为 1.55 ，对现患房颤的合并优势比为 2.80。但在 MR 分析中，FT4（每标准差）对新发或现患房颤的工具估计均无统计学意义，总体风险比为 0.84，优势比为 1.32。
促甲状腺激素与房颤关系：在观察性研究水平数据的随机效应荟萃分析中，促甲状腺激素对新发房颤的合并风险比为 1.00，对现患房颤的合并优势比为 0.99。
9.4.2 汇总水平分析
FT4 与房颤关系：参考范围内 FT4 每标准差的 IVW - RE 优势比为 1.01。加权中位数分析显示，FT4 每标准差的房颤优势比为 0.87 。
FT3:FT4 比值与房颤关系：基因预测的 FT3:FT4 比值每增加 1 个标准差，房颤风险增加，优势比为 1.33。
促甲状腺激素与房颤关系：基因预测的促甲状腺激素每增加 1 个标准差，与房颤呈负相关，IVW - RE 优势比为 0.88。
甲状腺功能减退与房颤关系：基因预测的甲状腺功能减退与房颤呈负相关，IVW - RE 优势比为 0.94。
甲状腺过氧化物酶抗体与房颤关系：基因预测的甲状腺过氧化物酶抗体每增加 1 个标准差，在 MR 分析中与房颤无显著关联。
甲状腺功能亢进与房颤关系：基因预测的甲状腺功能亢进与房颤相关，IVW - RE 优势比为 1.04，MR - Egger 优势比为 1.31，存在水平多效性（P = 0.045）。
那么，总结一下:
不好看懂是吧……翻译成中文：基因上增加的FT₃:FT₄比率和甲状腺功能亢进（但不包括参考范围内的FT₄）与增加的房颤（AF）相关，而参考范围内增加的甲状腺素和甲状腺功能减退与减少的房颤相关，这支持了涉及垂体-甲状腺-心脏轴的通路。
我们可以，发现观察性研究是支持FT4和房颤的有关系的，但是这种关联没有在MR上得到验证。
没有一个明确的结果就意味着研究失败了吗？当然不是，看似矛盾的结果反而给了我们更多的思考空间，促使我们对结果进行反思并尝试解释，这是一次研究很重要的组成部分，也是一篇好文章必须有的部分。
9.5 所以最后还要进行一些反思和研究讨论
9.5.1 因果关系的复杂性：研究结果表明，甲状腺功能与房颤之间的因果关系很复杂。FT4水平本身并未显示出直接因果关系，而FT3:FT4比值、促甲状腺激素水平以及甲状腺功能减退症和亢进症则显示出与房颤的因果关联。这是观察性研究无法得出的结果，也进一步拓宽了我们对它们之间的认识。
9.5.2 文章结果联系到过往研究中甲状腺-垂体-心脏轴的作用：这个研究支持甲状腺功能通过垂体-甲状腺-心脏轴影响房颤易感性的假设，但具体因果因素仍需进一步研究。
9.5.3 研究局限性：研究样本主要来自欧洲人群，结果的外推性可能有限。这里要联系到我们之前提到的人群异质性现象，也就是这个结果是在欧洲人数据中得出的，可能不适合亚洲人也不适合美洲人。
此外，研究中对FT3:FT4比值的工具变量较少，限制了对FT3作用的深入分析。
9.6一些关于读paper的题外话
其实除去英语看不懂导致读paper容易读破防之外，只要掌握文章的核心脉络，就可以大致看懂文章并学到知识。多读一些相应的研究，每个人都能偷偷地成为MR研究者。
考虑到咱小木屋本身的情况，大家可能看到这篇推文的时候可能还处于没读过多少paper的时期，（虽然写这篇推文的人读得也不一定多，悲）
但是这里有一些可能有帮助的tips：
看文章先看abstract。你可以通过这个部分基本判定文章是否内容符合自己研究学习的要求，同时大概get到文章的脉络，这将对把握整篇文章很有帮助。
看不懂英语不要紧，首先现在翻译足够发达（这里推荐翻译器用DeepL，比较准确，不是卖广告，虽然我倒也希望它给我一点广告费）
但是要注意的是尽可能不要总是依赖翻译器，否则读个十年八年的文章那也还是看不懂英语。对于同一个领域的研究paper来说，总有许多专有名词会反复出现，多看几遍往往就有印象了。
有一些读paper或者下载paper相关的软件可以了解一下，比如Zotero和Endnote，使用教程请移步b站大学
合理使用人工智能。这个时代怎么可能不用ai呢……问题是怎么用会好一点？答案或需要自己探索。不过现在别老用ai写文章，人机味冲的话很容易被看出来的。
略读十篇不如精读一篇。这句话是王师兄送的。走马观花的学习方式可能并不适合大多数人，囫囵吞枣往往导致事倍功半。完整地弄明白一篇paper，熟记于心深究细节甚至动手复现，或许是更好的选择。
十、R语言示例代码
那现在大家已经基本具备了复现一篇MR文章的能力了，万事俱备，只欠东风，不如动手试一试吧？阴暗地（不是）从网上下载两个GWAS数据，打开Rstudio，copy一下这份小编写了超级多注释的本地双变量孟德尔随机化（TSMR）代码——
GWAS网站的介绍和数据下载方式在这篇推文里没有啦，大家要不看看咱们小木屋其他推文呢，说不定到时候会出（
准备工作：提前装几个R包，装完一次之后就不需要这几行代码了（运行时需要把“#”删掉，下同）
## install.packages("devtools")##devtools::install_github("MRCIEU/TwoSampleMR") ##devtools::install_github("rondolab/MR-PRESSO",force = TRUE)##install.packages("vroom")#报错改为英文，方便阅读理解Sys.setenv(LC_MESSAGES="E")#查找和设置工作路径getwd()setwd("你的文件路径")#加载所需的R包library(TwoSampleMR)library(MRPRESSO)library(vroom)
第一步：数据的导入
#使用vroom 函数读取文件,col_names = TRUE表示文件的第一行包含列名(即数据的表头)，R 会把它作为数据框的一部分来读取。a<-vroom('GCST90274802.h.tsv.gz',col_names = TRUE)#如果你下载的文件是vcf格式，那么需要先解压vcf文件并使用如下代码#library(VariantAnnotation)#library(gwasvcf)#devtools::install_github("mrcieu/gwasglue")#library(gwasglue)#vcf<-readVcf("ieu-a-7.vcf")#a<-gwasvcf_to_TwoSampleMR(vcf)
第二步：数据相关性处理
##①提取有用的列+命名为标准名称----------colnames(a)##查看列名head(a)#查看原始数据a的前几行#每个下载的数据内容都有区别，我们可以用上述代码查看，后续实操中，我们要根据实际情况选择有用于后续分析的列保留，并给他们新的标准的名字，以便后续更加方便的调用分析a1 = a[,c(10,3,4,5,6,8,11)]#在这次的数据里，我们所需要的就是a中的这几列，分别对应的是"SNP"，"beta"，"se"，"eaf"，"effect_allele"，"other_allele"，"p"，"N"，提取出来赋给a1。当然，不同的数据需要提取出来的列是不一样的，需要根据自己的数据进行查看和判断。##②标准化列名："SNP"，"beta"，"se"，"eaf"，"effect_allele"，"other_allele"，"p"，"N"colnames(a1) = c("SNP","effect_allele","other_allele","beta","se","p","N") #给它们命名（这些具体的含义请看上面GWAS数据段的介绍）##③提取强相关的SNP-------b<-subset(a,p_value<5e-08) #这个时候就发现b就小很多啦##④数据转换成后缀为.exposure的标准化数据-(也就相当于正式赋予这个数据“暴露”的属性）--------colnames(b)b1<-format_data( b, type = "exposure", snps = NULL, #自动确定snp的列 header = TRUE, #数据有标题行 phenotype_col = "Phenotype", snp_col = "rsid", beta_col = "beta", se_col = "standard_error", effect_allele_col = "effect_allele", other_allele_col = "other_allele", pval_col = "p_value", eaf_col = "effect_allele_frequency")
第三步：数据独立性处理(clumping)
##⑤去除LD----------exp<-clump_data(b1,clump_kb = 10000,clump_r2 = 0.001#这里可能会会报错，是因为没有配置token,要用下面的代码配置#自己的api的申请方法详见：#https://mp.weixin.qq.com/s/nLiIHvzIPZnIg0COxXFM1g#“你自己的token”换成申请到的token字符串，注意要加上双引号Sys.setenv(opengwas_jwt="你自己的token")#检查一下有没有配置成功ieugwasr::get_opengwas_jwt()#成功的会提醒你:由于LD，Removing了若干个viriants###同理啦，读取结局数据---------library(vroom)c<-vroom('finngen_R10_E4_GRAVES_STRICT.gz',col_names = TRUE)##转换为标准化列名##标准化列名："SNP"，"beta"，"se"，"eaf"，"effect_allele"，"other_allele"，"p"，"N"colnames(c)out = c[,c(3,4,5,7,9,10,11)]colnames(out) = c("other_allele","effect_allele","SNP","p","beta","se","eaf")
第四步、数据的进一步处理
##3.暴露与结局取交集，确定哪些SNP在两者中共有-----------snps<-intersect(exp$SNP,out$SNP)#寻找两者共有的snpb<-as.data.frame(cbind(snps,1)) #找到共同的snp并且用“1"这个符号标注##4.将交集SNP的信息从暴露、结局中提取出来---------##%in%是啥意思【a1 %in% a2表示的是向量a1中，同时出现在向量a2中的元素】#即：提取 exp 数据框中，SNP列与 b$snps 中共同的SNP对应的行数据，并赋值给exp1exp1<-subset(exp,SNP %in% b$snps)colnames(out)#对outcome数据进行标准化格式转换，将数据转换成结局格式out1<-format_data( out,##拥有“标准化列名”的原始data的名称 type = "outcome", snps = b$snps, header = TRUE, phenotype_col = "Phenotype", snp_col = "SNP", beta_col = "beta", se_col = "se", eaf_col = "eaf", effect_allele_col = "effect_allele", other_allele_col = "other_allele", pval_col = "p",)
第五步：等位基因对齐--harmonise暴露【exp1】和结局【out1】
#目的是①协调不同向的snp（A-G和G-A）②去除回文序列③踢出不兼容的snp-----------dat<-harmonise_data(exposure_dat = exp1,outcome_dat = out1,action = 3)#dat<-dat[-1,]
第六步：正式MR分析
result <- mr(dat)result #这里可以看一下分析出来的结果OR <-generate_odds_ratios(result) #将β值（效应量）转换成OR值（包含其置信区间），利于后续研究#使用MR-PRESSO检验，以识别和纠正工具变工具变量的多重效应），并且检测数据是否存在失真。（分析时对象里至少包含3个以上工具变量）MRPEO <- mr_presso(BetaOutcome ="beta.outcome", BetaExposure = "beta.exposure", SdOutcome ="se.outcome", SdExposure = "se.exposure", OUTLIERtest = TRUE,DISTORTIONtest = TRUE, data = dat, NbDistribution = 1000, SignifThreshold = 0.05)#为每个SNP单独估算暴露对结局的因果效应#默认Wald比值res_single <- mr_singlesnp(dat)ORR <-generate_odds_ratios(res_single) #每个SNP暴露对结局的影响#如果不想使用默认的方法，首先可以查看一下当前支持的mr方法：mr_method_list()#选择你喜欢的mr方法，进行指定并运算，比如：result_add<-mr(dat,method_list=c("mr_ivw","mr_egger_regression","mr_weighted_median"))result_add
第七步：数据敏感性分析
#异质性检验het <- mr_heterogeneity(dat)het#多效性检验pleio <- mr_pleiotropy_test(dat)pleio#逐个剔除检验single <- mr_leaveoneout(dat)mr_leaveoneout_plot(single)
第八步：画图
#散点图mr_scatter_plot(result,dat)#森林图res_single <- mr_singlesnp(dat)mr_forest_plot(res_single)#漏斗图mr_funnel_plot(res_single)
有时候，在数据处理的过程中我们会遇到一些别样的问题，需要具体情况具体解决哦：
rsid缺失、GRCh37与GRCh38的互相转化：相关问题的解决##
以下部分比较难，新手可以先跳过哦
#缺失rsid的解决办法#假设我们已经拥有数据框data，包含有SNP的chr，所在位点pos（可能命名为position、base pair location，但本质一样）等信息#第二步 安装需要使用的包，注意可能需要使用到bioconductor::install，这需要我们提前install并library bioconductor这个包#注：此处我们针对的是GRCh37的版本，如果data中的SNP信息是GRCh38版本的，则需要更换包数据源，或采取本地方法，或提前先将38版本的信息转为37版本#第三种方法不推荐，因为可能引入缺失因素，导致工具变量减少library(BSgenome.Hsapiens.1000genomes.hs37d5)library(SNPlocs.Hsapiens.dbSNP144.GRCh37)library(MungeSumstats)library(data.table)##要将原始数据CHR,BP,A1,A2,P,EAF,BETA做好格式标准化library(dplyr)colnames(data)new_data <- data %>% rename(CHR = chromosome, BP = base_pair_location, A1 = effect_allele, A2 = other_allele , P = p_value , EAF =effect_allele_frequency , BETA = beta)reformatted <- format_sumstats( new_data, dbSNP = "144", ref_genome = "GRCh37", nThread = 4, # 多线程，根据你的电脑配置（4的话对电脑配置要求比较高，一般2或者3就差不多了） save_path = "你想存的文件的名字.tsv.gz", return_data = FALSE # 如果是TRUE，直接返回到R右边的环境)# 文件已生成到当前工作路径，名字是"你想存的文件的名字.tsv.gz"#除此之外，若能够获得GRh37/38版本的SNP信息，通常是拥有chr、pos、rsid三列的“对照表”类文档，那么可以使用下面的代码#这个代码对原始数据的要求更小，更加便利final <- merge(rawdata,GRCh37or38, by.x.= c("chr","pos"), by.y.= c("chromosome","position"))#注意参数根据具体文档内容进行调整即可
#GRCh37与GRCh38的互相转化
#GRCh37与GRCh38的互相转化#这个操作并不常见，但是若确有需要，可以按照下面的方法进行#注意，这个方法会有一定的缺失率，所以不到万不得已不要使用~#第一步 加载包install.packages("BiocManager")BiocManager::install("rtracklayer")library(rtracklayer)#第二步 下载两个版本分别的链文件，在这里获得：https://genome.ucsc.edu/cgi-bin/hgLiftOver#第三步 转化# 假设有 GRCh37 的 SNP 位点信息grch37_snps <- data.frame( chr = c("chr1", "chr2", "chr3"), pos = c(123456, 7891011, 11223344))# 创建 GRanges 对象library(GenomicRanges)grch37_gr <- GRanges( seqnames = grch37_snps$chr, ranges = IRanges(start = grch37_snps$pos, end = grch37_snps$pos))# 加载 liftOver 链文件（GRCh37 -> GRCh38）chain <- import.chain("hg19ToHg38.over.chain")# 使用 liftOver 转换坐标grch38_gr <- liftOver(grch37_gr, chain)# 提取转换后的位点信息grch38_snps <- data.frame( chr = seqnames(grch38_gr), pos = start(grch38_gr))# 查看转换后的 SNP 信息print(grch38_snps)#或者是：# 假设有 GRCh38 的 SNP 位点信息grch38_snps <- data.frame( chr = c("chr1", "chr2", "chr3"), pos = c(987654, 10111213, 22334455))# 创建 GRanges 对象grch38_gr <- GRanges( seqnames = grch38_snps$chr, ranges = IRanges(start = grch38_snps$pos, end = grch38_snps$pos))# 加载 liftOver 链文件（GRCh38 -> GRCh37）chain <- import.chain("hg38ToHg19.over.chain")# 使用 liftOver 转换坐标grch37_gr <- liftOver(grch38_gr, chain)# 提取转换后的位点信息grch37_snps <- data.frame( chr = seqnames(grch37_gr), pos = start(grch37_gr))# 查看转换后的 SNP 信息print(grch37_snps)
希望看到这里的朋友，即便还没有真的开始动手实操，也已经有了“我可以试着跑通MR”的自信啦。
我们最基础的东西似乎就讲完了，所以来看看MR的其他类型和应用，希望能对大家今后的研究有更多启发！
十一、MR的类别与应用
上面我们只是了解了一种最简单的MR，也就是TSMR，事实上MR还有许多变体，在实际研究中也可以大显身手，不妨再一起认识一下它们：
11.1 两样本孟德尔随机化TSMR
上文介绍的几乎都是TSMR的相关内容。为什么叫“两样本”呢？这是因为在MR分析中，遗传工具变量与暴露变量的关系来自一个样本，而遗传工具变量与结局变量的关系来自另一个独立样本。它是最基础的MR类型，通过工具变量只研究一个暴露和一个结果之间的因果关系。只有熟练掌握之后才能更好地进行其他变种的研究哦~
突击检查！MR-PRESSO的作用是什么，还记得吗？
不记得了要回头看看及时复习~
11.2 中介孟德尔随机化
有一天我们“突然”从之前的研究中发现：吸烟（暴露因素）会导致冠心病（结局）的发生，我们觉得很好奇，这样的因果关系，是否是通过吸烟致使的炎症反应，从而间接导致冠心病的发生呢？
想要知道“炎症反应”在暴露和结局之间扮演的角色，就需要用到中介MR。
（如图，C是总效应，加上中介变量后，A*B是间接效应，C'是直接效应，有一条公式：总效应=间接效应+直接效应。如果总效应、直接效应和间接效应方向都相同的情况下，我们还可以算出中介效应比例（间接效应/总效应））
方法（两步法MR Two-Step MR）：
第一步：确定存在C这个因果效应
第二步：通过常规MR分析步骤确定分别存在A、B两个因果关系
第三步：
当A、B、C均显著时，说明暴露到结局存在因果关联，AB可作为中介效应，C -（AB）作为直接效应；
若C不显著，但A和B显著，表明暴露到结局的关联完全由中介变量介导；
若C显著，但A和B至少有一个不显著，说明不存在该中介变量介导的中介效应。
*广义的中介MR还包含MVMR（多变量MR）
MVMR将多个表型作为暴露，共同探讨其与结局的关系。通过评估不同表型暴露的占比来确定中介效应。
11.3 药靶MR
首先，药物与机体生物大分子的结合部位即药物靶点。我们知道，常见的生物大分子有蛋白质、核酸等，它们所参与构成的离子通道、受体、酶、DNA，都可以是药物作用的靶点，故而这些生物大分子的相关特性也就可以作为暴露因素与相关的疾病进行MR分析，这就是药靶MR。
这个研究方法的特点：
常见的药靶MR专注于研究蛋白与其上游的基因，并采用eQTL与pQTL（可以简单理解成基因表达和蛋白表达的相关SNP位点）作为暴露因素与疾病作为结局进行MR分析。
相较于经典MR，药靶MR更加强调认准致病机制上某一个蛋白质或者其标志物，研究它在导致整个疾病通路过程中扮演的角色，并对其进行精准的干预调控；而这显然是在每一个基因座上都选一个SNP的传统MR所不能做到的。
11.4 本地、在线MR的联系与区别
11.4.1 联系
本地、在线MR的原理、核心方法都是一致的，在实操过程中，两者的结果可以相互补充和验证，以提升研究结果的可靠性和稳健性
11.4.2 区别
本地MR需要先行下载相关数据（如GWAS数据）到本地，在RStidio上分析时，暴露、结局的数据需要先使用代码读取进来；有时，由于读取进来的数据框并不是数据框的格式，还需要将数据整理成数据框格式，并对文件的列进行命名，才可以进行后续分析。
在线MR是通过网络直接连接到数据库，实时获取数据（例如利用API读取IEU OpenGwas 数据库的数据），利用extract_outcome(exposure)_data函数，提取到相应id的数据。
11.4.3 优缺点分析（选择时的注意事项）
本地MR不依赖网站和网络环境，相同时间内可以进行更大量的分析，是更加自由，灵活的分析选择；但是本地MR要求提前下载数据，下载的过程可能耗时较长，且对设备内存有要求，且导入、处理数据的过程都容易导致报错的出现从而妨碍后续进程。
在线MR提取数据方便简单，前期准备工作更加轻松，是更加方便、快捷的选择；然而在线MR的运行十分依赖网站和网络环境，如果出现网站崩溃，或网站给的token数的限制，将会使分析停滞不前，其分析的稳定性不如本地的MR
——需要注意一下，在线mr也往往需要申请token码，相当于准入凭证，然后将它导进代码里面在Rstudio跑一下，才能进行在线分析，而这个码是有使用时限的，过期了要重新申请，，，代码还是用下面这两行（刚刚本地进行LD的时候其实已经用过了）
#载入一下token(ieu在线数据需要，十四天有效) Sys.setenv("OPENGWAS_JWT"="你的token") #看token是不是生效了 ieugwasr::get_opengwas_jwt()
代码片段：可切换语言，无法单独设置文字格式
如果token生效了长这样↓↓这堆乱码是我的token码，不用在意
数据的导入
#使用vroom 函数读取文件,col_names = TRUE表示文件的第一行包含列名(即数据的表头)，R 会把它作为数据框的一部分来读取。a<-vroom('GCST90274802.h.tsv.gz',col_names = TRUE)#如果你下载的文件是vcf格式，那么需要先解压vcf文件并使用如下代码#library(VariantAnnotation)#library(gwasvcf)#devtools::install_github("mrcieu/gwasglue")#library(gwasglue)#vcf<-readVcf("ieu-a-7.vcf")#a<-gwasvcf_to_TwoSampleMR(vcf)
（这里的AD和depression都是疾病的名字，可以改的；在这个文档中看代码不太好看，建议复制粘贴到Rstudio中看）
##extract_instruments 函数：这是一个用于提取GWAS数据中工具变量的函数。工具变量是从GWAS研究中提取的与暴露显著相关的SNPs。##outcomes = "bbj-a-86"：指定要提取工具变量的GWAS研究ID。bbj-a-86是日本生物银行（BioBank Japan）的一个特定研究ID。##AD是一个数据框，包含提取的工具变量及其相关信息。AD <- extract_instruments(outcomes = "bbj-a-86", p1 = 5e-08, clump = TRUE, r2 = 0.001, kb = 10000)View(AD)##extract_outcome_data 函数：用于提取与工具变量相关的结果数据。结果数据通常来自另一个GWAS研究。##depression是一个数据框，包含与工具变量相关的抑郁症结果数据。depression <- extract_outcome_data(snps = AD$SNP, outcomes = "ukb-a-552", proxies = TRUE, rsq = 0.8, align_alleles = 1, palindromes = 1, opengwas_jwt = ieugwasr::get_opengwas_jwt(), splitsize = 10000, proxy_splitsize = 500)View(depression)##接下来进行harmonizedata<-harmonise_data(exposure_dat=AD,outcome_dat=depression,action=2)#检查结局的p值，p值小的要剔除掉View(data)#保存harmonization的数据write.csv(data,file = "data_harmonization.csv")##然后就可以正常对data进行mr分析了
不过嘛，跑在线有时候会碰到一些很难解决的问题……
（下图报错，显然就是在线MR链接崩了导致的无法在线读取GWAS数据，碰到这种问题只好暂时痛苦面具一下等它自己修好，或者换本地MR跑吧）
在最终进行MR分析的时候，要根据自己当时的需求和情况合理选择～
写在后面
到这里我们能想到的东西就讲的差不多了，实际应用中还会有各种奇奇怪怪的问题等着大家解决，这时候就到了发挥各位聪明才智的时候啦。我们把乱七八糟各种内容都塞了一点进来，希望刚进入小木屋的朋友们从这里开始对MR有一点粗浅的理解，少点迷茫，少走弯路。