分计算iv值_TwoSampleMR实战教程之提取IV在结局中的信息

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

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

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

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

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

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

library(TwoSampleMR)bmi_exp outcomes&＃61;&＃39;ieu-a-835&＃39;, clump&＃61;TRUE, r2&＃61;0.01, kb&＃61;5000,access_token &＃61; NULL )dim(bmi_exp)# [1] 80 15t2d_out snps&＃61;bmi_exp$SNP, outcomes&＃61;&＃39;ieu-a-26&＃39;, proxies &＃61; FALSE, maf_threshold &＃61; 0.01, access_token &＃61; NULL)dim(t2d_out)# [1] 80 16

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

snps&＃xff1a;它是一串以rs开头的SNP ID

outcomes&＃xff1a;它是outcome在MR base中的ID&＃xff1b;

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

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

access_token&＃xff1a;大陆用户必须设置成access_token&＃61;NULL。

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

米老鼠从DIAGRAM研究中下载了与&＃39;ieu-a-26&＃39;对应的完整GWAS数据然后提取IV&＃xff0c;代码如下&＃xff1a;

#install.packages(&＃39;data.table&＃39;) 安装data.table包library(data.table) # 加载R包t2d &＃39;DIAGRAMv3.2012DEC17.txt&＃39;,header&＃61;T) head(t2d) # 查看数据

t2d$phenotype &＃39;Type 2 diabetes&＃39; t2d$beta log(t2d$OR) t2d$se abs(t2d$beta/qnorm(t2d$P_VALUE/head(t2d) # 查看数据

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

由于原始的GWAS结果中没有phenotype、beta和se的信息&＃xff0c;因此米老鼠先将它读取到R中&＃xff0c;然后转换格式。米老鼠这里是先把原始的GWAS使用data.table包的fread()函数读到R中&＃xff0c;因为这个fread()函数读取大文件的速度非常快&＃xff0c;接着我再使用format_data()函数将该数据框转化成TwoSampleMR的格式&＃xff0c;关于format_data()函数的用法参见往期内容TwoSampleMR包实战教程之读取暴露文件。

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

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

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

最后&＃xff0c;米老鼠整理好了MR base里以IEU开头的outcome信息&＃xff0c;有兴趣的朋友可以私聊我获取。DIAGRAM consortium的GWAS数据库网站可以点击阅读原文后下载第一个文件即可&＃xff0c;当然也可以私聊米老鼠获取。