本文首发于“生信补给站”公众号 https://mp.weixin.qq.com/s/uIUU2CBxQLQgXPVTeLRUeQ
通过RNAseq|Lasso构建预后模型,绘制风险评分的KM 和 ROC曲线构建的预后模型KM显著,还需要验证其独立性?
本文就介绍下添加临床因素构建多因素COX模型 以及 森林图,诺莫图(列线图),校准曲线以及DCA决策曲线的绘制方法。
一 载入R包,数据
仍然使用之前处理过的TCGA的SKCM数据,以及生存数据和临床数据,梳理到一个RData中了,文末有数据获取方式
library(tidyverse) library(regplot) library(rms) library(survminer) library(survival) #载入数据 load("Expr_phe_cli_riskscore.RData") riskScore_cli2 <- riskScore_cli %>% inner_join(phe) head(riskScore_cli2)
二 多因素森林图
得到riskscore后,除了使用其他数据集验证外,还需要使用多因素COX模型验证其独立性 ,然后使用ggforest函数绘制森林图
library(survminer) #构建多因素COX multiCox <- coxph(Surv(OS.time, OS) ~ age + gender + tumor_stage + riskScore2, data = riskScore_cli2) multiCox #绘制森林图 multicox_ggforest <- ggforest(multiCox, #coxph得到的Cox回归结果 data = riskScore_cli2, #数据集 main = 'Hazard ratio of multi cox', #标题 cpositions = c(0.05, 0.15, 0.35), #前三列距离 fontsize = 0.8, #字体大小 refLabel = 'reference', #相对变量的数值标签,也可改为1 noDigits = 3 #保留HR值以及95%CI的小数位数 ) multicox_ggforest
可见添加年龄,性别以及分期矫正构建多因素COX后,riskscore仍然显著。更多参数可参考Forest plot(森林图) | Cox生存分析可视化
三 Nomogram 图
Nomogram也被称为诺莫图或者列线图,在期刊出现频率越来愈多,常用于评估肿瘤学和医学的预后情况,可将Cox回归的结果进行可视化。
1,nomogram绘制
#可以输入??datadist查看详细说明 dd=datadist(riskScore_cli2) options(datadist="dd") ## 构建COX比例风险模型 f <- psm(Surv(OS.time,OS) ~ age + gender + tumor_stage + riskScore2, data = riskScore_cli2,dist='lognormal') surv <- Survival(f) # 构建生存概率函数 ## time是以”天“为单位,此处绘制1年,3,5年的生存概率 nom <- nomogram(f, fun=list(function(x) surv(365, x), function(x) surv(1095, x), function(x) surv(1825, x) ), funlabel=c("1-year OS", "3-year OS", "5-year OS")) plot(nom, xfrac=.2)
2,regplot 绘制
使用regplot包绘制
library(regplot) Cox_nomo2 <- cph(Surv(OS.time, OS) ~ age + gender + tumor_stage + riskScore2, data = riskScore_cli2, x=T, y=T) regplot(Cox_nomo2, observation = riskScore_cli2[4,], #指定某一患者 interval ="confidence", title="Nomogram", plots=c("violin", "boxes"), clickable = T, failtime = c(365,1095,1825)) #设置随访时间1年、3年和5年
3,校准曲线
绘制1 ,3, 5的校准曲线,注意m值根据样本数量选择,体现在图中多少errbar ,大概选择1/3 或者1/4即可。
#1-year cox1 <- cph(Surv(OS.time,OS) ~ age + gender + tumor_stage + riskScore2, surv=T,x=T, y=T,time.inc = 1*365,data=riskScore_cli2) cal <- calibrate(cox1, cmethod="KM", method="boot", u=1*365, m= 100, B=1000) #3-year cox3 <- cph(Surv(OS.time,OS) ~ age + gender + tumor_stage + riskScore2, surv=T,x=T, y=T,time.inc = 1*365*3,data=riskScore_cli2) ca3 <- calibrate(cox3, cmethod="KM", method="boot", u=1*365*3, m= 100, B=1000) #5-year cox5 <- cph(Surv(OS.time,OS) ~ age + gender + tumor_stage + riskScore2, surv=T,x=T, y=T,time.inc = 1*365*5,data=riskScore_cli2) ca5 <- calibrate(cox5, cmethod="KM", method="boot", u=1*365*5, m= 100, B=1000)
使用plot 函数进行可视化即可,如果想1,3,5年的校准曲线绘制在一张图中,使用add = T 即可,但是要注意调节xlim 和 ylim 的范围 。
1)绘制单条
#绘制单条 plot(cal,lwd=2,lty=1,errbar.col="black", xlab ="Nomogram-Predicted Probability of 1-Year Survival", ylab="Actual 1-Year Survival",col="blue",sub=F)
2)add = T 绘制多条
plot(cal,lwd=2,lty=1,errbar.col="black", xlim = c(0.5,1) , ylim = c(0.4,1), xlab ="Nomogram-Predicted Probability of 1,3,5Year Survival", ylab="Actual 1,3,5Year Survival",col="blue",sub=F) plot(ca3, add = T , lwd=2,lty=1,errbar.col="black",col="red",sub=F) plot(ca5, add = T , lwd=2,lty=1,errbar.col="black",col="green",sub=F)
注意add = T 以及 调节xlim 和 ylim 的范围
四 DCA 决策曲线
使用yikeshu的ggDCA-R包进行绘制 ,GitHub - yikeshu0611/ggDCA
library(caret) #remotes::install_github('yikeshu0611/ggDCA') library(ggDCA) Cox_nomo2<-rms::cph(Surv(OS.time, OS) ~ age + gender + tumor_stage + riskScore2, data = riskScore_cli2) d_train <- dca(Cox_nomo2, times=1800) ggplot(d_train)
该包含有很多可以使用的场景,参数以及ggplot2的优化方案,详情见首发!决策曲线分析R包:ggDCA