本文参考文章：
1.使用R语言的cgdsr包获取TCGA数据
2.TCGA的28篇教程- 对TCGA数据库的任意癌症中任意基因做生存分析
3.TCGA之生存分析（1）基本概念与操作
4.TCGA之生存分析（2）
之前我写过利用R包TCGAbiolinks进行各种数据下载, 那么还有很多种其他的方法可以下载TCGA的数据。如下图：
NAR,2016 May 5; 44(8): e71.
这张图片比较了７种下载数据的方法，可以看到TCGAbiolinks几乎可以下载所有种类的数据。然而我看到好几个大牛们写的教程里，TCGA数据做生存分析的数据下载，都是利用上图中最后一种(CGDS)。我就专门查了一下为啥子，只查到一句：
cBioPortal是按照发表文章的方式来组织TCGA数据的。
既然这样，那就按照流程先走一遍。
（一）查看有多少不同的癌症数据集
> library(cgdsr) #没有这个包的：BiocManager::install("cgdsr")安装
> library(DT)
#Get list of cancer studies at server
#注意！！！坑来了：下一步如果你输入的是大部分教程里写的“"http://www.cbioportal.org/public-portal/"”，这里可能会下载不成功
#我试了很多遍，发现不能加后面的“public-portal”。不排除我自己电脑的问题，但是我用win10和Linux两种系统试了以后都是这个问题
> mycgds <- CGDS("http://www.cbioportal.org/")
#检查下载是否成功，如果是FAILED就是没成功。这里我如果上面加了网址后面的“public-portal”,这里显示的都是FAILED。
> test(mycgds)
getCancerStudies... OK
getCaseLists (1/2) ... OK
getCaseLists (2/2) ... OK
getGeneticProfiles (1/2) ... OK
getGeneticProfiles (2/2) ... OK
getClinicalData (1/1) ... OK
getProfileData (1/6) ... OK
getProfileData (2/6) ... OK
getProfileData (3/6) ... OK
getProfileData (4/6) ... OK
getProfileData (5/6) ... OK
getProfileData (6/6) ... OK
> all_TCGA_studies <- getCancerStudies(mycgds)
> DT::datatable(all_TCGA_studies)
这时你的Rstudio右侧会弹出这样的页面，有点像一个列表：
你也可以直接访问：http://www.cbioportal.org/datasets网站，查看每一个study的具体信息。
（二）查看你感兴趣的数据集
从上面的所有数据集里，选一个你感兴趣的数据集：
#search the dataset you are intested in
> hnscc2018 <- "hnsc_tcga_pan_can_atlas_2018"
查看你选择的这个数据集可用的样品列表：
> getCaseLists(mycgds,hnscc2018)[,c(1,2)]
case_list_id case_list_name
1 hnsc_tcga_pan_can_atlas_2018_all All samples
2 hnsc_tcga_pan_can_atlas_2018_3way_complete Complete samples
3 hnsc_tcga_pan_can_atlas_2018_cna Samples with CNA data
4 hnsc_tcga_pan_can_atlas_2018_log2CNA Samples with log2 copy-number data
5 hnsc_tcga_pan_can_atlas_2018_rna_seq_v2_mrna Samples with mRNA data (RNA Seq V2)
6 hnsc_tcga_pan_can_atlas_2018_cnaseq Samples with mutation and CNA data
7 hnsc_tcga_pan_can_atlas_2018_sequenced Samples with mutation data
查看这个数据集可以下载的基因组数据类别：
> getGeneticProfiles(mycgds,hnscc2018)[,c(1,2)]
genetic_profile_id genetic_profile_name
1 hnsc_tcga_pan_can_atlas_2018_gistic Putative copy-number alterations from GISTIC
2 hnsc_tcga_pan_can_atlas_2018_rna_seq_v2_mrna mRNA Expression, RSEM (Batch normalized from Illumina HiSeq_RNASeqV2)
3 hnsc_tcga_pan_can_atlas_2018_rna_seq_v2_mrna_median_Zscores mRNA expression z-scores relative to diploid samples (RNA Seq V2 RSEM)
4 hnsc_tcga_pan_can_atlas_2018_log2CNA Log2 copy-number values
5 hnsc_tcga_pan_can_atlas_2018_mutations Mutations
6 hnsc_tcga_pan_can_atlas_2018_fusion Fusions
7 hnsc_tcga_pan_can_atlas_2018_rna_seq_v2_mrna_median_all_sample_Zscores mRNA expression z-scores relative to all samples (log RNA Seq V2 RSEM)
选择你感兴趣的基因（列表）：
#choose a gene list you are interested in
> choose_genes <- c("TP53","MKI67")
（三）下载mRNA数据：
> mycaselist <- getCaseLists(mycgds,hnscc2018)[5,1]
> mygeneticprofile <- getGeneticProfiles(mycgds,hnscc2018)[2,1]
> expr <- getProfileData(mycgds,choose_genes,
mygeneticprofile,mycaselist)
> View(expr)
expr长这样：
（四）下载突变数据
> mycaselist <- getCaseLists(mycgds,hnscc2018)[7,1]
> mygeneticprofile <- getGeneticProfiles(mycgds,hnscc2018)[5,1]
> mut_df <- getProfileData(mycgds,choose_genes,
mygeneticprofile,mycaselist)
> View(mut_df)
处理一下突变数据表：
> mut_df <- apply(mut_df,2,as.factor)
＃去除“NAN”
> mut_df[mut_df == "NaN"] = ""
＃去除“NA”
> mut_df[is.na(mut_df)] = ""
> mut_df[mut_df != ''] = "MUT"
现在突变数据表长这样：
（五）下载拷贝数变异数据（CNA）
> mycaselist <- getCaseLists(mycgds,hnscc2018)[3,1]
> mygeneticprofile <- getGeneticProfiles(mycgds,hnscc2018)[1,1]
> cna <- getProfileData(mycgds,choose_genes,
mygeneticprofile,mycaselist)
> View(cna)
拷贝数变异数据：
处理一下拷贝数变异数据：
> rn <- rownames(cna)
> cna <- apply(cna,2,function(x)
as.character(factor(x, levels = c(-2:2),
labels = c("HOMDEL", "HETLOSS", "DIPLOID", "GAIN", "AMP"))))
> cna[is.na(cna)] = ""
> cna[cna == 'DIPLOID'] = ""
> rownames(cna)=rn
处理后的：
（六）下载临床数据
> mycaselist <- getCaseLists(mycgds,hnscc2018)[1,1]
> myclinicaldata <- getClinicalData(mycgds,mycaselist)
> View(myclinicaldata)
临床数据有92列：
你可以选择感兴趣的临床数据类型进行后续分析：
#这里我选了10列
>choose_columns=c('AGE','AJCC_PATHOLOGIC_TUMOR_STAGE','DFS_MONTHS','DFS_STATUS',
'GRADE','NEW_TUMOR_EVENT_AFTER_INITIAL_TREATMENT',
'OS_STATUS','OS_MONTHS','RADIATION_THERAPY','SEX'
)
> choose_clinicaldata <- myclinicaldata[,choose_columns]
关于上面选取的DFS/OS的month和status，大神的文章是这样说的（TCGA的28篇教程- 对TCGA数据库的任意癌症中任意基因做生存分析）：
无病生存期（Disease-free survival,DFS）的定义是指从随机化开始至疾病复发或由于疾病进展导致患者死亡的时间。该指标也常作为抗肿瘤药物III期临床试验的主要终点。某些情况下，DFS与OS相比，作为终点比较难以记录，因为它要求认真随访，及时发现疾病复发，而且肿瘤患者的死亡原因也很难确定。肿瘤患者常有合并症（如心血管病），这些合并症可能会干扰对DFS的判断。并且，肿瘤患者常死于医院外，不能常规进行尸检。
总生存期（Overall survival,OS) 的定义是指从随机化开始至因任何原因引起死亡的时间。该指标常常被认为是肿瘤临床试验中最佳的疗效终点。如果在生存期上有小幅度的提高，可以认为是有意义的临床受益证据。作为一个终点，生存期应每天进行评价，可通过在住院就诊时，通过与患者直接接触或者通过电话与患者交谈，这些相对比较容易记录。确认死亡的日期通常几乎没有困难，并且死亡的时间有其独立的因果关系。当记录至死亡之前的失访患者，通常截止到最后一次有记录的、与患者接触的时间。
综上，一般使用OS_MONTHS和OS_STATUS就可以来估计KM生存函数，画出最简单的生存曲线。
(七）做一个最简单的生存曲线
参考TCGA之生存分析（1）基本概念与操作：
使用survival包的survfit函数进行生存分析，survfit函数需要Surv函数创建的生存分析对象。如果只是简单观察某一组生存数据的生存曲线，而不进行分类、分层、cox回归等分析，则将formula的右侧参数写为1即可
> library(survival)
> dat <- choose_clinicaldata[choose_clinicaldata$OS_MONTHS >0,]
> table(dat$OS_STATUS)
DECEASED LIVING
219 303
#create the survival object, 'DECEASED' means 'occurrence'(估计KM生存曲线)
> my.surv <- Surv(dat$OS_MONTHS,dat$OS_STATUS=='DECEASED')
> kmfit1 <- survfit(my.surv~1)
> library("survminer")
> ggsurvplot(kmfit1,data = dat)
这是最简单的生存曲线，因为我们没有根据突变/拷贝数/基因表达/临床分期等，来画生存曲线。曲线最后趋于平稳，可能这些病人痊愈了，也可能是直到最后一次记录为止他们还健在。
（八）以“性别”分类的生存分析
> my.surv <- Surv(dat$OS_MONTHS,dat$OS_STATUS=='DECEASED')
> kmfit2 <- survfit(my.surv~dat$SEX) #这里只需要将上面的１换成dat$SEX
> ggsurvplot(kmfit2,data = dat)
看一下以性别分类的生存分析结果的显著性：
> survdiff(my.surv~dat$SEX)
Call:
survdiff(formula = my.surv ~ dat$SEX)
n=522, 1 observation deleted due to missingness.
N Observed Expected (O-E)^2/E (O-E)^2/V
dat$SEX=Female 141 71 56.3 3.86 5.25
dat$SEX=Male 381 148 162.7 1.34 5.25
Chisq= 5.2 on 1 degrees of freedom, p= 0.02
你还可以直接让ggsurvplot在画图时直接显示p值：
> ggsurvplot(kmfit2,data = dat,pval=T)
（九）根据肿瘤分期进行生存分析
> ggsurvplot(kmfit2,data = dat,pval=T)
> kmfit3 <- survfit(my.surv~dat$AJCC_PATHOLOGIC_TUMOR_STAGE)
> ggsurvplot(kmfit3,data = dat,pval=T)
（十）根据基因是否突变进行生存分析
在最开始，你定义了你感兴趣的基因，这时你可以利用你选的基因是否有突变来进行生存分析：
＃你会发现你的临床信息里的样品数和你的mut_df列表里的样品数不一样，你需要先取它们的交集
> length(intersect(rownames(choose_clinicaldata),rownames(mut_df)))
[1] 515
> dat2 <- cbind(choose_clinicaldata[rownames(mut_df),c('OS_STATUS','OS_MONTHS')],mut_df)
> dat2 <- dat2[dat2$OS_MONTHS > 0,]
> dat2 <- dat2[!is.na(dat$OS_STATUS),]
> dat2$OS_STATUS <- as.character(dat2$OS_STATUS)
> attach(dat2)
> View(dat2)
> my.surv <- Surv(OS_MONTHS,OS_STATUS=='DECEASED')
> kmfit4 <- survfit(my.surv~TP53,data = dat2)
> kmfit5 <- survfit(my.surv~MKI67,data = dat2)
> library("survminer")
> ggsurvplot(kmfit4,conf.int =F, pval = T,risk.table =T, ncensor.plot = TRUE,
title="Survival plot by TP53 Mutaiton")
> ggsurvplot(kmfit5,conf.int =F, pval = T,risk.table =T, ncensor.plot = TRUE,
title="Survival plot by MKI67 Mutaiton") #图不放了
> detach(dat2)
（十一）根据基因的拷贝数进行生存分析
> length(intersect(rownames(choose_clinicaldata),rownames(cna)))
[1] 523
> dat <- cbind(choose_clinicaldata[rownames(cna),c('OS_STATUS','OS_MONTHS')],cna)
> dat <- dat[dat$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')
> kmfit7 <- survfit(my.surv~TP53,data = dat)
> kmfit8 <- survfit(my.surv~MKI67,data = dat)
> ggsurvplot(kmfit7,conf.int =F, pval = T,risk.table =T, ncensor.plot = TRUE,
title="Survival plot by CNA")
（十二）根据基因表达情况进行生存分析
先看一下你选的基因在“LIVING”和“DECEASED"两组里的表达情况：
> library(survival)
> dat3=cbind(choose_clinicaldata[,c('OS_STATUS','OS_MONTHS')],expr[rownames(choose_clinicaldata),])
> dat3=dat3[dat3$OS_MONTHS > 0,]
> dat3=dat3[!is.na(dat3$OS_STATUS),]
> dat3$OS_STATUS=as.character(dat3$OS_STATUS)
> library(ggpubr)
> p <- ggboxplot(dat3, x="OS_STATUS", y="TP53", color = "OS_STATUS",
palette = "jco", add = "jitter")
> p+stat_compare_means(method = "t.test")
> dat3=dat3[!is.na(dat3$TP53),]
> dat3$TP53 = ifelse(dat3$TP53 > median(dat3$TP53),'high','low')
> attach(dat3)
> table(dat3$TP53)
high low
257 257
> library(survival)
> my.surv <- Surv(dat3$OS_MONTHS,dat3$OS_STATUS=='DECEASED')
> kmfit1 <- survfit(my.surv~TP53,data = dat3)
> ggsurvplot(kmfit1,conf.int =F, pval = T,risk.table =T, ncensor.plot = TRUE)
关于COX回归生存分析，可以参考本文开始提到的几篇文章。
网友评论