生存分析：将事件的结果（终点事件）和出现这一结果所经历的时间结合起来的一种统计分析方法。
生存分析的目的：
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多因素分析筛选的差异基因取交集
length(intersect(log_rank, intersect(unicox, multicox))) # 可以取交集用于后续分析，也可任选lr或cox之一进行分析。
