RNAseq|Lasso构建预后模型，绘制风险评分的KM 和 ROC曲线
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” 。
#默认的统计图，设定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）等 ，待补充。
◆ ◆ ◆ ◆ ◆
更多精心内容详见：精心整理（含图PLUS版）|R语言生信分析，可视化（R统计，ggplot2绘图，生信图形可视化汇总）