RNAseq-ML|CoxBoost生存分析完成预后模型变量筛选以及预测
RNAseq-ML|CoxBoost生存分析完成预后模型变量筛选以及预测
机器学习构建预后模型的文章很多，且越来越卷，动不动就是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] 54002，确定最优stepno
使用cv.CoxBoost函数确定最优stepno，取cv.res$optimal.step的值
#number of folds to be used for cross-validation
cv.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] 863，构建模型
使用上面得到的参数构建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 Model
1 train 0.6987871 CoxBoost
2 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
腾讯云开发者
