Log-rank批量生存分析
Cox 批量生存分析
作圖
1、Kaplan-Meier生存分析——KM plot
一句代碼即可實(shí)現(xiàn)
sfit <- survfit(Surv(time, event)~gender, data=meta)
image.png
meta——臨床信息表格
image.png
通過(guò)event和time就可以做出一張KM Plot
event里面，1代表已經(jīng)死了，0代表還活著
time里面以月為單位
meta表格里面的數(shù)據(jù)只要能分組就可以畫出KM 圖
用綠色框框出來(lái)的都是可以顯示在圖中的，示例如下：
image.png
2、COX回歸——風(fēng)險(xiǎn)評(píng)估
Cox回歸本質(zhì)上是一種回歸模型，它沒(méi)有直接使用生存時(shí)間，而是使用了風(fēng)險(xiǎn)
比（hazard ratio）作為因變量，該模型不用于估計(jì)生存率，而是用于因素分
析，也就是找到某一個(gè)危險(xiǎn)因素對(duì)結(jié)局事件發(fā)生的貢獻(xiàn)度。
? Cox回歸的重要統(tǒng)計(jì)指標(biāo)：風(fēng)險(xiǎn)比（hazard ratio）
? 當(dāng)HR>1時(shí)，說(shuō)明研究對(duì)象是一個(gè)危險(xiǎn)因素。
? 當(dāng)HR<1時(shí)，說(shuō)明研究對(duì)象是一個(gè)保護(hù)因素。
? 當(dāng)HR=1時(shí)，說(shuō)明研究對(duì)象對(duì)生存時(shí)間不起作用。
3、代碼實(shí)現(xiàn)
生存分析只需要tumor數(shù)據(jù)，不要normal，將其去掉，新表達(dá)矩陣數(shù)據(jù)命名為exprSet；
clinical信息需要進(jìn)一步整理，成為生存分析需要的格式，新臨床信息數(shù)據(jù)命名為meta。
由于不同癌癥的臨床信息表格組織形式不同，這里的代碼需要根據(jù)實(shí)際情況修改。
rm(list=ls())
options(stringsAsFactors = F)
load("TCGA-CHOLgdc.Rdata") #加載之前處理好的數(shù)據(jù)
library(stringr)
#clinical通常有幾十列，選出其中有用的幾列。
tmp = data.frame(colnames(clinical)) #變成data.frame有便于復(fù)制自己想要的列名
clinical = clinical[,c(
'bcr_patient_barcode', #barcode是病人的ID
'vital_status', #生存狀態(tài)
'days_to_death', #死亡日期
'days_to_last_followup', #最后隨訪的日期
'race_list', #人種
'days_to_birth', #年齡：這里其實(shí)不太對(duì)，應(yīng)該換成age_at_initial_pathologic_diagnosis，但是他沒(méi)有，只能用這個(gè)
'gender' , #性別
'stage_event' #生存狀態(tài)
)]
#其實(shí)days_to_last_followup應(yīng)該是找age_at_initial_pathologic_diagnosis，這表格里沒(méi)有，用days_to_birth計(jì)算一下年齡，暫且替代。
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'] #按照選列把腫瘤的樣本留下，把正常的樣本丟掉
#簡(jiǎn)化meta的列名
colnames(meta)=c('ID','event','death','last_followup','race','age','gender','stage')
#調(diào)整meta的ID順序與exprSet列名一致
meta=meta[match(str_sub(colnames(exprSet),1,12),meta$ID),] #讓病人ID和樣本ID一一匹配；1,12代表樣本的前12位，根據(jù)列名的前12位來(lái)調(diào)整meta中ID行的順序
all(meta$ID==str_sub(colnames(exprSet),1,12)) #檢查一下是否相等，返回值是TRUE就沒(méi)問(wèn)題
圖中有123，共3個(gè)問(wèn)題需要解決
1、這里死亡日期和隨訪日期應(yīng)該加起來(lái)才可以，并且這兩列有許多空值，都是NA，NA是沒(méi)有辦法參加計(jì)算的，會(huì)傳染，就是誰(shuí)加上NA最后就會(huì)等于NA，所以要把所有的NA變成0才能參與計(jì)算
2、年齡這一列有負(fù)值，年齡不可能是負(fù)的所以要變成正的
3、stage這一列我們只需要分期信息，也就是羅馬數(shù)字II IV這樣的數(shù)，其他的都要去掉
#整理生存分析的輸入數(shù)據(jù)----
#1.1由隨訪時(shí)間和死亡時(shí)間計(jì)算生存時(shí)間(月)
is.empty.chr = function(x){
ifelse(stringr::str_length(x)==0,T,F) #str_length(x)==0這一步就是在判斷長(zhǎng)度=0的就是空字符串，注意NA的長(zhǎng)度不是0，是1，而空字符串的長(zhǎng)度才是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 #這兩列加起來(lái)是天數(shù)，然后除以30就是月數(shù)
#1.2 根據(jù)生死定義event，活著是0，死的是1
meta$event=ifelse(meta$event=='Alive',0,1) #這里的Alive是因?yàn)楸砀裰幸彩沁@樣寫的，如果別的數(shù)據(jù)中全部是大寫也要換掉！
#1.3 年齡和年齡分組
meta$age=ceiling(abs(as.numeric(meta$age))/365) #ceiling是向上取整，周歲
meta$age_group=ifelse(meta$age>median(meta$age),'older','younger')
下面解決第3個(gè)問(wèn)題：只需要分期信息
分期信息
從圖中我們可以看到，這個(gè)字符串是可以按照空格來(lái)拆分的，然后取第二部分，但是要注意的是不同的數(shù)據(jù)可能這里的信息是不一樣的，這里的代碼需要自己去修改，有的是按照這樣來(lái)排列的，有的是只有stage I II 這樣子，所以具體情況具體對(duì)待，如果生搬硬套只會(huì)出錯(cuò)的！
#1.4 stage
library(stringr)
meta$stage
tmp = str_split(meta$stage,' ',simplify = T)[,2]
#先按照空格將字符串拆分出來(lái)，然后取第2部分
#simplify = T實(shí)現(xiàn)了把列表變成了一個(gè)向量或者矩陣，簡(jiǎn)化之后取第2列
str_count(tmp,"T") #數(shù)一下這個(gè)字符串有幾個(gè)T，str_count代表計(jì)數(shù)
str_locate(tmp,"T")[,1] #返回從第幾位到第幾位，T前面的那一位就是我們想要的分期結(jié)束的位置
tmp = str_sub(tmp,1,str_locate(tmp,"T")[,1]-1) #從T開(kāi)始的位置減去1就是我們要的分期
table(tmp)
meta$stage = tmp
統(tǒng)計(jì)一下分期
最后是人種，不需要改動(dòng)
#1.5 race 人種
table(meta$race)
save(meta,exprSet,group_list,cancer_type,file = paste0(cancer_type,"sur_model.Rdata"))
image.png
##################分割線：上面準(zhǔn)備好數(shù)據(jù)，下面就可以進(jìn)行生存分析了############
（其實(shí)最難的是準(zhǔn)備數(shù)據(jù)，后面的分析，可視化就沒(méi)有那么難了）
rm(list = ls())
load("TCGA-CHOLsur_model.Rdata")
library(survival)
library(survminer)
#性別
sfit <- survfit(Surv(time, event)~gender, data=meta) #這里用性別先畫一個(gè)為例，~gender可以換成其他的列，只要是有分組的就可以
ggsurvplot(sfit, conf.int=F, pval=TRUE)
image.png
還可以把它畫的更好看一點(diǎn)
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
如果是多個(gè)也可以組合到一起
#性別年齡
sfit1=survfit(Surv(time, event)~gender, data=meta)
sfit2=survfit(Surv(time, event)~age_group, data=meta)
splots <- list() #把上面兩個(gè)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這個(gè)函數(shù)實(shí)現(xiàn)了拼圖
dev.off()
image.png
單個(gè)基因
#單個(gè)基因
g = rownames(exprSet)[1] #rownames(exprSet)[1]就是第1個(gè)基因，以他為例
meta$gene = ifelse(as.integer(exprSet[g,]) > median(as.integer(exprSet[g,])),'high','low')
#上面一行代碼實(shí)現(xiàn)的是根據(jù)(exprSet[g,])是否大于中位數(shù)median來(lái)判斷hign和low，這樣就把病人分成了2組
sfit1=survfit(Surv(time, event)~gene, data=meta)
ggsurvplot(sfit1,pval =TRUE, data = meta, risk.table = TRUE)
單個(gè)基因
多個(gè)基因：和前面的拼圖思想是一樣的
#多個(gè)基因
gs=rownames(exprSet)[1:4] #取4個(gè)基因?yàn)槔?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批量生存分析
#這里是根據(jù)logrankfile把每個(gè)基因p值都挑出來(lái)，然后根據(jù)需求挑出那些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)
#可以根據(jù)names把這些p<0.05/p<0.01的基因挑出來(lái)
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
??數(shù)據(jù)分析太難了，是一條不歸路~~
聲明：以上代碼不是本人原創(chuàng)，只是生信技能樹(shù)學(xué)習(xí)數(shù)據(jù)挖掘班的筆記分享，所以不會(huì)答疑，因?yàn)槲乙膊恢拦??
