【原】RNAseq|Lasso構(gòu)建預(yù)后模型，繪制風(fēng)險評分的KM 和 ROC曲線

生信補(bǔ)給站 2023-02-14 發(fā)布于北京

展開全文

可以參考使用機(jī)器學(xué)習(xí)方法構(gòu)建預(yù)后模型的集大成者文獻(xiàn)，2010年NC的文章 Pan-cancer characterization of immune-related lncRNAs identifies potential oncogenic biomarkers。文中提供了示例代碼Data availability ：http://bio-bigdata./ImmLnc 。

本次先簡單的介紹一下通過lasso 的方式構(gòu)建預(yù)后模型的方法，文章末尾會簡單的介紹下構(gòu)建完成后還可以做哪些分析。

一載入R包，數(shù)據(jù)

仍然使用之前處理過的TCGA的SKCM數(shù)據(jù)，此外需要讀入生存數(shù)據(jù)和臨床數(shù)據(jù)

library(tidyverse)library(openxlsx)library("survival")library("survminer")
load("SKCM.uni-COX.RData")#讀取臨床數(shù)據(jù)surv <- read.table("TCGA-SKCM.survival.tsv",sep = "\t" , header = T,                   stringsAsFactors = FALSE ,                   check.names = FALSE)phe <- read.xlsx("cli.xlsx",sheet = 1)

二 Lasso構(gòu)建預(yù)后模型

1 ，Lasso的輸入數(shù)據(jù)

使用glmnet包進(jìn)行Lasso分析，首先構(gòu)建lasso的生存模型需要2個數(shù)據(jù)，一個是表達(dá)量的矩陣數(shù)據(jù)(x)，一個是隨訪數(shù)據(jù) (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]]) ，這里僅為示例。

其中第五列是基因的起始列，前面四列是隨訪信息，根據(jù)自己的實際數(shù)據(jù)情況修改。

2， lasso 模型以及交叉驗證

使用glmnet函數(shù)就可以一行代碼運(yùn)行l(wèi)asso模型，cv.glmnet函數(shù)進(jìn)行交叉驗證，注意生存數(shù)據(jù)時，family處為 “cox” 。

#默認(rèn)的統(tǒng)計圖，設(shè)定label = TRUE可以給曲線加上注釋，lasso <- glmnet(x, y, family = "cox", alpha = 1 , nlambda = 1000)plot(lasso)
#交叉驗證Lasso回歸#使用glmnet包中K折交叉驗證法進(jìn)行變量篩選，設(shè)置隨機(jī)種子數(shù)并定義10折交叉set.seed(123)#注 生存分析的時間不能是0fitCV <- cv.glmnet(x, y,                    family = "cox",                   type.measure = "deviance",                   nfolds = 10)plot(fitCV)

上圖的每一條線為一個基因；下圖的每一個豎為一個基因，兩條虛線分別為lambda.min 或者 lambda.1se的結(jié)果。

這就是文獻(xiàn)中常見的lasso結(jié)果圖，下一步就是提取lasso篩選出來的基因進(jìn)行多因素COX回歸分析。

3，構(gòu)建多因素COX模型

根據(jù)lambda.min 或者 lambda.1se 獲取篩選后的基因，然后構(gòu)建多因素COX模型。

lambda.min 篩選后的基因較多，lambda.1se相對較少，一般會比較兩種情況下的模型結(jié)果然后確定選擇哪一種。這里直接使用lambda.min結(jié)果進(jìn)行示例

1）獲取lasso篩選出的基因

#λ值重新建模,選擇lambda.minfitCV$lambda.mincoefficient <- coef(fitCV, s = "lambda.min")
#系數(shù)不等于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）構(gòu)建COX模型

#提取sig_gene_mult_cox基因，構(gòu)建預(yù)后模型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

然后輸出的結(jié)果主要就是為了計算risk score（每個樣本的風(fēng)險評分），名字可以根據(jù)前面的分析過程以及想說明的問題進(jìn)行命名（免疫風(fēng)險評分，銅死亡風(fēng)險評分等）。

使用predict函數(shù)計算風(fēng)險評分

#predict函數(shù)計算風(fēng)險評分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

這樣就得到了每個樣本的評分結(jié)果。

三 KM 以及 ROC可視化

得到riskscore后還需要再使用其他數(shù)據(jù)集（GEO ，文獻(xiàn)數(shù)據(jù)，自測數(shù)據(jù)等）進(jìn)行驗證，后續(xù)會涉及。

再就是進(jìn)行一些可視化來確定該riskscore的預(yù)后是否顯著，是否獨(dú)立，相比較當(dāng)前的預(yù)后因子（分期，年齡等）是否更好？

1，KM曲線

一般可以使用KM曲線來看某因素是否預(yù)后顯著。先將riskscore進(jìn)行二分類，常見的是按照中位數(shù)（median）分為高風(fēng)險組和低風(fēng)險組，也有按照1/4進(jìn)行區(qū)分，也可以使用最優(yōu)cutoff方式R生存分析|關(guān)心的變量KM曲線不顯著，還有救嗎？進(jìn)行高低分組。

######riskScore 二分繪制KM##########riskScore_cli <- riskScore %>%   inner_join(surv)#按照中位數(shù)分為高低風(fēng)險兩組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"), #標(biāo)簽                             legend.title="RiskScore",                             title="Overall survival", #標(biāo)題                             ylab="Cumulative survival (percentage)",xlab = " Time (Days)", #更改橫縱坐標(biāo)                             censor.shape = 124,censor.size = 2,conf.int = FALSE, #刪失點(diǎn)的形狀和大小                             break.x.by = 720#橫坐標(biāo)間隔)lasso_KM

更多參數(shù)設(shè)置詳見 R|生存分析 - KM曲線，必須擁有姓名和顏值

2，ROC曲線

ROC（Receiver Operating Characteristic Curve），主要是用來確定一個模型的閾值，同時在一定程度上也可以衡量這個模型的好壞。一般情況下該曲線都應(yīng)該處于(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)))