生存分析2-5-1/2

作者: 小胡同学ime | 来源:发表于2022-01-29 16:19 被阅读0次

生存分析2-5-1/2
生存分析（2）
R语言配色方案01-colortools
一文掌握Rmarkdown分析报表
R语言生存分析03-Cox比例风险模型
2-5-1 GridView
TCGA相关分析
六、生存分析
TCGA生存分析③
R语言ggplot2绘制箱线图

title: "生存分析前的数据整理"
output: html_document
editor_options:
chunk_output_type: console

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
knitr::opts_chunk$set(fig.width = 6,fig.height = 6,collapse = TRUE)
knitr::opts_chunk$set(message = FALSE)

生存分析只需要tumor数据，不要normal，将其去掉，新表达矩阵数据命名为exprSet；

clinical信息需要进一步整理，成为生存分析需要的格式，新临床信息数据命名为meta。

由于不同癌症的临床信息表格组织形式不同，这里的代码需要根据实际情况修改。

rm(list=ls())
options(stringsAsFactors = F)
load("TCGA-CHOL_gdc.Rdata")  #输入处理数据
library(stringr)

1.简化临床信息，选出要用的列

tmp = data.frame(colnames(clinical))
clinical = clinical[,c(
  'bcr_patient_barcode',
  'vital_status',    #生
  'days_to_death',    #随访到死亡的日期
  'days_to_last_followup',  #失访时间
  'race_list', 
  'days_to_birth',
  'gender' ,
  'stage_event'
)]   #选出所需关键信息数据

#其实days_to_last_followup应该是找age_at_initial_pathologic_diagnosis，这表格里没有，用days_to_birth计算一下年龄，暂且替代。
dim(clinical)
#rownames(clinical) <- NULL
rownames(clinical) <- clinical$bcr_patient_barcode
clinical[1:4,1:4]

meta = clinical  #改了个名
exprSet=exp[,Group=='tumor']  #生存分析时正常样本不需要分析，只分析tumor样本
exprSet = log2(edgeR::cpm(exprSet)+1)  #取cpm

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

#空着的值改为NA
meta[meta==""]=NA

2.实现表达矩阵与临床信息的匹配

有的病人会有两个或两个以上的肿瘤样本，就有重复。两种可行的办法：

（1）以病人为中心，对表达矩阵的列按照病人ID去重复，每个病人只保留一个样本。（本文档）

（2）以样本为中心，把meta里的病人ID替换成样本ID，这样同一个病人的两个样本就会有两行完全一致的临床信息。（zz.R）

# 以病人为中心，表达矩阵按病人ID去重复
k = !duplicated(str_sub(colnames(exprSet),1,12));table(k)
exprSet = exprSet[,k]

#调整meta的ID顺序与exprSet列名一致
meta=meta[match(str_sub(colnames(exprSet),1,12),meta$ID),]
identical(meta$ID,str_sub(colnames(exprSet),1,12))

3.整理生存分析的输入数据

#1.1由随访时间和死亡时间计算生存时间(月)
meta$time = ifelse(meta$event=="Alive",
                   meta$last_followup,
                   meta$death)
meta$time = as.numeric(meta$time)/30  #按月份计算生存期

#1.2 根据生死定义event，活着是0，死的是1
meta$event=ifelse(meta$event=='Alive',
                  0,
                  1)
table(meta$event)

#1.3 年龄和年龄分组
meta$age=ceiling(abs(as.numeric(meta$age))/365)  #ceiling向上去整数，floor向下去整数

meta$age_group=ifelse(meta$age>median(meta$age,na.rm = T),'older','younger')  #根据中位数进行取整数操作；na.rm = T设置不能让na值影响计算
table(meta$age_group)

#1.4 stage 
library(stringr) 
head(meta$stage)

a = str_extract_all(meta$stage,"I|V");head(a)
b = sapply(a,paste,collapse = "");head(b)
meta$stage = b

#1.5 race
table(meta$race)
save(meta,exprSet,cancer_type,file = paste0(cancer_type,"_sur_model.Rdata"))

title: "生存分析"
output: html_document
editor_options:
chunk_output_type: console

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
knitr::opts_chunk$set(fig.width = 8,fig.height = 6,collapse = TRUE)
knitr::opts_chunk$set(message = FALSE,warning = FALSE)

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)
ggsurvplot(sfit, conf.int=F, pval=TRUE)

image.png

计算年龄上两者有无区别，p=0.66，更高级的图如下，改变带入数据就可用

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

image.png

#性别年龄   #展示拼图
sfit1=survfit(Surv(time, event)~gender, data=meta)  #得到性别数据
sfit2=survfit(Surv(time, event)~age_group, data=meta)  #得到年龄数据
splots <- list()   #说明splots是列表
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)
#拼图函数arrange_ggsurvplots接受对象是列表
dev.off()

image.png

#单个基因  #根据基因表达量的高低进行分组，为患者做生存分析
g = rownames(exprSet)[1] #截选一个基因
meta$gene = ifelse(exprSet[g,]> median(exprSet[g,]),'high','low')   #新增加gene列，如果样本A的g基因大于所有所有样本的中位数，则返回high，否则low
sfit1=survfit(Surv(time, event)~gene, data=meta)
ggsurvplot(sfit1,pval =TRUE, data = meta, risk.table = TRUE)

image.png

#多个基因
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)

image.png

2.logrank批量生存分析

 #非常厉害，对表达矩阵的每个基因进行生存分析取P值
logrankfile = paste0(cancer_type,"_log_rank_p.Rdata")  #使用是改掉"_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) #分别计算出小于0.01及小于0.05的基因，缩小目标基因范围
table(log_rank_p<0.05)

image.png

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) #挑选<0.01的基因矩阵
table(cox_results$p<0.05)#挑选<0.05的基因矩阵

lr = names(log_rank_p)[log_rank_p<0.05];length(lr) 
cox = rownames(cox_results)[cox_results$p<0.05];length(cox)  
length(intersect(lr,cox))   #取交集