赞
踩
我要一键出图
- #install.packages("devtools")
- #devtools::install_github("mrcieu/gwasglue", force = TRUE)
-
- #if (!require("BiocManager", quietly = TRUE))
- # install.packages("BiocManager")
- #BiocManager::install("VariantAnnotation")
-
- #install.packages("remotes")
- #remotes::install_github("MRCIEU/TwoSampleMR")
-
- #install.packages("dplyr")
- #install.packages("tidyr")
- #install.packages("CMplot")
-
- road<-getwd()
- setwd(road) #设置工作目录
-
- #引用包
- library(VariantAnnotation)
- library(gwasglue)
- library(dplyr)
- library(tidyr)
- library(CMplot)
- library(TwoSampleMR)
-
- exposureName="BMI" #图形中展示暴露数据的名称
- outcomeName="Coronary heart disease" #图形中展示结局数据的名称
-
- inputFile="bbj-a-57.vcf.gz" #暴露数据输入文件(根据下载暴露数据的文件名称进行修改)
- outcomeFile="ukb-a-561.vcf.gz" #结局数据输入文件(改成已下载的结局数据文件)
-
-
- #读取输入文件, 并对输入文件进行格式转换
- vcfRT1 <- readVcf(inputFile)
- exposuredata=gwasvcf_to_TwoSampleMR(vcf=vcfRT1, type="exposure")
- #读取结局数据的vcf文件,并对数据进行格式转换
- vcfRT2 <- readVcf(outcomeFile)
- outcomeData=gwasvcf_to_TwoSampleMR(vcf=vcfRT2, type="outcome")
- #根据pvalue<5e-08对结果进行过滤
- outTab<-subset(exposuredata, pval.exposure<5e-08)
- write.csv(outTab, file="exposure.pvalue.csv", row.names=F)
-
- #准备绘制暴露变量的曼哈顿图的数据
- exposuredata=exposuredata[,c("SNP", "chr.exposure", "pos.exposure", "pval.exposure")]
- colnames(exposuredata)=c("SNP","CHR","BP","pvalue")
-
- #绘制线性的曼哈顿图
- CMplot(exposuredata, plot.type="m",
- LOG10=TRUE, threshold=5e-08, threshold.lwd=3, threshold.lty=1, signal.cex=0.2,
- chr.den.col=NULL, cex=0.2, bin.size=1e5, ylim=c(0,50),
- file="pdf", file.output=TRUE, width=15, height=9, verbose=TRUE)
-
- #绘制圈图
- CMplot(exposuredata, plot.type="c",
- LOG10=TRUE, threshold=5e-08, threshold.lwd=3, threshold.lty=1, signal.cex=0.2,
- chr.den.col=NULL, cex=0.2, bin.size=1e5, ylim=c(0,100),
- file="pdf", file.output=TRUE, width=7, height=7, verbose=TRUE)
-
-
-
- exposureFile="exposure.pvalue.csv" #输入文件
- #setwd(road) #设置孟德尔结果数据保存的工作目录
-
- #读取输入文件
- exposure_dat<-read_exposure_data(filename=exposureFile,
- sep = ",",
- snp_col = "SNP",
- beta_col = "beta.exposure",
- se_col = "se.exposure",
- effect_allele_col = "effect_allele.exposure",
- other_allele_col = "other_allele.exposure",
- eaf_col = "eaf.exposure",
- samplesize_col = "samplesize.exposure",
- clump = F)
-
- #去除连锁不平衡的SNP,一般标准为kb=10000,r2=0.001
- exposure_dat_clumped <- clump_data(exposure_dat, clump_kb=10000, clump_r2=0.001)
- write.csv(exposure_dat_clumped, file="exposure.LD.csv", row.names=F)
-
- inputFile="exposure.LD.csv" #输入文件
- #setwd(road) #设置工作目录
-
- #读取输入文件
- dat<-read.csv(inputFile, header=T, sep=",", check.names=F)
-
- #计算F检验值
- N=dat[1,"samplesize.exposure"] #获取样品的数目
- dat=transform(dat,R2=2*((beta.exposure)^2)*eaf.exposure*(1-eaf.exposure)) #计算R2
- dat=transform(dat,F=(N-2)*R2/(1-R2)) #计算F检验值
-
- #根据F值>10进行过滤, 删除弱工具变量
- Ffilter=10 #F值过滤条件,
- outTab=dat[dat$F>Ffilter,]
- write.csv(dat, "exposure.F.csv", row.names=F)
-
- library(MendelianRandomization) #引用包
- inputFile="exposure.F.csv" #输入文件
- #setwd(road) #可设置专属文件的独有工作目录
-
- #读取已经去除连锁不平衡、低F值的暴露因素文件
- dat=read.csv(inputFile, header=T, sep=",", check.names=F)
-
- #对SNP分组
- snpId=dat$SNP
- y=seq_along(snpId)
- chunks <- split(snpId, ceiling(y/100))
-
- #对分组进行循环,每次得到一个分组
- outTab=data.frame()
- for(i in names(chunks)){
- #混杂因素分析
- confounder=phenoscanner(
- snpquery = chunks[[i]],
- catalogue = "GWAS",
- pvalue = 1e-05,
- proxies = "None",
- r2 = 0.8,
- build = 37)
- outTab=rbind(outTab, confounder$results)
- }
- #输出SNP相关性状的表格
- write.csv(outTab, "confounder.result.csv", row.names=F)
-
- #输出去除混杂因素的结果
- delSnp=c("rs13078960", "rs2030323") #混杂SNP的名称(需修改)
- dat=dat[!dat$SNP %in% delSnp,]
- write.csv(dat, "exposure.confounder.csv", row.names=F)
-
- exposureFile="exposure.confounder.csv" #输入经各种过滤的暴露数据文件
-
- #读取暴露数据
- exposure_dat<-read_exposure_data(filename=exposureFile,
- sep = ",",
- snp_col = "SNP",
- beta_col = "beta.exposure",
- se_col = "se.exposure",
- effect_allele_col = "effect_allele.exposure",
- other_allele_col = "other_allele.exposure",
- eaf_col = "eaf.exposure",
- clump = F)
-
-
-
- #从结局数据中提取工具变量
- outcomeTab<-merge(exposure_dat, outcomeData, by.x="SNP", by.y="SNP")
- write.csv(outcomeTab[,-(2:ncol(exposure_dat))], file="outcome.csv")
-
- #读取整理好的结局数据
- outcome_dat<-read_outcome_data(snps=exposure_dat$SNP,
- filename="outcome.csv", sep = ",",
- snp_col = "SNP",
- beta_col = "beta.outcome",
- se_col = "se.outcome",
- effect_allele_col = "effect_allele.outcome",
- other_allele_col = "other_allele.outcome",
- pval_col = "pval.outcome",
- eaf_col = "eaf.outcome")
-
- #将暴露数据和结局数据合并
- exposure_dat$exposure=exposureName
- outcome_dat$outcome=outcomeName
- dat2<-harmonise_data(exposure_dat=exposure_dat,
- outcome_dat=outcome_dat)
-
- #输出用于孟德尔随机化的工具变量
- outTab=dat2[dat2$mr_keep=="TRUE",]
- write.csv(outTab, file="table.SNP.csv", row.names=F)
-
- #MR-PRESSO异常值检测(偏倚的SNP)
- presso=run_mr_presso(dat2)
- write.csv(presso[[1]]$`MR-PRESSO results`$`Outlier Test`, file="table.MR-PRESSO.csv")
-
- #孟德尔随机化分析
- mrResult=mr(dat)
- #选择孟德尔随机化的方法
- #mr_method_list()$obj
- #mrResult=mr(dat, method_list=c("mr_ivw", "mr_egger_regression", "mr_weighted_median", "mr_simple_mode", "mr_weighted_mode"))
- #对结果进行OR值的计算
- mrTab=generate_odds_ratios(mrResult)
- write.csv(mrTab, file="table.MRresult.csv", row.names=F)
-
- #异质性分析
- heterTab=mr_heterogeneity(dat)
- write.csv(heterTab, file="table.heterogeneity.csv", row.names=F)
-
- #多效性检验
- pleioTab=mr_pleiotropy_test(dat)
- write.csv(pleioTab, file="table.pleiotropy.csv", row.names=F)
-
- #绘制散点图
- pdf(file="pic.scatter_plot.pdf", width=7.5, height=7)
- mr_scatter_plot(mrResult, dat)
- dev.off()
-
- #森林图
- res_single=mr_singlesnp(dat) #得到每个工具变量对结局的影响
- pdf(file="pic.forest.pdf", width=7, height=6.5)
- mr_forest_plot(res_single)
- dev.off()
-
- #漏斗图
- pdf(file="pic.funnel_plot.pdf", width=7, height=6.5)
- mr_funnel_plot(singlesnp_results = res_single)
- dev.off()
-
- #留一法敏感性分析
- pdf(file="pic.leaveoneout.pdf", width=7, height=6.5)
- mr_leaveoneout_plot(leaveoneout_results = mr_leaveoneout(dat))
- dev.off()
-
- pFilter=1 #pvalue过滤条件
- #setwd(road) #设置工作目录
-
- ############定义森林图函数############
- bioForest=function(inputFile=null, forestFile=null, forestCol=null){
- #读取输入文件
- rt=read.csv(inputFile, header=T, sep=",", check.names=F)
- row.names(rt)=rt$method
- rt=rt[rt$pval<pFilter,]
- method <- rownames(rt)
- or <- sprintf("%.3f",rt$"or")
- orLow <- sprintf("%.3f",rt$"or_lci95")
- orHigh <- sprintf("%.3f",rt$"or_uci95")
- OR <- paste0(or,"(",orLow,"-",orHigh,")")
- pVal <- ifelse(rt$pval<0.001, "<0.001", sprintf("%.3f", rt$pval))
-
- #输出图形
- pdf(file=forestFile, width=7, height=4.6)
- n <- nrow(rt)
- nRow <- n+1
- ylim <- c(1,nRow)
- layout(matrix(c(1,2),nc=2),width=c(3.5,2))
-
- #绘制森林图左边的信息
- xlim = c(0,3)
- par(mar=c(4,2.5,2,1))
- plot(1,xlim=xlim,ylim=ylim,type="n",axes=F,xlab="",ylab="")
- text.cex=0.8
- text(0,n:1,method,adj=0,cex=text.cex)
- text(1.9,n:1,pVal,adj=1,cex=text.cex);text(1.9,n+1,'pvalue',cex=1,font=2,adj=1)
- text(3.1,n:1,OR,adj=1,cex=text.cex);text(2.7,n+1,'OR',cex=1,font=2,adj=1)
-
- #绘制右边的森林图
- par(mar=c(4,1,2,1),mgp=c(2,0.5,0))
- xlim = c(min(as.numeric(orLow)*0.975,as.numeric(orHigh)*0.975,0.9),max(as.numeric(orLow),as.numeric(orHigh))*1.025)
- plot(1,xlim=xlim,ylim=ylim,type="n",axes=F,ylab="",xaxs="i",xlab="OR")
- arrows(as.numeric(orLow),n:1,as.numeric(orHigh),n:1,angle=90,code=3,length=0.05,col="darkblue",lwd=3)
- abline(v=1, col="black", lty=2, lwd=2)
- boxcolor = ifelse(as.numeric(or)>1, forestCol, forestCol)
- points(as.numeric(or), n:1, pch = 15, col = boxcolor, cex=2)
- axis(1)
- dev.off()
- }
-
- #调用函数,绘制森林图
- bioForest(inputFile="table.MRresult.csv", forestFile="forest.pdf", forestCol="red")
`,_+=1})),x=t}let U=l,A=l;l<0?(a=0,i<=o&&(a=s.div(s.sub(i,d),2))):0==l?i<=o&&(a=s.div(s.sub(i,d),2)):(t=s.add(4,36),t=s.add(s.add(d,t),s.add(s.mul(_,38),36)),l>S.sub(i,T)&&(A=S.sub(i,T)));let O="",L=t?s.runtime.getURL("img/video-default.png"):"https://res.stayfork.app/scripts/BB8CD00276006365956C32A6556696AD/icon.png",D='
Copyright © 2003-2013 www.wpsshop.cn 版权所有,并保留所有权利。