survminer.jpg

survminer的基础用法指南

本文全部代码及示例数据领取：公众号后台回复surv领取。

生存分析是将事件的结果和出现这一结果所经历的时间结合起来分析的一类统计分析方法。不仅考虑事件是否出现，而且还考虑事件出现的时间长短，因此这类方法也被称为事件时间分析（time-to-event analysis）。surminer R包提供了便于生存分析和可视化的功能。

安装与加载

install.packages("survminer")
library("survminer")

ggsurvplot: 绘制生存曲线

#调用survival包
require("survival")
#survival包自带肺癌数据集：lung，查看数据样式
head(lung)

> head(lung)
  inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss
1    3  306      2  74   1       1       90       100     1175      NA
2    3  455      2  68   1       0       90        90     1225      15
3    3 1010      1  56   1       0       90        90       NA      15
4    5  210      2  57   1       1       90        60     1150      11
5    1  883      2  60   1       0      100        90       NA       0
6   12 1022      1  74   1       1       50        80      513       0

#lung数据集：NCCTG晚期肺癌患者的生存率。
inst # 机构代码；
time # 生存天数；
status # 生存状态，1为删失，2为死亡；
age # 年龄；
sex # 性别，1为男性，2为女性；
ph.ecog、ph.karno、pat.karno # 为病人和患者评分，这里用不到；
meal.cal # 进食时消耗的卡路里；
wt.loss # 最近6个月内的体重下降。

绘制无分组的生存曲线

#利用survival包的Surv函数对时间及生存状态进行拟合(Kaplan-Meier法)
fit1 <- survfit(Surv(time, status) ~ 1, data = lung)
#绘制生存曲线
ggsurvplot(fit, palette = "#2E9FDF") #颜色自定

single_sur

绘制两组的生存曲线

#以sex进行分组为例
#基础图形
fit <- survfit(Surv(time, status) ~ sex, data = lung)
ggsurvplot(fit, data = lung)

sex_sur

#更改删失点形状及大小，默认为"+", 可选"|"
ggsurvplot(fit, data = lung, censor.shape="|", censor.size = 4)

sex_sur_2

# Use brewer color palette "Dark2"
ggsurvplot(fit, linetype = "strata", 
           conf.int = TRUE, pval = TRUE,
           palette = "Dark2")

sex_sur_dark

# Use grey palette 
ggsurvplot(fit, linetype = "strata", 
           conf.int = TRUE, pval = TRUE,
           palette = "grey")

sex_sur_grey

#绘制累计风险曲线
ggsurvplot(fit, data = lung, 
           conf.int = TRUE, # 增加置信区间
           fun = "cumhaz") # 绘制累计风险曲线

累计风险曲线

#添加总患者生存曲线
ggsurvplot(fit, # 创建的拟合对象
           data = lung,  # 指定变量数据来源
           conf.int = TRUE, # 显示置信区间
           pval = TRUE, # 添加P值
           surv.median.line = "hv",  # 添加中位生存时间线
           add.all = TRUE) # 添加总患者生存曲线

添加总患者生存曲线

#个性化的生存曲线
ggsurvplot(
  fit,
  data = lung,
  size = 1,                 # 线条大小
  palette =
    c("#E7B800", "#2E9FDF"),# 分组颜色
  conf.int = TRUE,          # 是否添加置信区间
  pval = TRUE,              # 是否添加p值
  risk.table = TRUE,        # 是否添加风险表
  risk.table.col = "strata",# 分线表颜色
  legend.labs =
    c("Male", "Female"),    # 图例标签
  risk.table.height = 0.25, # 生存曲线图下所有生存表的高度，数值0-1之间
  ggtheme = theme_bw()      # 主题
)

Customized survival curves

#更个性化的生存曲线
ggsurvplot(
   fit,                    
   data = lung,            
   risk.table = TRUE,       
   pval = TRUE,           
   conf.int = TRUE,        
   xlim = c(0,500),         # xlim, ylim # 指定x轴和y轴的范围，如xlim = c(0,30), ylim = c(0,1)                       
   xlab = "Time in days",   # xlab, ylab 分别指x轴和y轴标签
   break.time.by = 100,     # 设定坐标轴刻度间距
   ggtheme = theme_light(), 
 risk.table.y.text.col = T, # 分线表y轴文字颜色.
  risk.table.y.text = FALSE # 风险表y轴展示条形图例不展示文字标签
)

Customized survival curves.png

#Uber定制生存曲线
ggsurv <- ggsurvplot(
           fit,                     
           data = lung,            
           risk.table = TRUE,      
           pval = TRUE,             
           conf.int = TRUE,                                 
           palette = c("#E7B800", "#2E9FDF"),
           xlim = c(0,500),      
           xlab = "Time in days",   
           break.time.by = 100,     
           ggtheme = theme_light(), 
          risk.table.y.text.col = T,
          risk.table.height = 0.25, 
          risk.table.y.text = FALSE,
          ncensor.plot = TRUE,      # 绘制删失事件图
          ncensor.plot.height = 0.25, 
          conf.int.style = "step",  # 设置置信区间的类型，有"ribbon"(默认),"step"两种
          surv.median.line = "hv",  # 添加中位生存时间线,可用值有”none”、”hv”、”h”、”v”；其中v绘制垂直线，h绘制水平线。
          legend.labs =
            c("Male", "Female")    
        )
ggsurv

Uber customized survival curves

#引入customize_labels 函数进行更高级的定制
customize_labels <- function (p, font.title = NULL,
                              font.subtitle = NULL, font.caption = NULL,
                              font.x = NULL, font.y = NULL, font.xtickslab = NULL, font.ytickslab = NULL)
{
  original.p <- p
  if(is.ggplot(original.p)) list.plots <- list(original.p)
  else if(is.list(original.p)) list.plots <- original.p
  else stop("Can't handle an object of class ", class (original.p))
  .set_font <- function(font){
    font <- ggpubr:::.parse_font(font)
    ggtext::element_markdown (size = font$size, face = font$face, colour = font$color)
  }
  for(i in 1:length(list.plots)){
    p <- list.plots[[i]]
    if(is.ggplot(p)){
      if (!is.null(font.title)) p <- p + theme(plot.title = .set_font(font.title))
      if (!is.null(font.subtitle)) p <- p + theme(plot.subtitle = .set_font(font.subtitle))
      if (!is.null(font.caption)) p <- p + theme(plot.caption = .set_font(font.caption))
      if (!is.null(font.x)) p <- p + theme(axis.title.x = .set_font(font.x))
      if (!is.null(font.y)) p <- p + theme(axis.title.y = .set_font(font.y))
      if (!is.null(font.xtickslab)) p <- p + theme(axis.text.x = .set_font(font.xtickslab))
      if (!is.null(font.ytickslab)) p <- p + theme(axis.text.y = .set_font(font.ytickslab))
      list.plots[[i]] <- p
    }
  }
  if(is.ggplot(original.p)) list.plots[[1]]
  else list.plots
}

#Uber platinum customized survival curves
# 更改标签
#生存曲线的标签
ggsurv$plot <- ggsurv$plot + labs(
  title    = "Survival curves", #主标题
  subtitle = "Based on Kaplan-Meier estimates", #副标题
  caption  = "created with survminer"  #说明
  )
ggsurv$plot

ggsurv$plot.png

#风险表标签 
ggsurv$table <- ggsurv$table + labs(
  title    = "Note the risk set sizes",
  subtitle = "and remember about censoring.",
  caption  = "source code: website.com"
  )
ggsurv$table

ggsurv$table.png

#删失图标签 
ggsurv$ncensor.plot <- ggsurv$ncensor.plot + labs(
  title    = "Number of censorings",
  subtitle = "over the time.",
  caption  = "source code: website.com"
  )
ggsurv$ncensor.plot

ggsurv$ncensor.plot.png

# 更改字体大小、类型和颜色
ggsurv <- customize_labels(
  ggsurv,
  font.title    = c(16, "bold", "darkblue"),#用长度为3的向量分别指定大小、类型、颜色
  font.subtitle = c(15, "bold.italic", "purple"),
  font.caption  = c(14, "plain", "orange"),
  font.x        = c(14, "bold.italic", "red"),
  font.y        = c(14, "bold.italic", "darkred"),
  font.xtickslab = c(12, "plain", "darkgreen")
)
ggsurv

ggsurv

绘制多组生存曲线图

#绘制多组生存曲线
#使用内置数据集colon
#查看数据格式
head(colon)

> head(colon)
  id study      rx sex age obstruct perfor adhere
1  1     1 Lev+5FU   1  43        0      0      0
2  1     1 Lev+5FU   1  43        0      0      0
3  2     1 Lev+5FU   1  63        0      0      0
4  2     1 Lev+5FU   1  63        0      0      0
5  3     1     Obs   0  71        0      0      1
6  3     1     Obs   0  71        0      0      1
  nodes status differ extent surg node4 time etype
1     5      1      2      3    0     1 1521     2
2     5      1      2      3    0     1  968     1
3     1      0      2      3    0     0 3087     2
4     1      0      2      3    0     0 3087     1
5     7      1      2      2    0     1  963     2
6     7      1      2      2    0     1  542     1

colon数据集：B/C期结肠癌辅助化疗治疗数据
d # 患者编号
study # 所有患者都是1
rx # 治疗方式，有三种：观察、Levamisole、Levamisole+5-FU
sex # 性别，男性为 1，女性为 0
age # 年龄
obstruct # 肿瘤是否阻塞结肠，有为1，无为0
perfor # 结肠是否穿孔，有为1，无为0
adhere # 肿瘤是否粘附附近器官，有为1，无为0
nodes # 检出淋巴结的数目
time # 至发生终点事件或删失的时间
status # 生存状态，1为发生终点事件，0为删失
differ # 肿瘤分化程度 (1=well, 2=moderate, 3=poor)
extent # 局部转移程度（1=粘膜下层，2=肌肉，3=浆膜，4=相邻结构）
surg # 从手术到登记的时间 (0=short, 1=long)
node4 #  超过4个阳性淋巴结，有为1，无为0
etype # 事件类型: 1=复发,2=死亡

#rx分组
fit2 <- survfit( Surv(time, status) ~ rx + adhere,
    data = colon )
ggsurvplot(fit2, pval = TRUE, 
           break.time.by = 800,
           risk.table = TRUE,
           risk.table.height = 0.5
           ) 

多组生存曲线图

参考

1. https://rpkgs.datanovia.com/survminer/
2. http://www.sthda.com/english/wiki/survminer-r-package-survival-data-analysis-and-visualization

往期文章：

ggcorrplot | 简单的相关性热图绘制

ggpubr|让数据可视化更加优雅

ggsci | 让你的配色Nature化