使用R语言的RTCGAToolbox包获取TCGA数据

作者: 因地制宜的生信达人 | 来源:发表于2018-07-01 20:47 被阅读249次

TCGA数据源

众所周知,TCGA数据库是目前最综合全面的癌症病人相关组学数据库，包括的测序数据有：

DNA Sequencing
miRNA Sequencing
Protein Expression
mRNA Sequencing
Total RNA Sequencing
Array-based Expression
DNA Methylation
Copy Number

知名的肿瘤研究机构都有着自己的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

其中FireBrowse被包装到R包RTCGAToolbox里面： http://bioconductor.org/packages/release/bioc/manuals/RTCGAToolbox/man/RTCGAToolbox.pdf

这里就介绍如何使用R语言的 RTCGAToolbox 包来获取任意TCGA数据吧。该包与2014年发表在plos one杂志；RTCGAToolbox: A New Tool for Exporting TCGA Firehose Data - PLOS

其实Firehose官方就提供过非常方便的命令行工具来根据他们的数据存放规则非常方便的获取数据，外网速度一般是10M/S，非常好用。

背景知识

TCGA上的数据量庞大，数据种类丰富，分析方法复杂，并不是所有人都能轻松下载、管理和分析这些数据。对于大部分研究人员来说，从如此海量的原始测序数据开始分析是不可行也是不必要的。实际上，我们可以下载经过预处理后的数据（pre-processed data），不仅数据量会小很多，分析起来也更快、更可靠。Broad institute开发的Firehose就能够提供这样的数据。

虽然Firehose为我们做好前期的处理工作，但在R里面还缺一个“搜索引擎”，所以RTCGAToolbox就应运而生。

RTCGAToolbox是Bioconductor上的一个软件包，它的作用就是查询、下载和组织TCGA Firehose的数据，还提供一些简单的数据分析和可视化工具。除此之外，下载好的数据也可以很方便的导入到Bioconductor的其他分析流程中。对于R用户来说，所有的TCGA数据分析工作（从数据下载一直到可视化图表）都可在一个pipeline中完成，能够极大地提高工作效率。

了解并获取FireBrowse的数据

#包下载
#source("https://bioconductor.org/biocLite.R")
#biocLite("RTCGAToolbox")
#加载包
library(RTCGAToolbox)
#哪些癌症数据可以下载
getFirehoseDatasets()

##  [1] "ACC"      "BLCA"     "BRCA"     "CESC"     "CHOL"     "COADREAD"
##  [7] "COAD"     "DLBC"     "ESCA"     "FPPP"     "GBMLGG"   "GBM"     
## [13] "HNSC"     "KICH"     "KIPAN"    "KIRC"     "KIRP"     "LAML"    
## [19] "LGG"      "LIHC"     "LUAD"     "LUSC"     "MESO"     "OV"      
## [25] "PAAD"     "PCPG"     "PRAD"     "READ"     "SARC"     "SKCM"    
## [31] "STAD"     "STES"     "TGCT"     "THCA"     "THYM"     "UCEC"    
## [37] "UCS"      "UVM"

#数据库中更新时间
getFirehoseRunningDates()

##  [1] "20160128" "20151101" "20150821" "20150601" "20150402" "20150204"
##  [7] "20141206" "20141017" "20140902" "20140715" "20140614" "20140518"
## [13] "20140416" "20140316" "20140215" "20140115" "20131210" "20131114"
## [19] "20131010" "20130923" "20130809" "20130715" "20130623" "20130606"
## [25] "20130523" "20130508" "20130421" "20130406" "20130326" "20130309"
## [31] "20130222" "20130203" "20130116" "20121221" "20121206" "20121114"
## [37] "20121102" "20121024" "20121020" "20121018" "20121004" "20120913"
## [43] "20120825" "20120804" "20120725" "20120707" "20120623" "20120606"
## [49] "20120525" "20120515" "20120425" "20120412" "20120321" "20120306"
## [55] "20120217" "20120124" "20120110" "20111230" "20111206" "20111128"
## [61] "20111115" "20111026"

getFirehoseAnalyzeDates()

##  [1] "20160128" "20150821" "20150402" "20141017" "20140715" "20140416"
##  [7] "20140115" "20130923" "20130523" "20130421" "20130326" "20130222"
## [13] "20130116" "20121221" "20121024" "20120913" "20120825" "20120725"
## [19] "20120623" "20120525" "20120425" "20120321" "20120217" "20120124"
## [25] "20111230" "20111128" "20111026" "20110921" "20110728" "20110525"
## [31] "20110421" "20110327" "20110217" "20110114" "20101223"

可以看到TCGA的各种癌症都在列表中了，这里用的是简称，比如BRCA就是乳腺癌。

而第二个不同的时间，指的是TCGA数据库在发展过程中样本量的增加，而FireBrowse是按照时间来定期运行程序处理数据的，所以一般来说用最新版的结果，就会涵盖TCGA里面的所有的样本了。

接下来下载所需要的数据，这里以乳腺癌为例，数据下载完后会直接放在你的工作目录，不同地方下载的速度不一样。

## 下载数据，需要选择癌症种类，数据分析时间，还有数据的种类
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)

可以看到同时把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)

了解从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.

可以看到包含了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），自动选择适合的工具
1. 拷贝数和基因表达量的相关性分析
1. 基因突变率分析
1. 生存分析
1. 可视化报告

没有下载表达矩阵，所以基因表达量的差异分析和相关性分析，针对表达信息的生存分析没办法做，以及针对差异分析结果的可视化报告都是无法运行的

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

看看生存情况

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)

1.png

detach(survData)

接下来可以根据各个基因的突变信息，拷贝数变异信息，以及其它临床信息把病人进行分组，进行上次分析检验，cox回归分析等等。

优缺点分析

两个优点：

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

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

使用R语言的RTCGAToolbox包获取TCGA数据

TCGA数据源

背景知识

了解并获取FireBrowse的数据

了解从FireBrowse下载到的S4对象

5大分析方法

优缺点分析

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读

TCGA数据挖掘