森林图

做完COX回归分析后，常用ggforest来画森林图，今天介绍forestplot绘制森林图

一、安装包并加载

install.packages("forestplot")
library(forestplot)

二、读取数据

dat = read.csv('ggforest.csv',header = T)
head(dat)
#      id   HR HR.95L HR.95H   pvalue
#1   MISP 1.30   1.16   1.46 6.08e-06
#2  FXYD3 1.37   1.21   1.55 1.15e-06
#3   PSG8 1.30   1.15   1.47 2.17e-05
#4   SOX3 1.23   1.12   1.36 3.67e-05
#5 B3GNT3 1.21   1.10   1.34 9.36e-05
#6   MLPH 1.30   1.16   1.45 5.61e-06

这个数据需要做一些处理：

• 行名取第一列的基因名
• 删除第一列

rownames(dat) = dat$id  #行名
dat = dat[,-1]  #删除第一列
#         HR HR.95L HR.95H   pvalue
#MISP   1.30   1.16   1.46 6.08e-06
#FXYD3  1.37   1.21   1.55 1.15e-06
#PSG8   1.30   1.15   1.47 2.17e-05
#SOX3   1.23   1.12   1.36 3.67e-05
#B3GNT3 1.21   1.10   1.34 9.36e-05
#MLPH   1.30   1.16   1.45 5.61e-06

三、画图

gene=rownames(dat)
hr=sprintf("%.3f",dat$HR)
hrLow=sprintf("%.3f",dat$HR.95L)
hrHigh=sprintf("%.3f",dat$HR.95H)
pVal=ifelse(dat$pvalue<0.001, "<0.001", sprintf("%.3f", dat$pvalue)) 
Hazard.ratio=paste0(hr,"(",hrLow,"-",hrHigh,")")

pdf(file="forestplot.pdf", width = 6, height =4.5)
n=nrow(dat)
nRow=n+1
ylim=c(1,nRow)
layout(matrix(c(1,2),nc=2),width=c(3,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="")
text.cex=0.8 
text(0,n:1,gene,adj=0,cex=text.cex)
text(1.5-0.5*0.2,n:1,pVal,adj=1,cex=text.cex);text(1.5-0.5*0.2,n+1,'pvalue',cex=text.cex,font=2,adj=1)
text(3,n:1,Hazard.ratio,adj=1,cex=text.cex);text(3,n+1,'Hazard ratio',cex=text.cex,font=2,adj=1,)

par(mar=c(4,1,1.5,1),mgp=c(2,0.5,0))
xlim = c(0,max(as.numeric(hrLow),as.numeric(hrHigh)))
plot(1,xlim=xlim,ylim=ylim,type="n",axes=F,ylab="",xaxs="i",xlab="Hazard ratio")
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, 'red', 'blue')
points(as.numeric(hr), n:1, pch = 15, col = boxcolor, cex=1.5)
axis(1)
dev.off()

1.png

部分代码参考：

https://www.jianshu.com/p/a051ce974196