一. 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./Encyclopedia/pages/TCGA_Barcode/ 接触和分析过TCGA数据的朋友肯定会经常处理TCGA barcode的前15位(有时12位),实际从上图可以看出TCGA的barcode设计总共有28位之多。每一个短横杠衔接的都是含不同意义的序列,如下图: 另外,将barcode的组成从层次结构(树)来看,是这样的: 可参考 https://docs.gdc./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./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 Quantification" #下载miRNA数据 data.type = "miRNA Quantification" #下载Copy Number Variation数据 data.type = "Copy Number Segment" (4)workflow.type 该数据类型有很多种,根据data.type的不同而不同,不同的数据类型,有其对应的参数可供选择。比如Gene 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 array, WXS,CGH array, VALIDATION, Gene array,WGS, MSI-Mono-Dinucleotide Assay, miRNA 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 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 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") 结果图 |
|