Closed ixxmu closed 2 years ago
众所周知,TCGA数据库是目前最综合全面的癌症病人相关组学数据库,包括的多组学数据有:
知名的肿瘤研究机构都有着自己的TCGA数据库探索工具,比如:
其中Broad研究所的就是FireBrowse,主页在:http://www.firebrowse.org/
这个网页工具当然是非常强大,不过咱们生信工程师喜欢的仍然是编程语言,所以有一个RTCGAToolbox的R包可以帮助我们通过代码来玩转它的网页工具,安装它超级简单:
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("RTCGAToolbox")
library(RTCGAToolbox)
# 查看哪些癌症数据可以下载
getFirehoseDatasets()
getFirehoseRunningDates( )
getFirehoseAnalyzeDates( )
# 可以看到工具在 20160128 就停止更新了
可以看到不同的癌症数据集,不同的更新时间,如下所示:
> getFirehoseDatasets()
[1] "ACC" "BLCA" "BRCA" "CESC" "CHOL" "COADREAD" "COAD" "DLBC"
[9] "ESCA" "FPPP" "GBMLGG" "GBM" "HNSC" "KICH" "KIPAN" "KIRC"
[17] "KIRP" "LAML" "LGG" "LIHC" "LUAD" "LUSC" "MESO" "OV"
[25] "PAAD" "PCPG" "PRAD" "READ" "SARC" "SKCM" "STAD" "STES"
[33] "TGCT" "THCA" "THYM" "UCEC" "UCS" "UVM"
这里用的是简称,比如BRCA就是乳腺癌。
而不同的时间,指的是TCGA数据库在发展过程中样本量的增加, 而FireBrowse是按照时间来定期运行程序处理数据的,所以一般来说用最新版的结果,就会涵盖TCGA里面的所有的样本了。
接下来下载所需要的数据,这里以乳腺癌为例,数据下载完后会直接放在你的工作目录,不同地方下载的速度不一样。
## 下载数据,需要选择癌症种类,数据分析时间,还有数据的种类
options(timeout=10000)
# 一般来说,我们会选择最新的数据,工具在 20160128 就停止更新了
brcaData = getFirehoseData (dataset="BRCA", runDate="20160128",
forceDownload = TRUE,
clinical=TRUE, Mutation=TRUE)
save(brcaData,file='brcaData.RTCGAToolbox.Rdata')
这里测试了临床信息和突变信息的数据下载,因为它们比较小,所以下载速度会很快,这里下载的数据包括:
trying URL 'http://gdac.broadinstitute.org/runs/stddata__2016_01_28/data/BRCA/20160128/gdac.broadinstitute.org_BRCA.Clinical_Pick_Tier1.Level_4.2016012800.0.0.tar.gz'
Content type 'application/x-gzip' length 164047 bytes (160 KB)
trying URL 'http://gdac.broadinstitute.org/runs/stddata__2016_01_28/data/BRCA/20160128/gdac.broadinstitute.org_BRCA.Mutation_Packager_Calls.Level_3.2016012800.0.0.tar.gz'
Content type 'application/x-gzip' length 10454503 bytes (10.0 MB)
trying URL 'http://gdac.broadinstitute.org/runs/analyses__2016_01_28/data/BRCA/20160128/gdac.broadinstitute.org_BRCA-TP.CopyNumber_Gistic2.Level_4.2016012800.0.0.tar.gz'
Content type 'application/x-gzip' length 53856803 bytes (51.4 MB)
如果你的网络对这个FireBrowse不友好,上面的代码可能会运行错误,或者好几个小时才能下载几十个M的小文件。
可以看到同时把CopyNumber_Gistic2数据下载给我了,而我想要的somatic mutation信息在 Mutation_Packager_Calls 里面,临床信息当然是必须的咯。
其实就是根据参数拼接了两个URL而已,原理非常简单,但是它有个好处就是,不仅仅是下载了数据,而且返回了包含这些数据的S4对象。
还有很多其它参数可以控制下载的数据类型,包括:
如果你的目标是探索其它ngs组学数据,当然是可以自己变化参数慢慢熟悉它们。
在R语言里面,S4对象是一个门槛,熟悉它基本上很多教程就可以无师自通了,比如单细胞流程的 seurat,我们的这个FireBrowse下载到的S4对象 也不例外。
load(file='brcaData.RTCGAToolbox.Rdata')
brcaData
## BRCA FirehoseData objectStandard run date: 20160128
## Analysis running date: 20160128
## Available data types:
## clinical: A data frame of phenotype data, dim: 1097 x 18
## GISTIC: A FirehoseGISTIC for copy number data
## Mutation: A data.frame, dim: 90490 x 67
## To export data, use the 'getData' function.
可以看到这个FireBrowse下载到的S4对象包含了3种数据,分别是临床信息,somatic的mutation,以及拷贝数变异信息。
这里需要用包定义好的函数来从S4对象里面获取数据,就是biocExtract函数:
biocExtract(object, type = c("clinical", "RNASeqGene", "miRNASeqGene",
"RNASeq2GeneNorm", "CNASNP", "CNVSNP", "CNASeq", "CNACGH", "Methylation",
"Mutation", "mRNAArray", "miRNAArray", "RPPAArray", "GISTIC", "GISTICA",
"GISTICT"))
这个函数可以提前几乎所有的信息,当然前提是我们确实下载了那些信息哦。
首先提取临床信息:
clinicData=biocExtract(brcaData,'clinical')
## working on: clinical
colnames(clinicData)
## [1] "Composite Element REF"
## [2] "years_to_birth"
## [3] "vital_status"
## [4] "days_to_death"
## [5] "days_to_last_followup"
## [6] "tumor_tissue_site"
## [7] "pathologic_stage"
## [8] "pathology_T_stage"
## [9] "pathology_N_stage"
## [10] "pathology_M_stage"
## [11] "gender"
## [12] "date_of_initial_pathologic_diagnosis"
## [13] "days_to_last_known_alive"
## [14] "radiation_therapy"
## [15] "histological_type"
## [16] "number_of_lymph_nodes"
## [17] "race"
## [18] "ethnicity"
DT::datatable(clinicData,
extensions = 'FixedColumns',
options = list(
#dom = 't',
scrollX = TRUE,
fixedColumns = TRUE
))
然后获取突变信息,如下所示:
mutationData = biocExtract(brcaData,'Mutation')
## working on: Mutation
length(mutationData@assays)
## [1] 993
class(mutationData@assays[[1]])
## [1] "GRanges"
## attr(,"package")
## [1] "GenomicRanges"
是一个GRanges 对象, 可以就按照 GRanges的操作手册来探索它。
RTCGAToolbox 提供了5个基本的数据分析工具:
对应的5个函数如下所示:
getMutationRate , Make a table for mutation rate of each gene in the cohort
getReport, Draws a circle plot into working directory
getSurvival ,Perform survival analysis based on gene expression data
getDiffExpressedGenes ,Perform differential gene expression analysis for mRNA expression data.
getCNGECorrelation, Perform correlation analysis between gene expression and copy number data
因为我们上面的示例代码并没有下载表达矩阵,所以基因表达量的差异分析和相关性分析,针对表达信息的生存分析没办法做,以及针对差异分析结果的可视化报告都是无法运行的。
getDiffExpressedGenes(brcaData)
corRes <- getCNGECorrelation(brcaData)
corRes
showResults(corRes[[1]])
可以运行的就是看看突变率,还有针对临床资料的生存分析了。
mt=getMutationRate(brcaData)
head(mt)
## Genes MutationRatio
## ACPP ACPP 0.006042296
## ALG13 ALG13 0.007049345
## AMY2A AMY2A 0.006042296
## B4GALT1 B4GALT1 0.004028197
## CARD6 CARD6 0.009063444
## CCDC114 CCDC114 0.005035247
tail(mt[order(mt$MutationRatio),])
## Genes MutationRatio
## GATA3 GATA3 0.09969789
## MUC16 MUC16 0.10070493
## CDH1 CDH1 0.11581067
## TTN TTN 0.19436052
## TP53 TP53 0.31117825
## PIK3CA PIK3CA 0.32628399
可以看到在我们的这个TCGA数据库的BRCA队列里面,突变频率比较高的基因是 TP53和PIK3CA,也确实是乳腺癌的明星基因。
看看这个TCGA数据库的BRCA队列生存情况:
head(clinicData[,1:4])
## Composite Element REF years_to_birth vital_status
## tcga.5l.aat0 value 42 0
## tcga.5l.aat1 value 63 0
## tcga.a1.a0sp value 40 0
## tcga.a2.a04v value 39 1
## tcga.a2.a04y value 53 0
## tcga.a2.a0cq value 62 0
## days_to_death
## tcga.5l.aat0 <NA>
## tcga.5l.aat1 <NA>
## tcga.a1.a0sp <NA>
## tcga.a2.a04v 1920
## tcga.a2.a04y <NA>
## tcga.a2.a0cq <NA>
survData <- data.frame(Samples=rownames(clinicData),
Time=as.numeric(clinicData[,4]),
Censor=as.numeric(clinicData[,3])
)
library(survival)
table(survData$Censor)
##
## 0 1
## 945 152
attach(survData)
my.surv <- Surv(Time,Censor)
kmfit1 <- survfit(my.surv~1)
plot(kmfit1)
可以看到,死亡这样的结局时间发生的概率并不多哦,这个乳腺癌的生存情况比较好,一般来说我们的生存分析是需要分组后去做检验,这样就知道我们的分组是否有统计学意义。我在生信技能树多次分享过生存分析的细节;
前面的示例代码里面,就可以根据 突变频率比较高的基因是 TP53和PIK3CA对 这个TCGA数据库的BRCA队列 进行分组,然后统计学检验,当然也可以联合两个基因突变再看生存效应啦。
两个优点:
最大的缺点就是只能下载全部的基因信息,这样下载速度肯定不能很快,而很多时候,我们只是对感兴趣的基因想探索一下而已。
另外,既然它是broad的FireBrowse包装盒,那我们当然可以直接使用broad的FireBrowse工具咯,命令行版本哈!
就需要大家熟练掌握R语言,比如把上面的基础绘图全部使用ggplot语法重新绘制,并且让它更美观,甚至惊艳。
以及需要掌握TCGA数据库及其背后的癌症数据集的背景知识了,这些都是需要时间积累的,不能一蹴而就。
强烈建议你推荐给身边的博士后以及年轻生物学PI,多一点数据认知,让他们的科研上一个台阶:
https://mp.weixin.qq.com/s/VvJ4Lqnc-YrUD3ZWHdSAdA