RNAseq-ML|CoxBoost生存分析完成预后模型变量筛选以及预测

RNAseq-ML|CoxBoost生存分析完成预后模型变量筛选以及预测 by 生信补给站

机器学习构建预后模型的文章很多，且越来越卷，动不动就是10种模型的101种组合，这个系列会逐一的介绍这些常用于预后模型变量筛选和模型构建的机器学习方法。

机器学习模型1-RNAseq|Lasso构建预后模型，绘制风险评分的KM 和 ROC曲线

机器学习模型2-RNAseq-ML|randomForestSRC完成随机森林生存分析-预后模型库+1

机器学习模型3-RNAseq-ML|弹性网络回归算法Enet（Elastic Net）完成预后模型变量筛选-模型库+2

本次介绍CoxBoost生存分析，一种用boosting做COX模型的方法，同样既可以用于变量筛选也可以用于直接预测。

友情提示最后的lapply方式批量得到多个数据集的模型Cindex结果，很实用。

一数据输入，处理

沿袭使用前面Lasso得到的SKCM.uni-COX2.RData数据（筛选过的单因素预后显著的基因），后面的更多机器学习的推文均会使用该数据，后台回复“机器学习”即可。

library(tidyverse)library(survival)library(survminer)library(CoxBoost)
load("SKCM.uni-COX2.RData")module_expr.cox <- module_expr.cox %>% select(- "_PATIENT") %>%   column_to_rownames("sample")
# 7:3 拆分set.seed(1234) #ind <- sample(nrow(module_expr.cox2),nrow(module_expr.cox2) * 0.7 )train <- module_expr.cox2[ind,]test <- module_expr.cox2[-ind,]## 确保训练集和验证集的基因一致 gene_com <- intersect(colnames(train) ,colnames(test))training <- train %>% select(gene_com)testing <- test %>% select(gene_com)
training[1:4,1:4]#                OS OS.time      TYRP1   IGKV4_1#TCGA-D3-A5GO-06A  0    4195  0.3807211 2.1135527#TCGA-D3-A1Q8-06A  1     854  0.5981395 7.5734985#TCGA-FW-A5DY-06A  0     587  0.2427982 8.5617523#TCGA-EE-A29B-06A  1    2588 10.3296179 0.4794816

添加seed可以保证training 和 testing 一致，可复现以及比较不同机器学习方法的好坏。

二构建CoxBoost模型

CoxBoost函数使用如下，可以看到除了随访时间和结局外，还有stepno 和 penalty两个参数需要确认。

CoxBoost(time=obs.time,status=status,x=x,         stepno=100,         penalty=100)

1，确定最优penalty

使用optimCoxBoostPenalty函数筛选当前数据的最优penalty ，将得到的pen$penalty 定为最终模型的参数

pen <- optimCoxBoostPenalty(training[,'OS.time'],                            training[,'OS'],                            as.matrix(training[,-c(1,2)]),                            trace=TRUE,                            #start.penalty=500,                            parallel = T)
#> pen$penalty#[1] 5400

2，确定最优stepno

使用cv.CoxBoost函数确定最优stepno，取cv.res$optimal.step的值

#number of folds to be used for cross-validationcv.res <- cv.CoxBoost(training[,'OS.time'],                      training[,'OS'],                      as.matrix(training[,-c(1,2)]),                      maxstepno=500,                      K=5,                      type="verweij",                      penalty= pen$penalty,                      multicore=1)
#> cv.res$optimal.step#[1] 86

3，构建模型

使用上面得到的参数构建CoxBoost模型

fit <- CoxBoost(training[,'OS.time'],                training[,'OS'],                as.matrix(training[,-c(1,2)]),                stepno=cv.res$optimal.step,                penalty=pen$penalty)summary(fit)plot(fit)

可以在summary中复制筛选后的基因，也可以使用以下方式获取

step.logplik<-predict(fit,newdata=as.matrix(training[,-c(1,2)]),                      newtime=training[,'OS.time'],                      newstatus=training[,'OS'],                      at.step=0:69,                      type="logplik")
print(fit$xnames[fit$coefficients[which.max(step.logplik),]!=0])

#[1] "CYTL1"    "KIT"      "CCDC8"    "ABCC2"    "IFITM1"   "CRABP2"   "BOK"      "TTYH2"    "CXCL11"  #[10] "PYROXD2"  "HSPA7"    "MT2A"     "NEAT1"    "DYNLT3"   "MIR4664"  "KIAA0040" "HPDL"     "UBA7"

和summary得到的结果是一样的，这样方便后续函数中自动化的变量筛选，验证。

三 CoxBoost模型验证

同样介绍两种预后模型的验证方式：一是筛选变量然后构建COX预后模型，二是直接预测，验证预后效果。

1，筛选变量构建COX模型

直接在矩阵文件中筛选上述的基因，然后构建COX模型，以及后续的一系列分析，参考前面即可。

2，直接预测结果预后

不筛选变量直接预测，直接预测得到Cindex值，这里介绍使用lapply的方式一次得到多个数据集的Cindex结果。

当多个数据集，多种构建模型的方法都如下输出时候，就是101组合机器学习的雏形了。

val_dd_list <- list(train=training,                    test1=testing)result <- data.frame()set.seed(1234)
rs <- lapply(val_dd_list,function(x){cbind(x[,1:2],                                           RS=as.numeric(predict(fit,                                                                 newdata=x[,-c(1,2)],                                                                 newtime=x[,1],                                                                  newstatus=x[,2],                                                                  type="lp")))})
cc <- data.frame(Cindex=sapply(rs,function(x){  as.numeric(summary(coxph(Surv(OS.time,OS)~RS,x))$concordance[1])}))%>%  rownames_to_column('ID')cc$Model <- paste0('CoxBoost')result <- rbind(result,cc)
result

     ID    Cindex    Model1 train 0.6987871 CoxBoost2 test1 0.6429533 CoxBoost

后续还可以做的更多分析这里就不赘述了，参考之前的推文。

参考资料:

[1] A machine learning framework develops a DNA replication stress model for predicting clinical outcomes and therapeutic vulnerability in primary prostate cancer

[2] Machine learning-based integration develops an immune-derived lncRNA signature for improving outcomes in colorectal cancer

◆ ◆ ◆ ◆ ◆

精心整理（含图PLUS版）|R语言生信分析，可视化（R统计，ggplot2绘图，生信图形可视化汇总）

RNAseq纯生信挖掘思路分享？不，主要是送你代码！（建议收藏）

ixxmu / mp_duty