TCGA
- library('survival')
- library('survminer')
生存分析需要三个 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',
- '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(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'), # 指定图例分组标签
- )
0条评论