RNAseq|Lasso构建预后模型，绘制风险评分的KM 和 ROC曲线
经过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” 。
lasso <- glmnet(x, y, family = "cox", alpha = 1 , nlambda = 1000)plot(lasso)set.seed(123)fitCV <- cv.glmnet(x, y,family = "cox",type.measure = "deviance",nfolds = 10)plot(fitCV)
' fill='%23FFFFFF'%3E%3Crect x='249' y='126' width='1' height='1'%3E%3C/rect%3E%3C/g%3E%3C/g%3E%3C/svg%3E)
上图的每一条线为一个基因；下图的每一个竖为一个基因，两条虚线分别为lambda.min 或者 lambda.1se的结果。
这就是文献中常见的lasso结果图，下一步就是提取lasso筛选出来的基因进行多因素COX回归分析。
3，构建多因素COX模型
根据lambda.min 或者 lambda.1se 获取筛选后的基因，然后构建多因素COX模型。
lambda.min 筛选后的基因较多，lambda.1se相对较少，一般会比较两种情况下的模型结果然后确定选择哪一种 。这里直接使用lambda.min结果进行示例
1）获取lasso筛选出的基因
λ值重新建模,选择lambda.minlambda.mincoefficient <- 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 %>%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）等 ，待补充。
◆ ◆ ◆ ◆ ◆
更多精心内容详见：精心整理（含图PLUS版）|R语言生信分析，可视化（R统计，ggplot2绘图，生信图形可视化汇总）