TwoSampleMR包实战教程之去除连锁不平衡(LD)
TwoSampleMR包实战教程之去除连锁不平衡(LD)
第一部分：背景知识
在做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就显得是一门技术活了！
第二部分：实战演练
在上期介绍读取暴露文件的内容中，我和大家使用的是”ieu-a-2”这个ID，它是一个与体质指数（body mass index，BMI）相关的GWAS结果。
在获取暴露信息时，米老鼠选取的是默认参数（clump=TRUE, r2=0.001, kb=10000），共获取了79个IV，具体代码如下：
library(TwoSampleMR)
bmi <-extract_instruments(outcomes='ieu-a-2',access_token = NULL) #用默认参数获取IV
dim(bmi) #查看IV个数
[1] 79 15 # 结果显示提取79个IV当然，我们也可以在提取数据的时候使用r2和p1这两个参数直接去除连锁不平衡，具体代码如下：
bmi <-extract_instruments(outcomes='ieu-a-2',
clump=TRUE, r2=0.01,
kb=5000,access_token = NULL
)
dim(bmi)
[1] 90 15 # 结果显示提取90个IV如果我们用的是自己的暴露文件该怎么办呢？
考虑到之前已经教过大家读取自己的表型文件了，因此这里米老鼠就直接以BMI数据为例介绍去除连锁不平衡。米老鼠这里先获取BMI的所有SNP（不进行clump），然后我们再利用clump_data()函数来去除连锁不平衡，具体代码如下：
bmi<-extract_instruments(outcomes='ieu-a-2',clump=FALSE, access_token = NULL)
dim(bmi)
[1] 2041 15 # 结果显示共提取2041个SNP
exp_dat <-clump_data(bmi,clump_r2=0.01,clump_kb=5000)
dim(exp_dat)
[1] 90 15 # 结果显示clump之后还剩90个IV这里第一步提取IV时我直接设置clump=FALSE，这样就可以把所有SNP都提取出来了。这里获取的bmi文件其实就相当于我们自己读入的文件，接下来我再使用clump_data()函数来去除连锁不平衡，这里我使用参数r2=0.01, kb=5000（个人觉得这个设置比较温和）。
在clump_data()函数中，它的clump_r2和clump_kb参数分别与extract_instruments()函数里的r2和kb一致，它有一个pop参数需要注意一下，这个参数是用来指定参考基因组的人种的，默认值为pop=’EUR’，也即认为暴露的GWAS是欧洲人群，如果咱们研究的是中国人，则可以设置pop=’EAS’。
最后，我们不难发现剩下的IV和使用如下代码的结果一致：
bmi <-extract_instruments(outcomes='ieu-a-2',
clump=TRUE, r2=0.01,
kb=5000,access_token = NULL
)关于去除连锁不平衡的方法暂时就先分享到这里，之前Ammy在简书上写过一篇关于使用TwoSampleMR时碰到的坑，大家可以点击“阅读原文”查看。