下载数据#载入必要的包
library('cgdsr')
library('tidyverse')
#获取乳腺癌的表达量数据
mycgds <- CGDS('http://www./')
getCancerStudies(mycgds) %>% filter(grepl('brca',.[,1])) %>%DT::datatable()
mycancerstudy <- 'brca_tcga'
#getProfileData用于获取特定的基因组数据（RNA、Mutation、CNV等等），
#需要先确定样本列表、基因组数据类型及特定的基因
getCaseLists(mycgds,mycancerstudy)[,c(1,2)] #查看当前可用的样本列表
getGeneticProfiles(mycgds,mycancerstudy)[,c(1,2)] #查看当前可以下载的基因组数据类别，根据需要的数据类别选择相应的值
choose_genes <- c('BRCA1','BRCA2') #关注两个基因
BRCA可用的样本列表如下： 
BRCA可以下载的遗传数据列表如下： 
##下载mRNA表达数据##
mycaselist <- getCaseLists(mycgds,mycancerstudy)[8,1]
mygeneticprofile <- getGeneticProfiles(mycgds,mycancerstudy)[6,1]
expr <- getProfileData(mycgds,choose_genes,
mygeneticprofile,mycaselist)
##下载获取临床数据##
mycaselist <- getCaseLists(mycgds,mycancerstudy)[1,1]
myclinicaldata <- getClinicalData(mycgds,mycaselist) #共有107列信息
choose_columns <- c('AJCC_METASTASIS_PATHOLOGIC_PM','AJCC_NODES_PATHOLOGIC_PN','AJCC_PATHOLOGIC_TUMOR_STAGE','AJCC_TUMOR_PATHOLOGIC_PT',
'AGE','SEX','VITAL_STATUS','OS_STATUS','OS_MONTHS',
'DFS_MONTHS','DFS_STATUS')
choose_clinicaldata <- myclinicaldata[,choose_columns] #只选择其中的11列用于后续分析
colnames(choose_clinicaldata)[1:4] <- c('M','N','STAGE','T')
简单生存分析使用survival包的survfit函数进行生存分析，survfit函数需要Surv函数创建的生存分析对象。如果只是简单观察某一组生存数据的生存曲线，而不进行分类、分层、cox回归等分析，则将formula的右侧参数写为1即可，如my.surv~1， Surv(dat$OS_MONTHS,dat$OS_STATUS=='DECEASED')~1等等。 使用survminer的ggsurvplot函数进行生存分析的可视化，使用方法很简单，将survfit函数的结果直接传入此函数即可，如果没有将数据集传给survfit函数，则在ggsurvplot函数中需要指明数据集（也就是参数data，如data=dat）。 ###简单生存分析
library(survival)
dat <- choose_clinicaldata %>% filter(OS_MONTHS > 0)
table(dat$OS_STATUS)
#DECEASED LIVING
# 154 935
my.surv <- Surv(dat$OS_MONTHS,dat$OS_STATUS=='DECEASED') #创建生存分析对象,DECEASED代表事件发生
kmfit1 <- survfit(my.surv~1)
library('survminer')
ggsurvplot(kmfit1,data = dat)
分类数据的生存分析以性别为分类变量，同简单分析的区别在于，将分类变量放置到formula的右侧，如my.surv~dat$SEX， Surv(dat$OS_MONTHS,dat$OS_STATUS=='DECEASED') ~dat$SEX。 #以性别为分类变量
my.surv <- Surv(dat$OS_MONTHS,dat$OS_STATUS=='DECEASED')
kmfit2 <- survfit(my.surv~dat$SEX)
ggsurvplot(kmfit2,data = dat)
还可以使用strata对分层变量进行控制： #还可以使用strata参数控制分层变量
survfit(my.surv~dat$SEX+strata(dat$AGE)) #AGE为分层变量
survfit(my.surv~dat$SEX+strata(dat$AGE,dat$M)) #AGE和M都是分层变量
ggsurvplot是可以自动添加p值的，使用pval=T参数，另外如前面所述，ggsurvplot会重新进行显著性计算，如果survfit函数没有指定data参数的话，那么需要在ggsurvplot中传入正确的data参数，比如： my.surv <- Surv(dat$OS_MONTHS,dat$OS_STATUS=='DECEASED')
kmfit3 <- survfit(my.surv~dat$STAGE)
ggsurvplot(kmfit3,data = dat,pval=T)
其实也可以使用survival包的survdiff函数手动进行显著性计算，同样的可以使用strata参数控制分层变量： survdiff(my.surv~dat$SEX)
survdiff(my.surv~dat$SEX+strata(dat$AGE,dat$M)) #还可以使用strata参数控制分层变量
cox回归生存分析survival包的coxph函数用于cox回归生存分析，详细的结果还可以使用summary函数进行查看。 绘图时，由于ggsurvplot函数需要survfit对象，因此需要将cox回归的结果传入survfit函数，也就是survfit(m, data = dat)，但是ggsurvplot在进行绘图时需要获取survfit的formula，并重新进行计算。如果此时直接传入survfit(m, data = dat)的结果的话，那么ggsurvplot就会错误读取到formula是m，而不是真正的formula：my.surv ~ AGE + SEX + N + T + M + STAGE，此时就会报错。所以对于cox回归的生存分析，使用ggsurvplot绘图时需要人工指定formula：fit$call$formula <- m$call$formula，如下： dat <- choose_clinicaldata[choose_clinicaldata$OS_MONTHS > 0,]
dat <- dat[!is.na(dat$OS_STATUS),]
dat$OS_STATUS <- as.character(dat$OS_STATUS)
attach(dat)
my.surv <- Surv( OS_MONTHS,OS_STATUS=='DECEASED')
m<-coxph(my.surv ~ AGE + SEX + N + T + M + STAGE, data = dat) #cox回归生存分析
#summary（m)
fit <- survfit(m, data = dat)
fit$call$formula <- m$call$formula #此步骤必不可少，否则ggsurvplot无法出图
ggsurvplot(fit, palette = '#2E9FDF', ggtheme = theme_minimal())
detach(dat)
参考资料 TCGA的28篇教程- 对TCGA数据库的任意癌症中任意基因做生存分析 【r<-统计|绘图】使用R进行生存分析——一文打尽https://www.jianshu.com/p/4ad9ba730719
|
