TCGA数据生存分析 |
您所在的位置：网站首页 › 生存语音包 › TCGA数据生存分析 |
|
本文参考文章： 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 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 DT::datatable(all_TCGA_studies)这时你的Rstudio右侧会弹出这样的页面，有点像一个列表： 你也可以直接访问：http://www.cbioportal.org/datasets网站，查看每一个study的具体信息。 （二）查看你感兴趣的数据集从上面的所有数据集里，选一个你感兴趣的数据集： #search the dataset you are intested in > hnscc2018 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 mycaselist mygeneticprofile expr View(expr)expr长这样： （四）下载突变数据 > mycaselist mygeneticprofile mut_df View(mut_df)处理一下突变数据表： > mut_df mut_df[mut_df == "NaN"] = "" ＃去除“NA” > mut_df[is.na(mut_df)] = "" > mut_df[mut_df != ''] = "MUT"现在突变数据表长这样： （五）下载拷贝数变异数据（CNA） > mycaselist mygeneticprofile cna View(cna)拷贝数变异数据： 处理一下拷贝数变异数据： > rn cna cna[is.na(cna)] = "" > cna[cna == 'DIPLOID'] = "" > rownames(cna)=rn处理后的： （六）下载临床数据 > mycaselist myclinicaldata 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 library(survival) > dat 0,] > table(dat$OS_STATUS) DECEASED LIVING 219 303 #create the survival object, 'DECEASED' means 'occurrence'(估计KM生存曲线) > my.surv kmfit1 library("survminer") > ggsurvplot(kmfit1,data = dat)这是最简单的生存曲线，因为我们没有根据突变/拷贝数/基因表达/临床分期等，来画生存曲线。曲线最后趋于平稳，可能这些病人痊愈了，也可能是直到最后一次记录为止他们还健在。 （八）以“性别”分类的生存分析 > my.surv kmfit2 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 ggsurvplot(kmfit3,data = dat,pval=T) （十）根据基因是否突变进行生存分析在最开始，你定义了你感兴趣的基因，这时你可以利用你选的基因是否有突变来进行生存分析： ＃你会发现你的临床信息里的样品数和你的mut_df列表里的样品数不一样，你需要先取它们的交集 > length(intersect(rownames(choose_clinicaldata),rownames(mut_df))) [1] 515 > dat2 dat2 0,] > dat2 dat2$OS_STATUS attach(dat2) > View(dat2) > my.surv kmfit4 kmfit5 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 dat 0,] > dat dat$OS_STATUS attach(dat) > my.surv kmfit7 kmfit8 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 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 kmfit1 ggsurvplot(kmfit1,conf.int =F, pval = T,risk.table =T, ncensor.plot = TRUE)关于COX回归生存分析，可以参考本文开始提到的几篇文章。 |
今日新闻 |
推荐新闻 |
| CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |