决策树与随机森林（8）—— 随机生存森林实战

刘花花的科研日记 医学僧的科研日记 2020-09-11 04:02
收录于合集
从MSK-IMPACT数据库下载数据，数据来源：
cbioportal.org/study/summary?id=msk_impact_2017
数据不详细介绍，简单说下，由59种肿瘤组成，共10945个样本。我们这里提取了临床信息和 Mutation signature 来进行建模。临床信息包括生存信息和远处转移信息。该文主要演示随机生存森林，Mutation signature 结合转移信息也可以构建随机森林，有兴趣者建议尝试下。
由于之前对 randomforest 进行详解过，所以这里就不用 randomForest 包进行随机森林的构建，而是使用另外一个包 randomForestSRC.

1. 数据清洗

library(survival)library(survminer)library(randomForestSRC)library(timeROC)library(data.table)library(tidyr)library(dplyr)library(tibble)library(caret)
## 临床数据清洗clin <- fread('data_clinical_patient.txt',data.table = F,skip = 4)ID <- fread('data_clinical_sample.txt',data.table = F,skip = 4)ID <- ID[,c(1,2,10)] ##选取ID信息和转移信息clin <- merge(ID,clin,by=1)clin <- clin[,c(2,3,5,7)]%>%na.omit()clin$VITAL_STATUS <- ifelse(clin$VITAL_STATUS=='ALIVE',0,1)clin$OS_MONTHS <- clin$OS_MONTHS/12 ##转化为年为单位的时间colnames(clin) <- c('ID','PM','OS','OS.time')clin <- clin[clin$OS.time!=0,] ##删除生存时间为0的样本length(unique(clin$ID)) ## 查看有无重复样本

图片

我们最后一共得到了8130个符合要求的样本
##Mutation signature 数据清洗ms <- fread('data_mutational_signature_contribution.txt',data.table = F)[,-c(1,3:5)]ms <- ms%>%column_to_rownames('NAME')%>%t()%>%as.data.frame()head(ms[,1:5])

图片

## merge clinical data and Mutation signature ## divide all data into trian and validation datadata <- merge(clin,ms,by.x=1,by.y = 0)head(data[,1:8])set.seed(2020911)size <- createDataPartition(data$OS,p=0.7,list = F)  # 随机选择70%的数据作为Train datatrain <- data[size,]validation <- data[-size,]
具有完整临床信息和 Mutation signature 信息的样本一共7268个，以 7：3 的比例分为训练集和验证集。


2. 模型构建

##构建随机生存森林rsf_t <- rfsrc(Surv(OS.time,OS)~.,data = train[,-c(1:2)],               ntree = 1000,nodesize = 15,##该值建议多调整                 splitrule = 'logrank',               importance = T,               proximity = T,               forest = T,               seed = 2020911)rsf_t

图片

plot(rsf_t)

图片

左图描述了 oob error rate 随着数的数量变化的趋势，可以看到，当 ntree > 200 左右的时候，oob error rate 趋于稳定。
右图则描述了 30 个 signature 的重要性，可以看出 signature 28 的重要性最强

3. KM+ROC 分析

## 获得每个样本的 riskscore，进一步进行下游分析score_t <- data.frame(train[,c(1,3,4)],Score=rsf_t$predicted)cut <- surv_cutpoint(score_t,'OS.time','OS','Score')cut

图片

plot(cut)

图片

Optimal Cutpoint = 252.56
## 生存分析cat <- surv_categorize(cut)fit <- survfit(Surv(OS.time,OS)~Score,cat)mytheme <- theme_survminer(font.legend = c(14,"plain", "black"),                           font.x = c(14,"plain", "black"),                           font.y = c(14,"plain", "black")) ## 自定义主题
ggsurvplot(fit,cat,           palette = 'jco',           size=1.3,           pval=T,           legend.labs=c("High","Low"),            legend.title='Score',           xlab="Time(years)",           ylab='Overall survival',           ggtheme = mytheme)

图片

## ROCcol <- c("#0073C2FF","firebrick1","orange") ## 自定义颜色tt <- timeROC(score_t$OS.time,score_t$OS,score_t$Score,              cause = 1,weighting = 'marginal',              times = seq(0.5,3,0.5),ROC = T,iid = T)tt$AUC

图片

plotAUCcurve(tt,conf.int = T,col = "tomato")

图片

plot(tt,time=1,title=FALSE,lwd=1.5,col=col[1])plot(tt,time=2,col=col[2],add=TRUE,title=FALSE,lwd=1.5)plot(tt,time=3,col=col[3],add=TRUE,title=FALSE,lwd=1.5)id <- c(paste0("1-Year AUC = ",round(tt$AUC[2],3)),        paste0("2-Year AUC = ",round(tt$AUC[4],3)),        paste0("3-Year AUC = ",round(tt$AUC[6],3)))legend("bottomright",id,       fill=col[1:3],       bty="o",cex=1,       border = NA)abline(0,1,lty=2,lwd=0.5)

图片

KM分析显示了 High risk 与不良的预后有着显著的关系；
ROC则显示了该模型在训练集中的准确性较高。

4. 模型验证

## 训练集rsf_v<- predict(rsf_t,newdata = validation,proximity=T)score_v <- data.frame(validation[,c(1,3,4)],Score=score_v$predicted)score_v$Group <- ifelse(score_v$Score>cut$cutpoint$cutpoint,'High','Low')
fit <- survfit(Surv(OS.time,OS)~Group,score_v)ggsurvplot(fit,           palette = 'jco',           size=1.3,           pval=T,           legend.labs=c("High","Low"),            legend.title='Score',           xlab="Time(years)",           ylab='Overall survival',           ggtheme = mytheme)

图片

vv <- timeROC(score_v$OS.time,score_v$OS,score_v$Score,              cause = 1,weighting = 'marginal',              times = seq(0.5,3,0.5),ROC = T,iid = T)vv$AUC

图片

可以看出，尽管KM也很显著，但是ROC则表明了该模型在训练集中的准确度很低。
在该文的第3-4张插图中，也表明错误率很高：44.8%。这其实已经说明了模型的泛化能力很弱。
我们可以通过调整参数，比如 ntree，nodesize，nsplit，mtry 等参数来优化模型，不过由于错误率太高，我估计优化后效果也不会很好。所以，该数据可能不太适合用随机生存森林来处理。

Summary：

使用随机生存森林的思路大致就是这样，当然也可以通过 rf 筛选重要变量后使用其它算法进行建模，方式很多样，重在尝试，以结果为导向，总有一个Model是让你满意的。
对于外部验证队列，可以直接使用TCGA和ICGC中提供的信息。
另外，临床信息中还有远处转移信息，可以基于转移信息和30个MS 构建模型，而且，MSK中不止有 MS 信息，还有 CNA，mutation，SV 等信息，都可以进行尝试下。
收录于合集 #决策树与随机森林
 8个
上一篇决策树与随机森林（7）—— randomForest package
喜欢此内容的人还喜欢

微信扫一扫
关注该公众号