R语言-使用RTCGAToolbox包获取TCGA数据(1)

RTCGAToolbox同样也是用于获取和处理TCGA数据的一个R包。很多癌症研究所都开发了自己的TCGA探索工具，比如Broad Institute开发了FireBrowse，而RTCGAToolbox就是用于获取和处理FireBrowse的R包。

这种R包的使用主流程一般为：下载数据 – 处理数据 – 可视化。下面依次来看本包在这三方面的安排：

下载数据

在下载数据之前，先安装RTCGAToolbox包，并载入所需工具包：

#安装RTCGAToolbox包
library("RTCGAToolbox")
#载入数据处理工具包
library("tidyverse")

RTCGAToolbox包中的getFirehoseData函数用于从Firehose下载TCGA数据。

函数的具体参数如下：

getFirehoseData(dataset, runDate = "20160128", gistic2Date = "20160128", 
        RNASeqGene = FALSE, clinical = TRUE,miRNASeqGene = FALSE, 
        RNASeq2GeneNorm = FALSE, CNASNP = FALSE,CNVSNP = FALSE, 
        CNASeq = FALSE, CNACGH = FALSE,    Methylation = FALSE, 
        Mutation = FALSE, mRNAArray = FALSE,miRNAArray = FALSE, 
        RPPAArray = FALSE, GISTIC = FALSE,
    RNAseqNorm = "raw_counts", RNAseq2Norm = "normalized_count",
    forceDownload = FALSE, destdir = ".", fileSizeLimit = 500,
    getUUIDs = FALSE, ...)

关键的参数为癌症类型dataset，时间runDate，GISTIC数据分析时间gistic2Data，以及要获取的数据类型（默认下载clinical数据）：RNAseq、clinical、miRNA、CNA、mutation等等，比如需要下载RNAseq数据，那么设置参数RNASeqGene =TRUE，需要CNA数据，那么设置参数CNASeq = TRUE。

dataset：癌症类型列表，使用函数getFirehoseDatasets获取；

runDate：数据分析时间，使用函数getFirehoseRunningDates获取，选最新就行；

gistic2Data：Gistic数据分析时间，使用函数getFirehoseAnalyzeDates获取，选最新就行；

其他可以下载的数据类型的含义为：

比如需要获取癌症ACC（Acute Myeloid Leukemia）在Firehose数据库中最新的数据，数据类别包含RNAseq、clinical、GISTIC以及mutation，如下：

FirehoseData <- getFirehoseData(dataset=getFirehoseDatasets()[1], 
                runDate=getFirehoseRunningDates()[1], gistic2Date=getFirehoseAnalyzeDates()[1], 
                RNASeqGene=TRUE,clinical=TRUE,GISTIC=TRUE,Mutation=TRUE)
FirehoseData
#ACC FirehoseData objectStandard run date: 20160128 
#Analysis running date: 20160128 
#Available data types: 
#  clinical: A data frame of phenotype data, dim:  92 x 18 
#  GISTIC: A FirehoseGISTIC for copy number data 
#  Mutation: A data.frame, dim:  20166 x 52 
#To export data, use the 'getData' function.

可以发现下载的数据只有clinical、GISTIC和Mutation，没有RNAseq。单独下载RNAseq数据，发现没有可用的数据类型，应该是ACC没有RNAseq数据。

getFirehoseData(dataset='ACC',runDate=getFirehoseRunningDates()[1], RNASeqGene=T,clinical=F)
#ACC FirehoseData objectStandard run date: 20160128 
#Analysis running date: NA 
#Available data types: 
#To export data, use the 'getData' function.

数据处理

提取数据

RTCGAToolbox里面有一个示例数据包RTCGASample，比较小，但是又包含了足够的数据：包括RNAseq、clinical、GISTIC以及Mutaion，以下就以它为例进行分析。

下载之后的TCGA数据是FirehoseData对象，可以使用三个函数（基本上是等价的）完成相应的数据提取：biocExtract、getData、selectType。

比如从示例数据RTCGASample中提取临床数据：

biocExtract(RTCGASample,type="clinical")%>%head
getData(RTCGASample,type= "clinical")%>%head
selectType(RTCGASample,"clinical")%>%head

结果输出是完全一样的。

type参数可用的选择有："clinical", "RNASeqGene", "miRNASeqGene", "RNASeq2GeneNorm", "CNASNP", "CNVSNP", "CNASeq", "CNACGH", "Methylation", "Mutation", "mRNAArray", "miRNAArray", "RPPAArray", "GISTIC", "GISTICA", "GISTICT", "GISTICP"。

数据处理

此外RTCGAToolbox还可以完成一些数据处理的工作：差异表达基因分析：getDiffExpressedGenes，基因表达和拷贝数变异的相关系数：getCNGECorrelation，突变率：getMutationRate。

差异表达基因分析如下：

t1 <- getDiffExpressedGenes(RTCGASample,logFC = 0.5) #绘图有热图等
#logFC默认值是2，筛选出来的基因太少，只有4个，此处为了叙述方便，改为0.5

同时会绘制热图：

t1的内容是一个S4类对象，可以使用@提取相关内容，比如t1[[1]]@Toptable就是常见的差异表达分析表（如下图），也可以使用函数showResults来提取数据：

拷贝数和基因表达的关系

getCNGECorrelation函数可以获得拷贝数变化和基因表达的相关关系：

corRes = getCNGECorrelation(RTCGASample,adj.pval = 1,raw.pval = 1)
corRes[[1]]@Correlations%>%head
#showResults(corRes[[1]])的输出结果是一样的，如下：
> corRes[[1]]@Correlations%>%head
  GeneSymbol         Cor adj.p.value    p.value
1      LYZL1 -0.14654634   0.9317555 0.36687874
2     SEPHS1 -0.37693816   0.6287624 0.01650501
3   FLJ45983  0.14150035   0.9428853 0.38378545
4     TUBAL3 -0.06913254   0.9712442 0.67165645
5      PRKCQ -0.18959020   0.8688328 0.24131149
6      WDR37 -0.20867359   0.8688328 0.19628861

getMutationRate函数的道理是一样的，不过多叙述了。

可视化

本软件有两个主要的可视化函数getReport和getSurvival。

getReport

按照预计，getReport可以直接绘制一个Circle图，包含差异基因的logFC变化、拷贝数变化及突变，然而实际上这个函数已经无法使用了。

原本的操作为：

require("RCircos")
require("Homo.sapiens")
# 基因坐标
Homo.sapiens %>%genes(columns="SYMBOL")%>% as.data.frame() %>% 
  mutate(SYMBOL_1=as.character(SYMBOL))%>%  dplyr::select(7,1,5,2,3) %>%   
  distinct(SYMBOL_1,.keep_all = T)%>% drop_na() %>%
  dplyr::rename(SYMBOL=SYMBOL_1)-> locations
rownames(locations) <- locations[,1]
# 差异基因分析
t1 <- getDiffExpressedGenes(RTCGASample,logFC = 1)
getReport(dataObject=RTCGASample,DGEResult1=t1[[1]],geneLocations=locations)

然后就会在工作目录下创建一个pdf文件，里面就是一个绘制好的Circle图。然而现在直接运行的话会报错：

Please use RCircos.Set.Core.Components() instead.

就是说getReport引用的RCircle包的一个函数修改了，但是getReport并没有更新。详情可以见https://www.jianshu.com/p/9b1559c3a5a7，里面叙述的很清楚。

那么是不是按照这篇文章修改一下getReport函数就可以了呢？好像由于RCircos又有更新，虽然这一次可以创建一个PDF文件，但是这个PDF文件无法打开，所以暂时还是放弃这个函数吧。

getSurvival

绘制KM生存曲线：

data(RTCGASample)

RTCGASample %>% getData("clinical") %>% select(3:5) %>% rownames_to_column()%>%
  mutate(Time=ifelse(is.na(daystolastfollowup),daystodeath,daystolastfollowup))%>%
  rename(Samples=rowname,Cencor=vitalstatus) %>%select(Samples,Time,Cencor) -> survData

getSurvival(dataObject=RTCGASample,geneSymbols=c("TAP2"),sampleTimeCensor=survData)
#TAP2基因是差异表达基因里面最显著且变化最大的基因。

结果如下，生存率和TAP2看着还有有关系的，分的很开。绘图比较原始，随便看看就行。