安装

第一次使用的时候需要安装

install.packages(‘devtools’)

# 安装最新版的TwoSampleMR包
devtools::install_github(‘MRCIEU/TwoSampleMR’)

# 安装指定版本的TwoSampleMR包
devtools::install_github(‘MRCIEU/TwoSampleMR@0.4.26’)

加载安装好的包

1	library(TwoSampleMR)

显示以下输出

TwoSampleMR version 0.5.6 
[>] New: Option to use non-European LD reference panels for clumping etc
[>] Some studies temporarily quarantined to verify effect allele
[>] See news(package='TwoSampleMR') and https://gwas.mrcieu.ac.uk for further details

读取暴露文件

有两种读取暴露文件的方法

使用TwoSampleMR获取MR base提供的数据
使用TwoSampleMR包读取本地文件

使用TwoSampleMR获取MR base提供的数据

因为MR base提供了很多数据，因此需要选择要使用的数据，也就是暴露的ID。假设我们选择的暴露是BMI，它的id是ieu-a-2

1 2	bmi <- extract_instruments(outcomes='ieu-a-2', access_token=NULL) # 获取暴露数据 head(bmi) # 查看暴露数据

以下的还不知道是什么意思

关于extract_instruments()的使用，有几个参数需要大家注意一下：

（1）第一个就是access_token这个参数，对于中国大陆地区的用户必须设置该参数为access_token=NULL，这样才能顺利获取数据，否则就需要开VPN获取谷歌授权。

（2）第二个是参数p1，它是用来指定暴露中SNP的p值的，它的默认值是p1=5e-8，因此只有p值小于5e-8的SNP才会提取出来。当然如果没有SNP小于5e-8的话，我们通常可以设置p1=1e-5，不过这个时候就需要认真评估弱工具变量偏倚了。

（3）第三个重要参数是clump以及与之相关的r2和kb，clump是一个逻辑型参数，只有clump=TRUE和clump=FALSE这两种情况。如果选择了参数值为clump=FALSE的话，那么r2和kb这两个参数就无效了，也即先不去除含有连锁不平衡的SNP。当clump=TRUE时，我们可以用r2和kb来确定去除连锁不平衡SNP的条件，具体内容我会在下期内容中进行详细讲解。默认情况下是clump=TRUE，r2=0.001和kb=10000。

使用TwoSampleMR包读取本地文件

可以看这篇推送。我自己就还没试过

去除连锁不平衡(LD)

我们可以选择在提取数据的时候去除连锁不平衡的IV，也可以先提取数据，再去除。效果是完全一样的。

背景知识

在做MR研究时，我们有一步非常重要，那就是去除存在连锁不平衡的IV。连锁不平衡主要是使用两个参数r2和kb来衡量：

（1）r2：它是0~1之间的数据，r2=1表示两个SNP间是完全的连锁不平衡关系，r2=0则表示两个SNP间是完全连锁平衡的，也即这两个SNP的分配是完全随机的。

（2）kb：它其实就是指考虑连锁不平衡的区域长度。在遗传学上我们认为在染色体上距离很近的遗传位点通常是“捆绑”在一起遗传给后代的，这也就导致距离很近的位点之间的r2会很大。在TwoSampleMR包中，我们去除连锁不平衡主要考虑的也是这两个参数。

举个例子，如果设置r2=0.001和kb=10000，那这就表示去掉在10000kb范围内与最显著SNP的r2大于0.001的SNP；而设置成r2=0.3和kb=1000，那这表示的就是去掉在1000kb范围内与最显著SNP的r2大于0.3的SNP。

从上面的解读中不难看出，随着r2的变小与kb的变大，被去除的存在连锁不平衡的SNP会越来越多，而最终剩下的IV会越来越少。在我之前的推送中曾以实例和大家探讨了IV数目对结果的影响（孟德尔随机化之高密度脂蛋白胆固醇（HDL-C）与心肌梗死的因果关系），一般IV个数越少，存在的混杂和多效性也就越少，但相应的统计效力不足；而IV数目的增多虽然能提高统计效力，但也会带来更多的偏倚。因此，调节好参数r2和kb就显得是一门技术活了！

在提取时去除LD

1	bmi <- extract_instruments(outcomes='ieu-a-2', clump=TRUE, r2=0.01, kb=5000, access_token=NULL)

先提取数据，再去除LD

1 2	bmi <- extract_instruments(outcomes='ieu-a-2', clump=FALSE, access_token=NULL) exp_dat <- clump_data(bmi, clump_r2=0.01, clump_kb=5000)

提取IV在结局中的信息

在读取完暴露文件并去除掉存在连锁不平衡的SNP后，我们接下来要做的一件事就是提取IV在结局中的信息，完成这一步主要有两种方法：

（1）利用TwoSampleMR获取MR base提供的结局信息

（2）读取自己结局的GWAS文件并提取相关信息

第一种方法使用起来非常简洁高效，可以批量读取多个结局文件，但是存在的问题是有的结局数据可能有问题（真的吗）；

第二种方法一次读取一个GWAS文件，如果批量处理的话可能会占用大量内存，得不偿失。接下来我将为大家详细介绍一下这两种方法，希望大家能明白这两种读取方法的差异。

读取自己结局的GWAS文件并提取相关信息

首先咱们先提取IV的信息并去除存在连锁不平衡的SNP，这里咱们还是以BMI作为暴露，但是ID号需要改成'ieu-a-835'，这主要是因为之前ID号’ieu-a-2’的GWAS是在混合人群中做的（也即把欧洲人、非洲人等不同人群合在一起做的GWAS），而’ieu-a-835’则是在欧洲人中做的。在之前的理论学习中，我曾和大家解释过人群的混杂会带来估计结果的偏倚，因此我们需要选择遗传背景一致的人群进行MR研究（如暴露和结局的GWAS都是在欧洲人群中进行的）。


library(TwoSampleMR)
bmi_exp <- extract_instruments(
 outcomes='ieu-a-835',
 clump=TRUE, r2=0.01,
 kb=5000,access_token = NULL
  )
dim(bmi_exp)
# [1] 80 15
t2d_out <- extract_outcome_data(
 snps=bmi_exp$SNP,
 outcomes='ieu-a-26',
 proxies = FALSE,
 maf_threshold = 0.01,
 access_token = NULL
)
dim(t2d_out)

这里我要和大家简单介绍一下extract_outcome_data()函数的关键参数：

snps：它是一串以rs开头的SNP ID

outcomes：它是outcome在MR base中的ID；

proxies：它表示是否使用代理SNP，默认值是TRUE，也即当一个SNP在outcome中找不到时可以使用与其存在强连锁不平衡的SNP信息来替代，我个人喜欢设置成FALSE。

maf_threshold：它表示的是SNP在outcome中的最小等位基因频率，默认值是0.3，不过大样本GWAS可以适当调低，我这里设置的是0.01。

access_token：大陆用户必须设置成access_token=NULL。

利用TwoSampleMR获取MR base提供的结局信息

见推送

计算并解读MR结果

mydata <- harmonise_data(
       exposure_dat=bmi_exp,
       outcome_dat=t2d_out,
       action= 2
       )

最后，咱们只要简单使用mr()函数即可，代码如下：

1 2	res <- mr(mydata) res

除了上述5种计算方法外，TwoSampleMR包还提供了很多计算方法，比如随机效应模型和固定效应模型等，感兴趣的朋友可以通过mr_method_list()函数来了解。

敏感性分析

（1）异质性检验（Heterogeneity test）：主要是检验各个IV之间的差异，如果不同IV之间的差异大，那么这些IV的异质性就大。

（2）多效性检验（Pleiotropy test）：主要检验多个IV是否存在水平多效性，常用MR Egger法的截距项表示，如果该截距项与0差异很大，说明存在水平多效性。除此以外，MR-PRESSO包也是常用的检验水平多效性的R包。

（3）逐个剔除检验（Leave-one-out sensitivity test）：主要是逐个剔除IV后计算剩下IV的MR结果，如果剔除某个IV后其它IV估计出来的MR结果和总结果差异很大，说明MR结果对该IV是敏感的。

coding-TwoSampleMR

安装