众所周知，TCGA数据库是目前最综合全面的癌症病人相关组学数据库，包括的多组学数据有：
- DNA Sequencing (WGS/WES)
- mRNA/miRNA Sequencing
- Protein Expression Array
- Array-based Expression
- DNA Methylation Array
- Copy Number Array
知名的肿瘤研究机构都有着自己的TCGA数据库探索工具，比如：
- Broad Institute FireBrowse portal, The Broad Institute
- cBioPortal for Cancer Genomics, Memorial Sloan-Kettering Cancer Center
- TCGA Batch Effects, MD Anderson Cancer Center
- Regulome Explorer, Institute for Systems Biology
- Next-Generation Clustered Heat Maps, MD Anderson Cancer Center
其中Broad研究所的就是FireBrowse，主页在：http://www.firebrowse.org/
这个网页工具当然是非常强大，不过咱们生信工程师喜欢的仍然是编程语言，所以有一个RTCGAToolbox的R包可以帮助我们通过代码来玩转它的网页工具，安装它超级简单：
代码语言：javascript复制if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("RTCGAToolbox")
RTCGAToolbox有什么
代码语言：javascript复制library(RTCGAToolbox)
# 查看哪些癌症数据可以下载
getFirehoseDatasets()
getFirehoseRunningDates( )
getFirehoseAnalyzeDates( )
# 可以看到工具在 20160128 就停止更新了
可以看到不同的癌症数据集，不同的更新时间，如下所示：
代码语言：javascript复制> 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里面的所有的样本了。
接下来下载所需要的数据，这里以乳腺癌为例，数据下载完后会直接放在你的工作目录，不同地方下载的速度不一样。
代码语言：javascript复制## 下载数据，需要选择癌症种类，数据分析时间，还有数据的种类
options(timeout=10000)
# 一般来说，我们会选择最新的数据，工具在 20160128 就停止更新了
brcaData = getFirehoseData (dataset="BRCA", runDate="20160128",
forceDownload = TRUE,
clinical=TRUE, Mutation=TRUE)
save(brcaData,file='brcaData.RTCGAToolbox.Rdata')
这里测试了临床信息和突变信息的数据下载，因为它们比较小，所以下载速度会很快，这里下载的数据包括：
代码语言：javascript复制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对象。
还有很多其它参数可以控制下载的数据类型，包括：
- clinical - Get the clinical data slot
- RNASeqGene - RNASeqGene
- RNASeq2GeneNorm - Normalized
- miRNASeqGene - micro RNA SeqGene
- CNASNP - Copy Number Alteration
- CNVSNP - Copy Number Variation
- CNASeq - Copy Number Alteration
- CNACGH - Copy Number Alteration
- Methylation - Methylation
- mRNAArray - Messenger RNA
- miRNAArray - micro RNA
- RPPAArray - Reverse Phase Protein Array
- Mutation - Mutations
- GISTICA - GISTIC v2 (’AllByGene’ only)
- GISTICT - GISTIC v2 (’ThresholdedByGene’ only)
- GISTIC - GISTIC v2 scores and probabilities (both)
如果你的目标是探索其它ngs组学数据，当然是可以自己变化参数慢慢熟悉它们。
了解从FireBrowse下载到的S4对象
在R语言里面，S4对象是一个门槛，熟悉它基本上很多教程就可以无师自通了，比如单细胞流程的 seurat，我们的这个FireBrowse下载到的S4对象 也不例外。
代码语言：javascript复制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函数：
代码语言：javascript复制biocExtract(object, type = c("clinical", "RNASeqGene", "miRNASeqGene",
"RNASeq2GeneNorm", "CNASNP", "CNVSNP", "CNASeq", "CNACGH", "Methylation",
"Mutation", "mRNAArray", "miRNAArray", "RPPAArray", "GISTIC", "GISTICA",
"GISTICT"))
这个函数可以提前几乎所有的信息，当然前提是我们确实下载了那些信息哦。
首先提取临床信息：
代码语言：javascript复制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
))
然后获取突变信息，如下所示：
代码语言：javascript复制mutationData = biocExtract(brcaData,'Mutation')
## working on: Mutation
length(mutationData@assays)
## [1] 993
class(mutationData@assays[[1]])
## [1] "GRanges"
## attr(,"package")
## [1] "GenomicRanges"
是一个GRanges 对象， 可以就按照 GRanges的操作手册来探索它。
5大分析方法
RTCGAToolbox 提供了5个基本的数据分析工具：
- 1. 差异表达分析（比较肿瘤和正常组织的基因表达量），根据不同的平台（RNA-Seq或Microarray），自动选择适合的工具
- 2. 拷贝数和基因表达量的相关性分析
- 3. 基因突变率分析
- 4. 生存分析
- 5. 可视化报告
对应的5个函数如下所示：
代码语言：javascript复制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
因为我们上面的示例代码并没有下载表达矩阵，所以基因表达量的差异分析和相关性分析，针对表达信息的生存分析没办法做，以及针对差异分析结果的可视化报告都是无法运行的。
代码语言：javascript复制getDiffExpressedGenes(brcaData)
corRes <- getCNGECorrelation(brcaData)
corRes
showResults(corRes[[1]])
可以运行的就是看看突变率，还有针对临床资料的生存分析了。
代码语言：javascript复制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队列生存情况：
代码语言：javascript复制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)
可以看到，死亡这样的结局时间发生的概率并不多哦，这个乳腺癌的生存情况比较好，一般来说我们的生存分析是需要分组后去做检验，这样就知道我们的分组是否有统计学意义。我在生信技能树多次分享过生存分析的细节；
- 人人都可以学会生存分析（学徒数据挖掘）
- 学徒数据挖掘之谁说生存分析一定要按照表达量中位值或者平均值分组呢？
- 基因表达量高低分组的cox和连续变量cox回归计算的HR值差异太大?
- 学徒作业-两个基因突变联合看生存效应
前面的示例代码里面，就可以根据 突变频率比较高的基因是 TP53和PIK3CA对 这个TCGA数据库的BRCA队列 进行分组，然后统计学检验，当然也可以联合两个基因突变再看生存效应啦。
优缺点分析
两个优点：
- 通过一个函数自动完成所有数据下载的工作（包括下载，解压，读入文件，删除压缩文件），极为方便
- 读入的TCGA数据被自动封装在一个S4的对象中，我们可以通过各种接口来轻松的访问它内部的数据，一个有条理的数据组织结构可以大大提高程序的可读性和可维护性
最大的缺点就是只能下载全部的基因信息，这样下载速度肯定不能很快，而很多时候，我们只是对感兴趣的基因想探索一下而已。
另外，既然它是broad的FireBrowse包装盒，那我们当然可以直接使用broad的FireBrowse工具咯，命令行版本哈！
更多玩法
就需要大家熟练掌握R语言，比如把上面的基础绘图全部使用ggplot语法重新绘制，并且让它更美观，甚至惊艳。
以及需要掌握TCGA数据库及其背后的癌症数据集的背景知识了，这些都是需要时间积累的，不能一蹴而就。
