读取暴露因素数据的方法

1. 在线读取获取数据

library("TwoSampleMR")
library("ieugwasr")
bmi_exp_dat <- extract_instruments(outcomes = 'ieu-a-2')
##过滤
bmi_exp_dat <- clump_data(bmi_exp_dat)

上面的extract_instruments即从https://gwas-api.mrcieu.ac.uk/这个API地址获取对应的数据。

具体查询地址是https://gwas.mrcieu.ac.uk/，从这里可以搜索各种GWAS数据集。

2. 自己创建一个暴露数据或者直接使用一个data.frame数据框

random_df <- data.frame(
  SNP = c("rs1", "rs2"),
  beta = c(1, 2),
  se = c(1, 2),
  effect_allele = c("A", "T")
)
random_df

random_exp_dat <- format_data(random_df, type = "exposure") ##规范为MR分析的标准格式
random_exp_dat

自己创建数据集的格式至少具有4列，格式如下：

head(random_df)
 SNP beta se effect_allele
1 rs1    1  1             A
2 rs2    2  2             T

format_data函数是把这个数据框转为MR分析需要的exposure格式
如果环境中直接有一个data.frame，则可以使用format_data把该数据框转为exposure格式。

3. 从本地读取文件，构建暴露文件

bmi2_file <- system.file("extdata/bmi.csv", package = "TwoSampleMR")
bmi_exp_dat <- read_exposure_data(
  filename = bmi2_file,
  sep = ",",
  snp_col = "rsid",
  beta_col = "effect",
  se_col = "SE",
  effect_allele_col = "a1",
  other_allele_col = "a2",
  eaf_col = "a1_freq",
  pval_col = "p-value",
  units_col = "Units",
  gene_col = "Gene",
  samplesize_col = "n"
)

这里使用的是TwoSampleMR包里自带的示例数据。
其中本地文件必须的列是 snp_col， beta_col, se_col,effect_allele_col, other_allele_col,eaf_col,pval_col

4. 从其他包加载数据

library(MRInstruments) #这是一个提供数据的包
data(gwas_catalog)
head(gwas_catalog)
bmi_gwas <-  subset(gwas_catalog, grepl("Speliotes", Author) & Phenotype == "Body mass index")
bmi_exp_dat <- format_data(bmi_gwas)

这里是选择了Speliotes作者提供的BMI数据集.

代谢组数据

data(metab_qtls) #全血中121种代谢物的qtl的数据集
head(metab_qtls)
ala_exp_dat <- format_metab_qtls(subset(metab_qtls, phenotype == "Ala")) #这里是提取丙氨酸的位点

蛋白组数据

data(proteomic_qtls) #全血中47种蛋白质代谢物的qtl数据集
head(proteomic_qtls)
#这里是提取ApoH蛋白的数据
apoh_exp_dat <-  format_proteomic_qtls(subset(proteomic_qtls, analyte == "ApoH"))

eQTL表达位点

# 32432 个基因在44 个组织的基因表达水平的eqtl
data(gtex_eqtl)
head(gtex_eqtl)
#此处是提取IRAK1BP1基因在皮下脂肪组织的位点
irak1bp1_exp_dat <- format_gtex_eqtl(subset(gtex_eqtl,gene_name == "IRAK1BP1" & tissue == "Adipose Subcutaneous"))

甲基化组数据

data(aries_mqtl)
head(aries_mqtl)
#此处是提取出生时cg25212131的甲基化水平的qtl位点
cg25212131_exp_dat <- format_aries_mqtl(subset(aries_mqtl, cpg == "cg25212131" & age == "Birth"))

#ieuGWAS数据库，这个即是第一条的在线读取方法
ieugwasr::api_status()
#获取该数据库所有可用的数据的信息
ao <- available_outcomes()
head(subset(ao, select = c(trait, id)))
bmi2014_exp_dat <- extract_instruments(outcomes = 'ieu-a-2')

这个很容易报错，因为国内访问ieu数据库不稳定。

确定好使用哪种方法获得暴露数据后，使用clump_data()函数对暴露数据的SNP进行过滤，以去除LD连锁不平衡和R2，Pval等值不符合要求的SNP.

---

读取结局因素的方法

1. 从ieu数据库获取结局文件

chd_out_dat <-  extract_outcome_data(snps=bmi_exp_dat$SNP,outcomes = "ieu-a-7") 

前面的snps参数是过滤提取包含这些SNP的数据。后面的outcomes是数据库中对应的某个冠心病的编号,后面也可以一次提供多个编号

2. 从本地读取GWAS的summary文件，生成结局文件

outcome_dat <- read_outcome_data(
    snps = bmi_exp_dat$SNP,
    filename = "gwas_summary.csv",
    sep = ",",
    snp_col = "rsid",
    beta_col = "effect",
    se_col = "SE",
    effect_allele_col = "a1",
    other_allele_col = "a2",
    eaf_col = "a1_freq",
    pval_col = "p-value",
    units_col = "Units",
    gene_col = "Gene",
    samplesize_col = "n"
)

gwas_summary.csv文件的内容格式如下：

rsid,effect,SE,a1,a2,a1_freq,p-value,Units,Gene,n
rs10767664,0.19,0.030612245,A,T,0.78,5.00E-26,kg/m2,BDNF,225238
rs13078807,0.1,0.020408163,G,A,0.2,4.00E-11,kg/m2,CADM2,221431
rs1514175,0.07,0.020408163,A,G,0.43,8.00E-14,kg/m2,TNNI3K,207641
rs1558902,0.39,0.020408163,A,T,0.42,5.00E-120,kg/m2,FTO,222476

3. 可以使用其他包辅助，从vcf文件获取暴露文件
   从https://gwas.mrcieu.ac.uk/下载vcf文件到本地

wget https://gwas.mrcieu.ac.uk/files/ieu-a-2/ieu-a-2.vcf.gz
wget https://gwas.mrcieu.ac.uk/files/ieu-a-2/ieu-a-2.vcf.gz.tbi
wget https://gwas.mrcieu.ac.uk/files/ieu-a-7/ieu-a-7.vcf.gz
wget https://gwas.mrcieu.ac.uk/files/ieu-a-7/ieu-a-7.vcf.gz.tbi

#0.5从本地读取数据进行MR分析
#参考地址 https://mrcieu.github.io/gwasglue/articles/mr.html#using-gwas-vcf-files
library(ieugwasr)
library(TwoSampleMR)
#devtools :: install_github("mrcieu/gwasglue")
library(gwasglue)
library(gwasvcf)
library(dplyr)

##这个是获取表型显著的GWAS的SNP的位点信息
expd1<- extract_instruments(outcomes = "ieu-a-2")

gwasvcf::set_bcftools()

#参数chromopos提供的应该是一个格式为chr:pos的向量，这个向量是某个GWAS表型显著的SNP位点
expd2 <- gwasvcf::query_gwas("ieu-a-2.vcf.gz",chrompos=paste0(expd1$chr.exposure, ":", expd1$pos.exposure))
expd2 <- gwasglue::gwasvcf_to_TwoSampleMR(expd2, type="exposure")

outd2 <- gwasvcf::query_gwas("ieu-a-7.vcf.gz", chrompos = paste0(expd1$chr.exposure, ":", expd1$pos.exposure))
outd2 <- gwasglue::gwasvcf_to_TwoSampleMR(outd2, "outcome")

dat2 <- TwoSampleMR::harmonise_data(expd2, outd2)
mr <- TwoSampleMR::mr(dat2)
mr

注意这个从vcf文件获取暴露文件和结局文件，需要提前准备好表型的GWAS结果位点信息，提供给参数chrompos.如果你和我一样无法访问ieu数据库，则无法使用上面示例的方法expd1<- extract_instruments(outcomes = "ieu-a-2")来获取显著的snp位点信息，此时就要自己想办法拿到GWAS的结果信息，获取对应的位点信息。