小程序开发外包 孟德尔当场化：下载IEU OPEN GWAS数据腹地处治
发布日期：2024-08-08 15:32 点击次数：125
💡专注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：
app图片
先找袒露数据的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.9extract_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、简书同名账号哔哩哔哩：阿越便是我 本站仅提供存储工作，总计本色均由用户发布，如发现存害或侵权本色，请点击举报。