生存分析2-5-1/2
2022-01-29 本文已影响0人
小胡同学ime
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)) #取交集
截屏2021-09-23 下午6.59.14.png
image.png
