測(cè)序上游分析系列：

mRNA-seq轉(zhuǎn)錄組二代測(cè)序從raw reads到表達(dá)矩陣：上中游分析pipeline
miRNA-seq小RNA高通量測(cè)序pipeline：從raw reads，鑒定已知miRNA-預(yù)測(cè)新miRNA冀泻，到表達(dá)矩陣【一】
miRNA-seq小RNA高通量測(cè)序pipeline：從raw reads常侣，鑒定已知miRNA-預(yù)測(cè)新miRNA，到表達(dá)矩陣【二】

其他文章系列：ggplot2作圖篇：

（1）基于ggplot2的RNA-seq轉(zhuǎn)錄組可視化：總述和分文目錄
（2）測(cè)序結(jié)果概覽：基因表達(dá)量rank瀑布圖弹渔，高密度表達(dá)相關(guān)性散點(diǎn)圖胳施，注釋特定基因及errorbar的表達(dá)相關(guān)性散點(diǎn)圖繪制
（3）單/多個(gè)基因在組間同時(shí)展示的多種選擇：分組小提琴圖、分組/分面柱狀圖肢专、單基因蜂群點(diǎn)圖拼圖的繪制
（4）配對(duì)樣本基因表達(dá)趨勢(shì)：優(yōu)化前后的散點(diǎn)連線圖+拼圖繪制

常見的TCGA數(shù)據(jù)挖掘辦法之一巾乳，是通過(guò)差異基因分析獲得差異表達(dá)基因您没，然后從中篩選出部分表達(dá)水平與患者生存相關(guān)聯(lián)的候選基因，對(duì)它們的表達(dá)水平進(jìn)行多因素cox回歸構(gòu)建風(fēng)險(xiǎn)模型胆绊，評(píng)估風(fēng)險(xiǎn)模型的預(yù)測(cè)能力（ROC曲線），并用Kaplan-Meier生存分析評(píng)估模型的風(fēng)險(xiǎn)評(píng)分是否有意義欧募。

本例仍使用之前教程選擇的的TCGA-LUSC白色人種肺鱗癌表達(dá)譜數(shù)據(jù)压状。
TCGA數(shù)據(jù)下載整合操作見前文：從TCGA數(shù)據(jù)庫(kù)下載并整合清洗高通量腫瘤表達(dá)譜-臨床性狀數(shù)據(jù)
DESeq2差異分析見前文：TCGA數(shù)據(jù)整合后進(jìn)行DESeq2差異表達(dá)分析和基于R的多種可視化

本文使用的數(shù)據(jù)為下載的HTSeq-counts數(shù)據(jù)，已經(jīng)通過(guò)了整合和刪拾數(shù)據(jù)的清洗跟继，并通過(guò)DESeq2完成了差異分析种冬。本文用到的變量和對(duì)象為：

(1) dds_DE: the object generated by DESeq2. It contains the results of differential expression such as symbols of DE-genes and normalized counts.
(2) condition_table_cancer: the data frame containing the TCGA_IDs, overall_survival and vital status of cancer patients.

需要的R包：DESeq2包（轉(zhuǎn)錄組差異分析），survival包（cox回歸）舔糖，survminer包（Kaplan-Meier Plot可視化）娱两，dplyr包（字符串處理），glmnet包（Lasso回歸）金吗，ggplot2包（數(shù)據(jù)可視化）十兢，GGally包（繪制相關(guān)性矩陣圖），rms包（計(jì)算VIF方差膨脹因子）摇庙，survivalROC包（繪制time dependent ROC曲線）旱物，plotROC包（繪制ROC曲線）。

install.packages(c('DESeq2', 'survival', 'survminer', 'dplyr', 'glmnet', 'ggplot2', 'GGally', 'rms', 'survivalROC', 'plotROC'))
library('DESeq2')
library('survival')
library('survminer')
library('dplyr')
library('glmnet')
library('ggplot2')
library('GGally')
library('rms')
library('survivalROC')
library('plotROC')

1. 獲取差異表達(dá)基因卫袒，以及差異基因的標(biāo)準(zhǔn)化counts矩陣宵呛。

我們假設(shè)現(xiàn)在DESeq2已經(jīng)做好差異分析，獲得了dds_DE object夕凝！

使用results函數(shù)獲取cancer組相對(duì)normal組差異表達(dá)p值<0.05的基因和相關(guān)信息宝穗。隨后選取符合|log2FoldChange|差異大于2（|FoldChange|大于4）和FDR<0.01的差異基因子集。resSigAll對(duì)象的rownames記錄了差異基因名信息码秉。

隨后獲取基因的normalized counts matrix并做vst轉(zhuǎn)換逮矛，僅提取差異基因的normalized counts進(jìn)行后續(xù)分析。

res_DE <- results(dds_DE, alpha=0.05,contrast=c("sample_type","cancer","normal"))
resOrdered <- res_DE[order(res_DE$padj),]
resSig <- subset(resOrdered,padj < 0.01)
resSigAll <- subset(resSig, (log2FoldChange < (-2)|log2FoldChange > 2))
rld <- vst(dds_DE, blind = T)
expr_norm <- assay(rld)
DESeq_norm_vst_for_survival <- expr_norm[resSigAll@rownames,]

2. 構(gòu)建回歸分析所需的信息表泡徙。

采用Survival R包對(duì)下載數(shù)據(jù)進(jìn)行回歸分析橱鹏，需要構(gòu)建一個(gè)至少包含樣本生存時(shí)間、刪失狀態(tài)的data frame堪藐。本例中莉兰，在數(shù)據(jù)整合建立condition_table時(shí)已經(jīng)獲得了一個(gè)cancer patients的子集condition_table_cancer，直接使用它建立信息表礁竞，并合并潛在預(yù)后預(yù)測(cè)因子：候選基因的normalized counts（t()轉(zhuǎn)換至行為樣本糖荒，列為基因）。

注：survival包認(rèn)定0為censored（無(wú)事件發(fā)生/刪失）模捂，1為diseased（有結(jié)局事件發(fā)生）捶朵。所以在本例中'Alive'轉(zhuǎn)換為0蜘矢，'Dead'轉(zhuǎn)換為1。
此外有基因名中帶有R不被認(rèn)可的dash '-'综看，所以將data frame中colnames的 '-' 轉(zhuǎn)換為下劃線 '_'品腹。

library('dplyr')
survival_cancer <- condition_table_cancer
#Alive=0, Dead=1
survival_cancer$censoring_status <- gsub(survival_cancer$censoring_status, pattern = 'Alive', replacement = '0') %>%
  gsub(pattern = 'Dead', replacement = '1') %>% as.numeric()
rownames(survival_cancer) <- survival_cancer$TCGA_IDs
survival_cancer <- cbind(survival_cancer, t(DESeq_norm_vst_for_survival)[survival_cancer$TCGA_IDs,])
#formula function in the next step will not recognize gene names like 'NKX1-2'. We have to change '-' to '_'.
colnames(survival_cancer) <- gsub(colnames(survival_cancer), pattern = '-', replacement = '_')

survival_cancer數(shù)據(jù)框的部分展示

3. 對(duì)差異基因進(jìn)行批量單因素COX回歸。

很多研究使用了先批量單因素COX回歸選出更少的候選基因红碑，隨后進(jìn)行多因素COX回歸建模的策略舞吭。雖然這個(gè)策略由于忽略了變量間的相互關(guān)系而非常有爭(zhēng)議，但是這里暫且不論策略對(duì)錯(cuò)析珊，我們先用代碼嘗試實(shí)現(xiàn)批量單因素COX回歸羡鸥。

單因素COX回歸建模依賴于survival包c(diǎn)oxph()函數(shù)，PH檢驗(yàn)依賴于cox.zph()函數(shù)忠寻。這里的思路是自建一個(gè)函數(shù)惧浴，然后輸入需批量操作的基因名vector，和已經(jīng)加入這些基因normalized counts的survival_cancer信息表奕剃。通過(guò)R特有的向量化操作lappy同時(shí)進(jìn)行批量計(jì)算各基因的風(fēng)險(xiǎn)比衷旅，并檢驗(yàn)變量是否符合cox等比例風(fēng)險(xiǎn)模型的前提PH假設(shè)，得到符合z值p值<0.05和符合PH假定（cox.zph檢驗(yàn)p值>0.05）的基因祭饭。最終輸出含有g(shù)ene名芜茵，單因素cox的β值，危險(xiǎn)比HR倡蝙，z值的p值九串，Wald檢驗(yàn)和似然比檢驗(yàn)p值的data frame。

注：gene名需要時(shí)刻注意寺鸥，將'-'轉(zhuǎn)換為'_'才能與survival_cancer匹配猪钮。函數(shù)中已鑲嵌這個(gè)轉(zhuǎn)換。

#filter potential useful sig genes using univariate cox regression.
uni_cox_in_bulk <- function(gene_list, survival_info_df){
  library('survival')
  gene_list <- gsub(gene_list, pattern = '-', replacement = '_')
  uni_cox <- function(single_gene){
    formula <- as.formula(paste0('Surv(overall_survival, censoring_status)~', single_gene))
    surv_uni_cox <- summary(coxph(formula, data = survival_cancer))
    ph_hypothesis_p <- cox.zph(coxph(formula, data = survival_cancer))$table[1,3]
    if (surv_uni_cox$coefficients[,5]<0.05 & ph_hypothesis_p>0.05){  #get the pvalue
      single_cox_report <- data.frame('uni_cox_sig_genes'=single_gene,
                                      'beta'=surv_uni_cox$coefficients[,1],
                                      'Hazard_Ratio'=exp(surv_uni_cox$coefficients[,1]),
                                      'z_pvalue'=surv_uni_cox$coefficients[,5],
                                      'Wald_pvalue'=as.numeric(surv_uni_cox$waldtest[3]),
                                      'Likelihood_pvalue'=as.numeric(surv_uni_cox$logtest[3]))
      single_cox_report
    }
  }
  uni_cox_list <- lapply(gene_list, uni_cox)
  do.call(rbind, uni_cox_list)
}
uni_cox_df <- uni_cox_in_bulk(gene_list = resSigAll@rownames, survival_info_df = survival_cancer)

最終得到279個(gè)有意義的基因胆建。輸出的data frame如圖：

批量單因素cox回歸的部分結(jié)果展示

Optional:

以上的代碼嵌套了兩次自定義函數(shù)烤低，雖然不是最優(yōu)化的做法，但較好地利用了R的向量化特性笆载。如下代碼的目的是等價(jià)的扑馁，也是批量單因素cox回歸，但放棄了向量化凉驻，在R中使用了for循環(huán)腻要，運(yùn)行速度明顯慢于上方代碼。

uni_cox_sig_genes <- c()
beta_co <- c()
HR <- c()
z_p <- c()
Wald_p <- c()
Likelihood_p <- c()
for (candidate_gene in gsub(resSigAll@rownames, pattern = '-', replacement = '_')){
  formula <- as.formula(paste0('Surv(overall_survival, censoring_status)~', candidate_gene))
  surv_uni_cox <- summary(coxph(formula, data = survival_cancer))
  ph_hypothesis_p <- cox.zph(coxph(formula, data = survival_cancer))$table[1,3]
  if (surv_uni_cox$coefficients[,5]<0.05 & ph_hypothesis_p>0.05){  #get the pvalue
    uni_cox_sig_genes <- append(uni_cox_sig_genes, candidate_gene)
    beta_co <- append(beta_co, surv_uni_cox$coefficients[,1])
    HR <- append(HR, exp(surv_uni_cox$coefficients[,1]))
    z_p <- append(z_p, surv_uni_cox$coefficients[,5])
    Wald_p <- append(Wald_p, as.numeric(surv_uni_cox$waldtest[3]))
    Likelihood_p <- append(Likelihood_p, as.numeric(surv_uni_cox$logtest[3]))
  }
}
uni_cox_df <- data.frame('uni_cox_sig_genes'=uni_cox_sig_genes, 'beta'=beta_co, 'Hazard_Ratio'=HR, 'z_pvalue'=z_p, 'Wald_pvalue'=Wald_p, 'Likelihood_pvalue'=Likelihood_p)

4. Lasso回歸篩選具有代表性的變量涝登。

單因素回歸篩選出的變量仍舊較多雄家，我們嘗試用Lasso回歸獲得較少有意義的變量。Lasso回歸的前提是變量數(shù)（這里是差異基因胀滚，n=2726）>樣本數(shù)（這里是cancer patients數(shù)趟济，n=344）乱投，回歸的目的是使deviance最小化，結(jié)果是將不重要變量的系數(shù)壓縮為0顷编，僅保留較少的重要變量系數(shù)不為0戚炫。

Lasso回歸依賴glmnet包。我個(gè)人下載這個(gè)包的時(shí)候遇到了R版本不匹配的問(wèn)題媳纬，如果有朋友也遇到無(wú)法install.packages的問(wèn)題嘹悼，可以右轉(zhuǎn)github下載，或網(wǎng)頁(yè)版CRAN手動(dòng)下載安裝层宫，這里不多說(shuō)。

注：glmnet需要輸入x的格式為矩陣（matrix）其监，y的格式為double萌腿，否則會(huì)報(bào)錯(cuò)。

#about glmnet: x should be in format of matrix, and time&status in y should be in double format.
x <- as.matrix(survival_cancer[,gsub(resSigAll@rownames, pattern = '-', replacement = '_')])
y <- survival_cancer[,c('overall_survival', 'censoring_status')]
names(y) <- c('time', 'status')
y$time <- as.double(y$time)
y$status <- as.double(y$status)
y <- as.matrix(survival::Surv(y$time, y$status))
lasso_fit <- cv.glmnet(x, y, family='cox', type.measure = 'deviance')
coefficient <- coef(lasso_fit, s=lasso_fit$lambda.min)
Active.Index <- which(as.numeric(coefficient) != 0)
active.coefficients <- as.numeric(coefficient)[Active.Index]
sig_gene_multi_cox <- rownames(coefficient)[Active.Index]

這一次我們僅僅獲得了4個(gè)candidates：

Lasso回歸非常簡(jiǎn)單粗暴又好用

5. 利用篩選出來(lái)的少數(shù)candidates建立多因素COX回歸模型抖苦。

話不多說(shuō)毁菱，依舊用survival包的cox.ph()建模，用coxzph()檢驗(yàn)每個(gè)因子的PH假設(shè)锌历。然后計(jì)算回歸模型中每個(gè)因子的VIF和相關(guān)系數(shù)贮庞，判斷因素間可能存在的共線性。通過(guò)PH假設(shè)和共線性檢驗(yàn)的變量再進(jìn)行第二次建模究西，隨后繪制COX回歸森林圖窗慎。

#perform the multi-variates cox regression using qualified genes.
formula_for_multivariate <- as.formula(paste0('Surv(overall_survival, censoring_status)~', paste(sig_gene_multi_cox, sep = '', collapse = '+')))
multi_variate_cox <- coxph(formula_for_multivariate, data = survival_cancer)
#check if variances are supported by PH hypothesis.
ph_hypo_multi <- cox.zph(multi_variate_cox)
#The last row of the table records the test results on the GLOBAL model. Delete it.
ph_hypo_table <- ph_hypo_multi$table[-nrow(ph_hypo_multi$table),]
#Remove variances not supported by ph hypothesis and perform the 2nd regression.
formula_for_multivariate <- as.formula(paste0('Surv(overall_survival, censoring_status)~', paste(rownames(ph_hypo_table)[ph_hypo_table[,3]>0.05], sep = '', collapse = '+')))
multi_variate_cox_2 <- coxph(formula_for_multivariate, data = survival_cancer)

到這里我們可以發(fā)現(xiàn)4個(gè)candidates中有1個(gè)未通過(guò)coxzph()的檢驗(yàn)，只有3個(gè)基因納入了第二次的建模卤材。

ADGRD1 failed the test and should be removed

關(guān)于檢測(cè)變量間的共線性遮斥，需要綜合考慮：變量間的相關(guān)系數(shù)<0.5和vif平方根<2均被提出是可行的檢驗(yàn)方法。

我們進(jìn)行相關(guān)系數(shù)扇丛，vif的計(jì)算术吗，并進(jìn)行相關(guān)性矩陣的可視化。

#check the co-linearity between samples.
correlation <- cor(survival_cancer[,rownames(ph_hypo_table)[ph_hypo_table[,3]>0.05]], method = 'pearson')
library('GGally')
ggpairs(survival_cancer[,rownames(ph_hypo_table)[ph_hypo_table[,3]>0.05]], 
        axisLabels = 'show')+
  theme_bw()+
  theme(panel.background = element_rect(colour = 'black', size=1, fill = 'white'), 
        panel.grid = element_blank())
library('rms')
vif <- rms::vif(multi_variate_cox_2)
#Some people said if the square root of VIF >2, they might be co-linear.
sqrt(vif) < 2

可以得到結(jié)果：變量vif均未達(dá)到共線性的程度帆精。不需進(jìn)一步剔除较屿。

變量vif均未達(dá)到共線性的程度

相關(guān)性矩陣可視化結(jié)果如下：

相關(guān)性矩陣的可視化

接下來(lái)用森林圖可視化這個(gè)模型。

ggforest(model = multi_variate_cox_2, data = survival_cancer, main = 'Hazard ratios of candidate genes', fontsize = 1)

可視化結(jié)果如下卓练，記錄了模型中各因子的危險(xiǎn)比HR隘蝎，置信區(qū)間，模型的AIC值昆庇，concordance index末贾。結(jié)果說(shuō)明模型中OR2W3為獨(dú)立的顯著危險(xiǎn)因素，CHEK2為獨(dú)立的顯著保護(hù)因素整吆，且整個(gè)cox模型也具有顯著意義：

森林圖展示

6. 多方面評(píng)價(jià)cox模型的預(yù)測(cè)能力拱撵。

（1）Concordance index：

模型的一致性系數(shù)反映了其預(yù)測(cè)能力辉川。C-index=0.5則完全隨機(jī)，=1為預(yù)測(cè)與實(shí)際完全一致拴测。>=0.9為高乓旗，0.7-0.9間為中等，0.5-0.7間為低集索。

C_index <- multi_variate_cox_2$concordance['concordance']
if(C_index >= 0.9){
  print('High accuracy')
}else{ 
  if(C_index < 0.9 & C_index >= 0.7){
    print('Medium accuracy')
      }else{
       print('Low accuracy')
      }
}

Anyway, just have fun...

C-index反映模型預(yù)測(cè)能力較差

（2）time-dependent ROC curve:

ROC曲線用于反映預(yù)測(cè)模型的準(zhǔn)確性和精確性屿愚，AUC曲線下面積越大，ROC曲線轉(zhuǎn)折點(diǎn)越靠近左上角务荆，說(shuō)明預(yù)測(cè)能力越好妆距。一般認(rèn)為AUC>=0.7可以夠看。
對(duì)于COX回歸函匕，可以使用time-dependent ROC曲線檢測(cè)模型預(yù)測(cè)特定時(shí)間預(yù)后的能力娱据。

COX建模的意義在于用模型計(jì)算出的risk score預(yù)測(cè)病人的生存狀態(tài)。首先用模型得到的系數(shù)計(jì)算每個(gè)樣本的risk score（gene1β1+gene2β2+gene3*β3）：

依舊是自建函數(shù)盅惜，輸入survival_cancer信息表中剩，candidate genes vector和多因素cox回歸對(duì)象，輸出包含每個(gè)樣本risk score的risk_score_table.

#calculate the risk score of each sample.
riskscore <- function(survival_cancer_df, candidate_genes_for_cox, cox_report) {
  library('dplyr')
  risk_score_table <- survival_cancer_df[,candidate_genes_for_cox]
  for(each_sig_gene in colnames(risk_score_table)){
    risk_score_table$each_sig_gene <- risk_score_table[,each_sig_gene]*(summary(cox_report)$coefficients[each_sig_gene,1])
  }
  risk_score_table <- cbind(risk_score_table, 'total_risk_score'=exp(rowSums(risk_score_table))) %>%
    cbind(survival_cancer_df[,c('TCGA_IDs','overall_survival','censoring_status')])
  risk_score_table <- risk_score_table[,c('TCGA_IDs','overall_survival','censoring_status', candidate_genes_for_cox, 'total_risk_score')]
  risk_score_table
}
candidate_genes_for_cox2 <- c(rownames(ph_hypo_table)[ph_hypo_table[,3]>0.05])
risk_score_table_multi_cox2 <- riskscore(survival_cancer, candidate_genes_for_cox2, multi_variate_cox_2)

risk_score_table的部分展示

這里我們自建函數(shù)抒寂，批量做3-5年多個(gè)時(shí)間點(diǎn)的ROC曲線结啼，比較出AUC最大的時(shí)點(diǎn)。

multi_ROC <- function(time_vector, risk_score_table){
  library('survivalROC')
  single_ROC <- function(single_time){
  for_ROC <- survivalROC(Stime = risk_score_table$overall_survival,
                         status = risk_score_table$censoring_status,
                         marker = risk_score_table$total_risk_score,
                         predict.time = single_time, method = 'KM')
  data.frame('True_positive'=for_ROC$TP, 'False_positive'=for_ROC$FP, 
             'Cut_values'=for_ROC$cut.values, 'Time_point'=rep(single_time, length(for_ROC$TP)),
             'AUC'=rep(for_ROC$AUC, length(for_ROC$TP)))
  }
  multi_ROC_list <- lapply(time_vector, single_ROC)
  do.call(rbind, multi_ROC_list)
}
#We evaluate 11 AUCs between 3-5 years.
for_multi_ROC <- multi_ROC(time_vector = c(365*seq(3,5,0.2)), risk_score_table = risk_score_table_multi_cox2)

輸出的for_multi_ROC dataframe中含有每個(gè)時(shí)間點(diǎn)的True positive和False positive隨不同cut off值變化的趨勢(shì)屈芜。

for_multi_ROC dataframe的部分展示

不同時(shí)間點(diǎn)ROC曲線的可視化：

#visualization of the ROC curves of multiple time points.
pROC<-ggplot(for_multi_ROC, aes(x = False_positive, y = True_positive, label = Cut_values, color = Time_point)) + 
  geom_roc(labels = F, stat = 'identity', n.cuts = 0) + 
  geom_abline(slope = 1, intercept = 0, color = 'red', linetype = 2)+
  theme_bw()+
  theme(panel.background = element_rect(colour = 'black', size=1, fill = 'white'), 
        panel.grid = element_blank())+
  annotate("text",x = 0.75, y = 0.15,
           label = paste("AUC max = ", round(AUC_max, 2), '\n', 'AUC max time = ', AUC_max_time, ' days', sep = ''))
pROC

不同時(shí)間點(diǎn)ROC曲線的比較郊愧，可見粉色曲線具有最高的預(yù)測(cè)效力

AUC max為0.63，可見其實(shí)預(yù)測(cè)能力都很差- -菜雞互啄罷了沸伏。只是案例操作糕珊，無(wú)妨。

(3) Kaplan-Meier生存分析毅糟。

Kaplan-Meier分析用于比較兩組之間的生存狀態(tài)红选。所以我們需要把計(jì)算出來(lái)的risk_score按照一個(gè)截點(diǎn)分為高/低風(fēng)險(xiǎn)組。

ROC曲線的另外一個(gè)功能就是可以用于risk score分組最佳截點(diǎn)的選擇姆另。

首先選取AUC最大的時(shí)間點(diǎn)喇肋，使用這個(gè)時(shí)間點(diǎn)的ROC曲線進(jìn)行Kaplan-Meier分析。

然后截取ROC曲線的轉(zhuǎn)折點(diǎn)迹辐，也就是真陽(yáng)性和假陽(yáng)性差值最大的點(diǎn)蝶防，作為risk score的cut-off值，高于這個(gè)值為高風(fēng)險(xiǎn)組明吩，低于這個(gè)值為低風(fēng)險(xiǎn)組间学，并把信息并入risk_score_table中。

AUC_max <- max(for_multi_ROC$AUC)
#maybe AUCs are identical in different time points. So select the last time point indicating longer survival.
AUC_max_time <- for_multi_ROC$Time_point[which(for_multi_ROC$AUC == AUC_max)]
AUC_max_time <- AUC_max_time[!duplicated(AUC_max_time)]
AUC_max_time <- AUC_max_time[length(AUC_max_time)]
for_multi_ROC$Time_point <- as.factor(for_multi_ROC$Time_point)
#find the optimal cutoff value within the ROC curve of the optimal time point.
optimal_time_ROC_df <- for_multi_ROC[which(for_multi_ROC$Time_point == AUC_max_time),]
cut.off <- optimal_time_ROC_df$Cut_values[which.max(optimal_time_ROC_df$True_positive-optimal_time_ROC_df$False_positive)]
high_low <- (risk_score_table_multi_cox2$total_risk_score > cut.off)
high_low[high_low == TRUE] <- 'high'
high_low[high_low == FALSE] <- 'low'
risk_score_table_multi_cox2 <- cbind(risk_score_table_multi_cox2, high_low)

此時(shí)，已經(jīng)將risk score這個(gè)連續(xù)變量轉(zhuǎn)變?yōu)榱?high'和'low'的二分類變量低葫。使用survival包建立KM分析對(duì)象详羡，然后使用survminer包可視化作KM-plot并進(jìn)行Log-rank分析。

注：為了最佳的可視化效果嘿悬，根據(jù)之前獲得的最佳預(yù)測(cè)時(shí)間（本例為1825天=5年）实柠，將Overall_survival超過(guò)預(yù)測(cè)時(shí)間的樣本編輯一下，將生存時(shí)間改為1825天善涨，生存狀態(tài)改為0窒盐。（可理解為觀察終點(diǎn)為5年，這些樣本在這個(gè)時(shí)點(diǎn)成功完成觀察而未發(fā)生事件）

#KM_plot generation.
library('survminer')
#first edit the status of patients with OS > AUC max time. (censoring status=0 (Alive), OS=365*5 days)
risk_score_table_multi_cox2$censoring_status[which(risk_score_table_multi_cox2$overall_survival > AUC_max_time)] <- 0
risk_score_table_multi_cox2$overall_survival[which(risk_score_table_multi_cox2$overall_survival > AUC_max_time)] <- AUC_max_time
fit_km <- survfit(Surv(overall_survival, censoring_status) ~high_low, data = risk_score_table_multi_cox2)     
ggsurvplot(fit_km, conf.int = F,pval = T,legend.title="total risk score",
           legend.labs=c(paste0('>',as.character(round(cut.off,2))), paste0('<=',as.character(round(cut.off,2)))), risk.table = T, 
           palette = c('red','blue'), surv.median.line = 'hv')

KM-plot說(shuō)明高低風(fēng)險(xiǎn)組的預(yù)后情況有顯著差異

至此時(shí)钢拧，一整套生存分析就完成了蟹漓！
寫這么多，我太難了...