TCGA学习03：生存分析

2020-05-08 本文已影响0人小贝学生信

TCGA学习01：数据下载与整理 - 简书
 TCGA学习02：差异分析 - 简书
 TCGA学习03：生存分析 - 简书
 TCGA学习04：建模预测-cox回归 - 简书
 TCGA学习04：建模预测-lasso回归 - 简书
 TCGA学习04：建模预测-随机森林&向量机 - 简书

第三步：生存分析

0、数据准备

（1）原始数据下载

因为生存分析考虑到样本多一点会好一些，因此使用R包TCGA-biolinks下载了一个五百多样本的数据来分析，的确挺方便的。

library(TCGAbiolinks)
TCGAbiolinks:::getGDCprojects()$project_id
cancer_type="TCGA-LUAD"
clinical=GDCquery_clinic(project=cancer_type,type="clinical")
dim(clinical)
clinical[1:4,1:4]
head(colnames(clinical))
query <- GDCquery(project = cancer_type, 
                  data.category = "Transcriptome Profiling", 
                  data.type = "miRNA Expression Quantification")
GDCdownload(query, method = "api", files.per.chunk = 50)
expdat <- GDCprepare(query = query)

（2）数据整理

表达矩阵

library(tibble)
rownames(expdat) <- NULL
#将第一列设置为行名
expdat <- column_to_rownames(expdat,var = "miRNA_ID")
#取出1，4，7...即所有的样本count列（其它列不需要）
exp = t(expdat[,seq(1,ncol(expdat),3)])
exp[1:3,1:3]
#修改列名
rownames(exp)=substr(rownames(exp),12,39)
exp[1:3,1:3]
dim(exp)
#表达矩阵一般行名为基因，列名为样本
exp=t(exp)
exp[1:3,1:3]
group_list=ifelse(as.numeric(substr(colnames(exp),14,15)) < 10,'tumor','normal')
table(group_list)
#仅取tumor样本，521个
exp_tumor=exp[,group_list=='tumor']

原始表达矩阵共567个样本

临床信息

meta = clinical
#同样以第一列作为列名
meta <- column_to_rownames(meta,var = "submitter_id")
colnames(meta)
#筛选我们感兴趣的临床信息
meta=meta[,colnames(meta) %in% c("vital_status",
                                 "days_to_last_follow_up",
                                 "days_to_death",
                                 "race",
                                 "gender",
                                 "age_at_index",
                                 "tumor_stage")]
dim(meta) #585个

head(colnames(exp_tumor))
head(rownames(meta))
#调整、筛选临床样本信息数量与顺序与表达矩阵相同
meta=meta[match(substr(colnames(exp_tumor),1,12),rownames(meta)),]

dim(exp_tumor)
head(colnames(exp_tumor))
dim(meta)
head(rownames(meta))

表达矩阵与临床信息样本数量、顺序一致

基于上述数据，生存分析用到的临床信息类型大致有包括一下六类

#1、计算生存时间
meta$days_to_death[is.na(meta$days_to_death)] <- 0   #缺失值标记为0
meta$days_to_last_follow_up[is.na(meta$days_to_last_follow_up)] <- 0
meta$days=as.numeric(meta[,2])+as.numeric(meta[,3])
#时间以月份记，保留两位小数
meta$time=round(meta$days/30,2)

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

#3 年龄分组(部分样本缺失，考虑可能的影响应该不大)
meta$age_at_index[is.na(meta$age_at_index)] <- 0
meta$age_at_index=as.numeric(meta$age_at_index)
meta$age_group=ifelse(meta$age_at_index>median(meta$age_at_index),'older','younger')
table(meta$age_group)

#4 癌症阶段
table(meta$tumor_stage)

#5 race 人种
table(meta$race)

#6 性别 gender
table(meta$gender)

关于TCGA临床数据的生存时间，还有一种说法是https://www.biostars.org/p/397150/ 。对于大部分病人，计算结果是一样的。

生存分析的一些概念，如生存时间，删失值可参考生存分析与R一文

save(exp_tumor,meta,file="tosur.RData")
#利用这两个文件接下来就可以开始生存分析了

1、初步了解生存分析图

rm(list=ls())
load("tosur.RData")
library(survival)
library(survminer)

sfit <- survfit(Surv(time, event)~gender, data=meta)
print(sfit)
ggsurvplot(sfit, conf.int=F, pval=TRUE)

生存曲线（survival curve）则是将每个时间点的生存率连接在一起的曲线，一般随访时间为X轴，生存率为Y轴；曲线平滑则说明高生存率，反之则低生存率；中位生存率（median survival time）越长，则说明预后较好

ggsurvplot(sfit, conf.int=F, pval=TRUE)

补充

Surv(time, event)根据生死状态，标记出哪些生存时间数据是删失值，用加号表示区分。
Surv(time, event)
survfit()用于计算生存曲线，属于Kaplan-Meier方法；~gender表示以性别分组绘图，注意此方法仅适用二分类变量进行分组。如果不想分组，只看总体，~1即可。
最后配合 ggsurvplot()函数可视化生存曲线；通过设置参数，可在生存曲线基础上，再绘制两张辅助图。

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

Number at risk图通过risk.table =TRUE设置，表示 the number of subjects at risk at time t. 我的理解是在该生存时间长度内，还存活的case；较直接得反映出两组的一个差异情况
Number of censoring图通过ncensor.plot = TRUE设置，表示the number of censored subjects, who exit the risk set, without an event, at time t. 即该时间段内的删失值。但是我绘制出来的和别人的不一样，也读不太懂，不知道是不是和R的warning提示有关，存疑。
warning
别人的图..
三图合一
对性别和年龄生存分析拼图

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)

性别和年龄的生存分析

2、基因表达的生存分析

思路：设定某一基因的表达基准，将临床样本分为高表达high与低表达low两组，再进行生存分析

(1) 直接指定感兴趣基因，可以是1个或多个

exprSet=exp_tumor  #套流程
g = rownames(exprSet)[150] # 随便选一个
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)

2.1

(2) 指定多个基因拼图

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

(3) 批量计算所有表达矩阵的生存分析p值，从而挑选出“显著基因”

log-rank方法--survdiff()

mySurv=with(meta,Surv(time, event))
head(exprSet)
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)
})
#上述如果返回报错说只有一组，需要将一些低表达数据删除点
#expr = expr[apply(expr, 1, function(x) {
#sum(x > 1) > 9  #过滤掉低表达基因
#}), ]  
log_rank_p=sort(log_rank_p)

image.png

此外花花老师还介绍了cox批量生存分析的方法，因为运行时报错，不知道咋回事，先就记录下代码吧。之后有机会再回过头看看这三种方法的原理。

#cox批量生存分析----
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)
table(cox_results[,4]<0.01)
table(cox_results[,4]<0.05)

Error

参考资料
1、生存分析（SurvivalAnalysis）fanyucai新浪博客
2、R语言生存分析入门 - 简书
3、R语言-Survival analysis（生存分析） | 生信笔记
4、R语言绘制生存曲线估计|生存分析|如何R作生存曲线图人工智能大数据部落-CSDN博客
5、生存分析与R_人工智能_走在码农路上的医学狗-CSDN博客

当然还是最要感谢花花老师的教程~~

闲聊几句，可能要等下半年才开学了，诶~~~