TCGA_生存分析

library("survival")
library("survminer")

生存分析需要三个 vector，在一个dataframe中：

1. 生存时间，以mouths或者days作单位；
2. 结局，"Dead"或者"Alive"，"Alive"是截尾数据，"Dead"是完全数据；
3. 分组信息。

Age_old vs young

一、读入数据

data_cl <- read.csv(file = "Results/测序or临床数据下载/data_cl.csv", header=T, row.names=2,check.names=FALSE)
t_needed=c("vital_status",
           "days_to_last_follow_up",
           "days_to_death",
           "age_at_index",
           "tumor_grade")
meta=data_cl[,t_needed] #筛选需要的临床信息
meta=meta[meta$vital_status %in% c('Alive','Dead'),] #排除结局为"Not Reported"的Sample

View(data_cl)

data_cl$vital_status 中有两个为 Not Reported，需要排除

二、整理数据

1 计算生存时间

meta$days_to_last_follow_up[is.na(meta$days_to_last_follow_up)] = 0 #is.na()用于返回是否为缺失值
meta$days_to_death[is.na(meta$days_to_death)] = 0   
meta$days<-ifelse(meta$vital_status=='Alive',meta$days_to_last_follow_up,meta$days_to_death)

meta$mouth=round(meta$days/30,2) #以month为单位，保留两位小数

2 添加age_group列（分组数据）

meta$age_group = ifelse(meta$age_at_index>median(meta$age_at_index),'old','young')

 三、分析

Surv() 函数输出带有截尾信息的生存时间数据；

survfit() 函数根据生存时间数据、分组信息，并基于”K-M法“输出拟合数据

survData = Surv(time = meta$mouth,                 #生存时间数据
                event = meta$vital_status=='Dead') #判断结局，完全数据/截尾数据

KMfit <- survfit(survData ~ meta$age_group)  # ~ 后是指定的分组

head(survData)        有"+"为截尾数据

[1]  13.20+ 116.73   73.27+  50.57+  15.43+  23.67

survData[1:10,1:2]        survData有2个列维度，"status==0"为截尾数据

         time status
 [1,]  13.20      0
 [2,] 116.73      1
 [3,]  73.27      0
 [4,]  50.57      0
 [5,]  15.43      0
 [6,]  23.67      0
 [7,]  64.90      0
 [8,]  77.47      0
 [9,]  87.60      0
[10,]   3.23      0

四、拟合生存曲线

ggsurvplot()参数详细介绍

ggsurvplot(KMfit,                     #拟合对象
           data = meta,               #变量数据来源
           pval = TRUE,               #P值
           surv.median.line = "hv",   #中位生存时间线
           risk.table = TRUE,         #风险表
           xlab = "Follow up time(m)",  #x轴标签
           break.x.by = 10,          #x轴刻度间距
           #legend = c(0.8,0.75),     #图例位置
           #legend.title = "",        #图例标题
           #legend.labs = c("old", "young"),  #图例分组标签
           ) 

 Gene Expression_high vs low

一、读入数据（临床数据+表达数据）

meta 中留下有结局的Sample；

exp 中合并重复的列名（一个Sample可能有多个组织样本，这里采用取平均的方式去重）

data_cl <- read.csv(file = "Results/测序or临床数据下载/data_cl.csv", header=T, row.names=2,check.names=FALSE)
t_needed=c("vital_status",
           "days_to_last_follow_up",
           "days_to_death",
           "age_at_index",
           "tumor_grade")
meta=data_cl[,t_needed] #筛选需要的临床信息
meta=meta[meta$vital_status %in% c('Alive','Dead'),] #排除结局为"Not Reported"的Sample

exp <- read.csv(file = "Results/测序or临床数据下载/dataFilt.csv", header=T, row.names=1,check.names=FALSE)
colnames(exp)=str_sub(colnames(exp),1,12)
colmeans=function(x){
  exp_m=as.matrix(x)
  exp_t=t(exp_m)
  exp_t=limma::avereps(exp_t)
  t(exp_t)
}
exp=colmeans(exp) #取平均去除重复列名的Sample

View(data_cl)

data_cl$vital_status 中有两个为 Not Reported，需要排除

二、整理数据

1 计算生存时间

meta$days_to_last_follow_up[is.na(meta$days_to_last_follow_up)] = 0 #is.na()用于返回是否为缺失值
meta$days_to_death[is.na(meta$days_to_death)] = 0   
meta$days<-ifelse(meta$vital_status=='Alive',meta$days_to_last_follow_up,meta$days_to_death)

meta$mouth=round(meta$days/30,2) #以month为单位，保留两位小数

2 筛选meta中有表达信息的Sample

t_index = rownames(meta)[rownames(meta) %in% colnames(exp)]
meta=meta[t_index,]

三、分析，拟合生存曲线

1 确定基因和高低表达

Gene = "CHAC1"
meta$Expression_level = ifelse(exp[Gene,rownames(meta)]>median(exp[Gene,]),'high','low')

2 作图 

survData = Surv(time=meta$mouth, #月份数据
                  event=meta$vital_status=='Dead') #判断哪些是截尾数据
KMfit <- survfit(survData ~ meta$Expression_level)  #~后指定分组

ggsurvplot(KMfit, # 创建的拟合对象
           data = meta,  # 指定变量数据来源
           pval = TRUE, # 添加P值
           surv.median.line = "hv",  # 添加中位生存时间线
           risk.table = TRUE, # 添加风险表
           ncensor.plot = FALSE, #？？图
           xlab = "Follow up time(m)", # 指定x轴标签
           break.x.by = 10, # 设置x轴刻度间距
           palette = c("#E7B800", "#2E9FDF"),
           #legend = c(0.8,0.75), # 指定图例位置
           legend.title = Gene, # 设置图例标题
           #legend.labs = c("old", "young"), # 指定图例分组标签
)