通过R包RTCGAToolbox链接FireBrowse来探索TCGA等公共数据

生信技能树

发布于 2022-07-26 10:05:25

3750

发布于 2022-07-26 10:05:25

文章被收录于专栏：生信技能树

众所周知，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包可以帮助我们通过代码来玩转它的网页工具，安装它超级简单：

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("RTCGAToolbox")

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对象。

还有很多其它参数可以控制下载的数据类型，包括：

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对象也不例外。

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的操作手册来探索它。

5大分析方法

RTCGAToolbox 提供了5个基本的数据分析工具：

\1. 差异表达分析（比较肿瘤和正常组织的基因表达量），根据不同的平台（RNA-Seq或Microarray），自动选择适合的工具
\2. 拷贝数和基因表达量的相关性分析
\3. 基因突变率分析
\4. 生存分析
\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队列进行分组，然后统计学检验，当然也可以联合两个基因突变再看生存效应啦。

优缺点分析

两个优点：

通过一个函数自动完成所有数据下载的工作（包括下载，解压，读入文件，删除压缩文件），极为方便
读入的TCGA数据被自动封装在一个S4的对象中，我们可以通过各种接口来轻松的访问它内部的数据，一个有条理的数据组织结构可以大大提高程序的可读性和可维护性

最大的缺点就是只能下载全部的基因信息，这样下载速度肯定不能很快，而很多时候，我们只是对感兴趣的基因想探索一下而已。

另外，既然它是broad的FireBrowse包装盒，那我们当然可以直接使用broad的FireBrowse工具咯，命令行版本哈！