下載數(shù)據(jù)#載入必要的包
library('cgdsr')
library('tidyverse')
#獲取乳腺癌的表達(dá)量數(shù)據(jù)
mycgds <- CGDS('http://www./')
getCancerStudies(mycgds) %>% filter(grepl('brca',.[,1])) %>%DT::datatable()
mycancerstudy <- 'brca_tcga'
#getProfileData用于獲取特定的基因組數(shù)據(jù)（RNA、Mutation,、CNV等等）,，
#需要先確定樣本列表、基因組數(shù)據(jù)類型及特定的基因
getCaseLists(mycgds,mycancerstudy)[,c(1,2)] #查看當(dāng)前可用的樣本列表
getGeneticProfiles(mycgds,mycancerstudy)[,c(1,2)] #查看當(dāng)前可以下載的基因組數(shù)據(jù)類別,，根據(jù)需要的數(shù)據(jù)類別選擇相應(yīng)的值
choose_genes <- c('BRCA1','BRCA2') #關(guān)注兩個(gè)基因
BRCA可用的樣本列表如下： 
BRCA可以下載的遺傳數(shù)據(jù)列表如下： 
##下載mRNA表達(dá)數(shù)據(jù)##
mycaselist <- getCaseLists(mycgds,mycancerstudy)[8,1]
mygeneticprofile <- getGeneticProfiles(mycgds,mycancerstudy)[6,1]
expr <- getProfileData(mycgds,choose_genes,
mygeneticprofile,mycaselist)
##下載獲取臨床數(shù)據(jù)##
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列用于后續(xù)分析
colnames(choose_clinicaldata)[1:4] <- c('M','N','STAGE','T')
簡(jiǎn)單生存分析使用survival包的survfit函數(shù)進(jìn)行生存分析,，survfit函數(shù)需要Surv函數(shù)創(chuàng)建的生存分析對(duì)象。如果只是簡(jiǎn)單觀察某一組生存數(shù)據(jù)的生存曲線,，而不進(jìn)行分類,、分層、cox回歸等分析,，則將formula的右側(cè)參數(shù)寫(xiě)為1即可,，如my.surv~1， Surv(dat$OS_MONTHS,dat$OS_STATUS=='DECEASED')~1等等,。 使用survminer的ggsurvplot函數(shù)進(jìn)行生存分析的可視化,，使用方法很簡(jiǎn)單，將survfit函數(shù)的結(jié)果直接傳入此函數(shù)即可,，如果沒(méi)有將數(shù)據(jù)集傳給survfit函數(shù)，則在ggsurvplot函數(shù)中需要指明數(shù)據(jù)集（也就是參數(shù)data,，如data=dat）,。 ###簡(jiǎn)單生存分析
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') #創(chuàng)建生存分析對(duì)象,DECEASED代表事件發(fā)生
kmfit1 <- survfit(my.surv~1)
library('survminer')
ggsurvplot(kmfit1,data = dat)
分類數(shù)據(jù)的生存分析以性別為分類變量，同簡(jiǎn)單分析的區(qū)別在于,，將分類變量放置到formula的右側(cè),，如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對(duì)分層變量進(jìn)行控制： #還可以使用strata參數(shù)控制分層變量
survfit(my.surv~dat$SEX+strata(dat$AGE)) #AGE為分層變量
survfit(my.surv~dat$SEX+strata(dat$AGE,dat$M)) #AGE和M都是分層變量
ggsurvplot是可以自動(dòng)添加p值的,，使用pval=T參數(shù)，另外如前面所述,，ggsurvplot會(huì)重新進(jìn)行顯著性計(jì)算,，如果survfit函數(shù)沒(méi)有指定data參數(shù)的話，那么需要在ggsurvplot中傳入正確的data參數(shù),，比如： my.surv <- Surv(dat$OS_MONTHS,dat$OS_STATUS=='DECEASED')
kmfit3 <- survfit(my.surv~dat$STAGE)
ggsurvplot(kmfit3,data = dat,pval=T)
其實(shí)也可以使用survival包的survdiff函數(shù)手動(dòng)進(jìn)行顯著性計(jì)算,，同樣的可以使用strata參數(shù)控制分層變量： survdiff(my.surv~dat$SEX)
survdiff(my.surv~dat$SEX+strata(dat$AGE,dat$M)) #還可以使用strata參數(shù)控制分層變量
cox回歸生存分析survival包的coxph函數(shù)用于cox回歸生存分析，詳細(xì)的結(jié)果還可以使用summary函數(shù)進(jìn)行查看,。 繪圖時(shí),，由于ggsurvplot函數(shù)需要survfit對(duì)象,，因此需要將cox回歸的結(jié)果傳入survfit函數(shù)，也就是survfit(m, data = dat),，但是ggsurvplot在進(jìn)行繪圖時(shí)需要獲取survfit的formula,，并重新進(jìn)行計(jì)算。如果此時(shí)直接傳入survfit(m, data = dat)的結(jié)果的話,，那么ggsurvplot就會(huì)錯(cuò)誤讀取到formula是m,，而不是真正的formula：my.surv ~ AGE + SEX + N + T + M + STAGE，此時(shí)就會(huì)報(bào)錯(cuò),。所以對(duì)于cox回歸的生存分析,，使用ggsurvplot繪圖時(shí)需要人工指定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無(wú)法出圖
ggsurvplot(fit, palette = '#2E9FDF', ggtheme = theme_minimal())
detach(dat)
參考資料 TCGA的28篇教程- 對(duì)TCGA數(shù)據(jù)庫(kù)的任意癌癥中任意基因做生存分析 【r<-統(tǒng)計(jì)|繪圖】使用R進(jìn)行生存分析——一文打盡https://www.jianshu.com/p/4ad9ba730719
|
