master

分支 (1)

管理

管理

master

lover
/
森林图美化.R

#1.载入包
library(tableone)
library(survival)
library(readxl)
##画森林图的包
library(forestplot)
library(stringr)
#2.清理工作环境
rm(list = ls())
#3.读入数据
aa <- read_excel("C:/Users/99405/Desktop/肖方2024.1.10.xlsx")
#4.查看数据前6行
head(aa)
#5.查看数据数据性质
str(aa)
#6.查看结局，1=复发，0=未复发


#设置因子顺序
aa$铁蛋白<-factor(aa$铁蛋白,levels=c("≤1000ug/L",">1000ug/L"))
aa$ECOG <- factor(aa$ECOG)
#  一、构建森林图的基本表
#一-1.cox多因素回归分析
mul_cox<-coxph(Surv(Survival_time,是否死亡==1)~年龄+性别+MNC+CD34+疾病类型+GROUP+ECOG+血型匹配 + 预处理方案 + 铁蛋白,data=aa)
#一-2 multi1：提取：变量+HR+95%CI+95%CI
mul_cox1 <- summary(mul_cox)
colnames(mul_cox1$conf.int)
multi1<-as.data.frame(round(mul_cox1$conf.int[, c(1, 3, 4)], 2))
#一-3、multi2：提取：HR(95%CI)和P
multi2<-ShowRegTable(mul_cox,
                     exp=TRUE,
                     digits=2,
                     pDigits =3,
                     printToggle = TRUE,
                     quote=FALSE,
                     ciFun=confint)
#一-4.将两次提取结果合并成表；取名result
result <-cbind(multi1,multi2);result
#一-5.行名转为表格第一列，并给予命名"Characteristics"
result<-tibble::rownames_to_column(result, var = "Characteristics");result


# 二、森林图的基本图形fig1
fig1<- forestplot(result[,c(1,5,6)], #告诉函数，合成的表格result的第1，5，6列还是显示数字
                  mean=result[,2],   #告诉函数，表格第2列为HR，它要变成森林图的小方块
                  lower=result[,3],  #告诉函数表格第3列为95%CI，
                  upper=result[,4],  #表格第5列为95%CI，它俩要化作线段，穿过方块
                  zero=1,            #告诉函数，零线或参考线为HR=1即x轴的垂直线
                  boxsize=0.3,       #设置小黑块的大小
                  graph.pos=2)       #森林图应插在图形第2列
fig1

# 三、森林图图形修饰
# 三-1、显示所有变量
#1.删除部分变量名，只保留亚变量即："ER_Positive"变为"Positive"
result$Characteristics<-str_remove(result$Characteristics,"性别|疾病类型|GROUP|ECOG|血型匹配|预处理方案|铁蛋白")

#2. 给参考变量插入空行
#2-1.这步代码不用改
ins <- function(x) {c(x, rep(NA, ncol(result)-1))}
##2-2：插入空行，形成一个新表
for(i in 5:6) {result[, i] = as.character(result[, i])}
result<-rbind(c("Characteristics", NA, NA, NA, "HR(95%CI)","p"),
              result[1, ],
              ins("Gender"),
              ins("female"),
              result[2:4, ],
              ins("疾病类型"),
              ins("SAA"),
              result[5, ],
              ins("GROUP"),
              ins("Firstline Haplo-HSCT"),
              result[6, ],
              ins("ECOG"),
              ins("0"),
              result[7, ],
              ins("血型匹配"),
              ins("次"),
              result[8:10, ],
              ins("预处理方案"),
              ins("BU/CY+ATG"),
              result[11, ],
              ins("铁蛋白"),
              ins("≤1000ug/L"),
              result[12:nrow(result), ],
              c(NA, NA, NA, NA, NA,NA)
)
for(i in 2:4) {result[, i] = as.numeric(result[, i])}

fig2 <- forestplot(result[,c(1,5,6)],
                   mean=result[,2],
                   lower=result[,3],
                   upper=result[,4],
                   zero=1,
                   boxsize=0.5,
                   graph.pos=2)
fig2
# 三-2、显示所有亚组的患者数
#1-1
myVars <- c("年龄","性别","MNC","CD34","疾病类型","GROUP","ECOG","血型匹配","预处理方案","铁蛋白")
#1-2
catVars <-  c("性别","疾病类型","GROUP","ECOG","血型匹配","预处理方案","铁蛋白")
#1-3
table1<- print(CreateTableOne(vars=myVars,
                              data = aa,
                              factorVars = catVars),
               showAllLevels=TRUE)

#2. 在基线表table1里插入空行，使它的行数和变量跟result一致
N<-rbind(c(NA,NA),
         table1[2, ],
         c(NA, NA),
         table1[3:6,],
         c(NA,NA),
         table1[7:8,],
         c(NA,NA),
         table1[9:10,],
         c(NA,NA),
         table1[11:12,],
         c(NA,NA),
         table1[13:16,],
         c(NA, NA),
         table1[17:18,],
         c(NA,NA),
         table1[19:20,],
         c(NA, NA))
N<-N[,-1]
N<-data.frame(N)
#3.把N表和result表合在一起
result1<-cbind(result,N)
#调顺序。变为:变量-N-HR......顺序
result1<-result1[,c(1,7,2:6)]

#4.优化第一行。第一行行名中加入"Number(%)"
for(i in 2:7) {result1[, i] = as.character(result1[, i])}
result1<-rbind(c("Characteristics","Number (%)",NA,NA,NA,"HR (95%CI)","P.value"),
               result1[2:nrow(result1),])
for(i in 3:5) {result1[, i] = as.numeric(result1[, i])}


#画图fig-3，注：因为多了一列，所以要注意改代码数字
fig3<-forestplot(result1[,c(1,2,6,7)],
                 mean=result1[,3],
                 lower=result1[,4],
                 upper=result1[,5],
                 zero=1,
                 boxsize=0.4,
                 graph.pos=3)
fig3
#三-3、图形细节优化

#Characteristics下面加一个空行,美化图形
result1<-rbind(result1[1,],c(NA,NA,NA,NA,NA,NA,NA),
               result1[2:nrow(result1),])

fig3_1<-forestplot(result1[,c(1,2,6,7)],
                   mean=result1[,3],
                   lower=result1[,4],
                   upper=result1[,5],
                   zero=1,
                   boxsize=0.3,
                   graph.pos= "right" ,
                   hrzl_lines=list("1" = gpar(lty=1,lwd=2), #表格上框线
                                   "3" = gpar(lty=2), #第二条虚线
                                   "30"= gpar(lwd=2,lty=1,columns=c(1:4)) ),#表格下框线，column设置其覆盖的列
                   graphwidth = unit(0.3,"npc"), #森林图占整张图的比例 这里是0.3即30%
                   xlab="HR", #森林图坐标轴标签
                   xticks=c(0,1,10,20) , #森林图坐标轴上的刻度
                   is.summary=c(T,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F),#设置有无加粗行，这样即只有第一行表头字体加粗
                   txt_gp=fpTxtGp(
                     label=gpar(cex=1), #表格字符的文字大小
                     ticks=gpar(cex=1), #森林图刻度的文字大小
                     xlab=gpar(cex=1.5), #森林图标题的文字大小
                     title=gpar(cex=3)),  #图标题文字大小
                   lwd.zero=1, # 森林图参考线的线宽
                   lwd.ci=1.5, # 森林图置信区间的线宽
                   lwd.xaxis=2, #森林图x轴线宽
                   lty.ci=1.5,  #森林图置信区间线型
                   ci.vertices =T, #森林图置信区间左右端点是否显示竖直线条
                   ci.vertices.height=0.4, #森林图置信区间左右端点线条长度
                   clip=c(0,20), #森林图HR坐标轴显示的范围
                   ineheight=unit(8, 'mm'), # 每个研究条目的高度
                   line.margin=unit(8, 'mm'),# 线的边距
                   colgap=unit(8, 'mm'),# 表格列之间的间隔。
                   fn.ci_norm="fpDrawDiamondCI",  #绘制置信区间的函数应该不用改
                   title="多因素Cox回归森林图",
                   col=fpColors(box ='#021eaa',
                                lines ='#021eaa',
                                zero = "black")) # 框、线和参考线的颜色

fig3_1