本文首发于“生信补给站”公众号 https://mp.weixin.qq.com/s/OmUXjkdySWmOeEkvrMPblQ
经过RNAseq|批量单因素生存分析 + 绘制森林图分析后得到了预后显著的基因集。后续的常见做法是通过机器学习（lasso，随机森林，SVM等）方法进行变量（基因）筛选，然后构建预后模型。
可以参考使用机器学习方法构建预后模型的集大成者文献，2010年NC的文章 Pan-cancer characterization of immune-related lncRNAs identifies potential oncogenic biomarkers。文中提供了示例代码Data availability ：http://bio-bigdata.hrbmu.edu.cn/ImmLnc 。
本次先简单的介绍一下通过lasso 的方式构建预后模型的方法，文章末尾会简单的介绍下构建完成后还可以做哪些分析。
一 载入R包，数据
仍然使用之前处理过的TCGA的SKCM数据，此外需要读入生存数据和临床数据
library(tidyverse) library(openxlsx) library("survival") library("survminer") load("SKCM.uni-COX.RData") #读取临床数据 surv <- read.table("TCGA-SKCM.survival.tsv",sep = "\t" , header = T, stringsAsFactors = FALSE , check.names = FALSE) phe <- read.xlsx("cli.xlsx",sheet = 1)
二 Lasso构建预后模型
1 ，Lasso的输入数据
使用glmnet包进行Lasso分析，首先构建lasso的生存模型需要2个数据，一个是表达量的矩阵数据(x)，一个是随访数据 (y)
library(glmnet) DEG_met_expr.lasso <- module_expr.cox %>% select(sample ,OS,OS.time,KM_sig$Variable) %>% column_to_rownames("sample") DEG_met_expr.lasso[1:2,1:5] # OS OS.time TYRP1 IGKV4_1 IGHG1 #TCGA-EE-A2GJ-06A 1 2270 9.736880 1.011539 5.668523 #TCGA-EE-A2GI-06A 0 1482 7.836613 8.484382 11.246364 #lasso需要矩阵类型，使用as.matrix #x <- as.matrix(DEG_met_expr.lasso[,5:dim(DEG_met_expr.lasso)[2]]) # x <- as.matrix(DEG_met_expr.lasso[,5:100]) y <- Surv(DEG_met_expr.lasso$OS.time,DEG_met_expr.lasso$OS==1)
注：正常情况下x <- as.matrix(data[,5:dim(data)[2]]) ，这里仅为示例 。
其中第五列是基因的起始列 ，前面四列是随访信息，根据自己的实际数据情况修改。
2， lasso 模型以及交叉验证
使用glmnet函数就可以一行代码运行lasso模型，cv.glmnet函数进行交叉验证，注意生存数据时，family处为 “cox” 。
#默认的统计图，设定label = TRUE可以给曲线加上注释， lasso <- glmnet(x, y, family = "cox", alpha = 1 , nlambda = 1000) plot(lasso) #交叉验证Lasso回归 #使用glmnet包中K折交叉验证法进行变量筛选，设置随机种子数并定义10折交叉 set.seed(123) #注 生存分析的时间不能是0 fitCV <- cv.glmnet(x, y, family = "cox", type.measure = "deviance", nfolds = 10) plot(fitCV)
上图的每一条线为一个基因；下图的每一个竖为一个基因，两条虚线分别为lambda.min 或者 lambda.1se的结果。
这就是文献中常见的lasso结果图，下一步就是提取lasso筛选出来的基因进行多因素COX回归分析。
3，构建多因素COX模型
根据lambda.min 或者 lambda.1se 获取筛选后的基因，然后构建多因素COX模型。
lambda.min 筛选后的基因较多，lambda.1se相对较少，一般会比较两种情况下的模型结果然后确定选择哪一种 。这里直接使用lambda.min结果进行示例
1）获取lasso筛选出的基因
#λ值重新建模,选择lambda.min fitCV$lambda.min coefficient <- coef(fitCV, s = "lambda.min") #系数不等于0的为纳入的变量（基因） Active.index <- which(as.numeric(coefficient) != 0) Active.coefficient <- as.numeric(coefficient)[Active.index] sig_gene_mult_cox <- rownames(coefficient)[Active.index] #查看具体哪些基因 sig_gene_mult_cox #[1] "SFRP1" "LOXL4" "BCAN" "BAALC" "CXCL10" "KIT" "KCNN4" "MZB1"
2）构建COX模型
#提取sig_gene_mult_cox基因，构建预后模型 DEG_met_expr.lasso_cox <- DEG_met_expr.lasso %>% dplyr::select(OS,OS.time,sig_gene_mult_cox) multiCox <- coxph(Surv(OS.time, OS) ~ ., data = DEG_met_expr.lasso_cox) summary(multiCox)
4，计算risk score
然后输出的结果主要就是为了计算risk score（每个样本的风险评分），名字可以根据前面的分析过程以及想说明的问题进行命名（免疫风险评分，铜死亡风险评分等）。
使用predict函数计算风险评分
#predict函数计算风险评分 riskScore=predict(multiCox,type="risk",newdata=DEG_met_expr.lasso_cox) riskScore<-as.data.frame(riskScore) riskScore$sample <- rownames(riskScore) head(riskScore,2) # riskScore sample #TCGA-EE-A2GJ-06A 0.6582442 TCGA-EE-A2GJ-06A #TCGA-EE-A2GI-06A 0.6881212 TCGA-EE-A2GI-06A
这样就得到了每个样本的评分结果。
三 KM 以及 ROC可视化
得到riskscore后还需要再使用其他数据集（GEO ，文献数据，自测数据等）进行验证，后续会涉及。
再就是进行一些可视化来确定该riskscore的预后是否显著，是否独立，相比较当前的预后因子（分期，年龄等）是否更好？
1，KM曲线
一般可以使用KM曲线来看 某因素 是否预后显著 。先将riskscore进行二分类，常见的是按照中位数（median）分为高风险组和低风险组，也有按照1/4进行区分，也可以使用最优cutoff方式R生存分析|关心的变量KM曲线不显著，还有救吗？进行高低分组。
######riskScore 二分绘制KM########## riskScore_cli <- riskScore %>% inner_join(surv) #按照中位数分为高低风险两组 riskScore_cli$riskScore2 <- ifelse(riskScore_cli$riskScore > median(riskScore_cli$riskScore), "High","Low") #KM分析 fit <- survfit(Surv(OS.time, as.numeric(OS)) ~ riskScore2, data=riskScore_cli) lasso_KM <- ggsurvplot(fit, data = riskScore_cli, 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_KM
更多参数设置详见 R|生存分析 - KM曲线 ，必须拥有姓名和颜值
2，ROC曲线
ROC（Receiver Operating Characteristic Curve），主要是用来确定一个模型的阈值，同时在一定程度上也可以衡量这个模型的好坏。 一般情况下该曲线都应该处于(0, 0)和(1, 1)连线的上方（如果在下方改变marker的方向）。使用ROC 曲线可以比较直观的展示模型的好坏，处于ROC 曲线下方的那部分面积的大小越大越好，也就是Area Under roc Curve(AUC)值。
绘制ROC曲线的方式很多种，这里使用timeROC绘制 1年，3年和5年的ROC曲线
library(timeROC) with(riskScore_cli, ROC_riskscore <<- timeROC(T = OS.time, delta = OS, marker = riskScore, cause = 1, weighting = "marginal", times = c(365,1080,1800), ROC = TRUE, iid = TRUE) ) plot(ROC_riskscore, time = 365, col = "red", add = F,title = "") plot(ROC_riskscore, time = 1080, col = "blue", add = T) plot(ROC_riskscore, time = 1800, col = "purple", add = T) legend("bottomright",c("1-Year","3-Year","5-Year"),col=c("red","blue","purple"),lty=1,lwd=2) text(0.5,0.2,paste("1-Year AUC = ",round(ROC_riskscore$AUC[1],3))) text(0.5,0.15,paste("3-Year AUC = ",round(ROC_riskscore$AUC[2],3))) text(0.5,0.1,paste("5-Year AUC = ",round(ROC_riskscore$AUC[3],3)))
更多可参考R|timeROC-分析
（1）可以看风险高低两组间的免疫浸润程度是否差异RNAseq|免疫浸润也杀疯了，cibersoert？xCELL？ESTIMATE？你常用哪一个？
（2）可以和临床指标一起构建多因素COX模型，查看该riskscore的独立性Forest plot（森林图） | Cox生存分析可视化
（3）可以看风险高低两组间的差异情况，进而富集分析或者GSEA，GSVA等分析clusterProfiler|GSEA富集分析及可视化
（4）还可以看高低两组间的 药物反应（IC50） 和 免疫反应（IPS，TIDE）等 ，待补充。
