六、生存分析

1.KM-plot

rm(list = ls())
load("TCGA-CHOL_sur_model.Rdata")
library(survival)
library(survminer)
#年龄
sfit <- survfit(Surv(time, event)~age_group, data=meta)#更换age_group可以换分组画图
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)

2.logrank批量生存分析

logrankfile = paste0(cancer_type,"_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)
  })#对每个基因都计算p值
  log_rank_p=sort(log_ra        nk_p)
  save(log_rank_p,file = logrankfile)
}
load(logrankfile)
table(log_rank_p<0.01) 
#> 
#> FALSE  TRUE 
#> 30213   135
table(log_rank_p<0.05) 
#> 
#> FALSE  TRUE 
#> 29445   903

3.cox批量生存分析

coxfile = paste0(cancer_type,"_cox.Rdata")
if(!file.exists(coxfile)){
  cox_results <-apply(exprSet , 1 , function(gene){
  #gene= exprSet[1,]
  meta$gene = gene
  #可直接使用连续型变量
  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 
#> 30231   117
table(cox_results$p<0.05)
#> 
#> FALSE  TRUE 
#> 29572   776

lr = names(log_rank_p)[log_rank_p<0.05];length(lr)#
#> [1] 903
cox = rownames(cox_results)[cox_results$p<0.05];length(cox)
#> [1] 776
length(intersect(lr,cox))
#> [1] 189

*生信技能树课程笔记