TwoSampleMR包实战教程之读取暴露文件

在上期内容中，我和大家简单介绍的TwoSampleMR这个R包的主要功能和安装方法，今天我将为大家讲解使用该包进行孟德尔随机化研究的第一步------读取暴露文件。

在米老鼠的实践中，通常有两种读取暴露文件的方法：

（1）第一种是直接使用TwoSampleMR包提供的MR base数据库提供的GWAS数据，这个方法要求网络状态良好；

（2）第二种是使用自己的GWAS数据并将其读入到TwoSampleMR。

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

在使用这个方法之前，你首先需要知道这个暴露的ID号。这里我们以体质指数（body mass index，BMI）作为暴露，该暴露的ID是‘ieu-a-2’，以此获取它的相关数据。

library(TwoSampleMR) #加载R包
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。

2. 使用TwoSampleMR包读取本地文件

这里的本地文件就是你自己的GWAS文件，假设米老鼠的数据如下图所示，那么我们该如何读入TwoSampleMR中进行孟德尔随机化分析呢？

这里我们需要使用read_exposure_data()这个函数或者format_data()这个函数，这里我先介绍read_exposure_data()函数。

exp_dat <- read_exposure_data(
   filename = 'MICAD_gwas.txt',
   clump = FALSE,
       sep= "\t",
   snp_col = "SNPID",
    beta_col = "log_OR",
    se_col = "se",
    effect_allele_col ="effect_allele",
   other_allele_col = "other_allele",
   eaf_col = "effect_allele_freq",
   pval_col = "Pvalue"
)
head(exp_dat)

接下来我将和大家详细介绍各个参数：

参数filename就是用来指定文件名的，clump用来指定是否去除连锁不平衡（默认先不去除），sep用于指定文件的分割符（该文件使用tab分割），snp_col用于指定SNP所在的列名，beta_col用于指定beta值所在的列名，effect_allele_col用于指定效应等位基因列名，other_allele_col用于指定其它等位基因列名，eaf_col用于指定效应等位基因频率的列名，pval_col用于指定P值的列名。

如果我们已经把数据读入R中并希望将其转化为TwoSampleMR的格式，这时候我们需要使用format_data()，我们需要注意到format_data()的输入参数是R语言的数据框，而read_exposure_data()的输入参数是一个本地文件名，这是两者最主要的区别。

micad <-read.table('MICAD_gwas.txt',header=T)
exp_dat <- format_data(
 micad,
 type='exposure',
 snp_col = "SNPID",
  beta_col = "log_OR",
  se_col = "se",
  effect_allele_col ="effect_allele",
 other_allele_col = "other_allele",
 eaf_col = "effect_allele_freq",
 pval_col = "Pvalue"
  )

这里大部分参数和read_exposure_data()一致，不过第一个参数是R语言的数据框，第二个是指定将数据框转化成exposure格式的，type也可以设置成type=’outcome’，用于将数据转化为outcome形式。

上面就是我们常用的读取暴露文件的方法，希望能给大家带来帮助。