TwoSampleMR实战教程之提取IV在结局中的信息

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

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

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

第一种方法使用起来非常简洁高效，可以批量读取多个结局文件，但是存在的问题是有的结局数据可能有问题（米老鼠做研究的过程确认过）；第二种方法一次读取一个GWAS文件，如果批量处理的话可能会占用大量内存，得不偿失。接下来我将为大家详细介绍一下这两种方法，希望大家能明白这两种读取方法的差异。

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

首先咱们先提取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 15t2d_out <- extract_outcome_data( snps=bmi_exp$SNP, outcomes='ieu-a-26', proxies = FALSE, maf_threshold = 0.01, access_token = NULL)dim(t2d_out)# [1] 80 16

这里我要和大家简单介绍一下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。

2. 从自己的GWAS结果中提取IV在结局中的信息

米老鼠从DIAGRAM研究中下载了与'ieu-a-26'对应的完整GWAS数据然后提取IV，代码如下：

#install.packages('data.table') 安装data.table包library(data.table) # 加载R包t2d <-fread('DIAGRAMv3.2012DEC17.txt',header=T) # 使用fread()快速读取大文件head(t2d) # 查看数据

t2d$phenotype <- 'Type 2 diabetes' # 添加phenotype列t2d$beta <- log(t2d$OR) # 转换OR值为betat2d$se <-abs(t2d$beta/qnorm(t2d$P_VALUE/2,lower.tail=F)) # 计算sehead(t2d) # 查看数据

t2d_out <- format_data( dat=t2d, type = "outcome", snps = bmi_exp$SNP, header = TRUE, phenotype_col = "phenotype", snp_col = "SNP", beta_col = "beta", se_col = "se", effect_allele_col = "RISK_ALLELE", other_allele_col = "OTHER_ALLELE", pval_col = "P_VALUE", ncase_col = "N_CASES", ncontrol_col = "N_CONTROLS", chr_col = "CHROMOSOME", pos_col = "POSITION")head(t2d_out)

由于原始的GWAS结果中没有phenotype、beta和se的信息，因此米老鼠先将它读取到R中，然后转换格式。米老鼠这里是先把原始的GWAS使用data.table包的fread()函数读到R中，因为这个fread()函数读取大文件的速度非常快，接着我再使用format_data()函数将该数据框转化成TwoSampleMR的格式。

当然，如果各位小伙伴们的文件里的信息很完整，可以考虑使用read_outcome_data()函数，它的用法与read_expsoure_data()类似，具体可以使用?read_outcome_data查询。

大家可以自行比对一下本文中的两个方法得出的t2d_out结果，米老鼠看了一下，他们的beta和p值是一致，但se有微小差异，这是由于精确位数不同导致，不会对结果产生实质影响。

另外，今天的内容里还介绍了OR和beta的转换以及se的计算，这个非常重要，希望大家掌握！

打赏