孟德尔随机化：下载IEU OPEN GWAS数据本地处理

发布日期：2024-07-22 03:08    点击次数：80

💡专注R语言在🩺生物医学中的使用

设为“星标”，精彩不错过

之前的推文详细介绍了TwoSampleMR这个R包的用法，但是这个流程是在线获取数据的，有些朋友可能会因为网络问题没法获取成功，直接卡死在第一步！哎，感叹下国内的科研环境真是难！

所以今天记录一下如何从IEU下载原始数据，然后在本地进行处理，获得和在线数据一样的结果。

为了和在线的结果一致，暴露因素还是选择BMI，它的ID是：ieu-a-2，结局还是选择CHD(冠心病)，它的ID是：ieu-a-7。

本文目录：

IEU数据下载

准备R包

读取暴露数据

根据P值过滤

去除连锁不平衡

读取结局数据

进行MR分析

IEU数据下载

打开以下网址：https://gwas.mrcieu.ac.uk/ ，点击datasets：

图片

先找暴露数据的SNP。

在GWAS ID后面的框里输入ieu-a-2，然后点击Filter，在出来的结果中点击ieu-a-2这个结果：

图片

就可以来到下载界面，点击Download VCF即可下载暴露数据的完整的GWAS数据：

图片

下载结局数据也是完全一样的流程，就不重复演示了。我这里为了和前面的在线获取进行比较，结局也选择的是ieu-a-7。

两个数据都下载后是这样的：

图片

这样一个暴露的数据，一个结局的数据，都是VCF格式的，我们就下载好了。下面就是读取数据，进行处理，看看能不能得到和在线获取的数据一样的结果。

这两个数据还是蛮大的，内存小的电脑可能会直接卡死。

准备R包

读取数据及后续进行处理需要首先安装以下几个R包：

# 记得先改镜像#options(BioC_mirror="https://mirrors.nju.edu.cn/bioconductor/")BiocManager::install("VariantAnnotation")# 对网络有要求哦remotes::install_github("MRCIEU/TwoSampleMR")remotes::install_github("mrcieu/gwasvcf")remotes::install_github("mrcieu/gwasglue")

加载一下：

library(gwasvcf)library(gwasglue)library(VariantAnnotation)library(TwoSampleMR)

读取暴露数据

使用VariantAnnotation包进行读取，直接读，非常简单，这是官方的教程：

# 注意路径不要写错vcf_exposure <- VariantAnnotation::readVcf("./ieu/ieu-a-2.vcf.gz")class(vcf_exposure)## [1] "CollapsedVCF"## attr(,"package")## [1] "VariantAnnotation"vcf_exposure## class: CollapsedVCF ## dim: 2554285 1 ## rowRanges(vcf):##   GRanges with 5 metadata columns: paramRangeID, REF, ALT, QUAL, FILTER## info(vcf):##   DataFrame with 1 column: AF## info(header(vcf)):##       Number Type  Description     ##    AF A      Float Allele Frequency## geno(vcf):##   List of length 9: ES, SE, LP, AF, SS, EZ, SI, NC, ID## geno(header(vcf)):##       Number Type   Description                                                ##    ES A      Float  Effect size estimate relative to the alternative allele    ##    SE A      Float  Standard error of effect size estimate                     ##    LP A      Float  -log10 p-value for effect estimate                         ##    AF A      Float  Alternate allele frequency in the association study        ##    SS A      Float  Sample size used to estimate genetic effect                ##    EZ A      Float  Z-score provided if it was used to derive the EFFECT and...##    SI A      Float  Accuracy score of summary data imputation                  ##    NC A      Float  Number of cases used to estimate genetic effect            ##    ID 1      String Study variant identifier

这样就读取好了，这个函数的读取的结果默认其实是一个S4对象，关于这个对象我下次再详细探索（和TCGAbiolinks下载的SummarisedExperiment对象很像）。我们可以对这个对象进行一些探索，但是这个和今天的主题没有关系，所以就不介绍了，直接下一步。

根据P值过滤

下面是根据P值进行过滤，你可以先使用query_gwas进行过滤，再通过gwasvcf_to_TwoSampleMR转成TwoSampleMR需要的格式，如下所示：

exposure_pval_filter <- query_gwas(vcf = vcf_exposure, pval = 5e-8)exposure_pval_filter <- gwasvcf_to_TwoSampleMR(vcf = exposure_pval_filter)colnames(exposure_pval_filter)##  [1] "chr.exposure"           "pos.exposure"           "other_allele.exposure" ##  [4] "effect_allele.exposure" "beta.exposure"          "se.exposure"           ##  [7] "pval.exposure"          "eaf.exposure"           "samplesize.exposure"   ## [10] "ncase.exposure"         "SNP"                    "ncontrol.exposure"     ## [13] "exposure"               "mr_keep.exposure"       "pval_origin.exposure"  ## [16] "id.exposure"dim(exposure_pval_filter)## [1] 2041   16head(exposure_pval_filter)##   chr.exposure pos.exposure other_allele.exposure effect_allele.exposure## 1            1     47684677                     T                      G## 2            1     47690438                     G                      T## 3            1     49376820                     T                      C## 4            1     49438005                     G                      A## 5            1     49589847                     A                      G## 6            1     49828663                     A                      G##   beta.exposure se.exposure pval.exposure eaf.exposure samplesize.exposure## 1       -0.0168      0.0030  2.181976e-08       0.5333              339152## 2        0.0165      0.0030  3.416017e-08       0.4746              338820## 3        0.0209      0.0035  2.371974e-09       0.2250              339110## 4       -0.0225      0.0031  3.773115e-13       0.6333              338453## 5       -0.0227      0.0031  2.123244e-13       0.5833              318585## 6       -0.0213      0.0032  2.710816e-11       0.6667              339129##   ncase.exposure       SNP ncontrol.exposure exposure mr_keep.exposure## 1             NA  rs977747                NA  ieu-a-2             TRUE## 2             NA rs2984618                NA  ieu-a-2             TRUE## 3             NA rs1343424                NA  ieu-a-2             TRUE## 4             NA rs3127553                NA  ieu-a-2             TRUE## 5             NA  rs657452                NA  ieu-a-2             TRUE## 6             NA rs7531656                NA  ieu-a-2             TRUE##   pval_origin.exposure id.exposure## 1             reported      HSY8sy## 2             reported      HSY8sy## 3             reported      HSY8sy## 4             reported      HSY8sy## 5             reported      HSY8sy## 6             reported      HSY8sy

过滤后只剩下2041个了。

或者也可以先转成TwoSampleMR需要的格式，再使用subset等函数自己过滤下。我选择第1种方法。

去除连锁不平衡

再下面就是去除连锁不平衡的数据。extract_instruments默认的clump_kb=10000，clump_r2=0.001，所以我这里也选择这个条件。

# 这一步也是联网进行的exposure_clumped <- clump_data(exposure_pval_filter,                                clump_kb=10000, clump_r2=0.001,                               pop = "EUR")# 参考人种选欧洲人，默认值# 保存下#save(exposure_clumped,file = "./ieu/exposure_clumped.rdata")dim(exposure_clumped)## [1] 78 16head(exposure_clumped)##     chr.exposure pos.exposure other_allele.exposure effect_allele.exposure## 1              1     47684677                     T                      G## 5              1     49589847                     A                      G## 118            1     72837239                     T                      C## 238            1     78048331                     T                      C## 276            1     96924097                     C                      T## 302            1    110082886                     C                      T##     beta.exposure se.exposure pval.exposure eaf.exposure samplesize.exposure## 1         -0.0168      0.0030  2.181976e-08       0.5333              339152## 5         -0.0227      0.0031  2.123244e-13       0.5833              318585## 118        0.0331      0.0030  1.880182e-28       0.6083              338123## 238        0.0201      0.0031  4.567726e-11       0.4250              339065## 276        0.0221      0.0030  1.433838e-13       0.5750              337797## 302        0.0659      0.0087  5.059411e-14       0.0339              313621##     ncase.exposure        SNP ncontrol.exposure exposure mr_keep.exposure## 1               NA   rs977747                NA  ieu-a-2             TRUE## 5               NA   rs657452                NA  ieu-a-2             TRUE## 118             NA  rs7531118                NA  ieu-a-2             TRUE## 238             NA rs17381664                NA  ieu-a-2             TRUE## 276             NA rs11165643                NA  ieu-a-2             TRUE## 302             NA  rs7550711                NA  ieu-a-2             TRUE##     pval_origin.exposure id.exposure## 1               reported      HSY8sy## 5               reported      HSY8sy## 118             reported      HSY8sy## 238             reported      HSY8sy## 276             reported      HSY8sy## 302             reported      HSY8sy

这个exposure_clumped就是最终的暴露数据了，最后剩下78个SNP，这个数据和在线版基本一样（理论上应该是完全一模一样的，但是不知道为啥这里少了一个，你可以换其他数据试试看），而且通过和各列进行比较，发现共有的数据都是一样的。

但是这一步其实也是联网进行的，所以如果你要本地LD的话，可以使用ieugwasr::ld_clump这个函数实现。

首先需要下载LD的参考数据集到本地，下载地址是：http://fileserve.mrcieu.ac.uk/ld/1kg.v3.tgz，下载后解压（全部解压后是7G+的大小），得到以下数据：

图片

我们这个数据参考人种是欧洲人，所以会用到EUR.bed、EUR.bim、EUR.fam这3个文件。

然后再下载一个plink.exe，下载地址是：https://github.com/explodecomputer/genetics.binaRies/raw/master/binaries/Windows/plink.exe

然后就可以运行以下代码实现本地化的LD了：

test <- ld_clump(  # 只能是这3列，且列名必须是下面这样的  dat = dplyr::tibble(rsid=exposure_pval_filter$SNP,                      pval=exposure_pval_filter$pval.exposure,                      id=exposure_pval_filter$exposure),  clump_kb=10000, clump_r2=0.001,clump_p = 1,  bfile = "ieu/1kg/EUR", # 参考数据集，路径不要写错，只要写前缀EUR即可  plink_bin = "ieu/plink.exe" # plink.exe的存放路径不要写错  )dim(test)## [1] 78  3head(test)## # A tibble: 6 × 3##   rsid           pval id     ##   <chr>         <dbl> <chr>  ## 1 rs977747   2.18e- 8 ieu-a-2## 2 rs657452   2.12e-13 ieu-a-2## 3 rs7531118  1.88e-28 ieu-a-2## 4 rs17381664 4.57e-11 ieu-a-2## 5 rs11165643 1.43e-13 ieu-a-2## 6 rs7550711  5.06e-14 ieu-a-2

结果也是78个，和在线版的是一模一样的，但是结果只有3列，但是只要只选择上面78个SNP就好了：

library(dplyr)exposure_local_ld <- exposure_pval_filter %>%   filter(SNP %in% test$rsid)dim(exposure_local_ld)## [1] 78 16head(exposure_local_ld)##   chr.exposure pos.exposure other_allele.exposure effect_allele.exposure## 1            1     47684677                     T                      G## 2            1     49589847                     A                      G## 3            1     72837239                     T                      C## 4            1     78048331                     T                      C## 5            1     96924097                     C                      T## 6            1    110082886                     C                      T##   beta.exposure se.exposure pval.exposure eaf.exposure samplesize.exposure## 1       -0.0168      0.0030  2.181976e-08       0.5333              339152## 2       -0.0227      0.0031  2.123244e-13       0.5833              318585## 3        0.0331      0.0030  1.880182e-28       0.6083              338123## 4        0.0201      0.0031  4.567726e-11       0.4250              339065## 5        0.0221      0.0030  1.433838e-13       0.5750              337797## 6        0.0659      0.0087  5.059411e-14       0.0339              313621##   ncase.exposure        SNP ncontrol.exposure exposure mr_keep.exposure## 1             NA   rs977747                NA  ieu-a-2             TRUE## 2             NA   rs657452                NA  ieu-a-2             TRUE## 3             NA  rs7531118                NA  ieu-a-2             TRUE## 4             NA rs17381664                NA  ieu-a-2             TRUE## 5             NA rs11165643                NA  ieu-a-2             TRUE## 6             NA  rs7550711                NA  ieu-a-2             TRUE##   pval_origin.exposure id.exposure## 1             reported      HSY8sy## 2             reported      HSY8sy## 3             reported      HSY8sy## 4             reported      HSY8sy## 5             reported      HSY8sy## 6             reported      HSY8sy

读取结局数据

和读取暴露数据完全一样的步骤，就不多做说明了。

# 释放下内存rm(list = ls());gc()##            used  (Mb) gc trigger   (Mb)  max used   (Mb)## Ncells 10747449 574.0   57439848 3067.7  43933768 2346.4## Vcells 19903045 151.9  118845657  906.8 148489984 1132.9vcf_outcome <- readVcf("./ieu/ieu-a-7.vcf.gz")outcome_origin <- gwasvcf_to_TwoSampleMR(vcf = vcf_outcome)rm(vcf_outcome);gc()##             used   (Mb) gc trigger   (Mb)  max used   (Mb)## Ncells  19231085 1027.1  130298919 6958.8 162873648 8698.4## Vcells 167142346 1275.2  788577285 6016.4 870562547 6641.9

extract_outcome_data是从暴露数据的SNP中挑选结局相关的数据，相当于取个交集，所以我们下面也先对暴露数据和结局数据取交集，获取共同的SNP。

load(file = "./ieu/exposure_clumped.rdata")# 取交集common_snp <- merge(exposure_clumped,outcome_origin,by="SNP")$SNPlength(common_snp)## [1] 78

然后对结局数据使用format_data函数修改格式，提取结局数据：

outcome_data <- format_data(dat = outcome_origin,                      type = "outcome", # 类型是outcome                      snps = common_snp, # 共有的SNP                      snp_col = "SNP",                      beta_col = "beta.exposure",                      se_col = "se.exposure",                      eaf_col = "eaf.exposure",                      effect_allele_col = "effect_allele.exposure",                      other_allele_col = "other_allele.exposure",                      pval_col = "pval.exposure",                      samplesize_col = "samplesize.exposure",                      ncase_col = "ncase.exposure",                      ncontrol_col = "ncontrol.exposure",                      id_col = "id.exposure"                      )

这个outcome_data就是我们最终的结局数据了，和在线版获取的也是基本一致哦：

dim(outcome_data)## [1] 78 14head(outcome_data)##   other_allele.outcome effect_allele.outcome beta.outcome se.outcome## 1                    T                     G    -0.013896  0.0096146## 2                    A                     G    -0.007670  0.0093844## 3                    T                     C     0.017675  0.0097331## 4                    T                     C     0.016188  0.0102099## 5                    C                     T     0.005150  0.0092844## 6                    C                     T    -0.048354  0.0305635##   pval.outcome eaf.outcome samplesize.outcome ncase.outcome        SNP## 1   0.14837509    0.550761             184305            NA   rs977747## 2   0.41374714    0.566124             184305            NA   rs657452## 3   0.06937452    0.480633             184305            NA  rs7531118## 4   0.11284907    0.349349             184305            NA rs17381664## 5   0.57910458    0.553284             184305            NA rs11165643## 6   0.11363078    0.027650             184305            NA  rs7550711##   ncontrol.outcome id.outcome outcome mr_keep.outcome pval_origin.outcome## 1               NA     FUNiWE outcome            TRUE            reported## 2               NA     FUNiWE outcome            TRUE            reported## 3               NA     FUNiWE outcome            TRUE            reported## 4               NA     FUNiWE outcome            TRUE            reported## 5               NA     FUNiWE outcome            TRUE            reported## 6               NA     FUNiWE outcome            TRUE            reported

最后修改下暴露因素和结局的名称：

outcome_data$id.outcome <- "CHD"exposure_clumped$id.exposure <- "BMI"

进行MR分析

暴露数据和结局数据都有了，下面就是协调下数据，然后就可以进行MR分析了，后面的步骤就和之前一模一样了！

dat_har <- harmonise_data(exposure_clumped, outcome_data)res <- mr(dat_har)

在往后的分析就不重复了。

联系我们，关注我们

免费QQ交流群1：613637742免费QQ交流群2：608720452公众号消息界面关于作者获取联系方式知乎、CSDN、简书同名账号哔哩哔哩：阿越就是我 本站仅提供存储服务，所有内容均由用户发布，如发现有害或侵权内容，请点击举报。