生存分析需要三个 vector,在一个dataframe中:
- 生存时间,以mouths或者days作单位;
- 结局,'Dead'或者'Alive','Alive'是截尾数据,'Dead'是完全数据;
- 分组信息。
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', 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()参数详细介绍
surv.median.line = 'hv', #中位生存时间线 xlab = 'Follow up time(m)', #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', 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) exp_t=limma::avereps(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)]
三、分析,拟合生存曲线
1 确定基因和高低表达
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, # 创建的拟合对象 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'), # 指定图例分组标签
|