森林图的绘制

res=read.table('PUD_PTSD_MR.txt',header = T,sep = '\t')
res= generate_odds_ratios(res)
write.table(res,file='PUD_PTSD.txt',sep = '\t',row.names = F,quote = F)
library(forestplot)
inputFile="PUD_PTSD.txt"
outFile="PUD_PTSD_forest.pdf"
rt=read.table(inputFile,header=T,sep="\t",row.names=5,check.names=F)
gene=rownames(rt)#读取基因列
hr=sprintf("%.3f",rt$"or")#获取HR列取小数点后3位
hrLow=sprintf("%.3f",rt$"or_lci95")#获取95%置信区间取HR.95L列小数点后3位
hrHigh=sprintf("%.3f",rt$"or_uci95")#获取95%置信区间HR.95H列小数点后3位
pVal=ifelse(rt$pval<0.001, "<0.001", sprintf("%.3f", rt$pval)) #获取P值
Hazard.ratio=paste0(hr,"(",hrLow,"-",hrHigh,")")#合并为一个数据集
#输出图形
pdf(file=outFile, width = 8, height =4.5)
n=nrow(rt)
nRow=n+1
ylim=c(1,nRow)
layout(matrix(c(1,2),nc=2),width=c(4,2))
#森林图左边的基因信息
xlim = c(0,3)
par(mar=c(4,2,1.5,1.5))
plot(1,xlim=xlim,ylim=ylim,type="n",axes=F,xlab="",ylab="")#axes=F表示禁止生成坐标轴
text.cex=0.8 #放大0.8倍
text(0,n:1,gene,adj=0,cex=text.cex,family='serif')#显示基因列信息
text(1.5-0.5*0.2,n:1,pVal,adj=1,cex=text.cex,family='serif');text(1.5-0.5*0.2,n+1,'pvalue',cex=text.cex,font=2,adj=1,family='serif')#显示pvalue列信息
text(3,n:1,Hazard.ratio,adj=1,cex=text.cex,family='serif');text(3,n+1,'OR(95%)',cex=text.cex,font=2,adj=1,family='serif')#显示计算出的Hazard ratio列信息
par(mar=c(4,1,1.5,1),mgp=c(2,0.5,0))
xlim = c(0.75,max(as.numeric(hrLow),as.numeric(hrHigh)))
plot(1,xlim=xlim,ylim=ylim,type="n",axes=F,ylab="",xaxs="i",xlab="Odd Ratio(95%)")#设置x轴的标题
arrows(as.numeric(hrLow),n:1,as.numeric(hrHigh),n:1,angle=90,code=3,length=0.03,col="darkblue",lwd=2.5)
abline(v=1,col="black",lty=2,lwd=2) #添加中线，设置中线的位置，颜色，类型，宽度
boxcolor = ifelse(as.numeric(hr) > 1, 'green', 'black')#设置中线的取值
points(as.numeric(hr), n:1, pch = 15, col = boxcolor, cex=1.5)#pch表示点的样式，设置点的大小，颜色
axis(1)
dev.off()

dat <- subset(dat, dat$palindromic !="TRUE")
data=select(dat,c('beta.exposure','se.exposure','pval.exposure','beta.outcome','se.outcome','pval.outcome'))
mr_presso(BetaOutcome = "beta.outcome", BetaExposure = "beta.exposure", SdOutcome = 'se.outcome', SdExposure = "se.exposure", OUTLIERtest = TRUE, DISTORTIONtest = TRUE, data = data, NbDistribution = 1000,  SignifThreshold = 0.05)

set.seed(100)
varlist <- with(X, sample(snp, size=1000000, replace=FALSE))
params <- est_cause_params(X, varlist)
params$rho
r2_thresh = 0.01
pval_thresh = 1e-3
X_clump <- X %>%
  rename(rsid = snp,
         pval = p1) %>%
  ieugwasr::ld_clump(dat = .,
                     clump_r2 = r2_thresh,
                     clump_p = pval_thresh,
                     plink_bin = genetics.binaRies::get_plink_binary(), 
                     pop = "EUR",
                     bfile = './EUR/EUR')
top_vars <- X_clump$rsid
res <- cause(X=X, variants = top_vars, param_ests = params)
res$elpd
summary(res, ci_size=0.95)$tab
summary(res)$p
summary(res)$tab
plot(res)
elpd_table <- recompute_elpd_table(res)
pnorm(res$elpd$z[3])
write.table(res,file = 'res')