生存分析曲线作图代码

### 生存曲线作图全代码教程
### (•̤̀ᵕ•̤́๑)ᵒᵏᵎᵎᵎᵎ
######################################################################################
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值
### 但中位值进行的分组可能并没有生存意义,这时候可以换种分组方式
### 科研界也没有规定一定要以中位值进行分组
生物信息学

GEO探针ID转换教程(代码+注释)

2020-8-13 21:33:59

生物信息学

TCGA临床数据下载与整理

2020-8-13 21:34:23

加入Q群
0 条回复 A文章作者 M管理员
    暂无讨论,说说你的看法吧
个人中心
今日签到
有新私信 私信列表
搜索