### 生存曲线作图全代码教程 ### (•̤̀ᵕ•̤́๑)ᵒᵏᵎᵎᵎᵎ ###################################################################################### rm(list = ls()) ### 安装CRAN上的R包："survival" if(length(getOption("CRAN"))==0) options(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/") # 设置镜像 if(!require("survival")) install.packages("survival",update = F,ask = F) # 正式安装 ### 安装Bioconductor上的R包："survminer" if(length(getOption("CRAN"))==0) options(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/") # 设置镜像 if(!require("BiocManager")) install.packages("BiocManager",update = F,ask = F) # 安装Bioconductor包 if(length(getOption("BioC_mirror"))==0) options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/") # 再设置BioC镜像 if(!require("survminer")) BiocManager::install("survminer",update = F,ask = F) # 正式安装 ###################################################################################### ### 第一步：准备数据集 ### 调用survival包中的"lung"数据集 data(lung) ### 查看lung数据集长啥样 head(lung) ### 长这样： # inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss # 3 306 2 74 1 1 90 100 1175 NA # 3 455 2 68 1 0 90 90 1225 15 # 3 1010 1 56 1 0 90 90 NA 15 # 5 210 2 57 1 1 90 60 1150 11 # 1 883 2 60 1 0 100 90 NA 0 # 12 1022 1 74 1 1 50 80 513 0 # time：时间；status：病人状态(2为Dead，1为Censored)；inst：组织编码；其它：各种变量 ### 在我们准备的数据集中，必须要有三列信息：时间(time)、事件(event)和预测变量(variable) # 在我们调用的lung数据集中，时间信息是"time"，事件信息是"status"，预测变量我们这里选择"sex" # 接下来评估性别对生存的影响 # 这里的"sex"是分类变量，如果预测变量为连续变量怎么办呢？比如基因表达量、预后评分或ALT含量等 # 一般我们需要对这种连续变量离散化来进行KM分析，常见的离散化方法有以下四种： # 中位值 # ROC的约登指数 # 耶鲁大学出的X-tile软件(官网可下载) # survminer包的surv_cutpoint()函数 ### 补充完这一点，我们接着预测lung数据集中sex变量对病人生存的影响 ### 第二步：计算生存曲线survfit() fit <- survfit(Surv(time,status) ~ sex,data = lung) ### 如果想用其它变量计算，把sex换下就好啦 ###################################################################################### ### 第三步：可视化生存曲线 ### 1.低配版生存曲线 surv <- ggsurvplot(fit,data = lung) surv ### 2.附上关键信息 surv <- ggsurvplot( fit, # survfit函数返回的对象 data = lung, # 自己的数据集 pval = TRUE, # 是否显示P值 pval.method = TRUE, # 是否显示P值的检验方法 legend.title = 'Sex', # 自定义图例的标题 legend.labs = c('Male','Female') # 自定义分组变量的名字 ) surv ### 3.高配版生存曲线 surv <- ggsurvplot( fit, # survfit函数返回的对象 data = lung, # 自己的数据集 pval = TRUE, # 是否显示P值 pval.method = TRUE, # 是否显示P值的检验方法 legend.title='Sex', # 自定义图例的标题 legend.labs = c('Male','Female'), # 自定义分组变量的名字 conf.int = TRUE, # 是否显示置信区间 risk.table = TRUE, # 是否增加risk table risk.table.col = "strata", # 可以根据分组改变risk table的字体颜色 linetype = "strata", # 可以根据分组改变分组曲线的样式 surv.median.line = "hv", # 指出中位生存率（可加#阻断） ggtheme = theme_bw(), # 设置主题 palette = c("#1F78B4","#E31A1C") # 设置曲线的颜色 ) surv ### 4.加点操作 surv <- ggsurvplot( fit, # survfit函数返回的对象 data = lung, # 自己的数据集 pval = TRUE, # 是否显示P值 pval.method = TRUE, # 是否显示P值的检验方法 legend.title = 'Sex', # 自定义图例的标题 legend.labs = c('Male','Female'), # 自定义分组变量的名字 conf.int = TRUE, # 是否显示置信区间 cont.int.style = 'step', # 自定义置信区间的style xlab = 'Time in days', # 自定义x轴的label break.time.by = 200, # x轴上的时间以多少隔开(如果年为单位的话，一般设置为1) ggtheme = theme_classic(), # 设置主题(自定义ggplot2中的主题) risk.table = T, # 所有可用的value是：TRUE、FALSE、'absolute'、'percentage'、'nrisk_cumcensor'、'nrisk_cumevents' risk.table.y.text.col = T, # risk table左侧是否用对应分层变量的颜色注释 risk.table.y.text = F, # 是否显示分层变量的文本：如果不展示，则用颜色条进行展示 # surv.median.line = "hv", # 指出中位生存率（可加#阻断） palette = c("#1F78B4","#E31A1C") # 自定义曲线的颜色 ) surv ### 5.加上病人的删失详细信息 surv <- ggsurvplot( fit, # survfit函数返回的对象 data = lung, # 自己的数据集 pval = TRUE, # 是否显示P值 pval.method = TRUE, # 是否显示P值的检验方法 legend.title = 'Sex', # 自定义图例的标题 legend.labs = c('Male','Female'), # 自定义分组变量的名字 conf.int = TRUE, # 是否显示置信区间 cont.int.style = 'step', # 自定义置信区间的style xlab = 'Time in days', # 自定义x轴的label break.time.by = 200, # x轴上的时间以多少隔开(如果年为单位的话，一般设置为1) ggtheme = theme_classic(), # 设置主题(自定义ggplot2中的主题) risk.table = T, # 所有可用的value是：TRUE、FALSE、'absolute'、'percentage'、'nrisk_cumcensor'、'nrisk_cumevents' risk.table.y.text.col = T, # risk table左侧是否用对应分层变量的颜色注释 risk.table.y.text = F, # 是否显示分层变量的文本：如果不展示，则用颜色条进行展示 # surv.median.line = "hv", # 指出中位生存率（可加#阻断） palette = c("#1F78B4","#E31A1C"), # 自定义曲线的颜色 ncensor.plot = TRUE # 是否显示censor table：显示在时间t，删失对象的数目 ) surv ### 6.顶配版生存曲线 ### 前面我们生成了surv，这里直接用它，顶配版其实就是改改字体，加点标题或注释 surv ### 改变生存曲线的label surv$plot <- surv$plot + labs( title = "", # 自定义生存曲线的主标题 subtitle = "", # 自定义生存曲线的副标题 caption = "" # 自定义生存曲线的说明：出现在生存曲线右下角 ) surv ### 改变Risk Table的label surv$table <- surv$table + labs( title = "Number at risk", # risk table的主标题 subtitle = "", # risk table的次标题 caption = "" # risk table的说明 ) surv ### 这里title、subtitle和capton就不设置了，大家知道有这个功能就行，一般不用 ### 下面是删失table的label if(T){ surv$ncensor.plot <- surv$ncensor.plot + labs( title = "Number of censorings", # 删失table的主标题 subtitle = "over the time.", # 删失table的次标题 caption = "꒰⑅•ᴗ•⑅꒱" # 删失table的说明 )} surv surv <- ggpar( surv, font.title = c(16, "bold", "darkblue"), # 16号字体，粗体，darkblue色 font.subtitle = c(15, "bold.italic", "purple"), # 15号字体，粗斜体，紫色 font.caption = c(14, "plain", "orange"), # 14号字体，plain体，橘黄色 font.x = c(14, "bold.italic", "red"), font.y = c(14, "bold.italic", "darkred"), font.xtickslab = c(12, "plain", "darkgreen"), # 这里是改正坐标轴上的字体大小、样式和颜色 legend = "top" # 图例出现在哪个地方，这里是上方 ) surv ### 到这里已经过于花哨了◍'ㅅ'◍ ###################################################################################### ### 以上代码及注释基本可以满足需求 ### 但在某些论文中，生存曲线不仅仅要有P值和P值检验方法，还要有预测变量的HR(危险比) ### 那这个HR怎么生成呢？一句代码解决问题： summary(coxph(Surv(time, status) ~ sex, data = lung)) ### 结果如下： # Call: # coxph(formula = Surv(time, status) ~ sex, data = lung) # n= 228, number of events= 165 # # coef exp(coef) se(coef) z Pr(>|z|) # sex -0.5310 0.5880 0.1672 -3.176 0.00149 ** # --- # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 # # exp(coef) exp(-coef) lower .95 upper .95 # sex 0.588 1.701 0.4237 0.816 # # Concordance= 0.579 (se = 0.021 ) # Likelihood ratio test= 10.63 on 1 df, p=0.001 # Wald test = 10.09 on 1 df, p=0.001 # Score (logrank) test = 10.33 on 1 df, p=0.001 ### 以上结果中，exp(coef)是HR，lower .95是HR的左置信区间，upper .95是HR的右置信区间 ### 有了HR，我们该如何注释到生存曲线上？ ### 主要看pval的设置那一行： logrank_p <- round(surv_pvalue(fit)$pval,3) # 计算出logrank检验的P值 HR<- round(summary(coxph(Surv(time, status) ~ sex, data = lung))$coefficients[2],3) # 计算HR ggsurvplot( fit, # survfit函数返回的对象 data = lung, # 自己的数据集 pval = paste0('logranknp = ',logrank_p,'nHR = ',HR),# 是否显示P值 pval.method = TRUE, # 是否显示P值的检验方法 legend.title = 'Sex', # 自定义图例的标题 legend.labs = c('Male','Female'), # 自定义分组变量的名字 conf.int = TRUE, # 是否显示置信区间 cont.int.style = 'step', # 自定义置信区间的style xlab = 'Time in days', # 自定义x轴的label break.time.by = 200, # x轴上的时间以多少隔开(如果年为单位的话，一般设置为1) ggtheme = theme_classic(), # 设置主题(自定义ggplot2中的主题) risk.table = T, # 所有可用的value是：TRUE、FALSE、'absolute'、'percentage'、'nrisk_cumcensor'、'nrisk_cumevents' risk.table.y.text.col = T, # risk table左侧是否用对应分层变量的颜色注释 risk.table.y.text = F, # 是否显示分层变量的文本：如果不展示，则用颜色条进行展示 # surv.median.line = "hv", # 指出中位生存率（可加#阻断） palette = c("#1F78B4","#E31A1C"), # 自定义曲线的颜色 ) ### 以上是生存曲线作图的详细代码 ### 在实际应用过程中，要注意所选变量的P值小于0.05 ### 当我们的变量是连续变量的时候，需要选择合适的cutoff值，可以把中位值作为cutoff值 ### 但中位值进行的分组可能并没有生存意义，这时候可以换种分组方式 ### 科研界也没有规定一定要以中位值进行分组
暂无讨论，说说你的看法吧