Log-rank批量生存分析
Cox 批量生存分析
作图
1、Kaplan-Meier生存分析——KM plot
一句代码即可实现
sfit <- survfit(Surv(time, event)~gender, data=meta)
image.png
meta——临床信息表格
image.png
通过event和time就可以做出一张KM Plot
event里面，1代表已经死了，0代表还活着
time里面以月为单位
meta表格里面的数据只要能分组就可以画出KM 图
用绿色框框出来的都是可以显示在图中的，示例如下：
image.png
2、COX回归——风险评估
Cox回归本质上是一种回归模型，它没有直接使用生存时间，而是使用了风险
比（hazard ratio）作为因变量，该模型不用于估计生存率，而是用于因素分
析，也就是找到某一个危险因素对结局事件发生的贡献度。
• Cox回归的重要统计指标：风险比（hazard ratio）
• 当HR>1时，说明研究对象是一个危险因素。
• 当HR<1时，说明研究对象是一个保护因素。
• 当HR=1时，说明研究对象对生存时间不起作用。
3、代码实现
生存分析只需要tumor数据，不要normal，将其去掉，新表达矩阵数据命名为exprSet；
clinical信息需要进一步整理，成为生存分析需要的格式，新临床信息数据命名为meta。
由于不同癌症的临床信息表格组织形式不同，这里的代码需要根据实际情况修改。
rm(list=ls())
options(stringsAsFactors = F)
load("TCGA-CHOLgdc.Rdata") #加载之前处理好的数据
library(stringr)
#clinical通常有几十列，选出其中有用的几列。
tmp = data.frame(colnames(clinical)) #变成data.frame有便于复制自己想要的列名
clinical = clinical[,c(
'bcr_patient_barcode', #barcode是病人的ID
'vital_status', #生存状态
'days_to_death', #死亡日期
'days_to_last_followup', #最后随访的日期
'race_list', #人种
'days_to_birth', #年龄：这里其实不太对，应该换成age_at_initial_pathologic_diagnosis，但是他没有，只能用这个
'gender' , #性别
'stage_event' #生存状态
)]
#其实days_to_last_followup应该是找age_at_initial_pathologic_diagnosis，这表格里没有，用days_to_birth计算一下年龄，暂且替代。
dim(clinical)
#48 8
#rownames(clinical) <- NULL
rownames(clinical) <- clinical$bcr_patient_barcode #让病人ID等于行名
clinical[1:4,1:4]
meta = clinical
exprSet=exp[,group_list=='tumor'] #按照选列把肿瘤的样本留下，把正常的样本丢掉
#简化meta的列名
colnames(meta)=c('ID','event','death','last_followup','race','age','gender','stage')
#调整meta的ID顺序与exprSet列名一致
meta=meta[match(str_sub(colnames(exprSet),1,12),meta$ID),] #让病人ID和样本ID一一匹配；1,12代表样本的前12位，根据列名的前12位来调整meta中ID行的顺序
all(meta$ID==str_sub(colnames(exprSet),1,12)) #检查一下是否相等，返回值是TRUE就没问题
图中有123，共3个问题需要解决
1、这里死亡日期和随访日期应该加起来才可以，并且这两列有许多空值，都是NA，NA是没有办法参加计算的，会传染，就是谁加上NA最后就会等于NA，所以要把所有的NA变成0才能参与计算
2、年龄这一列有负值，年龄不可能是负的所以要变成正的
3、stage这一列我们只需要分期信息，也就是罗马数字II IV这样的数，其他的都要去掉
#整理生存分析的输入数据----
#1.1由随访时间和死亡时间计算生存时间(月)
is.empty.chr = function(x){
ifelse(stringr::str_length(x)==0,T,F) #str_length(x)==0这一步就是在判断长度=0的就是空字符串，注意NA的长度不是0，是1，而空字符串的长度才是0
}
is.empty.chr(meta[1,4])
meta[,3][is.empty.chr(meta[,3])]=0 #把第3列空字符串改为0
meta[,4][is.empty.chr(meta[,4])]=0 #把第4列空字符串改为0
meta$time=(as.numeric(meta[,3])+as.numeric(meta[,4]))/30 #这两列加起来是天数，然后除以30就是月数
#1.2 根据生死定义event，活着是0，死的是1
meta$event=ifelse(meta$event=='Alive',0,1) #这里的Alive是因为表格中也是这样写的，如果别的数据中全部是大写也要换掉！
#1.3 年龄和年龄分组
meta$age=ceiling(abs(as.numeric(meta$age))/365) #ceiling是向上取整，周岁
meta$age_group=ifelse(meta$age>median(meta$age),'older','younger')
下面解决第3个问题：只需要分期信息
分期信息
从图中我们可以看到，这个字符串是可以按照空格来拆分的，然后取第二部分，但是要注意的是不同的数据可能这里的信息是不一样的，这里的代码需要自己去修改，有的是按照这样来排列的，有的是只有stage I II 这样子，所以具体情况具体对待，如果生搬硬套只会出错的！
#1.4 stage
library(stringr)
meta$stage
tmp = str_split(meta$stage,' ',simplify = T)[,2]
#先按照空格将字符串拆分出来，然后取第2部分
#simplify = T实现了把列表变成了一个向量或者矩阵，简化之后取第2列
str_count(tmp,"T") #数一下这个字符串有几个T，str_count代表计数
str_locate(tmp,"T")[,1] #返回从第几位到第几位，T前面的那一位就是我们想要的分期结束的位置
tmp = str_sub(tmp,1,str_locate(tmp,"T")[,1]-1) #从T开始的位置减去1就是我们要的分期
table(tmp)
meta$stage = tmp
统计一下分期
最后是人种，不需要改动
#1.5 race 人种
table(meta$race)
save(meta,exprSet,group_list,cancer_type,file = paste0(cancer_type,"sur_model.Rdata"))
image.png
##################分割线：上面准备好数据，下面就可以进行生存分析了############
（其实最难的是准备数据，后面的分析，可视化就没有那么难了）
rm(list = ls())
load("TCGA-CHOLsur_model.Rdata")
library(survival)
library(survminer)
#性别
sfit <- survfit(Surv(time, event)~gender, data=meta) #这里用性别先画一个为例，~gender可以换成其他的列，只要是有分组的就可以
ggsurvplot(sfit, conf.int=F, pval=TRUE)
image.png
还可以把它画的更好看一点
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() #把上面两个gender和age_group存到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] #rownames(exprSet)[1]就是第1个基因，以他为例
meta$gene = ifelse(as.integer(exprSet[g,]) > median(as.integer(exprSet[g,])),'high','low')
#上面一行代码实现的是根据(exprSet[g,])是否大于中位数median来判断hign和low，这样就把病人分成了2组
sfit1=survfit(Surv(time, event)~gene, data=meta)
ggsurvplot(sfit1,pval =TRUE, data = meta, risk.table = TRUE)
单个基因
多个基因：和前面的拼图思想是一样的
#多个基因
gs=rownames(exprSet)[1:4] #取4个基因为例
splots <- lapply(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)
image.png
logrank批量生存分析
### 2.logrank批量生存分析
#这里是根据logrankfile把每个基因p值都挑出来，然后根据需求挑出那些p<0.05或者p<0.01的基因
logrankfile = paste0(cancer_type,"log_rank_p.Rdata")
if(!file.exists(logrankfile)){
mySurv=with(meta,Surv(time, event))
log_rank_p <- apply(exprSet , 1 , function(gene){
# gene=exprSet[1,]
meta$group=ifelse(gene>median(gene),'high','low')
data.survdiff=survdiff(mySurv~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)
table(log_rank_p<0.05)
#可以根据names把这些p<0.05/p<0.01的基因挑出来
lr = log_rank_p[log_rank_p<0.05]
names(lr) #就可以看到p<0.05有哪些基因
image.png
image.png
#############分割线：下面是COX回归############
coxfile = paste0(cancer_type,"cox.Rdata")
if(!file.exists(coxfile)){
cox_results <-apply(exprSet , 1 , function(gene){
# gene= exprSet[1,]
group=ifelse(gene>median(gene),'high','low')
survival_dat <- data.frame(group=group,stage=meta$stage,age=meta$age,
gender=meta$gender,
stringsAsFactors = F)
m=coxph(mySurv ~ gender + age + stage + group, data = survival_dat)
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['grouplow',])
})
cox_results=t(cox_results)
save(cox_results,file = coxfile)
}
load(coxfile)
table(cox_results[,4]<0.01)
table(cox_results[,4]<0.05)
lr = names(log_rank_p)[log_rank_p<0.01]
cox = rownames(cox_results)[cox_results[,4]<0.01]
length(intersect(lr,cox))
#76
image.png
😔数据分析太难了，是一条不归路~~
声明：以上代码不是本人原创，只是生信技能树学习数据挖掘班的笔记分享，所以不会答疑，因为我也不知道哈哈~😄
