【孟德尔随机化】森林图怎么画？

生信菜鸟团

发布于 2023-09-19 15:46:30

2.1K10

代码可运行

文章被收录于专栏：生信菜鸟团生信菜鸟团

运行总次数：0

代码可运行

本来想复现上周的那篇文献导读，但是找了很久都没找到他的原始数据，就……很绝望。。。所以给大家补一个森林图的画图方式~

我们常在MR文献中看到的森林图咋画呢？

今天简单来画一下：

假如原始表格长这样

View(outresult)

out_multi <- outresult
# out_multi[,2:ncol(out_multi)] <- as.numeric(unlist(out_multi[,2:ncol(out_multi)]))
hz <- paste(round(out_multi$HR,3),
            "(",round(out_multi$L95CI,3),
            "-",round(out_multi$H95CI,3),")",sep = "")


tabletext <- cbind(c(NA,"CELL",out_multi$CELL),
                   c(NA,"Coefficient",round(out_multi$coef,3)),
                   c(NA,"P value",ifelse(out_multi$pvalue<0.001,"P < 0.001",round(out_multi$pvalue,3))),
                   c(NA,"Hazard Ratio(95% CI)",hz))

png("../rawdata/based_unicox_Hazard_Ratio.png", width=1000, height=600)
pforest <- forestplot(labeltext=tabletext, 
           graph.pos=3,  #为Pvalue箱线图所在的位置
           col=fpColors(box="#D55E00", lines="#CC79A7", zero = "gray50"),
           mean=c(NA,NA,out_multi$HR),
           lower=c(NA,NA,out_multi$L95CI), #95%置信区间下限
           upper=c(NA,NA,out_multi$H95CI), #95%置信区间上限
           xlab="Unicox−analytic Hazard Ratio of immune cell population",
           boxsize=0.3,lwd.ci=2,   #箱子大小，线的宽度
           ci.vertices.height = 0.08,ci.vertices=TRUE, #置信区间用线宽、高、型
           zero=1,lwd.zero=1,      #zero线宽 基准线的位置
           colgap=unit(5,"mm"),    #列间隙
           xticks = c(0.5, 1,4), #横坐标刻度
           lwd.xaxis=1,            #X轴线宽
           lineheight = unit(1,"cm"), #固定行高
           graphwidth = unit(.3,"npc"), #图在表中的宽度比例
           cex=0.9, fn.ci_norm = fpDrawCircleCI, #误差条显示方式
           hrzl_lines=list("2" = gpar(lwd=2, col="black"),
                           "3" = gpar(lwd=2, col="black"), #第三行顶部加黑线，引号内数字标记行位置
                           "16" = gpar(lwd=2, col="black")),#最后一行底部加黑线,"29"中数字为nrow(tabletext)+1
           mar=unit(rep(0.1, times = 4), "mm"),#图形页边距
           # #fpTxtGp函数中的cex参数设置各个组件的大小
           txt_gp=fpTxtGp(label=gpar(cex=1),
                          ticks=gpar(cex=1.5),
                          xlab=gpar(cex = 1.25),
                          title=gpar(cex = 1.2))
           )
pforest
dev.copy(png, "../Hazard_Ratio.png")
dev.off()