专栏首页生物信息云TCGAbiolinks包下载TCGA数据进行表达差异分析-乳腺癌案例

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

一. 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")

结果图

本文分享自微信公众号 - MedBioInfoCloud(MedBioInfoCloud),作者:DoubleHelix

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2019-05-28

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • R语言基础绘图教程——第2章:散点图

    plot函数中,x和y分别表示所绘图形的横坐标和纵坐标;函数中的...为附加的参数。

    DoubleHelix
  • 基因芯片数据挖掘分析表达差异基因

    基因芯片(genechip)(又称DNA芯片、生物芯片)的原型是80年代中期提出的。基因芯片的测序原理是杂交测序方法,即通过与一组已知序列的核酸探针杂交进行核酸...

    DoubleHelix
  • R语言数据分析与挖掘(第六章):主成分分析(1)——主成分分析概论

    在许多领域的研究与应用中,往往需要对反映事物的多个变量进行大量的观测,收集大量数据以便进行分析寻找规律。多变量大样本无疑会为研究和应用提供了丰富的信息,但也在一...

    DoubleHelix
  • AI一分钟 | 特斯拉再融46亿;腾讯AI Lab宣布开源多标签图像数据集

    近日,《证券日报》记者登录上海市工商行政管理局官网发现,特斯拉(上海)有限公司的注册资本已由 1 亿元增至 46.7 亿元,这意味着马斯克凭借特斯拉这匾金字招牌...

    AI科技大本营
  • 通过Pandas实现快速别致的数据分析

    在您选择和准备数据进行建模之前,您需要事先了解一些基础内容。

    Giacin
  • PHP 多任务协程处理

    上周 有幸和同事一起在 SilverStripe 分享最近的工作事宜。今天我计划分享 PHP 异步编程,不过由于上周我聊过 ReactPHP;我决定讨论一些不一...

    柳公子
  • 【压缩率3000%】上交大ICCV:精度保证下的新型深度网络压缩框架

    【新智元导读】上海交通大学人工智能实验室的研究人员提出了一种新的方法,能够在保证网络模型精度的前提下对深度网络进行压缩。相关论文已被ICCV 2017接收,由上...

    新智元
  • 平均年薪30万!人才缺口500万!人工智能工程师为什么这么火?

    在一线城市,年薪10万仅能饱腹,就算熬夜加班苦干10年达到100万,相信也所剩无几。

    Leetcode名企之路
  • 中止请求和超时 跨域的HTTP请求 认证方式 JSONP

    作为同源策略的一部分,XMLHttpRequest对象可以发起HTTP请求,由于同源的影响,导致必须是同源的,

    mySoul
  • KnockoutJS语法

      假设我们的页面输入区域有一个div用来展示一件物品的名字,同时有一个输入框用来编辑这件物品的名字

    javascript.shop

扫码关注云+社区

领取腾讯云代金券