RNAseq-ML|randomForestSRC完成随机森林生存分析-预后模型库+1
RNAseq-ML|randomForestSRC完成随机森林生存分析-预后模型库+1
机器学习构建预后模型的文章很多，且越来越卷，动不动就是10种模型的101种组合，这个系列会逐一的介绍这些常用于预后模型变量筛选和模型构建的机器学习方法。
作者代码公开在github上了，GitHub - Zaoqu-Liu/IRLS: Machine learning-based integrative analysis develops an immune-derived lncRNA signature for improving clinical outcomes in colorectal cancer 可以自行下载 或者 后台回复 “机器学习”获取下载好的。前面介绍过了RNAseq|Lasso构建预后模型，绘制风险评分的KM 和 ROC曲线，本次介绍使用randomForestSRC完成随机森林的生存分析。
一 数据输入，处理
沿袭使用前面Lasso得到的SKCM.uni-COX.RData数据（筛选过的单因素预后显著的基因），后面的更多机器学习的推文均会使用该数据
#载入R包
library(tidyverse)
library(openxlsx)
library("survival")
library("survminer")
library(randomForestSRC)
load("SKCM.uni-COX.RData")
module_expr.cox2 <- module_expr.cox %>% select(- "_PATIENT") %>%
column_to_rownames("sample")
module_expr.cox2[1:6,1:6]
将数据处理成以上形式，含有 随访时间 + 生存状态 + 基因表达信息。1，数据集拆分正常情况下是TCGA构建模型，然后在GEO中进行验证。这里仅为示例，直接按照7：3的比例将TCGA数据拆分为训练集和验证集（后面会介绍更多拆分方法）
# 7:3 拆分
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-EE-A2MM-06A 1 5107 1.38460143 5.2408878
#TCGA-EE-A2GE-06A 0 5286 0.04187911 10.1611678
#TCGA-ER-A194-01A 1 1354 9.56901508 0.3122559
#TCGA-EB-A44R-06A 1 315 0.06131739 7.3046339注：训练集和验证集的基因一致，不然可能存在无法验证的情况。
二 构建随机森林生存模型
1，rfsrc函数构建RSF 生存模型
注意设置随机种子seed，方便以后复现；此外nodesize 值可以多设置几个尝试
fit <- rfsrc(Surv(OS.time,OS)~.,data = training,
ntree = 1000,
nodesize = 10,
splitrule = 'logrank',
importance = T,
proximity = T,
forest = T,
seed = 1234)
可以看到该模型含有320样本，537个基因。
2，重要性变量
使用随机森林生存分析进行变量筛选，主要依据的就是每个基因的重要性值 ，该数据在fit$importance中，这里示例查看TOP20 的基因注意：这里的重要性基因不会得到文献中常提到的基因前面的系数，系数可以通过将重要基因进行多因素COX生存分析得到。
importance_gene <- data.frame(fit$importance) %>%
rownames_to_column("gene") %>%
arrange(- fit.importance) %>%
head(20)
importance_gene
（1）使用plot函数直接可视化
plot(fit,10)
（2）使用ggplot2绘制柱形图使用reorder函数进行排序
ggplot(data=importance_gene, aes(x = reorder(gene, fit.importance),
y=fit.importance,fill=gene)) +
geom_bar(stat="identity") +
theme_classic() +
theme(legend.position = 'none') +
coord_flip()
这样就不会拥挤在一起，且可自定义颜色。
三 RSF模型验证
这里面介绍2种验证方式，第一种起到了和Lasso一样的筛选基因的作用，第二种是直接验证。
1，使用RSF得到的重要基因构建COX模型
（1）在上面的importance_gene文件中，根据fit.importance设置阈值，然后选出候选基因 或者
（2）在上面的importance_gene文件中，直接选择TOP多少的基因作为候选基因。
然后将候选基因构建多因素COX模型，这样就可以得到文献中常见的基因系数。
注：这里的阈值和TOP没有固定的cutoff ，结果导向即可。
2，RSF模型直接验证集预测
直接使用验证集验证模型，得到每个样本的系数，然后可以使用生存分析得到Cindex以及KM曲线等。
（1）C-index
fit.p <- predict(fit, as.data.frame(testing))
testing$RSF_p <- as.vector(fit.p$predicted)
#计算C index
testing_surv <- coxph(Surv(OS.time, OS) ~ fit.p$predicted,data = testing)
summary(testing_surv)$concordance C se(C)
0.64523954 0.03881865（2）KM曲线
testing$RSF_score <- ifelse(testing$RSF_p > median(testing$RSF_p),"High","Low")
fit <- survfit(Surv(OS.time, as.numeric(OS)) ~ RSF_score, data=testing)
ggsurvplot(fit, data = testing,
pval = T,
risk.table = T,
surv.median.line = "hv", #添加中位生存曲线
palette=c("red", "blue"), #更改线的颜色
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 = FALSE, #删失点的形状和大小
break.x.by = 720#横坐标间隔
)
这样就完成了随机森林生存模型筛选变量或者预测的介绍，Lasso之外可以多一种尝试了。
参考资料:
[1] Getting starting with the randomForestSRC R-package for random forest analysis of regression, classification, survival and more • Fast Unified Random Forests with randomForestSRC
[2] Machine learning-based integration develops an immune-derived lncRNA signature for improving outcomes in colorectal cancer
◆ ◆ ◆ ◆ ◆
精心整理（含图PLUS版）|R语言生信分析，可视化（R统计，ggplot2绘图，生信图形可视化汇总）