TCGAbiolinks包下载TCGA数据进行表达差异分析-乳腺癌案例

DoubleHelix

发布于 2019-08-07 10:33:51

15.9K1

发布于 2019-08-07 10:33:51

一. TCGAbiolinks介绍与安装

TCGAbiolinks -一个用于TCGA数据综合分析的R/BioConductor软件包，能够通过GDC Application Programming Interface (API)访问 National Cancer Institute (NCI) Genomic Data Commons (GDC) ，来搜索、下载和准备相关数据，以便在R中进行分析。

1.TCGAbiolinks包的安装

devtools::install_github(repo = "BioinformaticsFMRP/TCGAbiolinks")

也可以通过下面代码安装

# 当前R的版本是"3.5"，对应的TCGAbiolinks版本是"3.7" or "3.8"
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("TCGAbiolinks", version = "3.8")

安装成功后就是加载包了，如果加载还需要某些包，就先安装相应的包：

library(TCGAbiolinks)
library(dplyr)
library(DT)

2.参考文献

文章中使用了该包，请引用：

Colaprico A, Silva TC, Olsen C, Garofano L, Cava C, Garolini D, Sabedot T, Malta TM, Pagnotta SM, Castiglioni I, Ceccarelli M, Bontempi G and Noushmehr H. "TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data." Nucleic acids research (2015): gkv1507.

此外，如果您使用ELMER分析，请引用：

Yao, L., Shen, H., Laird, P. W., Farnham, P. J., & Berman, B. P. "Inferring regulatory element landscapes and transcription factor networks from cancer methylomes." Genome Biol 16 (2015): 105.

Yao, Lijing, Benjamin P. Berman, and Peggy J. Farnham. "Demystifying the secret mission of enhancers: linking distal regulatory elements to target genes." Critical reviews in biochemistry and molecular biology 50.6 (2015): 550-573.

3.TCGA条码（barcode）

TCGA条码由一组标识符组成。每个都专门标识一个TCGA数据元素。有关元数据标识符如何组成条形码的说明，请参见下图：

更多TCGA barcode信息：

https://docs.gdc.cancer.gov/Encyclopedia/pages/TCGA_Barcode/

接触和分析过TCGA数据的朋友肯定会经常处理TCGA barcode的前15位（有时12位），实际从上图可以看出TCGA的barcode设计总共有28位之多。每一个短横杠衔接的都是含不同意义的序列，如下图：

另外，将barcode的组成从层次结构（树）来看，是这样的：

可参考 https://docs.gdc.cancer.gov/Encyclopedia/pages/images/TCGA-TCGAbarcode-080518-1750-4378.pdf

二、数据下载准备

首先，默认读者已经掌握了TCGA数据库网页版的使用，如果对该数据库的使用还有问题，请参考文章：TCGA数据库使用教程。

然后，我们解释不同的下载方法和SummarizedExperiment对象，它是TCGAbiolinks中使用的默认数据结构。

1.TCGAbiolinks可以使用两种方法下载GDC数据：

client：此方法创建MANIFEST文件并使用 GDC Data Transfer Tool下载数据。此方法更可靠，但与api方法相比可能更慢。

api：此方法使用 GDC Application Programming Interface (API)下载数据。这将创建一个MANIFEST文件，并且下载的数据将是一个格式为tar.gz的压缩文件。如果文件的大小和数量太大，这个tar.gz文件会太大导致下载失败的可能性提高。为了解决这个问题，我们将使用files.per.chunk功能将文件拆分成多个小文件，例如，如果chunks.per.download等于10，我们将每个tar.gz分为10个小文件下载。

2.SummarizedExperiment对象

使用 SummarizedExperiment package，我们可以从SummarizedExperiment对象中提取三个主要的数据矩阵

colData(data)：获得样本信息的矩阵，包括了从对应TCGA论文中获得的临床数据以及肿瘤亚型信息

assay(data)：获得Assay信息的矩阵，就是每一个样本中基因的表达量

rowRanges(data)：获得特征（一般是指基因）信息的矩阵，包括特征的元数据，例如基因所在基因组范围

3.Summarized Experiment：注释信息

使用GDCprepare函数时，会调用一个参数SummarizedExperiment，该参数决定了输出类型为Summarized Experiment（默认选项）或数据框。为了创建一个Summarized Experiment对象，我们需要使用最新的基因组注释文件进行数据注释。比如：1）对于legacy数据（与hg19对齐的数据），TCGAbiolinks正在使用GRCh37.p13进行注释；2）对于harmonized数据（与hg38对齐的数据），TCGAbiolinks正在使用GRCh38.p7 (May 2017)进行注释。不幸的是，在GRCh38.p7 这样的注释文件更新后，比如一些基因缩写名称的改变/删除、更改基因坐标等。这可能会导致一些TCGA数据的丢失。例如，如果基因被删除，我们就不能再映射它了，那么在SummarizedExperiment中这些信息会丢失。如果设置SummarizedExperiment为FALSE，您将获得未修改的数据，并需要您自己去注释。此外，DNA甲基化数据并没有更新。但是可以在这里找到最新的元数据： http://zwdzwd.github.io/InfiniumAnnotation

4.GDCquery()参数解析

关于TCGAbiolinks包中的函数很多，这里重点介绍GDCquery()这个函数，他是下载数据中很重要的一个函数。其他函数可查看帮助文档，不知道有什么函数可参考文章：TCGAbiolinks包帮助文档（英文版）

（1）project

可以使用TCGAbiolinks:::getGDCprojects()$project_id)得到各个癌种的项目id，总共有45个ID值。

（2）data.category

data.category总共有8种

可以使用TCGAbiolinks:::getProjectSummary(project)查看project中有哪些数据类型，如查询"TCGA-LIHC"，有7种数据类型，case_count为病人数，file_count为对应的文件数。小编要下载表达谱，所以设置data.category="Transcriptome Profiling"

（3）data.type

筛选要下载的文件的数据类型。没有命令可以查看data.type里都有哪些数据类型，但是根据官网连接，查看data.type有12种，但设置参数的时候不代表所有的project和data.category都对应12种。先在官网查看后再设置。

#下载rna-seq的counts数据
data.type = "Gene Expression Quantification"
#下载miRNA数据
data.type = "miRNA Expression Quantification"
#下载Copy Number Variation数据
data.type = "Copy Number Segment"

（4）workflow.type

该数据类型有很多种，根据data.type的不同而不同，不同的数据类型，有其对应的参数可供选择。比如Gene Expression Quantification数据类型下workflow.type 有4种类型分别为：

HTSeq - FPKM-UQ：FPKM上四分位数标准化值

HTSeq - FPKM：FPKM值/表达量值

HTSeq - Counts：原始count数

STAR - Counts

具体可在GDC官网查看

（5）legacy

这个参数主要是因为TCGA数据有两个入口可以下载，GDC Legacy Archive 和 GDC Data Portal，区别主要是注释参考基因组版本不同分别是：GDC Legacy Archive（hg19和GDC Data Portal（hg38）。参数默认为FALSE，下载GDC Data Portal（hg38）。这里建议是，下载转录组层面的数据使用hg38，下载DNA层面的数据使用hg19，因为比如做SNP分析的时候很多数据库没有hg38版本的数据，都是hg19的。

（6）access

数据开放和不开放，有两个参数:controlled, open。

（7）platform

这里涉及到的平台种类非常多，可以具体去官网看每种数据都有什么平台的可以下载。这个参数可以省略不设置。

（8）file.type

主要是在GDC Legacy Archive下载数据的时候使用，可以参考官网说明。在GDC Data Portal下载数据，该参数省略不设置。

（9）barcode

A list of barcodes to filter the files to download。可以根据这个参数设置只下载某个样本等。如：

barcode = c("TCGA-14-0736-02A-01R-2005-01", "TCGA-06-0211-02A-02R-2005-01"）

（10）experimental.strategy

两个下载入口参数选择

GDC Data Portal：WXS, RNA-Seq, miRNA-Seq, Genotyping Array.

Legacy: WXS, RNA-Seq, miRNA-Seq, Genotyping Array, DNA-Seq, Methylation array, Protein expression array, WXS,CGH array, VALIDATION, Gene expression array,WGS, MSI-Mono-Dinucleotide Assay, miRNA expression array, Mixed strategies, AMPLICON, Exon array, Total RNA-Seq, Capillary sequencing, Bisulfite-Seq

（11）sample.type

A sample type to filter the files to download，可以对样本类型进行过滤下载。这里我要下载所有样本类型数据，不设置。部分值选择如下（全部可以查看官网）：如sample.type = "Recurrent Solid Tumor"

四. 数据下载实例

基因表达数据的下载

我们以乳腺癌(BRCA) 数据集的下载和分析作为案例进行讲解。

（1）加载包

# 可先设置好工作目录，
# 第一次学习最后是一个空文件夹，没运行一行代码，查看产生的数据和文件。
# 加载响应的包，默认已经安装好TCGAbiolinks包
library(TCGAbiolinks)
library(plyr)
library(limma)
library(biomaRt)
library(SummarizedExperiment)

（2）数据下载

query <- GDCquery(project = "TCGA-BRCA",
                  data.category = "Transcriptome Profiling",
                  data.type = "Gene Expression Quantification", 
                  workflow.type = "HTSeq - Counts")
# 从query中获取结果表，它可以选择带有cols参数的列，并使用rows参数返回若干行。
# 1222个barcode
samplesDown <- getResults(query,cols=c("cases"))

# samplesDown中筛选出TP（肿瘤）样本的barcode
#1102个barcode
dataSmTP <- TCGAquery_SampleTypes(barcode = samplesDown,
                                  typesample = "TP")
# samplesDown中筛选出NT(正常组织)样本的barcode
#113个barcode
dataSmNT <- TCGAquery_SampleTypes(barcode = samplesDown,
                                  typesample = "NT")

###选择100个正常组织和100个肿瘤组织样本作为洋酒对象
dataSmTP_short <- dataSmTP[1:100]
dataSmNT_short <- dataSmNT[1:100]

queryDown <- GDCquery(project = "TCGA-BRCA", 
                      data.category = "Transcriptome Profiling",
                      data.type = "Gene Expression Quantification", 
                      workflow.type = "HTSeq - Counts", 
                      barcode = c(dataSmTP_short, dataSmNT_short))
# 下载数据，默认存放位置为当前工作目录下的GDCdata文件夹中。
GDCdownload(query = queryDown)

（3）数据处理

#读取下载的数据并将其准备到R对象中，在工作目录生成brca_case1.rda文件，
#同时还产生Human_genes__GRCh38_p12_.rda文件
dataPrep1 <- GDCprepare(query = queryDown, save = TRUE, save.filename = "brca_case1.rda")

# 去除dataPrep1中的异常值，dataPrep数据中含有肿瘤组织和正常组织的数据
dataPrep <- TCGAanalyze_Preprocessing(object = dataPrep1, 
                                      cor.cut = 0.6,
                                      datatype = "HTSeq - Counts")
write.csv(dataPrep,file = "dataPrep.csv",quote = FALSE)
###############Tumor purity filtering###########

###vector containing all TCGA barcodes that hhave 60% tumor purity or more
# TCGAtumor_purity使用来自5种方法的5个估计值作为阈值对TCGA样本进行过滤
# 这5个值是estimate, absolute, lump, ihc, cpe，这里设置cpe=0.6
purityDATA <- TCGAtumor_purity(colnames(dataPrep1), 0, 0, 0, 0, 0.6)
# filtered 为被过滤的数据， pure_barcodes是我们要的数据
Purity.BRCA<-purityDATA$pure_barcodes

################DEA with Molecular subtypes############

####Molecular subtypes:

####Filtering data so all samples have a pam50 subtype for BRCA:

##diff contains TCGA samples that have an available molecular subtype
###Also Applying Tumor purity filtering by making a set difference
#documentation for TCGA_MolecularSubtype is available
# 将TCGAtumor_purity方法纯化后的数据与肿瘤组织中的数据(已去除异常值)的分子亚型数据进行setdiff运算
# 也就是将在Purity.BRCA肿瘤组织中都具有分子亚型的barcode返回，总共83个。
diff<-setdiff(Purity.BRCA,TCGA_MolecularSubtype(colnames(dataPrep[,dataSmTP_short]))$filtered)
write.csv(diff,file = "diff.csv",quote = FALSE)

# 获取所有diff中的数据，83个barcode的表达数据
dataFilt.brca.cancer<-dataPrep[,diff]
write.csv(dataFilt.brca.cancer,file = "dataFilt.brca.cancer.csv",quote = FALSE)

#从dataPrep中把正常的barcode的表达数据取出来，100个barcode表达数据。
dataFilt.brca.normal<-dataPrep[,dataSmNT_short]

#将正常组织和肿瘤组织的数据合并，183个barcode表达数据。
dataFilt.brca<-cbind(dataFilt.brca.normal, dataFilt.brca.cancer)
write.csv(dataFilt.brca,file = "dataFilt.brca.csv",quote = FALSE)

mol_subtypes<-c(rep("Normal", 100), TCGA_MolecularSubtype(colnames(dataFilt.brca.cancer))$subtypes$subtype)
write.csv(mol_subtypes,file = "mol_subtypes.csv",quote = FALSE)
mol_subtypes<-make.names(mol_subtypes)
#dataFilt.brca数据中的ensembl基因序号更换为基因名称
rownames(dataFilt.brca)<-rowData(dataPrep1)$external_gene_name

#将数据进行标准化
dataNorm.brca <- TCGAanalyze_Normalization(tabDF = dataFilt.brca,
                                           geneInfo = geneInfo,
                                           method = "gcContent")

# 将标准化后的数据再过滤，得到最终的数据
dataFilt.brca.final <- TCGAanalyze_Filtering(tabDF = dataNorm.brca,
                                             method = "quantile", 
                                             qnt.cut =  0.25)
write.csv(dataFilt.brca.final,file = "dataFilt.brca.final.csv",quote = FALSE)

（4）表达差异分析

DEG.brca.edgeR <- TCGAanalyze_DEA(MAT=dataFilt.brca.final,
                            pipeline="edgeR",
                            batch.factors = c("TSS"),
                            Cond1type = "Normal",
                            Cond2type = "Tumor",
                            fdr.cut = 0.01 ,
                            logFC.cut = 1,
                            voom = FALSE,
                            method = "glmLRT", ClinicalDF = data.frame(),
                            Condtypes = mol_subtypes,
                            contrast.formula = "Mycontrast = (BRCA.LumA+BRCA.LumB)/2 -Normal")
#in this DEA we use Normal as a reference, thus genes with LogFC > 1 are up regulated in the subtypes
#with respect to the normal samples and viceversa for the LogFC < -1.

write.csv(DEG.brca.edgeR, "DEGsBRCA_edgeR_091018.csv", quote = FALSE)
#3404 genes identified

DEG.brca.edgeR.filt.TSS<-DEG.brca.edgeR$Mycontrast[which(abs(DEG.brca.edgeR$Mycontrast$logFC) >= 6), ]
#57 genes
TCGAVisualize_volcano(DEG.brca.edgeR$Mycontrast$logFC, DEG.brca.edgeR$Mycontrast$FDR,
                      filename = "LuminalABvsNormal_FC6.TSS.edgeR.pdf", xlab = "logFC",
                      names = rownames(DEG.brca.edgeR$Mycontrast), show.names = "highlighted",
                      x.cut = 1, y.cut = 0.01, 
                      highlight = rownames(DEG.brca.edgeR$Mycontrast)[which(abs(DEG.brca.edgeR$Mycontrast$logFC) >= 6)],
                      highlight.color = "orange")

结果图