孟德尔立时化:下载IEU OPEN GWAS数据土产货不停
💡专注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、简书同名账号哔哩哔哩:阿越即是我 本站仅提供存储就业,悉数骨子均由用户发布,如发现存害或侵权骨子,请点击举报。