R中用GDCRNATools包下载TCGA数据
2020-03-09 本文已影响0人
养猪场小老板
用GDCRNATools下载TCGA数据，以TCGA-STAD为例下载RNAseq
1）数据下载，gdcRNADownload（）函数
###########用GDCRNATools下载TCGA数据###########
#设置镜像
rm(list=ls())
options("repos"="https://mirrors.ustc.edu.cn/CRAN/")
#意思是，如果没有安装BiocManager包，那就安装BiocManager
if(!require("BiocManager")) install.packages("BiocManager",update = F,ask = F)
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
#安装GDCRNATools
BiocManager::install("GDCRNATools")
#加载GDCRNATools
library(GDCRNATools)
#下载方式一：指定gdcRNADownload（）函数中的参数data.type和project.id
#gdcRNADownload（）函数将下载 HTSeq Counts data只需通过指定参数data.type='RNAseq'和project.id。
#gdcRNADownload（）函数将下载BCGSC miRNA Profiling data 只需通过指定参数 data.type='miRNAs'和project.id。
#下载TCGA-STAD数据#
gdcRNADownload(project.id = 'TCGA-STAD', #参数项目id
data.type = 'RNAseq', #也可改为miRNAs
write.manifest = FALSE,
method = 'gdc-client', #先下载gdc-client，然后再下载病人的数据
directory = 'GDCRNATOOLS/TCGA-STAD/RNAseq')#保存路径
# 下载临床数据#
#在gdcClinicalDownload（）函数中指定project.id下载.xml格式的临床信息
gdcClinicalDownload(project.id = 'TCGA-STAD',
write.manifest = FALSE,
directory = 'GDCRNATOOLS/TCGA-STAD/Clinical')
#################################################################################
#################################################################################
#下载方式二：指定gdcRNADownload（）函数中的参数manifest
####### 下载 RNAseq 数据 #######
gdcRNADownload(manifest = 'TCGA-STAD/TCGA-STAD.RNAseq.gdc_manifest.2019-11-23T14-40-52.txt',
directory = 'TCGA-STAD/RNAseq')
#######下载 miRNAs 数据 #######
gdcRNADownload(manifest = 'TCGA-STAD/TCGA-STAD.miRNAs.gdc_manifest.2019-11-22T15-36-57.txt',
directory = 'TCGA-STAD/miRNAs')
####### 下载临床数据 #######
gdcRNADownload(manifest = 'TCGA-STAD/TCGA-STAD.Clinical.gdc_manifest.2019-11-23T14-42-01.txt',
directory = 'TCGA-STAD/Clinical')
2） 数据处理
-
2.1）gdcParseMetadata（）函数
#第一种方式： 通过提供 metadata file解析```
####### RNAseq #######
metaMatrix.RNA <- gdcParseMetadata(metafile='TCGA-PRAD/TCGA-PRAD.RNAseq.metadata.2019-11-23T17-23-59.json')#该.json文件可在官网下载
####### miRNAs #######
metaMatrix.MIR <- gdcParseMetadata(metafile='TCGA-PRAD/TCGA-PRAD.miRNAs.metadata.2017-11-23T17-33-55.json')
#第二种方式：指定project.id 和 data.type
######## RNAseq #######
metaMatrix.RNA <- gdcParseMetadata(project.id = 'TCGA-PRAD',
data.type = 'RNAseq',
write.meta = FALSE)
####### miRNAs #######
metaMatrix.MIR <- gdcParseMetadata(project.id = 'TCGA-PRAD',
data.type = 'miRNAs',
write.meta = FALSE)
-
2.2） 过滤样本，两个函数gdcFilterDuplicate（过滤重复的样本）和gdcFilterSampleType（过滤非肿瘤非正常组织）
# a）过滤重复的样本gdcFilterDuplicate()
metaMatrix.RNA <- gdcFilterDuplicate(metaMatrix.RNA)# RNAseq
metaMatrix.MIR <- gdcFilterDuplicate(metaMatrix.MIR)# miRNAs
# b)过滤非肿瘤非正常组织gdcFilterSampleType()
metaMatrix.RNA <- gdcFilterSampleType(metaMatrix.RNA)#RNA
metaMatrix.MIR <- gdcFilterSampleType(metaMatrix.MIR)#miRNAs
-
2.3） 融合数据，两个函数gdcRNAMerge(原始数据)和gdcClinicalMerge(临床数据)
##gdcRNAMerge（）将RNAseq的原始计数数据合并到单个表达矩阵中，其中行为Ensembl id，列为sample。
##gdcClinicalMerge（）将临床数据合并到数据框中，其中行是sample ID，列是临床特征。 如果key.info = TRUE，则仅报告那些最常用的临床特征，否则，将报告所有临床信息。
#### 合并 RNAseq/miRNAs数据
rnaMatrix <- gdcRNAMerge(metadata = metaMatrix.RNA,
path = 'TCGA-PRAD/RNAseq/',
data.type = 'RNAseq')#RNAseq
#or
mirMatrix <- gdcRNAMerge(metadata = metaMatrix.MIR,
path = 'TCGA-PRAD/miRNAs/',
data.type = 'miRNAs')#miRNAs
####合并临床数据
clinicalDa <- gdcClinicalMerge(path = 'TCGA-PRAD/Clinical/', key.info = TRUE)
-
2.4) TMM标准化和Voom转换，gdcVoomNormalization（）函数
#低表达基因（超过一半的样本中logcpm <1）默认情况下会被滤除。 通过在gdcVoomNormalization（）中设置filter = TRUE可以保留所有基因。
####### RNAseq #######
rnaExpr <- gdcVoomNormalization(counts = rnaMatrix, filter = FALSE)
####### miRNAs #######
mirExpr <- gdcVoomNormalization(counts = mirMatrix, filter = FALSE)
3) 差异基因表达分析
-
3.1)差异分析，函数gdcDEAnalysis( )
#gdcDEAnalysis()参数method中有三种 limma，edgeR和DESeq2
#识别用户定义的任何两组之间的差异表达基因（DEG）或miRNA
DEGAll <- gdcDEAnalysis(counts = rnaMatrix,
group = metaMatrix.RNA$sample_type,
comparison = 'PrimaryTumor-SolidTissueNormal',
method = 'limma')
-
3.2)报告DE基因/ miRNAs，gdcDEReport（）函数
#通过在gdcDEReport（）中设置geneType参数，可以分别报告All DEGs, DE long non-coding genes, DE protein coding genes 和 DE miRNAs。报告中的Gene symbols 和 biotypes是Ensembl 90注释的。
### All DEGs
deALL <- gdcDEReport(deg = DEGAll, gene.type = 'all')
#### DE long-noncoding
deLNC <- gdcDEReport(deg = DEGAll, gene.type = 'long_non_coding')
#### DE protein coding genes
dePC <- gdcDEReport(deg = DEGAll, gene.type = 'protein_coding')
-
3.3） DEG 可视化
#火山图和条形图可分别通过gdcVolcanoPlot（）和gdcBarPlot（）函数可视化DE分析结果。
#DEG表达式矩阵上的聚类分析结果可以通过gdcHeatmap（）函数分析和绘制。
# a)火山图 Volcano plot
gdcVolcanoPlot(DEGAll)
# b)条形图 Barplot
gdcBarPlot(deg = deALL, angle = 45, data.type = 'RNAseq')
# c) 热图 Heatmap
degName = rownames(deALL)
gdcHeatmap(deg.id = degName, metadata = metaMatrix.RNA, rna.expr = rnaExpr)
