TCGA的28篇教程- 使用R语言的RTCGAToolbox包获取TCGA数据
前些天被TCGA的终结新闻刷屏，但是一直比较忙，还没来得及仔细研读，但是笔记本躺着的一些TCGA教程快发霉了，借此契机好好整理一下吧，预计28篇教程！
——jimmy
往期目录如下：
TCGA数据源
背景知识
了解并获取FireBrowse的数据
了解从FireBrowse下载到的S4对象
5大分析方法
优缺点分析
众所周知,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: clinicalcolnames(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: Mutationlength(mutationData@assays)## [1] 993class(mutationData@assays[[1]])## [1] "GRanges"
## attr(,"package")
## [1] "GenomicRanges"
对于 GRanges 对象，就按照 GRanges的操作手册来即可
5大分析方法
RTCGAToolbox 提供了5个基本的数据分析工具：
1. 差异表达分析（比较肿瘤和正常组织的基因表达量），根据不同的平台（RNA-Seq或Microarray），自动选择适合的工具
2. 拷贝数和基因表达量的相关性分析
3. 基因突变率分析
4. 生存分析
5. 可视化报告
没有下载表达矩阵，所以基因表达量的差异分析和相关性分析，针对表达信息的生存分析没办法做，以及针对差异分析结果的可视化报告都是无法运行的
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.005035247tail(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 152attach(survData)
my.surv <- Surv(Time,Censor)
kmfit1 <- survfit(my.surv~1)
plot(kmfit1)
detach(survData)
接下来可以根据各个基因的突变信息，拷贝数变异信息，以及其它临床信息把病人进行分组，进行上次分析检验，cox回归分析等等。
优缺点分析
两个优点：
1. 通过一个函数自动完成所有数据下载的工作（包括下载，解压，读入文件，删除压缩文件），极为方便
1. 读入的TCGA数据被自动封装在一个S4的对象中，我们可以通过各种接口来轻松的访问它内部的数据，一个有条理的数据组织结构可以大大提高程序的可读性和可维护性
最大的缺点就是只能下载全部的基因信息，这样下载速度肯定不能很快，而很多时候，我们只是对感兴趣的基因想探索一下而已。
后面的2~5期我们就会讲一下如何探索感兴趣的基因哈！
那么你当然可以直接使用broad的FireBrowse工具咯，命令行版本哈！
点击下面的阅读原文可以直达我以前在生信技能树发布的命令行工具教程，就不占用正文篇幅了哈！
纽约 时间 比加 州时间早三个小时，
New York is 3 hours ahead of California
但加州时间并没有变慢。
but it does not make California slow.
有人22岁就毕业了，
Someone graduated at the age of 22,
但等了五年才找到稳定的工作！
but waited 5 years before securing a good job!
有人25岁就当上CEO，
Someone became a CEO at 25,
却在50岁去世。
and died at 50.
也有人直到50岁才当上CEO，
While another became a CEO at 50,
然后活到90岁。
and lived to 90 years.
有人单身，
Someone is still single,
同时也有人已婚，
while someone else got married,
也有人又恢复单身了。
someone is single again.
欧 巴马 55岁就退休，
Obama retires at 55,
川普70岁才开始当总统 。
but Trump starts at 70.
世上每个人本来就有自己的发展时区。
Absolutely everyone in this world works based on their Time Zone.
身边有些人看似走在你前面，
People around you might seem to go ahead of you,
也有人看似走在你后面。
some might seem to be behind you.
但其实每个人在自己的时区有自己的步程。
But everyone is running their own RACE, in their own TIME.
不用嫉妒或嘲笑他们。
Don’t envy them or mock them.
他们都在自己的时区里，你也是！
They are in their TIME ZONE, and you are in yours!
生命就是等待正确的行动时机。
Life is about waiting for the right moment to act.
所以，放轻松。
So, RELAX.
你没有落后。
You’re not LATE.
你没有领先。
You’re not EARLY.
在命运为你安排的属于自己的时区里，一切都准时。
You are very much ON TIME, and in your TIME ZONE Destiny set up for you.