本文学习生存分析代码，KM-plot可视化、log-rank test得到p值、批量单因素cox得到HR值，本篇代码几乎不用改动
1.准备输入数据
rm(list = ls())
proj = "TCGA-KIRC"
load(paste0(proj,"_sur_model.Rdata"))
ls()
exprSet[1:4,1:4]
str(meta)
2.KM-plot
# 简单版本和进阶版本
library(survival)
library(survminer)
sfit <- survfit(Surv(time, event)~gender, data=meta)
ggsurvplot(sfit,pval=TRUE)
ggsurvplot(sfit,
palette = "jco",
risk.table =TRUE,
pval =TRUE,
conf.int =TRUE)
注意上面代码中的gender，这个位置必须是非连续型的比如stage，而不能是连续型的age否则如图1所示
图1
实在想用age这种连续型信息，该怎么作KM分析？例如年龄、基因？那就以某个值为依据，划分为两组
连续型数据的离散化
# 年龄
group = ifelse(meta$age>median(meta$age,na.rm = T),"older","younger")#要去除NA值，否则算不出中位数是多少
table(group)
sfit=survfit(Surv(time, event)~group, data=meta)# group不在meta里面，但也能直接用
ggsurvplot(sfit,pval =TRUE, data = meta, risk.table = TRUE)
# 基因
g = rownames(exprSet)[1];g
meta$gene = ifelse(exprSet[g,]> median(exprSet[g,]),'high','low')
table(meta$gene)
sfit=survfit(Surv(time, event)~gene, data=meta)
ggsurvplot(sfit,pval =TRUE, data = meta, risk.table = TRUE)
3.log-rank test
KM的p值是log-rank test得出的，可以批量操作
source("KM_cox_function.R")#计算的过程在这个函数里面，已经简化了
logrankfile = paste0(proj,"_log_rank_p.Rdata")
if(!file.exists(logrankfile)){
log_rank_p <- apply(exprSet , 1 , geneKM)
log_rank_p=sort(log_rank_p)
save(log_rank_p,file = logrankfile)
}#计算耗时较长，如果计算过一遍了，即已经有“TCGA-KIRC_log_rank_p.Rdata”这个文件了，就不会重新运行了
load(logrankfile)#计算的结果即log_rank_p
table(log_rank_p<0.01)
table(log_rank_p<0.05)
图2
4.Cox回归
Cox回归本质上是一种回归模型，它没有直接使用生存时间，而是使用了风险比（hazard ratio）作为因变量，该模型不用于估计生存率，而是用于因素分析，也就是找到某一个因素对结局事件发生的贡献度
比较重要的统计指标，除了上面提到的p值，还有风险比
图3
详细介绍可查看该文章https://www.sohu.com/a/280105039_743978
离散型和连续型都可以
KM-plot需要分组，而cox离散型和连续型数据都可以
但更倾向于用连续型数据，因为这样会分辨率更高，颗粒度更大。而且如果简单的分为两个组，无法保证被分到一组的就没有区别都没有。而分完组以后，即使有区别，其也无法保留在log_rank_test里了
