热点资讯
你的位置：制作物联网软件物联网软件的价格 > 物联网软件开发价格 > 制作物联网软件物联网软件的价格 孟德尔立时化：下载IEU OPEN GWAS数据腹地处治
制作物联网软件物联网软件的价格 孟德尔立时化：下载IEU OPEN GWAS数据腹地处治
发布日期：2024-11-06 09:16 点击次数：56
💡专注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，物联网app开发这个数据和在线版基本一样（表面上应该是统融合模一样的，但是不知谈为啥这里少了一个，你可以换其他数据碰红运），况且通过和各列进行比拟，发现共有的数据皆是一样的。
但是这一步其实亦然联网进行的，是以要是你要腹地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分析了，背面的流弊就和之前一模一样了！
号码频次：在第182期历史同期开奖中，号码0-9出现频次分别为：号码7出现3个，号码0、6出现4个，号码2、5出现6个，号码1出现7个，号码3、9出现8个，号码4、8出现10个，今年同期绝杀一码7，独胆看好3。
dat_har <- harmonise_data(exposure_clumped, outcome_data)res <- mr(dat_har)
在往后的分析就不肖似了。
有计划咱们制作物联网软件物联网软件的价格，姿色咱们
免费QQ通常群1：613637742免费QQ通常群2：608720452公众号音书界濒临于作家获取有计划方式知乎、CSDN、简书同名账号哔哩哔哩：阿越即是我 本站仅提供存储干事，系数履行均由用户发布，如发现存害或侵权履行，请点击举报。