1. Kaplan-Meier生存分析

Kaplan Meier是一種單因素生存分析数苫，它可用于研究1個因素對于生存時間的影響笋颤。
參考：生存分析之Kaplan-Meier曲線都告訴我們什么

思路

一般來說逼龟，我們做生存分析亿胸，會有（P<0.05）和（P>0.05）兩種結(jié)果。KM plot在生物醫(yī)學(xué)中很常見蟀苛，主要用來做預(yù)后分析益咬，比如可以根據(jù)表達量把病人分成兩組，然后比較哪組病人預(yù)后好帜平，進而可以得出基因表達量高低與病人預(yù)后好壞相關(guān)性的結(jié)論幽告。(但凡是能把病人分成兩個組的梅鹦，都可以拿來做生存分析，根據(jù)基因表達量评腺，人種帘瞭，性別，stage等等來分蒿讥，都是可以的。)
畫KM plot時抛腕，有時候會比較糾結(jié)怎樣對病人進行分組芋绸，如何來設(shè)置分組的cutoff。一般來說常見的幾種設(shè)置cutoff值得思路如下：
1：大多數(shù)情況下担敌，根據(jù)表達量從低到高對樣本進行排序摔敛，取前50%為低表達，后50%為高表達全封，然后畫KM plot马昙。
2：還有一些文章也會將樣本表達量均分為三組或者四組。
3：一些文章也會選一些其它的cutoff刹悴，比如前1/3和后2/3行楞，前25%和后25%（中間50%的數(shù)據(jù)去掉）。
參考：https://zhuanlan.zhihu.com/p/86307194#:~:text=KM%20plot

2. 結(jié)果解讀

左圖：KM-plot的橫坐標(biāo)是生存時間（可以是OS土匀，也可以是PFS等）子房，縱坐標(biāo)是生存率。起點是隨訪開始的時間就轧，曲線下降代表患者死亡证杭，曲線上的+號表示刪失（仍然存活的病人的最后隨訪時間）。p值是log-rank test的p值妒御，衡量的繪制的曲線之間生存率有沒有顯著差別解愤。
右上：左圖的美化版，曲線的陰影表示置信區(qū)間乎莉。
右中：每個時間節(jié)點有多少人在隨訪中
右下：刪失表格送讲，右上圖中每出現(xiàn)一個+，也就是刪失梦鉴，右下對應(yīng)的位置就會有一個條圖李茫。縱軸是在該時間點刪失的數(shù)目肥橙。

3. 數(shù)據(jù)整理??

生存分析需要的輸入數(shù)據(jù)有兩個：表達矩陣數(shù)據(jù)和臨床信息數(shù)據(jù)
生存分析只需要tumor數(shù)據(jù)魄宏，不要normal，將其去掉存筏，新表達矩陣數(shù)據(jù)命名為exprSet宠互；
clinical信息需要進一步整理味榛，成為生存分析需要的格式，新臨床信息數(shù)據(jù)命名為meta予跌。
由于不同癌癥的臨床信息表格組織形式不同搏色，這里的代碼需要根據(jù)實際情況修改。

rm(list=ls())
options(stringsAsFactors = F)
load("TCGA-CHOL_gdc.Rdata")
load("TCGA-CHOL_symbol_exp.Rdata")
library(stringr)

3.1 整理表達矩陣與臨床信息

• 整理表達矩陣

exprSet=exp[,Group=='tumor'] #只留下腫瘤樣本

#前面的過濾標(biāo)準(zhǔn)比較寬松券册，這里可以再次過濾
k = apply(exprSet,1, function(x){sum(x==0)<=0.25*ncol(exprSet)});table(k)
# k
# FALSE  TRUE 
#  4631 25521

exprSet = log2(edgeR::cpm(exprSet[k,])+1) #標(biāo)準(zhǔn)化
dim(exprSet) #25521個基因频轿，36個樣本
# [1] 25521    36

• 整理臨床信息

tmp = data.frame(colnames(clinical)) # 方便搜索列名
meta = clinical[,c(
  'bcr_patient_barcode',
  'vital_status',
  'days_to_death',
  'days_to_last_followup',
  'race_list',
  'days_to_birth',
  'gender' ,
  'stage_event'
)] # 挑出其中有用的列，重新生成數(shù)據(jù)框
dim(meta)
# [1] 48  8
head(meta)
#   bcr_patient_barcode vital_status days_to_death days_to_last_followup race_list days_to_birth gender        stage_event
# 1        TCGA-W5-AA2Q        Alive                                  50     WHITE        -25069   MALE 6thStage IIT2bN0M0
# 2        TCGA-W5-AA2O         Dead           640                           WHITE        -21118   MALE   6thStage IT1N0M0
# 3        TCGA-W5-AA39         Dead           170                           WHITE        -29594   MALE  7thStage IIT2N0M0
# 4        TCGA-W6-AA0S        Alive                                 352     WHITE        -16951 FEMALE   7thStage IT1N0MX
# 5        TCGA-3X-AAVE        Alive                                 203     ASIAN        -21943   MALE  7thStage IIT2N0M0
# 6        TCGA-W5-AA2I         Dead          1939                           WHITE        -24388   MALE   6thStage IT1N0M0

#其實days_to_last_followup應(yīng)該是找age_at_initial_pathologic_diagnosis烁焙，這表格里沒有航邢，用days_to_birth計算一下年齡，暫且替代骄蝇。

rownames(meta) <- meta$bcr_patient_barcode
meta[1:4,1:4]
#>              bcr_patient_barcode vital_status days_to_death
#> TCGA-W5-AA2Q        TCGA-W5-AA2Q        Alive              
#> TCGA-W5-AA2O        TCGA-W5-AA2O         Dead           640
#> TCGA-W5-AA39        TCGA-W5-AA39         Dead           170
#> TCGA-W6-AA0S        TCGA-W6-AA0S        Alive              
#>              days_to_last_followup
#> TCGA-W5-AA2Q                    50
#> TCGA-W5-AA2O                      
#> TCGA-W5-AA39                      
#> TCGA-W6-AA0S                   352

#簡化meta的列名
colnames(meta)=c('ID','event','death','last_followup','race','age','gender','stage')

#空著的值改為NA
meta[meta==""]=NA

3.2 實現(xiàn)表達矩陣與臨床信息的匹配

有的病人會有兩個或兩個以上的腫瘤樣本膳殷，就有重復(fù)。兩種可行的辦法：
（1）以病人為中心九火，對表達矩陣的列按照病人ID去重復(fù)赚窃，每個病人只保留一個樣本。（本文檔）
（2）以樣本為中心岔激，把meta里的病人ID替換成樣本ID勒极，這樣同一個病人的兩個樣本就會有兩行完全一致的臨床信息。（zz.R）

# 以病人為中心鹦倚，表達矩陣按病人ID去重復(fù)
k = !duplicated(str_sub(colnames(exprSet),1,12));table(k)
# k
# TRUE 
#   36
exprSet = exprSet[,k]

#調(diào)整meta的ID順序與exprSet列名一致
meta=meta[match(str_sub(colnames(exprSet),1,12),meta$ID),]
identical(meta$ID,str_sub(colnames(exprSet),1,12))
# [1] TRUE

3.3 整理生存分析的輸入數(shù)據(jù)

#1.1由隨訪時間和死亡時間計算生存時間(月)
table(meta$event) 
#
# Alive  Dead 
#    20    16
meta$time = ifelse(meta$event=="Alive", meta$last_followup, meta$death)
meta$time = as.numeric(meta$time)/30

#1.2 根據(jù)生死定義event河质，活著是0，死的是1
meta$event=ifelse(meta$event=='Alive', 0, 1)
table(meta$event)
# 
#  0  1 
# 20 16

#1.3 年齡和年齡分組
meta$age=ceiling(abs(as.numeric(meta$age))/365)

meta$age_group=ifelse(meta$age>median(meta$age,na.rm = T),'older','younger')
table(meta$age_group)
# 
#   older younger 
#      18      18

#1.4 stage 
library(stringr) 
head(meta$stage)
# [1] "6thStage IVT3N0M1"  "7thStage IIIT3N0M0" "7thStage IIT2bNXMX"
# [4] "7thStage IT1N0M0"   "7thStage IIT2N0M0"  "7thStage IT1N0M0"

a = str_extract_all(meta$stage,"I|V");head(a) #提取分期信息
b = sapply(a,paste,collapse = "");head(b)
#> [1] "IV"  "III" "II"  "I"   "II"  "I"
meta$stage = b

# 去掉生存信息不全或者生存時間小于0.1(月)的病人震叙，樣本納排標(biāo)準(zhǔn)不唯一松嘶，且差別很大
k1 = meta$time>=0.1;table(k1)
# k1
# FALSE  TRUE 
#     1    35
k2 = !(is.na(meta$time)|is.na(meta$event));table(k2)
# k2
# TRUE 
#   36
exprSet = exprSet[,k1&k2]
meta = meta[k1&k2,]

save(meta,exprSet,proj,file = paste0(proj,"_sur_model.Rdata"))

整理好的表達矩陣

整理好的臨床信息矩陣

4. 繪制生存曲線

# 格式：
sfit <- survfit(Surv(time,event)~group,data= meta)
# time指總生存期
# event是終點事件
# group是分組鸠踪，是臨床信息中的一列
# data是臨床信息表格

rm(list = ls())
load("TCGA-CHOL_sur_model.Rdata")
library(survival)
library(survminer)
#年齡
sfit <- survfit(Surv(time, event)~age_group, data=meta)
ggsurvplot(sfit, conf.int=F, pval=TRUE)

ggsurvplot(sfit,palette = c("#E7B800", "#2E9FDF"),
           risk.table =TRUE,pval =TRUE,
           conf.int =TRUE,xlab ="Time in months", 
           ggtheme =theme_light(), 
           ncensor.plot = TRUE)

#性別年齡
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)

dev.off()
#> null device 
#>           1
#單個基因
g = rownames(exprSet)[1]
meta$gene = ifelse(exprSet[g,]> median(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 <- lapply(gs, function(g){
  meta$gene=ifelse(exprSet[g,] > median(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)

5. 批量生存分析

參考：https://www.sohu.com/a/280105039_743978

5.1 logrank批量生存分析

logrankfile = paste0(proj,"_log_rank_p.Rdata")
if(!file.exists(logrankfile)){
  log_rank_p <- apply(exprSet , 1 , function(gene){
    # gene=exprSet[1,]
    meta$group=ifelse(gene>median(gene),'high','low')  
    data.survdiff=survdiff(Surv(time, event)~group,data=meta)
    p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)
    return(p.val)
  })
  log_rank_p=sort(log_rank_p)
  save(log_rank_p,file = logrankfile)
}
load(logrankfile)
table(log_rank_p<0.01) 
#> 
#> FALSE  TRUE 
#> 25430    91
table(log_rank_p<0.05) 
#> 
#> FALSE  TRUE 
#> 24859   662

# 挑一個p值小的基因來做KM_plot
g = names(log_rank_p)[1]
meta$gene = ifelse(exprSet[g,]> median(exprSet[g,]),'high','low')
sfit1=survfit(Surv(time, event)~gene, data=meta)
ggsurvplot(sfit1,pval =TRUE, data = meta, risk.table = TRUE)

5.2 cox批量生存分析

Cox 回歸本質(zhì)上是一種回歸模型 馏予，它沒有直接使用生存時間雳刺，而是使用了風(fēng)險比（ hazard ratio ）作為因變量，該模型不用于估計生存率划址，而是用于因素分析也就是 找到某一個危險因素對結(jié)局事件發(fā)生的貢獻度扔嵌。

Cox 回歸的重要統(tǒng)計指標(biāo)： 風(fēng)險比（ hazard ratio）
? 當(dāng) HR>1 時，說明研究對象是一個危險因素夺颤。
? 當(dāng) HR<1 時痢缎，說明研究對象是一個保護因素。
? 當(dāng) HR=1 時世澜，說明研究對象對生存時間不起作用独旷。

coxfile = paste0(proj,"_cox.Rdata")
if(!file.exists(coxfile)){
  cox_results <-apply(exprSet , 1 , function(gene){
  #gene= exprSet[1,]
  meta$gene = gene
  #可直接使用連續(xù)型變量
  m = coxph(Surv(time, event) ~ gene, data =  meta)
  #也可使用二分類變量
  #meta$group=ifelse(gene>median(gene),'high','low') 
  #m=coxph(Surv(time, event) ~ group, data =  meta)
  beta <- coef(m)
  se <- sqrt(diag(vcov(m)))
  HR <- exp(beta)
  HRse <- HR * se
  
  #summary(m)
  tmp <- round(cbind(coef = beta, 
                     se = se, z = beta/se, 
                     p = 1 - pchisq((beta/se)^2, 1),
                     HR = HR, HRse = HRse,
                     HRz = (HR - 1) / HRse, 
                     HRp = 1 - pchisq(((HR - 1)/HRse)^2, 1),
                     HRCILL = exp(beta - qnorm(.975, 0, 1) * se),
                     HRCIUL = exp(beta + qnorm(.975, 0, 1) * se)), 3)
  
  return(tmp['gene',]) 
  #return(tmp['grouplow',])#二分類變量
})
  cox_results=as.data.frame(t(cox_results))
  save(cox_results,file = coxfile)
}
load(coxfile)
table(cox_results$p<0.01)
#> 
#> FALSE  TRUE 
#> 25439    82
table(cox_results$p<0.05)
#> 
#> FALSE  TRUE 
#> 24979   542

lr = names(log_rank_p)[log_rank_p<0.05];length(lr)
#> [1] 662
cox = rownames(cox_results)[cox_results$p<0.05];length(cox)
#> [1] 542
length(intersect(lr,cox))
#> [1] 132

代碼來自2021生信技能樹數(shù)據(jù)挖掘課