RNAseq-ML | SuperPC 算法构建预后模型并预测

生信补给站

发布于 2024-03-26 14:17:28

2530

发布于 2024-03-26 14:17:28

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

机器学习构建预后模型的文章很多，且越来越卷，动不动就是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模型

1，数据准备

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#横坐标间隔
)