TwoSampleMR实战教程之提取IV在结局中的信息
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 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)
# [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值为beta
t2d$se <-abs(t2d$beta/qnorm(t2d$P_VALUE/2,lower.tail=F)) # 计算se
head(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的格式，关于format_data()函数的用法参见往期内容TwoSampleMR包实战教程之读取暴露文件。
当然，如果各位小伙伴们的文件里的信息很完整，可以考虑使用read_outcome_data()函数，它的用法与read_expsoure_data()类似，具体可以使用?read_outcome_data查询。
大家可以自行比对一下本文中的两个方法得出的t2d_out结果，米老鼠看了一下，他们的beta和p值是一致，但se有微小差异，这是由于精确位数不同导致，不会对结果产生实质影响。
另外，今天的内容里还介绍了OR和beta的转换以及se的计算，这个非常重要，希望大家掌握！
最后，米老鼠整理好了MR base里以IEU开头的outcome信息，有兴趣的朋友可以私聊我获取。DIAGRAM consortium的GWAS数据库网站可以点击阅读原文后下载第一个文件即可，当然也可以私聊米老鼠获取。