机器学习构建预后模型的文章很多,且越来越卷,动不动就是10种模型的101种组合,这个系列会逐一的介绍这些常用于预后模型变量筛选和模型构建的机器学习方法。 本次介绍SuperPC 算法完成生存分析并做预测。 沿袭使用前面Lasso得到的SKCM.uni-COX2.RData数据(筛选过的单因素预后显著的基因),后面的更多机器学习的推文均会使用该数据,后台回复“机器学习”即可。 library(tidyverse) library(openxlsx) library("survival") library("survminer") 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 一致,可复现以及可以比较不同机器学习方法的好坏。此外注意确保训练集和验证集的基因一致 。 SuperPC 模型稍微有些区别的是需要将数据处理成list形式,生存分析模型的y 和 censoring.status分别对应时间和状态;注意x 和 featurenames要去掉OS 和 OS.time ,也就是本例中的第1,2列,需要根据自己的实际情况进行修改。set.seed(1234) data <- list(x=t(training[,-c(1,2)]), y=training$OS.time, censoring.status=training$OS, featurenames=colnames(training)[-c(1,2)])
2,构建模型构建模型的时候注意type 选择survival ,s0.perc 参数需要的是0-1之间的值,默认时候为0.5 。 #between 0 and 1: the percentile of standard deviation values added to the denominator. Default is 0.5 (the median) fit <- superpc.train(data = data, type = 'survival', s0.perc = 0.5) #default summary(fit)
3,交叉验证
使用n.fold确定几折交叉验证,其他参数在不熟悉的情况下可以使用默认参数 cv.fit <- superpc.cv(fit,data, n.threshold = 20,#default ,Number of thresholds to consider n.fold = 5, #Number of cross-validation folds n.components=3, min.features=5,#default max.features=nrow(data$x), compute.fullcv= TRUE, compute.preval=TRUE)
superpc.plotcv(cv.fit)
1,验证集验证 将验证数据同样处理成list的形式,使用superpc.predict 函数进行预测 test <- list(x=t(testing[,-c(1,2)]), y=testing$OS.time, censoring.status=testing$OS, featurenames=colnames(testing)[-c(1,2)])
pred_superpc <- superpc.predict(fit,data, test, threshold = cv.fit$thresholds[which.max(cv.fit[["scor"]][1,])], n.components = 1)
pred_superpc <- as.numeric(pred_superpc$v.pred) #输出index结果 cindex_superpc = 1-rcorr.cens(pred_superpc ,Surv(testing$OS.time, testing$OS))[[1]] print(cindex_superpc)
[1] 0.6143049
2,结合临床数据进行预后KM曲线绘制 #结合 testing$pred_superpc <- pred_superpc
#KM曲线 testing$superpc_score <- ifelse(testing$pred_superpc > median(testing$pred_superpc),"High","Low")
fit <- survfit(Surv(OS.time, as.numeric(OS)) ~ superpc_score, data=testing)
ggsurvplot(fit, data = testing, pval = T, risk.table = T, surv.median.line = "hv", #添加中位生存曲线 legend.labs=c("High risk","Low risk"), #标签 legend.title="RiskScore", title="Overall survival", #标题 ylab="Cumulative survival (percentage)",xlab = " Time (Days)", #更改横纵坐标 censor.shape = 124,censor.size = 2,conf.int = TRUE, #删失点的形状和大小 break.x.by = 720#横坐标间隔 )
RNAseq纯生信挖掘思路分享?不,主要是送你代码!(建议收藏)
|