--- title: 'Example data for TCGA Workflow: Analyze cancer genomics and epigenomics data using Bioconductor packages' author: Tiago C. Silva, Antonio Colaprico, Catharina Olsen, Fulvio D’Angelo, Gianluca Bontempi Michele Ceccarelli , and Houtan Noushmehr date: "`r Sys.Date()`" output: BiocStyle::html_document vignette: > %\VignetteIndexEntry{'Example data for TCGA Workflow: Analyze cancer genomics and epigenomics data using Bioconductor packages'} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- # Data Introduction This package provides a dataset for those wishing to try out the [TCGA Workflow: Analyze cancer genomics and epigenomics data using Bioconductor packages](https://f1000research.com/articles/5-1542/v2) [@10.12688/f1000research.8923.2]. The data in this package are a subset of the TCGA data for LGG (Lower grade glioma) and GBM (Glioblastoma multiforme) samples. # Loading the data ```{r, fig.show='hold'} library(TCGAWorkflowData) data("elmerExample") data("GBMnocnvhg19") data("GBMIllumina_HiSeq") data("LGGIllumina_HiSeq") data("mafMutect2LGGGBM") data("markersMatrix") data("histoneMarks") data("biogrid") data("genes_GR") ``` # Data creation The following commands were used to create the data included with this package. ## Genes information Download gene information for hg19 using TCGAbiolinks, which uses biomart parckage. ```{r , eval = FALSE, message=FALSE,warning=FALSE, include=TRUE} library(GenomicRanges) library(TCGAbiolinks) ############################## ## Recurrent CNV annotation ## ############################## # Get gene information from GENCODE using biomart genes <- TCGAbiolinks:::get.GRCh.bioMart(genome = "hg19") genes <- genes[genes$external_gene_id != "" & genes$chromosome_name %in% c(1:22,"X","Y"),] genes[genes$chromosome_name == "X", "chromosome_name"] <- 23 genes[genes$chromosome_name == "Y", "chromosome_name"] <- 24 genes$chromosome_name <- sapply(genes$chromosome_name,as.integer) genes <- genes[order(genes$start_position),] genes <- genes[order(genes$chromosome_name),] genes <- genes[,c("external_gene_id", "chromosome_name", "start_position","end_position")] colnames(genes) <- c("GeneSymbol","Chr","Start","End") genes_GR <- makeGRangesFromDataFrame(genes,keep.extra.columns = TRUE) save(genes_GR,genes,file = "genes_GR.rda", compress = "xz") ``` ## GISTIC results Download and save a subset of GBM GISTIC results from GDAC firehose. ```{R, eval=FALSE, message=FALSE, warning=FALSE, include=TRUE} library(RTCGAToolbox) # Download GISTIC results lastAnalyseDate <- getFirehoseAnalyzeDates(1) gistic <- getFirehoseData("GBM",gistic2_Date = lastAnalyseDate, GISTIC = TRUE) # get GISTIC results gistic.allbygene <- getData(gistic,type = "GISTIC", platform = "AllByGene") gistic.thresholedbygene <- getData(gistic,type = "GISTIC", platform = "ThresholdedByGene") save(gistic.allbygene,gistic.thresholedbygene,file = "GBMGistic.rda", compress = "xz") ``` ## Copy number variations (CNVs) The following code will download segmented CNV from SNP array (Affymetrix Genome-Wide Human SNP Array 6.0) for 20 Glioblastoma multiforme (GBM) samples. ```{r , eval=FALSE, include=TRUE, results='asis'} library(TCGAbiolinks) query.gbm.nocnv <- GDCquery(project = "TCGA-GBM", data.category = "Copy number variation", legacy = TRUE, file.type = "nocnv_hg19.seg", sample.type = c("Primary solid Tumor")) # to reduce time we will select only 20 samples query.gbm.nocnv$results[[1]] <- query.gbm.nocnv$results[[1]][1:20,] GDCdownload(query.gbm.nocnv, chunks.per.download = 100) gbm.nocnv <- GDCprepare(query.gbm.nocnv, save = TRUE, save.filename = "GBMnocnvhg19.rda") ``` ## Gene expression data The following code will download 20 LGG (Lower grade glioma) and 20 GBM (Glioblastoma multiforme) samples that have gene expression data. The Gene expression data is the raw expression signal for expression of a gene. ```{r , eval=FALSE, include=TRUE, results='asis'} query <- GDCquery(project = "TCGA-GBM", data.category = "Gene expression", data.type = "Gene expression quantification", platform = "Illumina HiSeq", file.type = "results", sample.type = c("Primary solid Tumor"), legacy = TRUE) # We will use only 20 samples to make the example faster query$results[[1]] <- query$results[[1]][1:20,] GDCdownload(query) gbm.exp <- GDCprepare(query, save = TRUE, summarizedExperiment = TRUE, save.filename = "GBMIllumina_HiSeq.rda") query <- GDCquery(project = "TCGA-LGG", data.category = "Gene expression", data.type = "Gene expression quantification", platform = "Illumina HiSeq", file.type = "results", sample.type = c("Primary solid Tumor"), legacy = TRUE) # We will use only 20 samples to make the example faster query$results[[1]] <- query$results[[1]][1:20,] GDCdownload(query) lgg.exp <- GDCprepare(query, save = TRUE, summarizedExperiment = TRUE, save.filename = "LGGIllumina_HiSeq.rda") ``` ## DNA methylation and Gene expression data The following code will select 10 LGG (Lower grade glioma) and 10 GBM (Glioblastoma multiforme) samples that have both DNA methylation and gene expression data. This objects will be then prepared to the format accept by the Biocondcutor package `r Biocpkg("ELMER")`([link])(http://bioconductor.org/packages/release/bioc/html/ELMER.html). The DNA methylation will have only probes in chromossome 9 in order to make the example of the workflow faster. For a real analysis, all chromossomes should be used. The Gene expression data is the normalized results for expression of a gene. ```{r , eval=FALSE, include=TRUE, results='asis'} #----------- 8.3 Identification of Regulatory Enhancers ------- library(TCGAbiolinks) # Samples: primary solid tumor w/ DNA methylation and gene expression matched_met_exp <- function(project, n = NULL){ # get primary solid tumor samples: DNA methylation message("Download DNA methylation information") met450k <- GDCquery(project = project, data.category = "DNA methylation", platform = "Illumina Human Methylation 450", legacy = TRUE, sample.type = c("Primary solid Tumor")) met450k.tp <- met450k$results[[1]]$cases # get primary solid tumor samples: RNAseq message("Download gene expression information") exp <- GDCquery(project = project, data.category = "Gene expression", data.type = "Gene expression quantification", platform = "Illumina HiSeq", file.type = "normalized_results", sample.type = c("Primary solid Tumor"), legacy = TRUE) exp.tp <- exp$results[[1]]$cases # Get patients with samples in both platforms patients <- unique(substr(exp.tp,1,15)[substr(exp.tp,1,12) %in% substr(met450k.tp,1,12)] ) if(!is.null(n)) patients <- patients[1:n] # get only n samples return(patients) } lgg.samples <- matched_met_exp("TCGA-LGG", n = 10) gbm.samples <- matched_met_exp("TCGA-GBM", n = 10) samples <- c(lgg.samples,gbm.samples) #----------------------------------- # 1 - Methylation # ---------------------------------- query.met <- GDCquery(project = c("TCGA-LGG","TCGA-GBM"), data.category = "DNA methylation", platform = "Illumina Human Methylation 450", legacy = TRUE, barcode = samples) GDCdownload(query.met) met <- GDCprepare(query.met, save = FALSE) met <- subset(met,subset = as.character(GenomicRanges::seqnames(met)) %in% c("chr9")) #----------------------------------- # 2 - Expression # ---------------------------------- query.exp <- GDCquery(project = c("TCGA-LGG","TCGA-GBM"), data.category = "Gene expression", data.type = "Gene expression quantification", platform = "Illumina HiSeq", file.type = "normalized_results", legacy = TRUE, barcode = samples) GDCdownload(query.exp) exp <- GDCprepare(query.exp, save = FALSE) save(exp, met, gbm.samples, lgg.samples, file = "elmerExample.rda", compress = "xz") ``` ## Mutation data The following code will download Mutation annotation files (aligned against the genoem of reference hg38) for LGG and GBM samples and merge them into one single single file. The GDC Somatic Mutation Calling Workflow used is the mutect2. For more information please check [GDC](https://docs.gdc.cancer.gov/Data/Bioinformatics_Pipelines/DNA_Seq_Variant_Calling_Pipeline/). ```{r , eval=FALSE, include=TRUE, results='asis'} library(TCGAbiolinks) LGGmut <- GDCquery_Maf(tumor = "LGG", pipelines = "mutect2") GBMmut <- GDCquery_Maf(tumor = "GBM", pipelines = "mutect2") mut <- plyr::rbind.fill(LGGmut, GBMmut) save(mut, LGGmut, GBMmut,file = "mafMutect2LGGGBM.rda") ``` ## Probes meta file from broadinstitute website for Copy Number Variation Analysis (CNV) analysis ```{r , eval=FALSE, include=TRUE, results='asis'} gdac.root <- "ftp://ftp.broadinstitute.org/pub/GISTIC2.0/hg19_support/" file <- paste0(gdac.root, "genome.info.6.0_hg19.na31_minus_frequent_nan_probes_sorted_2.1.txt") # Retrieve probes meta file from broadinstitute website if(!file.exists(basename(file))) downloader::download(file, basename(file)) # -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=--=-==--=--==---=-=-=-=-=-=-=-=-=-=-=-=-=-=-=--=-==--=--==-- # For hg38 analysis please take a look on: # https://gdc.cancer.gov/about-data/data-harmonization-and-generation/gdc-reference-files # File: SNP6 GRCh38 Liftover Probeset File for Copy Number Variation Analysis # -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=--=-==--=--==---=-=-=-=-=-=-=-=-=-=-=-=-=-=-=--=-==--=--==-- markersMatrix <- readr::read_tsv(basename(file), col_names = FALSE, col_types = "ccn", progress = FALSE) save(markersMatrix, file = "markersMatrix.rda", compress = "xz") ``` ## Biogrid data Download biogrid information. ```{r , eval=FALSE, include=TRUE, results='asis'} ### read biogrid info ### Check last version in https://thebiogrid.org/download.php file <- "http://thebiogrid.org/downloads/archives/Release%20Archive/BIOGRID-3.4.146/BIOGRID-ALL-3.4.146.tab2.zip" if(!file.exists(gsub("zip","txt",basename(file)))){ downloader::download(file,basename(file)) unzip(basename(file),junkpaths =TRUE) } tmp.biogrid <- read.csv(gsub("zip","txt",basename(file)), header=TRUE, sep="\t", stringsAsFactors=FALSE) save(tmp.biogrid, file = "biogrid.rda", compress = "xz") ``` ## GISTIC2.0 auxiliary data Download hg19_support information. ```{r , eval=FALSE, include=TRUE, results='asis'} file <- "ftp://ftp.broadinstitute.org/pub/GISTIC2.0/hg19_support/CNV.hg19.bypos.111213.txt" if(!file.exists(basename(file))) downloader::download(file, basename(file)) commonCNV <- readr::read_tsv(basename(file), progress = FALSE) commonCNV <- as.data.frame(commonCNV) save(commonCNV,file = "CNV.hg19.bypos.111213.rda",compress = "xz") ``` ## Histone marks The code below was used to download histone marks specific for brain tissue using the AnnotationHub package that can access the Roadmap database. ```{r results='hide', eval=FALSE, echo=FALSE, message=FALSE,warning=FALSE} library(ChIPseeker) library(AnnotationHub) library(pbapply) library(ggplot2) #------------------ Working with ChipSeq data --------------- # Step 1: download histone marks for a brain and non-brain samples. #------------------------------------------------------------ # loading annotation hub database ah = AnnotationHub() # Searching for brain consolidated epigenomes in the roadmap database bpChipEpi_brain <- query(ah , c("EpigenomeRoadMap", "narrowPeak", "chip", "consolidated","brain","E068")) # Get chip-seq data histone.marks <- pblapply(names(bpChipEpi_brain), function(x) {ah[[x]]}) names(histone.marks) <- names(bpChipEpi_brain) save(histone.marks, file = "histoneMarks.rda", compress = "xz") ``` # Session info ```{r sessionInfo, results='asis', echo=FALSE} pander::pander(sessionInfo(), compact = FALSE) ```