RNAseq|构建预后模型后你还需要这些图，森林图，诺莫图，校准曲线，DCA决策曲线

生信补给站

发布于 2023-08-25 10:08:22

1.2K0

发布于 2023-08-25 10:08:22

文章被收录于专栏：生信补给站

通过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)