GDCRNATools--一个R包就能解决TCGA数据处理和可视化！
导语
GUIDE ╲
GDCRNATools是一个易于使用的用于整合GDC中lncRNA、mRNA和miRNA数据的R/Bioconductor软件包。
背景介绍
TCGA数据库作为癌症研究的首选公共数据库，整合了各种癌症的多组学数据。基因组数据共享数据库（GDC）维护着来自美国国家癌症研究所（NCI）计划的标准化基因组，临床和样本数据，包括TCGA和TARGET，它也接受来自非NCI支持的癌症研究计划的高质量数据集，例如来自Foundation Medicine的基因组数据。
GDCRNATools是一个R软件包，它提供了一个易于使用且全面的方法，用于下载，分析和可视化GDC中的RNA表达数据，重点在于解读癌症中与lncRNA-mRNA相关的ceRNA调控网络。同时还可以进行多种分析，包括基于limma、edgeR和DESeq2识别差异基因，基于clusterProfile和DO包进行的功能富集分析（GO、KEGG、DO），以及单因素生存分析和相关性分析等等。除了一些常规的可视化功能，如火山图、条形图、K-M图，还开发了一些简单的shiny应用，使用户可以在本地网页上的可视化结果。
GDCRNATools的安装
##使用bioconductor安装R包
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("GDCRNATools")
library(GDCRNATools)GDCRNATools的使用
在GDCRNATools中，内置了一些功能，供用户下载和处理GDC数据。用户还可以使用自己的数据，这些数据来源于UCSC Xena GDC hub，TCGAbiolinks或TCGA-Assembler等。
在这里，我们使用一个小的数据集来进行ceRNAs网络分析的最基本步骤。
01
数据准备工作
1、 HTSeq-Counts data的标准化
library(DT)
### 下载RNAcounts数据
data(rnaCounts)
### 下载miRNAcounts数据
data(mirCounts)
####### 标准化函数gdcVoomNormalization
rnaExpr <- gdcVoomNormalization(counts = rnaCounts, filter = FALSE)
mirExpr <- gdcVoomNormalization(counts = mirCounts, filter = FALSE)2、 Parse metadata
###生成胆管癌的metadata矩阵
metaMatrix.RNA <- gdcParseMetadata(project.id = 'TCGA-CHOL',
data.type = 'RNAseq',
write.meta = FALSE)
###对矩阵样本进行筛选
metaMatrix.RNA <- gdcFilterDuplicate(metaMatrix.RNA)
metaMatrix.RNA <- gdcFilterSampleType(metaMatrix.RNA)
##展示前6行
datatable(as.data.frame(metaMatrix.RNA[1:6,]), extensions = 'Scroller',
options = list(scrollX = TRUE, deferRender = TRUE, scroller = TRUE))
02
ceRNA网络分析
1、差异表达基因的识别
##输入counts数据，按sample_type对样本分类,选择计算差异表达的方法
DEGAll <- gdcDEAnalysis(counts = rnaCounts,
group = metaMatrix.RNA$sample_type,
comparison = 'PrimaryTumor-SolidTissueNormal',
method = 'limma')
##展示筛选结果
datatable(as.data.frame(DEGAll),
options = list(scrollX = TRUE, pageLength = 5))
### 筛选差异表达的lncRNA
deLNC <- gdcDEReport(deg = DEGAll, gene.type = 'long_non_coding')
### 筛选差异表达的编码RNA
dePC <- gdcDEReport(deg = DEGAll, gene.type = 'protein_coding')2、对于差异基因进行ceRNA网络分析
##生成ceRNA网络，RNA靶标来源starbase
ceOutput <- gdcCEAnalysis(lnc = rownames(deLNC),
pc = rownames(dePC),
lnc.targets = 'starBase',
pc.targets = 'starBase',
rna.expr = rnaExpr,
mir.expr = mirExpr)
##展示结果
datatable(as.data.frame(ceOutput),
options = list(scrollX = TRUE, pageLength = 5))
3、将ceRNA网络输出到cytoscape
##筛选p值显著的结果生成网络
ceOutput2 <- ceOutput[ceOutput$hyperPValue<0.01
& ceOutput$corPValue<0.01 & ceOutput$regSim != 0,]
### 生成边的信息
edges <- gdcExportNetwork(ceNetwork = ceOutput2, net = 'edges')
datatable(as.data.frame(edges),
options = list(scrollX = TRUE, pageLength = 5))
### 生成节点的信息
nodes <- gdcExportNetwork(ceNetwork = ceOutput2, net = 'nodes')
datatable(as.data.frame(nodes),
options = list(scrollX = TRUE, pageLength = 5))将节点和边的信息保存下来就可以输入到cytoscape画图啦！
数据分析实例
在这部分，我们将使用TCGA-CHOL整个数据集为例，详细说明GDCRNATools的使用方法。
01
RNA和miRNA表达数据的下载
在GDCRNATools中，提供了gdc-client方法通过指定data.type和project.id参数自动下载数据。
直接运行如下代码可能会报错，是因为R包内的gdc-client版本过低，我们需要自己下载新版本的gdc-client，解压后替换掉原来的gdc-client_v1.3.0,下载地址是：https://gdc.cancer.gov/access-data/gdc-data-transfer-tool
project <- 'TCGA-CHOL'
rnadir <- paste(project, 'RNAseq', sep='/')
mirdir <- paste(project, 'miRNAs', sep='/')
####### 下载RNAseq data,directory参数要写具体的下载路径（英文目录）
gdcRNADownload(project.id = 'TCGA-CHOL',
data.type = 'RNAseq',
write.manifest = FALSE,
method = 'gdc-client',
directory = rnadir)
####### 下载mature miRNA data
gdcRNADownload(project.id = 'TCGA-CHOL',
data.type = 'miRNAs',
write.manifest = FALSE,
method = 'gdc-client',
directory = mirdir)
####### 下载 clinical data
clinicaldir <- paste(project, 'Clinical', sep='/')
gdcClinicalDownload(project.id = 'TCGA-CHOL',
write.manifest = FALSE,
method = 'gdc-client',
directory = mirdir)02
metadata的处理
在gdcParseMetadata（）函数中指定project.id和data.type来获取清单文件中的数据信息来解析metadata。对患者的组织和基本临床信息（例如年龄，分期和性别等）进行数据分析。通过gdcFilterDuplicate（）去掉重复样本，使用gdcFilterSampleType（）过滤掉既不是原发性肿瘤（01）也不是实体组织的正常样本（11）。
#### Parse RNAseq metadata
metaMatrix.RNA <- gdcParseMetadata(project.id = 'TCGA-CHOL',
data.type = 'RNAseq',
write.meta = FALSE)
#### 去除 RNAseq metadata的重复样本
metaMatrix.RNA <- gdcFilterDuplicate(metaMatrix.RNA)
#### 去除 RNAseq metadata既不是原发性肿瘤也不是实体组织的正常样本
metaMatrix.RNA <- gdcFilterSampleType(metaMatrix.RNA)
#### Parse miRNAs metadata #######
metaMatrix.MIR <- gdcParseMetadata(project.id = 'TCGA-CHOL',
data.type = 'miRNAs',
write.meta = FALSE)
#### 去除miRNAs metadata的重复样本
metaMatrix.MIR <- gdcFilterDuplicate(metaMatrix.MIR)
#### 去除 miRNAs metadata既不是原发性肿瘤也不是实体组织的正常样本
metaMatrix.MIR <- gdcFilterSampleType(metaMatrix.MIR)03
合并原始counts数据
gdcRNAMerge（）将RNAseq的原始数据合并到单个表达矩阵中，其中行为Ensembl id，列为样本。如果不同样本的数据位于单独的文件夹中，可以指定organized = FALSE，否则，指定organized = TRUE。
####### 合并 RNAseq data
rnaCounts <- gdcRNAMerge(metadata = metaMatrix.RNA,
path = 'C:\\Users\\Administrator\\Desktop\\rnaseq', # the folder in which the data stored
organized = FALSE, # if the data are in separate folders
data.type = 'RNAseq')
####### 合并 miRNAs data
mirCounts <- gdcRNAMerge(metadata = metaMatrix.MIR,
path = 'C:\\Users\\Administrator\\Desktop\\mirna', # the folder in which the data stored
organized = FALSE, # if the data are in separate folders
data.type = 'miRNAs')
####### 合并 clinical data
clinicalDa <- gdcClinicalMerge(path = 'C:\\Users\\Administrator\\Desktop\\clinicaldir', key.info = TRUE)
clinicalDa[1:6,5:10]04
TMM标准化和Voom转换
通过运行gdcVoomNormalization（）函数，原始数据将通过edgeR中的TMM方法进行标准化，并通过limma中提供的voom方法进行进一步转换。默认情况下，低表达基因（超过一半的样本中logcpm <1）将被滤除。通过在gdcVoomNormalization（）中设置filter = TRUE可以保留所有基因。
### RNAseq data的标准化
rnaExpr <- gdcVoomNormalization(counts = rnaCounts, filter = FALSE)
### miRNAs data的标准化
mirExpr <- gdcVoomNormalization(counts = mirCounts, filter = FALSE)05
差异表达分析和可视化
##差异表达
DEGAll <- gdcDEAnalysis(counts = rnaCounts,
group = metaMatrix.RNA$sample_type,
comparison = 'PrimaryTumor-SolidTissueNormal',
method = 'limma')
data(DEGAll)
### All DEGs
deALL <- gdcDEReport(deg = DEGAll, gene.type = 'all')
### 筛选差异lncRNA
deLNC <- gdcDEReport(deg = DEGAll, gene.type = 'long_non_coding')
### 筛选差异编码基因
dePC <- gdcDEReport(deg = DEGAll, gene.type = 'protein_coding')火山图绘制
gdcVolcanoPlot(DEGAll)
发现只有半边结果，按照网上的参考改了一下函数，从新写了一个函数gdcVolcanoPlot2，这样画出来的结果就好看多啦！
gdcVolcanoPlot2<-function (deg.all, fc = 2, pval = 0.01, dotsize=0.8)
{
geneList <- deg.all
geneList$threshold <- c()
geneList$threshold[geneList$logFC > log(fc, 2) & geneList$FDR <
pval] <- 1
geneList$threshold[geneList$logFC >= -log(fc, 2) & geneList$logFC <=
log(fc, 2) | geneList$FDR >= pval] <- 2
geneList$threshold[geneList$logFC < -log(fc, 2) & geneList$FDR <
pval] <- 3
geneList$threshold <- as.factor(geneList$threshold)
lim <- max(max(geneList$logFC), abs(min(geneList$logFC))) +
0.5
volcano <- ggplot(data = geneList, aes(x = logFC,
y = -log10(FDR)))
volcano + geom_point(aes(color = threshold), alpha = 1,
size = dotsize) + xlab("log2(Fold Change)") + ylab("-log10(FDR)") +
scale_colour_manual(values = c("red", "black", "green3")) + xlim(c(-lim, lim)) +
geom_vline(xintercept = c(-log(fc, 2), log(fc, 2)), color = "darkgreen",
linetype = 3) + geom_hline(yintercept = -log(pval,
10), color = "darkgreen", linetype = 3) + theme_bw() +
theme(axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.border = element_rect(colour = "black"),
panel.background = element_blank()) + theme(legend.position = "none") +
theme(axis.text = element_text(size = 14), axis.title = element_text(size = 16))
}
##使用新函数绘图
gdcVolcanoPlot2(DEGAll, pval = 0.01)
上下调的柱状图绘制：
gdcBarPlot(deg = deALL, angle = 45, data.type = 'RNAseq')
热图的绘制：
degName = rownames(deALL)
gdcHeatmap(deg.id = degName, metadata = metaMatrix.RNA, rna.expr = rnaExpr)
06
ceRNA网络分析
gdcCEAnalysis（）还可以获取用户提供的miRNA-mRNA和miRNA-lncRNA相互作用数据集，而不仅仅是来源于starbase数据库
##使用来源于starbase的靶标数据集
ceOutput <- gdcCEAnalysis(lnc = rownames(deLNC),
pc = rownames(dePC),
lnc.targets = 'starBase',
pc.targets = 'starBase',
rna.expr = rnaExpr,
mir.expr = mirExpr)
##使用自行输入的靶标数据集
### 加载 miRNA-lncRNA互作数据
data(lncTarget)
### 加载miRNA-mRNA互作数据
data(pcTarget)
pcTarget[1:3]
ceOutput <- gdcCEAnalysis(lnc = rownames(deLNC),
pc = rownames(dePC),
lnc.targets = lncTarget,
pc.targets = pcTarget,
rna.expr = rnaExpr,
mir.expr = mirExpr)
###将节点和边的信息输出后，使用cytoscape绘图
ceOutput2 <- ceOutput[ceOutput$hyperPValue<0.01 &
ceOutput$corPValue<0.01 & ceOutput$regSim != 0,]
edges <- gdcExportNetwork(ceNetwork = ceOutput2, net = 'edges')
nodes <- gdcExportNetwork(ceNetwork = ceOutput2, net = 'nodes')
write.table(edges, file='edges.txt', sep='\t', quote=F)
write.table(nodes, file='nodes.txt', sep='\t', quote=F)
07
共表达分析
这部分可以使用gdcCorPlot（）函数绘制相关性图，也可以使用开发的shinyCorPlot（）函数制作一个交互式的本地网页，自行输入基因名来绘图。
gdcCorPlot(gene1 = 'ENSG00000251165',
gene2 = 'ENSG00000091831',
rna.expr = rnaExpr,
metadata = metaMatrix.RNA)
shinyCorPlot(gene1 = rownames(deLNC),
gene2 = rownames(dePC),
rna.expr = rnaExpr,
metadata = metaMatrix.RNA)
08
K-M曲线的绘制
K-M曲线的绘制同样也有两种方式，通过gdcKMPlot（）绘制基本的生存曲线，还有通过shinyKMPlot（）绘制交互式的网页。
gdcKMPlot(gene = 'ENSG00000136193',
rna.expr = rnaExpr,
metadata = metaMatrix.RNA,
sep = 'median')
shinyKMPlot(gene = rownames(deALL), rna.expr = rnaExpr,
metadata = metaMatrix.RNA)
09
功能富集分析
enrichOutput <- gdcEnrichAnalysis(gene = rownames(deALL), simplify = TRUE)
data(enrichOutput)
gdcEnrichPlot(enrichOutput, type = 'bar', category = 'GO', num.terms = 10)
gdcEnrichPlot(enrichOutput, type='bubble', category='GO', num.terms = 10)
文章参考：
http://bioconductor.org/packages/devel/bioc/vignettes/GDCRNATools/inst/doc/GDCRNATools.htm
小编总结
GDCRNATools作为集成了对GDC数据下载、处理和可视化的R包，可以说是功能十分全面了，可能是目前使用最方便的TCGA工具了，希望小编的介绍会对大家有所帮助，让小伙伴们多多了解GDCRNATools，为大家的科研工作助力！
END
腾讯云开发者
