TCGA数据挖掘（4）：生存分析

image.png

生存分析

image.png event：0代表或者，1：代表阳性结局；time：活着的时间；年龄是不符合生存曲线的

COX回归

HR>1,表示与病人的危险正相关

生存分析画图代码

生存分析前要进行数据的整理

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)

#clinical通常有几十列，选出其中有用的几列。
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']

#简化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),]
all(meta$ID==str_sub(colnames(exprSet),1,12))

#整理生存分析的输入数据----

#1.1由随访时间和死亡时间计算生存时间(月)

meta[,3][meta[,3]==""]=0
meta[,4][meta[,4]==""]=0
meta$time=(as.numeric(meta[,3])+as.numeric(meta[,4]))/30

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

#1.3 年龄和年龄分组
meta$age=ceiling(abs(as.numeric(meta$age))/365)
meta$age_group=ifelse(meta$age>median(meta$age),'older','younger')

#1.4 stage 
library(stringr) 
meta$stage
tmp = str_split(meta$stage,' ',simplify = T)[,2]
str_count(tmp,"T")
str_locate(tmp,"T")[,1]
tmp = str_sub(tmp,1,str_locate(tmp,"T")[,1]-1)
table(tmp)
meta$stage = tmp

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

生存模型

生存模型的目标：进一步缩小基因的范围 训练集构建模型，测试集评估模型