TCGA代码分析流程 - 3.2 生存分析：log-rank，Cox回归

生存分析：将事件的结果（终点事件）和出现这一结果所经历的时间结合起来的一种统计分析方法。

生存分析的目的：
1.生存率比较：估计处理组和对照组n年的生存率和中位生存期。
2.生存曲线比较：比较处理组和对照组的生存率是否有差别。
3.影响因素分析：分析变量与生存结局/事件的关系。
4.生存预测：根据变量预测患者n年的生存率。

从生存分析的方法上看，一般可以分为三类：
1.参数法：知道生存时间的分布模型，然后根据数据来估计模型参数，最后以分布模型来计算生存率。如Weibull回归、lognormal回归等。一般不用，因为目前认为我们不知道生存时间符合什么模型。
2.非参数法：不需要生存时间分布，根据样本统计量来估计生存率。常见方法Kaplan-Meier分析（乘积极限法）、寿命法（Life Tables）。而对于Kaplan-Meier法来说，其中的p值我们常用log-rank检验（对数秩检验）和Wilcoxon检验去求。
3.半参数法：也不需要生存时间的分布，但最终是通过模型来评估影响生存率的因素。最为常见的是Cox比例风险回归模型。

从生存分析的因素上看，有单因素分析和多因素分析：
1.单因素生存分析：描述生存过程。常用的单因素分析方法有乘积极限法（Kaplan-Meier分析）和寿命表法（Life Tables）。
2.多因素生存分析：分析生存过程的影响因素。常用的多因素生存分析方法有Cox比例风险回归模型。

在简单生存分析中，由于仅考虑单个影响因素（且为分类型或顺序型变量），故采取的是直接绘制生存曲线并作Log-Rank检验来判断影响因素和生存概率的方法。
而在Cox回归分析中，需要同时考虑多个影响因素（可为分类/顺序型变量，也可为数值型变量），故而绘制生存曲线的方式显然不合适，此时就需通过建模的方法来进一步分析。

单个变量的Cox回归和K-M法结果不一致时，此时我们还是应该选择Cox的结果，因为参数检验效力高于非参数检验。
多个变量的Cox回归中单个变量的显著性与K-M法不一致，此时我们应该选择cox的结果，因为K-M发只考虑单个变量，而Cox考虑了多个变量。

描述性统计分析：Kaplan-Meier分析：非参数估计，不要求总体的分布形式，直观地表现出两组或多组的生存率或死亡率，适合在文章中展示。
统计推断：组间比较：log-rank检验，两组之间生存率的差异是否显著。

Kaplan-Meier analysis 单因素生存分析绘图 KM-Plot

rm(list = ls())
library(survival)
library(survminer)
library(pbapply)
setwd("D:/BioInformationAnalyze/TCGAdata/CRC")
load("TCGA-COAD_04.lncRNA_SurvData.Rdata")  # 生存分析输入数据，包含tpm形式的表达矩阵和临床信息

需要的文件有基因表达文件（exprSet）和生存数据文件（meta），文件格式如下：

年龄

sfit <- survfit(Surv(time, event) ~ age_group, data = meta)
              # Surv()函数创建生存数据对象（主要输入生存时间和状态逻辑值）
      # survfit()函数对生存数据对象拟合生存函数，创建KM(Kaplan-Meier)生存曲线
                   # time表示总生存期
                         # event表示终点事件，1死亡，0存活。
                                # ~age_group表示以年龄分组绘图，注意此方法仅适用二分类变量进行分组。如果不想分组，只看总体，~1即可
# 如果想看整体的sfit结果，可以summary(sfit) 或 summary(sfit)$table，
# 想看每个时间点的更好的方式则是:surv_summary(sfit)：
                               # time:曲线上的时间点。
                               # n.risk:时间t时面临风险的受试者数量.
                               # n.event:在时间t发生的事件数。
                               # n.censor:在时间t时退出风险集且没有事件的受审查对象的数量。
                               # surv:生存概率的估计。
                               # std.err:生存的标准误差。
                               # upper:置信区间的上限。
                               # lower:置信区间的下限。
                               # strata:表示生存曲线的分层。如果分层不为 NULL，则结果中存在多条曲线。分层（因子）的水平是曲线的标签。
# 根据上述sfit结果，我们可以将KM生存数据进行可视化，主要使用Survminer包的ggsurvplot()函数:
ggsurvplot(sfit, conf.int = FALSE, pval = TRUE)
# 通过设置参数，可在生存曲线基础上，再绘制两张辅助图:
ggsurvplot(sfit, 
           palette = c("#E7B800", "#2E9FDF"),
           pval = TRUE, conf.int = TRUE,  # 显示p值和置信区间
           xlab = "Time in months",  # x轴标签
           linetype = "strata",  # 线条类型
         # fun = "event",  # 绘制cumulative events图，即事件发生累计概率图
           ggtheme = theme_light(),  # 更改 ggplot2 主题
           surv.median.line = "hv",  # 指定中位生存期
           risk.table = TRUE,  # Number at risk图，表示在该生存时间长度内，还存活的case；较直接得反映出两组的一个差异情况。
         # risk.table.col = "strata",  # 按分组改变risk table的颜色
           ncensor.plot = TRUE)  # Number of censoring图，即该时间段内的删失值。

性别、年龄

sfit1 <- survfit(Surv(time, event) ~ gender, data = meta)
sfit2 <- survfit(Surv(time, event) ~ age_group, data = meta)
splots <- list()
splots[[1]] <- ggsurvplot(sfit1, pval = TRUE, data = meta, risk.table = TRUE)
splots[[2]] <- ggsurvplot(sfit2, pval = TRUE, data = meta, risk.table = TRUE)
arrange_ggsurvplots(splots, print = TRUE, ncol = 2, nrow = 1, risk.table.height = 0.4)  # 内置拼图函数

单个基因

设定某一基因的表达基准，将临床样本分为高表达high与低表达low两组，再进行生存分析。这里用表达矩阵的tpm值进行计算。

g <- rownames(exprSet)[2]  # 挑选基因，可修改成自己感兴趣的基因
meta$gene <- ifelse(as.integer(exprSet[g,]) > median(as.integer(exprSet[g,])), "high", "low")
sfit1 <- survfit(Surv(time, event) ~ gene, data = meta)
ggsurvplot(sfit1, pval = TRUE, data = meta, risk.table = TRUE)

多个基因

gs <- rownames(exprSet)[1:4]  # 挑选基因，可修改成自己感兴趣的基因
splots <- pblapply(gs, function(g){
  meta$gene <- ifelse(as.integer(exprSet[g,]) > median(as.integer(exprSet[g,])), "high", "low")
  sfit1 <- survfit(Surv(time, event) ~ gene, data = meta)
  ggsurvplot(sfit1, pval = TRUE, data = meta, risk.table = TRUE)
}) 
arrange_ggsurvplots(splots, print = TRUE,  
                    ncol = 2, nrow = 2, risk.table.height = 0.4)  # 内置拼图函数

log-rank批量生存分析

Log-Rank test：对不同组的生存率进行假设检验，是无参数检验，近似于卡方检验，零假设是组间没有差异。

由于log-rank检验需要对自变量进行分组，无法使用连续型变量，而有时候在不知道连续性变量的最佳截断值或者自变量很多的时候，只能按中位数进行分组，所以结果并没有Cox回归用连续型变量得出的结果好。所以一般先用Cox或Lasso将变量降维后，再找最佳截断值来画生存曲线。

Log-Rank检验属单因素分析方法，应用条件是除比较因素外，影响生存率的其他因素组件均衡可比，否则应矫正各因素的影响，可采用Cox比例风险回归模型。
用apply()函数批量计算所有表达矩阵中基因的生存分析p值，从而挑选出“显著基因”。这里用表达矩阵的tpm值进行计算

mySurv <- with(meta, Surv(time, event == "1"))  # 创建生存对象，event == 1表示“1为结局时间发生“
log_rank_p <- pbapply(exprSet, 1, function(gene){
  # gene <- exprSet[1,]
  meta$group <- ifelse(gene > median(gene), "high", "low") 
  data.survdiff <- survdiff(mySurv ~ group, data = meta)
    # survdiff()用于检验两条或多条生存曲线之间是否存在差异，
      # 或者针对已知替代项检验单条曲线，也适用于分组资料以及多组间生存曲线的比较；
  p.val <- 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)  # 计算log-rank检验的p值
  return(p.val)
})
log_rank_p <- sort(log_rank_p)  # 取出每一个基因生存分析的P值，形成表
table(log_rank_p < 0.01)  # 哪些是P<0.01的
table(log_rank_p < 0.05)  # 哪些是P<0.05的
 
log_rank <- names(log_rank_p)[log_rank_p < 0.05]  # log-rank检验筛选的差异基因

log-rank批量生存分析结果如下：

KM曲线的优势：可以直接从图中看出组间差异。
KM曲线的局限性：1.无法得知暴露因素对结局的影响是非暴露组多少倍；2.没有控制其他变量。

要解决以上局限性，需引入Cox比例风险回归模型，进行多变量Cox回归分析，得到风险差异的具体估计值HR。
Cox比例风险回归分析的优势在于可以分析分类变量与数值变量，并且将生存分析的范围由单因素拓展到多因素的分析。
Cox比例风险回归模型可以同时分析几个因素的生存。另外，统计模型提供每个因素的效果大小。

Cox比例风险回归

Cox比例风险回归模型可以分析多个因素对生存事件的影响，而且允许有删失。是生存分析中重要的多因素分析方法。
Cox回归本质上是一种回归模型，它没有直接使用生存时间，而是使用了风险比（hazard ratio）作为因变量。
该模型不用于估计生存率，而是用于因素分析，也就是找到某一个危险因素对结局事件发生的贡献度。

Cox回归的重要统计指标：风险比（hazard ratio，HR），表示变量X变化一个单位引起风险的变化。
当HR > 1时，说明研究对象是一个危险因素。
当HR < 1时，说明研究对象是一个保护因素。
当HR = 1时，说明研究对象对生存时间不起作用。

Cox批量单因素生存分析 Univariate Cox Analysis

mySurv <- with(meta, Surv(time, event == "1"))  # 创建生存对象，event = 1表示“1为结局事件发生”
unicox_results <- pbapply(exprSet, 1, function(gene){  # 对表达矩阵exprSet中的每一行（基因）执行函数function
  # gene = as.numeric(exprSet[1,]) # 举个例子
  group <- ifelse(gene > median(gene), "high","low")  # 标记不同样本中该基因的上下调，将表达量转换为二分类变量
  survival_dat <- data.frame(group = group,
                             stage = meta$stage,
                             age = meta$age,
                             gender = meta$gender,
                             gene = gene,  # 基因表达量
                             stringsAsFactors = F)  # 选择用于生存分析的数据
  # 拟合Cox比例风险回归模型 m ：
  m <- coxph(mySurv ~ gender + age + stage + group, data = survival_dat)  # 表达量用二分类变量group
  # m <- coxph(mySurv ~ group, data = survival_dat)  # 表达量用二分类变量group
  # m <- coxph(mySurv ~ gene, data =  survival_dat)  # 表达量也可直接使用连续型变量gene
  # sm <- summary(m)  # 看拟合出的Cox比例风险回归模型的参数
    # n：总人数。
    # number of events：发生结局事件的人数，即死亡人数
    # coef(coefficients,回归系数)：结果为正(负)，表示男性（genderMALE）死亡风险高（低）于女性，或者说有较低（高）的生存概率。
    # exp(coef)：也称为风险比（HR）。HR = 1: 无影响，HR < 1: 降低风险，HR > 1: 增加风险。
    # se(coef)：coef的标准误
    # z：Wald检验的统计量，z = coef/se(coef)，是否具有统计学意义看p值
    # Pr(>|z|)：p值
    # low.95与upper.95表示HR值的95%置信度区间CI，如果95%CI跨越了1，一般就不认为该因素对生存有显著影响。
  beta <- coef(m)
  se <- sqrt(diag(vcov(m)))
  tmp <- round(cbind(coef = beta,  # 模型的偏回归系数β
                     se = se,  # 偏回归系数β的标准误
                     z = beta/se,  # Wald检验的统计量
                     p = 1 - pchisq((beta/se)^2, 1),  # p值，差异的统计学意义
                     HR = exp(beta),  # 风险值，计算公式为exp(coef)
                     HRCILL = exp(beta - qnorm(.975, 0, 1) * se),  # HR的95%置信区间低点
                     HRCIUL = exp(beta + qnorm(.975, 0, 1) * se)), 3)  # HR的95%置信区间高点
  return(tmp["grouplow",])  # 分类变量
  # return(tmp["gene",])  # 连续型变量
})
unicox_results <- as.data.frame(t(unicox_results))
table(unicox_results[,4] < 0.01)  # p值<0.01的有多少
table(unicox_results[,4] < 0.05)  # p值<0.05的有多少
 
unicox <- rownames(unicox_results)[unicox_results[,4] < 0.05]  # Cox单因素分析筛选的差异基因

Cox批量单因素生存分析结果如下：

Cox多因素生存分析 Multivariable Cox Analysis

mySurv <- with(meta, Surv(time, event == "1")) # 创建生存对象
survival_dat <- as.data.frame(t(exprSet))
# 拟合Cox比例风险回归模型 m ：
m <- coxph(mySurv ~ LINC02185 + AL391845.2 + LINC00858, data = survival_dat) # 表达量用连续型变量
# sm <- summary(m) # 看拟合出的Cox比例风险回归模型的参数
  # n：总人数。
  # number of events：发生结局事件的人数，即死亡人数
  # coef(coefficients,回归系数)：结果为正(负)，表示男性（genderMALE）死亡风险高（低）于女性，或者说有较低（高）的生存概率。
  # exp(coef)：也称为风险比（HR）。HR = 1: 无影响，HR < 1: 降低风险，HR > 1: 增加风险。
  # se(coef)：coef的标准误
  # z：Wald检验的统计量，z = coef/se(coef)，是否具有统计学意义看p值
  # Pr(>|z|)：p值
  # low.95与upper.95表示HR值的95%置信度区间CI，如果95%CI跨越了1，一般就不认为该因素对生存有显著影响。
beta <- coef(m)
se <- sqrt(diag(vcov(m)))
multicox_results <- round(cbind(coef = beta,  # 模型的偏回归系数β
                                se = se,  # 偏回归系数β的标准误
                                z = beta/se,  # Wald检验的统计量
                                p = 1 - pchisq((beta/se)^2, 1),  # p值，差异的统计学意义
                                HR = exp(beta),  # 风险值，计算公式为exp(coef)
                                HRCILL = exp(beta - qnorm(.975, 0, 1) * se),  # HR的95%置信区间低点
                                HRCIUL = exp(beta + qnorm(.975, 0, 1) * se)), 3)  # HR的95%置信区间高点
table(multicox_results[,4] < 0.01)  # p值<0.01的有多少
table(multicox_results[,4] < 0.05)  # p值<0.05的有多少
 
multicox <- rownames(multicox_results)[multicox_results[,4] < 0.05]  # Cox多因素分析筛选的差异基因