生存分析图表美化
2019中秋节我们的全国巡讲本来应该是去贵州或者兰州，但无奈报名人数太少，临时调到广州。中秋节专场吸引到了几位中山大学临床统计方向博士研究生，其中一位小伙伴把她以前的统计学技能全部R语言化了，会出一系列笔记。
做生存分析时，想模仿做这样的图。
也就是想做K-M曲线，带有risk-table，并且把5-year OS （或者HRs） 加在图的空白区域。
下面的例子试试。
准备和检查数据
# package load
library(magrittr) # 用pipeline符号。
library(survival)
library(survminer) # 用来做km曲线。
library(ggsci) # 当作颜色库来选颜色。生信技能树学到的小妙招。
library(dplyr) # @ data_frame()
library(gridExtra) # @ tableGrob()
library(flextable)# Data preparation and checking
data("colon")
head(colon) %>% flextable() %>% autofit()
|
id |
study |
rx |
sex |
age |
obstruct |
perfor |
adhere |
nodes |
status |
differ |
extent |
surg |
node4 |
time |
etype |
|
1 |
1 |
Lev+5FU |
1 |
43 |
0 |
0 |
0 |
5 |
1 |
2 |
3 |
0 |
1 |
1,521 |
2 |
|
1 |
1 |
Lev+5FU |
1 |
43 |
0 |
0 |
0 |
5 |
1 |
2 |
3 |
0 |
1 |
968 |
1 |
|
2 |
1 |
Lev+5FU |
1 |
63 |
0 |
0 |
0 |
1 |
0 |
2 |
3 |
0 |
0 |
3,087 |
2 |
|
2 |
1 |
Lev+5FU |
1 |
63 |
0 |
0 |
0 |
1 |
0 |
2 |
3 |
0 |
0 |
3,087 |
1 |
|
3 |
1 |
Obs |
0 |
71 |
0 |
0 |
1 |
7 |
1 |
2 |
2 |
0 |
1 |
963 |
2 |
|
3 |
1 |
Obs |
0 |
71 |
0 |
0 |
1 |
7 |
1 |
2 |
2 |
0 |
1 |
542 |
1 |
## outcome
summary(colon$time)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 8 566 1855 1538 2331 3329colon$time %>% density %>% plot
table(colon$status)##
## 0 1
## 938 920## exposure
table(colon$rx)##
## Obs Lev Lev+5FU
## 630 620 608
第一步：画带risk table的K-M曲线。
p=ggsurvplot(survfit(Surv(time,status==1)~rx,data = colon),
conf.int = T, # 标注95% CI。
surv.scale = "percent",
palette = pal_nejm("default")(3), # 选择ggsci包里边的颜色。
ylab = "Survival proportion",
font.x = c(18),
font.y = c(18),
censor.shape="|",
censor.size = 2,
legend.title = "",
legend.labs=levels(colon$rx), # legend 需要改下，不然多了rx字样。
font.legend = c(12, "plain"), # plot legend font size, color.
risk.table = T, # 显示risk table。
risk.table.title = "No. at risk",
risk.table.fontsize = 3, # table number font size
risk.table.height = 0.3,
surv.plot.height = 0.6,
title = "",
pval = T, # log-rank test得到的p value。
pval.size = 4,
pval.coord = c(2, .1) # p value在图中的位置。
) +
labs(x='') # 去掉km曲线的x 坐标的label，因为与risk table的x label是重复的。
p$table <- p$table +
theme(plot.title = element_text(size = 10))+labs(x='Time since randomization') # 调整risk table的字体大小.ggsurvplot 得到的table与plot均可用ggplot2里边的参数来调整。非常方便。
p
第二步：做5-year OS 表格。
# 得到5-year OS 数据list
a=summary(survfit(Surv(time,status==1) ~ rx, data = colon), times = 5*365)
# 提取我们要的 5-year OS和95% CI,并调整好格式。
mytable=data_frame(Group=a[["strata"]],surve5=a[["surv"]],lower=a[["lower"]],upper=a[["upper"]]) %>% mutate(Group=levels(colon$rx),surve5=paste0(round(100*surve5,2),'%'),lower=paste0(round(lower*100,2),'%'),upper=paste0(round(upper*100,2),'%'),OS_5y=paste0(surve5,' (',lower,'-',upper,')')) %>% .[,c(1,5)]
mytable %>%flextable() %>% autofit()
|
Group |
OS_5y |
|
Obs |
48.75% (44.98%-52.84%) |
|
Lev |
49.76% (45.95%-53.87%) |
|
Lev+5FU |
62.44% (58.68%-66.43%) |
第三步：将K-M 曲线和5-year OS表格整到一起。
p$plot=p$plot+annotation_custom(tableGrob(mytable,rows = NULL),xmin=50, xmax=2000, ymin=0.15, ymax=0.5)
p
Done!
尽力了，还是不是很漂亮，图与表的位置和尺寸再根据需要图片大小来微调吧。
5-year OS 表格的背景不会去掉，若有知道的小伙伴，欢迎留言交流。