TCGA Workflow: Analyze cancer genomics and epigenomics data using Bioconductor packages

Biotechnological advances in sequencing have led to an explosion of publicly available data via large international consortia such as The Cancer Genome Atlas (TCGA), The Encyclopedia of DNA Elements (ENCODE), and The NIH Roadmap Epigenomics Mapping Consortium (Roadmap). These projects have provided unprecedented opportunities to interrogate the epigenome of cultured cancer cell lines as well as normal and tumor tissues with high genomic resolution. The Bioconductor project offers more than 1,000 open-source software and statistical packages to analyze high-throughput genomic data. However, most packages are designed for specific data types (e.g. expression, epigenetics, genomics) and there is no one comprehensive tool that provides a complete integrative analysis of the resources and data provided by all three public projects. A need to create an integration of these different analyses was recently proposed. In this workflow, we provide a series of biologically focused integrative analyses of different molecular data. We describe how to download, process and prepare TCGA data and by harnessing several key Bioconductor packages, we describe how to extract biologically meaningful genomic and epigenomic data. Using Roadmap and ENCODE data, we provide a work plan to identify biologically relevant functional epigenomic elements associated with cancer. To illustrate our workflow, we analyzed two types of brain tumors: low-grade glioma (LGG) versus high-grade glioma (glioblastoma multiform or GBM). This workflow introduces the following Bioconductor packages: AnnotationHub, ChIPSeeker, ComplexHeatmap, pathview, ELMER, GAIA, MINET, RTCGAToolbox, TCGAbiolinks.


Introduction
Cancer is a complex genetic disease spanning multiple molecular events such as point mutations, structural variations, translocations and activation of epigenetic and transcriptional signatures and networks. The effects of these events take place at different spatial and temporal scales with interlayer communications and feedback mechanisms creating a highly complex dynamic system. To gain insight into the biology of tumors most of the research in cancer genomics is aimed at the integration of the observations at multiple molecular scales and the analysis of their interplay. Even if many tumors share similar recurrent genomic events, their relationships with the observed phenotype are often not understood. For example, although we know that the majority of the most aggressive form of brain tumors such as glioma harbor the mutation of a single gene (IDH), the mechanistic explanation of the activation of its characteristic epigenetic and transcriptional signatures are still far to be well characterized. Moreover, network-based strategies have recently emerged as an effective framework for the discovery functional disease drivers that act as main regulators of cancer phenotypes.
Indeed, recent technological developments allowed the deposition of large amounts of genomic and epigenomic data, such as gene expression, DNA methylation, and genomic localization of transcription factors, into freely available public international consortia like The Cancer Genome Atlas (TCGA), The Encyclopedia of DNA Elements (ENCODE), and The NIH Roadmap Epigenomics Mapping Consortium (Roadmap) 1 . An overview of the three consortia is described below: • The Cancer Genome Atlas (TCGA): The TCGA consortium, which is a National Institute of Health (NIH) initiative, makes publicly available molecular and clinical information for more than 30 types of human cancers including exome (variant analysis), single nucleotide polymorphism (SNP), DNA methylation, transcriptome (mRNA), microRNA (miRNA) and proteome. Sample types available at TCGA are: primary solid tumors, recurrent solid tumors, blood derived normal and tumor, metastatic, and solid tissue normal 2 .

• The Encyclopedia of DNA Elements (ENCODE): Found in 2003 by the National Human Genome
Research Institute (NHGRI), the project aims to build a comprehensive list of functional elements that have an active role in the genome, including regulatory elements that govern gene expression. Biosamples includes immortalized cell lines, tissues, primary cells and stem cells 3 .
• The NIH Roadmap Epigenomics Mapping Consortium: This was launched with the goal of producing a public resource of human epigenomic data in order to analyze biology and disease-oriented research. Roadmap maps DNA methylation, histone modifications, chromatin accessibility, and small RNA transcripts in stem cells and primary ex vivo tissues 4,5 .
Briefly, these three consortia provide large scale epigenomic data onto a variety of microarrays and next-generation sequencing (NGS) platforms. Each consortium encompasses specific types of biological information on specific type of tissue or cell and when analyzed together, it provides an invaluable opportunity for research laboratories to better understand the developmental progression of normal cells to cancer state at the molecular level and importantly, correlate these phenotypes with tissue of origins.
Although there exists a wealth of possibilities 6 in accessing cancer associated data, Bioconductor represent the most comprehensive set of open source, updated and integrated professional tools for the statistical analysis of large scale genomic data. Thus, we propose our workflow within Bioconductor to describe how to download, process, analyze and integrate cancer data to understand specific cancer-related specific questions. However, there is no tool that solves

Amendments from Version 1
In this revised version of our workflow, we made the following major changes: -The introduction now includes the GDC NCI data portal (https://gdc.nci.nih.gov/). -The codes to acquire TCGA data were rewritten in order to use the new version of the TCGAbiolinks package. TCGAbiolinks was entirely redesigned to query, download and prepare data from the GDC NCI data portal (https://gdc.nci.nih.gov/) instead of the inactive DCC TCGA data portal (https://tcga-data.nci.nih.gov).
Minor changes includes improvements based on the referee's comments, such as: -The inclusion of a paragraph to introduce GISTIC data -The inclusion of a paragraph explaining the differences between the open (TCGA level 3 and 4 data) and controlled data (TCGA level 1 and 2 data) and pointing to sources that might help the user request access to controlled data. -Hyperlinks and references were corrected.
-Improvements in the text by removing a few redundancies.
-We used data aligned to reference genome hg19 in all steps of the workflow.

See referee reports
REVISED the issue of integration in a comprehensive sequence and mutation information, epigenomic state and gene expression within the context of gene regulatory networks to identify oncogenic drivers and characterize altered pathways during cancer progression. Therefore, our workflow presents several Bioconductor packages to work with genomic and epigenomics data.
TCGAquery_subtype section in TCGAbiolinks vignette), and clinical information. Listing 1 shows how to use these functions to download DNA methylation and gene expression data from the GDC legacy database and 2 shows how to download copy number variation from harmonized data portal. Other examples, that access the harmonized data can be found in the TCGAbiolinks vignette.

Listing 2. Downloading TCGA copy number variation data from GDC harmonized database with TCGAbiolinks
If a summarizedExperiment object was chosen, the data can be accessed with three different accessors: assay for the data information, rowRanges to gets the range of values in each row and colData to get the sample information (patient, batch, sample type, etc) 9,10 . An example is shown in Listing 3.

Listing 3. summarizedExperiment accessors
The clinical data can be obtained using TCGAbiolinks through two methods. The first one will download only the indexed GDC clinical data which includes diagnoses (vital status, days to death, age at diagnosis, days to last follow up, days to recurrence), treatments (days to treatment, treatment id, therapeutic agents, treatment intent type), demographic (gender, race, ethnicity) and exposures (cigarettes per day, weight, height, alcohol history) information. This indexed clinical data can be obtained using the function GDCquery_clinical which can be used as described in Listing 4. This function has two arguments project ("TCGA-GBM","TARGETAML", etc) and type ("Clinical" or "Biospecimen"). The second method will download the xml files with all clinical data for the patient and retrieve the desired information from it. This will give access to all clinical data available which includes patient (tumor tissue site, histological type, gender, vital status, days to birth, days to last follow up, etc), drug (days to drug therapy start, days to drug therapy end, therapy types, drug name), radiation (days to radiation therapy start, days to radiation therapy end, radiation type, radiation dosage), new tumor event (days to new tumor event after initial treatment, new neoplasm event type, additional pharmaceutical therapy), follow up (primary therapy outcome success, follow up treatment success, vital status, days to last follow up, date of form completion), stage event (pathologic stage, tnm categories), admin (batch number, project code, disease code, Biospecimen Core Resource).

Listing 4. Downloading clinical data with TCGAbiolinks
Mutation information is stored in two types of Mutation Annotation Format (MAF): Protected and Somatic (or Public) MAF files, which are derived from the GDC annotated VCF files. Annotated VCF files often have variants reported on multiple transcripts whereas the protected MAF (*protected.maf) only reports the most critically affected one and the Somatic MAFs (*somatic.maf) are further processed to remove low quality and potential germline variants. To download Somatic MAFs data using TCGAbiolinks, GDCquery_maf function is provided (see Listing 5).
1 gbm.subtypes <-TCGAquery_subtype(tumor = "gbm") 2 brca.subtypes <-TCGAquery_subtype(tumor = "brca") Listing 6. Accessing subtype information retrieved from TCGA papers Downloading data from Broad TCGA GDAC. The Bioconductor package RTCGAToolbox 25 provides access to Firehose Level 3 and 4 data through the function getFirehoseData. The following arguments allows users to select the version and tumor type of interest: • dataset -Tumor to download. A complete list of possibilities can be viewed with getFirehoseDatasets function. • runDate -Stddata run dates. Dates can be viewed with getFirehoseRunningDates function. • gistic2_Date -Analyze run dates. Dates can viewed with getFirehoseAnalyzeDates function.

Listing 7. Downloading TCGA data files with RTCGAtoolbox
Finally, using RTCGAToolbox the user can retrieve CNV level 4 data, including the amplified or deleted genes identified by GISTIC which rates each segment based on the frequency of occurrence combined with the amplitude of aberration, using a permutation test to assess the statistical significance. Among GISTIC results there are two tables that can be accessed by RTCGAToolbox: • A gene-level table of copy number values for all samples. The copy number values in the table are in units of (copy number -2), so that no amplification or deletion is 0, genes with amplifications have positive values, and genes with deletions are negative values. The data are converted from marker level to gene level using the extreme method: a gene is assigned the greatest amplification or the least deletion value among the markers it covers.
• A gene-level table of discrete amplification and deletion indicators for all samples. A table value of 0 means no amplification or deletion above the threshold (diploid normal copy). Amplifications are positive numbers: 1 means amplification above the amplification threshold (low-level gain, 1 extra copy); 2 means amplifications larger to the arm level amplifications observed for the sample (high-level amplification, 2 or more extra copies). Deletions are represented by negative table values: -1 represents deletion beyond the threshold (possibly a heterozygous deletion); -2 means deletions greater than the minimum arm-level deletion observed for the sample (possibly a homozygous deletion).
More details about the GISTIC algorithm and its use are described in 26-28. (see Listing 8).

Genomic analysis
Copy number variations (CNVs) have a critical role in cancer development and progression. A chromosomal segment can be deleted or amplified as a result of genomic rearrangements, such as deletions, duplications, insertions and translocations. CNVs are genomic regions greater than 1 kb with an alteration of copy number between two conditions (e.g., Tumor versus Normal).
TCGA collects copy number data and allows the CNV profiling of cancer. Tumor and paired-normal DNA samples were analyzed for CNV detection using microarray and sequencing-based technologies. Level 3 processed data are the aberrant regions along the genome resulting from CNV segmentation, and they are available for all copy number technologies.
In this section, we will show how to analyze CNV level 3 data from TCGA to identify recurrent alterations in cancer genome. We analyzed GBM and LGG segmented CNV from SNP array (Affymetrix Genome-Wide Human SNP Array 6.0).
Pre-Processing Data. The only CNV platform available for both LGG and GBM in TCGA is "Affymetrix Genome-Wide Human SNP Array 6.0". Using TCGAbiolinks, we queried for CNV SNP6 level 3 data for primary solid tumor samples in the legacy database. Data for selected samples were downloaded and prepared in two separate rse objects (RangedSummarizedExperiment).

Listing 10. Recurrent CNV identification in cancer with GAIA
Recurrent amplifications and deletions were identified for both LGG ( Figure 1a) and GBM (Figure 1b), and represented in chromosomal overview plots by a statistical score (-log 10 corrected p-value for amplifications and log 10 corrected p-value for deletions). Genomic regions identified as significantly altered in copy number (corrected p-value < 10 -4 ) were then annotated to report amplified and deleted genes potentially related with cancer.

Gene annotation of recurrent CNV.
The aberrant recurrent genomic regions in cancer, as identified by GAIA, have to be annotated to verify which genes are significantly amplified or deleted. Using biomaRt we retrieved the genomic ranges of all human genes and we compared them with significant aberrant regions to select full length genes. An example of the result is shown in Table 1.

Overview of genomic alterations by circos plot
Genomic alterations in cancer, including CNV and mutations, can be represented in an effective overview plot named circos. We used circlize CRAN package to represent significant CNV (resulting from GAIA analysis) and recurrent   LGG. Circos plot can illustrate molecular alterations genome-wide or only in one or more selected chromosomes. The Figure 3 shows the resulting circos plot for all chromosomes, while the Figure 4 shows the plot for only the chromosome 17.

Missense_Mutation Frame_Shift_Del Nonsense_Mutation Frame_Shift_Ins
Transcriptomic analysis Pre-Processing Data. The LGG and GBM data used for following transcriptomic analysis were downloaded using TCGAbiolinks. We downloaded only primary solid tumor (TP) samples, which resulted in 516 LGG samples and 156 GBM samples, then prepared it in two separate rse objects (RangedSummarizedExperiment) saving them as an R object with a file name including both the name of the cancer and the name of the platform used for gene expression data (see Listing 14).

Missense_Mutation Frame_Shift_Del Nonsense_Mutation Frame_Shift_Ins
To pre-process the data, first, we searched for possible outliers using the TCGAanalyze_Preprocessing function, which performs an Array Array Intensity correlation AAIC (lines 14-17 and 26-29 of Listing 15). The Array-array intensity correlation plot (AAIC) is a re-adaptation of the function correlationPlot from the R/Bioconductor affyQCReport package 31 that shows an heat map of the array-array Spearman/Pearson rank correlation coefficients. The arrays are ordered using the phenotypic data (if available) in order to place arrays with similar samples adjacent to each other. Self-self correlations are on the diagonal and by definition have a correlation coefficient of 1.0. Data from similar tissues or treatments will tend to have higher coefficients. This plot is useful for detecting outliers, failed hybridizations, or mistracked samples.
In our example, we defined a square symmetric matrix of Pearson correlation among all samples in each cancer type (LGG or GBM). This matrix found 0 samples with low correlation (cor.cut = 0.6) that can be identified as possible outliers.
Second, using the TCGAanalyze_Normalization function, which encompasses the functions of the EDASeq package, we normalized mRNA transcripts.
EA: enrichment analysis. In order to understand the underlying biological process of DEGs we performed an enrichment analysis using TCGAanalyze_EA_complete function (see Listing 16).

Listing 16. Enrichment analysis
TCGAanalyze_EAbarplot outputs a bar chart as shown in Figure 5 with the number of genes for the main categories of three ontologies (i.e., GO:biological process, GO:cellular component, and GO:molecular function).
The Figure 5 shows canonical pathways significantly over-represented (enriched) by the DEGs. The most statistically significant canonical pathways identified in the DEGs are ranked according to their p-value corrected FDR (-Log10) (colored bars) and the ratio of list genes found in each pathway over the total number of genes in that pathway (ratio, red line).

PEA: Pathways enrichment analysis.
To verify if the genes found have a specific role in a pathway, the Bioconductor package pathview 34 can be used. Listing 17 shows an example how to use it. It can receive, for example, a named vector of gene with the expression level, the pathway.id which can be found in KEGG database, the species ('hsa' for Homo sapiens) and the limits for the gene expression (see Figure 6).

Inference of gene regulatory networks.
Starting with the set of differentially expressed genes, we infer gene regulatory networks using the following state-of-the art inference algorithms: ARACNE 35 , CLR 36 , MRNET 37 and C3NET 38 . These methods are based on mutual inference and use different heuristics to infer the edges in the network. These methods have been made available via Bioconductor/CRAN packages (MINET 39 and c3net 38 , respectively).
Many gene regulatory interactions have been experimentally validated and published. These 'known' interactions can be accessed using different tools and databases such as BioGrid 40 or GeneMANIA 41 . However, this knowledge is far from complete and in most cases only contains a small subset of the real interactome. The quality assessment of the inferred networks can be carried out by comparing the inferred interactions to those that have been validated. This comparison results in a confusion matrix as presented in Table 2. The performance of an algorithm can then be summarized using ROC (false positive rate versus true positive rate) or PR (precision versus recall) curves.
A weakness of this type of comparison is that an edge that is not present in the set of known interactions can either mean that an experimental validation has been tried and did not show any regulatory mechanism or (more likely) has not yet been attempted.
In the following, we ran the nce on i) the 2,901 differentially expressed genes identified in Section "Transcriptomic analysis".

Retrieving known interactions
As previously stated, different sources for protein-protein interactions are available (e.g. I2D, BioGrid database). In this example, we obtained a set of known interactions from the BioGrid database, but the users can chose their preferred database. There are 3,941 unique interactions between the 2,901 differentially expressed genes.

Using differentially expressed genes from TCGAbiolinks workflow
We start this analysis by inferring two gene regulatory networks (the corresponding number of edges are presented in Table 3) for the GBM data set and one gene set for the LGG data.

Epigenetic analysis
The DNA methylation is an important component in numerous cellular processes, such as embryonic development, genomic imprinting, X-chromosome inactivation, and preservation of chromosome stability 42 .
In mammals DNA methylation is found sparsely but globally, distributed in definite CpG sequences throughout the entire genome; however, there is an exception. CpG islands (CGIs) which are short interspersed DNA sequences that are enriched for GC. These islands are normally found in sites of transcription initiation and their methylation can lead to gene silencing 43 .
Thus, the investigation of the DNA methylation is crucial to understanding regulatory gene networks in cancer as the DNA methylation represses transcription 44 . Therefore, the DMR (Differentially Methylation Region) detection can help us investigate regulatory gene networks.
This section describes the analysis of DNA methylation using the Bioconductor package TCGAbiolinks 8 . For this analysis, and due to the time required to perform it, we selected only 10 LGG samples and 10 GBM samples that have both DNA methylation data from Infinium HumanMethylation450 and gene expression from Illumina HiSeq 2000 RNA Sequencing Version 2 analysis (lines 1-56 of the Listing 18 describes how to make the data acquisition). We started by checking the mean DNA methylation of different groups of samples, then performed a DMR in which we search for regions of possible biological significance, (e.g., regions that are methylated in one group and unmethylated in the other). After finding these regions, they can be visualized using heatmaps.

Visualizing the mean DNA methylation of each patient.
It should be highlighted that some pre-processing of the DNA methylation data was done. The DNA methylation data from the 450k platform has three types of probes cg (CpG loci), ch (non-CpG loci) and rs (SNP assay). The last type of probe can be used for sample identification and tracking and should be excluded for differential methylation analysis according to the ilumina manual. Therefore, the rs probes were removed (see Listing 18 lines 68). Also, probes in chromosomes X, Y were removed to eliminate potential artifacts originating from the presence of a different proportion of males and females 45 . The last pre-processing steps were to remove probes with at least one NA value (see Listing 18 lines 65).
After this pre-processing step and using the function TCGAvisualize_meanMethylation function, we can look at the mean DNA methylation of each patient in each group. It receives as argument a summarizedExperiment object with the DNA methylation data, and the arguments groupCol and subgroupCol which should be two columns from the sample information matrix of the summarizedExperiment object (accessed by the colData function) (see Listing 18 lines 70-74). Searching for differentially methylated CpG sites. The next step is to define differentially methylated CpG sites between the two groups. This can be done using the TCGAanalyze_DMR function (see Listing 19). The DNA methylation data (level 3) is presented in the form of beta-values that uses a scale ranging from 0.0 (probes completely unmethylated) up to 1.0 (probes completely methylated).
To find these differentially methylated CpG sites, first, the function calculates the difference between the mean DNA methylation (mean of the beta-values) of each group for each probe. Second, it test for differential expression between two groups using the Wilcoxon test adjusting by the Benjamini-Hochberg method. Arguments of TCGAanalyze_DMR was set to require a minimum absolute beta-values difference of 0.25 and an adjusted p-value of less than 10 −2 . After these tests, a volcano plot (x-axis: difference of mean DNA methylation, y-axis: statistical significance) is created to help users identify the differentially methylated CpG sites and return the object with the results in the rowRanges. Figure 9 shows the volcano plot produced by Listing 19. This plot aids the user in selecting relevant thresholds, as we search for candidate biological DMRs.

Listing 19. Finding differentially methylated CpG sites
To visualize the level of DNA methylation of these probes across all samples, we use heatmaps that can be generate by the bioconductor package complexHeatmap 30 . To create a heatmap using the complexHeatmap package, the user should provide at least one matrix with the DNA methylation levels. Also, annotation layers can be added and placed at the bottom, top, left side and right side of the heatmap to provide additional metadata description. The Listing 20 shows the code to produce the heatmap of a DNA methylation datum (Figure 10).

Listing 20. Creating heatmaps for DNA methylation using ComplexHeatmap
Motif analysis. Motif discovery is the attempt to extract small sequence signals hidden within largely non-functional intergenic sequences. These small sequence nucleotide signals (6-15 bp) might have a biological significance as they can be used to control the expression of genes. These sequences are called Regulatory motifs. The bioconductor package rGADEM 46,47 provides an efficient de novo motif discovery algorithm for large-scale genomic sequence data.
The user may be interested in looking for unique signatures in the regions defined by 'differentially methylated' to identify candidate transcription factors that could bind to these elements affected by the accumulation or absence of DNA methylation. For this analysis we use a sequence of 100 bases before and after the probe location (See lines 6-8 in the Listing 21). An object will be returned which contains all relevant information about your motif analysis (i.e., sequence consensus, pwm, chromosome, p-value, etc).

Integrative (Epigenomic & Transcriptomic) analysis
Recent studies have shown that providing a deep integrative analysis can aid researchers in identifying and extracting biological insight from high through put data 42, 53,54 . In this section, we will introduce a Bioconductor package called ELMER to identify regulatory enhancers using gene expression + DNA methylation data + motif analysis. In addition, we show how to integrate the results from the previous sections with important epigenomic data derived from both the ENCODE and Roadmap.

Integration of DNA methylation & gene expression data.
After finding differentially methylated CpG sites, one might ask whether nearby genes also undergo a change in expression either an increase or a decrease. DNA methylation at promoters of genes has been shown to be associated with silencing of the respective gene. The starburst plot is proposed to combine information from two volcano plots and is applicable for studies of DNA methylation and gene expression 55 . Even though, being desirable that both gene expression and DNA methylation data are from the same samples, the starburst plot can be applied as a meta-analysis tool, combining data from different samples 56 . We used the TCGAvisualize_starburst function to create a starburst plot. The log 10 (FDR-corrected P value) for DNA methylation is plotted on the x axis, and for gene expression on the y axis, for each gene. The horizontal black dashed line shows the FDR-adjusted P value of 10 -2 for the expression data and the vertical ones shows the FDR-adjusted P value of 10 -2 for the DNA methylation data. The starburst plot for the Listing 23 is shown in Figure 13. While the argument met.p.cut and exp.p.cut controls the black dashed lines, the arguments diffmean.cut and logFC.cut will be used to highlight the genes that respects these parameters (circled genes in Figure 13). For the example below we set higher p.cuts trying to get the most significant list of pair gene/probes. But for the next sections we will use exp.p.cut = 0.01 and logFC.cut = 1 as the previous sections.

Listing 23. Starburst plot for comparison of DNA methylation and gene expression
ChIP-seq analysis. ChIP-seq is used primarily to determine how transcription factors and other chromatin-associated proteins influence phenotype-affecting mechanisms. Determining how proteins interact with DNA to regulate gene expression is essential for fully understanding many biological processes and disease states. The aim is to explore significant overlap datasets for inferring co-regulation or transcription factor complex for further investigation. A summary of the association of each histone mark is shown in Table 4. Besides, ChIP-seq data exists in the ROADMAP database and can be obtained through the AnnotationHub package 57 or from Roadmap web portal. The Table 5 shows the description for all the Roadmap files that are available through AnnotationHub. After obtaining the ChIP-seq data, we can then identify overlapping regions with the regions identified in the starburst plot. The narrowPeak files are the ones selected for this step. For a complete pipeline with Chip-seq data, bioconductor provides excellent tutorials to work with ChIP-seq and we encourage our readers to review the following article 66 .
The first step shown in Listing 24 is to download the chip-seq data. The function query received as argument the annotationHub database (ah) and a list of keywords to be used for searching the data, EpigenomeRoadmap is selecting the roadmap database, consolidated is selecting only the consolidate epigenomes, brain is selecting the brain samples, E068 is one of the epigenomes for brain (keywords can be seen in the summary table) and narrowPeak is selecting the type of file. The data downloaded is a processed data from an integrative Analysis of 111 reference human epigenomes 67 .

Listing 24. Download chip-seq data
The Chipseeker package 68 implements functions that uses Chip-seq data to retrieve the nearest genes around the peak, to annotate genomic region of the peak, among others. Also, it provides several visualization functions to summarize the coverage of the peak, average profile and heatmap of peaks binding to TSS regions, genomic annotation, distance to TSS and overlap of peaks or genes.
After downloading the histone marks (see Listing 24), it is useful to verify the average profile of peaks binding to hypomethylated and hypermethylated regions, which will help the user understand better the regions found. Listing 25 shows an example of code to plot the average profile. Figure 14 shows the result.
To help the user better understand the regions found in the DMR analysis, we downloaded histone marks specific for brain tissue using the AnnotationHub package that can access the Roadmap database (Listing 24). Next, the Chipseeker was used to visualize how histone modifications are enriched onto hypomethylated and hypermethylated regions, (Listing 25). The enrichment heatmap and the average profile of peaks binding to those regions is shown in Figure 14 and Figure 15 respectively.
The hypomethylated and hypermethylated regions are enriched for H3K4me3, H3K9ac, H3K27ac and H3K4me1 which indicates regions of enhancers, promoters and increased activation of genomic elements. However, these regions are associated neither with transcribed regions nor Polycomb repression as the H3K36me3 and H3K27me3 heatmaps does not show an enrichment nearby the position 0, and the average profile also does not show a peak at position 0.

Listing 25. Average profile plot
To annotate the location of a given peak in terms of genomic features, annotatePeak assigns peaks to genomic annotation in "annotation" column of the output, which includes whether a peak is in the TSS, Exon, 5' UTR, 3' UTR, Intronic or Intergenic location.  H3K27me3, H3K36me3, H3K4me1, H3K4me3, H3K9ac, and H3K9me3. The figure indicates that the differentially methylated regions overlaps regions of enhancers, promoters and increased activation of genomic elements.

Identification of Regulatory Enhancers.
Recently, many studies suggest that enhancers play a major role as regulators of cell-specific phenotypes leading to alteration in transcriptomes realated to diseases [69][70][71][72] . In order to investigate regulatory enhancers that can be located at long distances upstream or downstream of target genes Bioconductor offer the Enhancer Linking by Methylation/Expression Relationship (ELMER) package. This package is designed to combine DNA methylation and gene expression data from human tissues to infer multi-level cis-regulatory networks. It uses DNA methylation to identify enhancers and correlates their state with expression of nearby genes to identify one or more transcriptional targets. Transcription factor (TF) binding site analysis of enhancers is coupled with expression analysis of all TFs to infer upstream regulators. This package can be easily applied to TCGA public available cancer data sets and custom DNA methylation and gene expression data sets 73 .
2. Identify distal enhancer probes with significantly different DNA methylation level in control group and experiment group.
3. Identify putative target genes for differentially methylated distal enhancer probes.
4. Identify enriched motifs for the distal enhancer probes which are significantly differentially methylated and linked to putative target gene.
5. Identify regulatory TFs whose expression associate with DNA methylation at motifs.

H3K36me3
This section shows how to use ELMER to analyze TCGA data using as example LGG and GBM samples.

Preparing the data for ELMER package.
After downloading the data with TCGAbiolinks package, some steps are still required to use TCGA data with ELMER. These steps can be done with the function TCGAprepare_elmer. This function for the DNA methylation data will remove probes with NA values in more than 20% samples and remove the annotation data, for RNA expression data it will take the log2(expression + 1) of the expression matrix in order to linearize the relation between DNA methylation and expression. Also, it will prepare the row names of the matrix as required by the package.

Listing 27. Preparing TCGA data for ELMER's mee object
Finally, the ELMER input is a mee object that contains a DNA methylation matrix, an gene expression matrix, a probe information GRanges, the gene information GRanges and a data frame summarizing the data. It should be highlighted that samples without both the gene expression and DNA methylation data will be removed from the mee object.
By default the function fetch.mee that is used to create the mee object will separate the samples into two groups, the control group (normal samples) and the experiment group (tumor samples), but the user can relabel the samples to compare different groups. For the next sections, we will work with both the experimental group (LGG) and control group (GBM).

Listing 29. Running ELMER analysis
When ELMER Identifies the enriched motifs for the distal enhancer probes which are significantly differentially methylated and linked to putative target gene, it will plot the Odds Ratio (x axis) for the each motifs found. The list of motifs for the hyper direction (probes hypermethylated in LGG group compared to the GBM group) is found in the Figure 17. We selected motifs that had a minimum incidence of 10 in the given probes set and the smallest lower boundary of 95% confidence interval for Odds Ratio of 1.1. These both values are the default from the ELMER package.
The analysis found 16 enriched motifs for the hyper direction and no enriched motifs for the hypo direction. After finding the enriched motifs, ELMER identifies regulatory transcription factors (TFs) whose expression is associated with DNA methylation at motifs. ELMER automatically creates a TF ranking plot for each enriched motif. This plot shows the TF ranking plots based on the association score (−log(Pvalue)) between TF expression and DNA methylation of the motif. We can see in Figure 18 that the top three TFs that are associated with that FOX motif are FOXD3, HMGA2 and HIST1H2BH.
The output of this step is a data frame with the following columns: 1. motif: the names of motif.
2. top.potential.TF: the highest ranking upstream TFs which are known recognized the motif.
3. potential.TFs: TFs which are within top 5% list and are known recognized the motif. top5percent: all TFs which are within top 5% list considered candidate upstream regulators Also, for each motif we can take a look at the three most relevant TFs. For example, Figure 19 shows the average DNA methylation level of sites with the FOX motif plotted against the expression of the TFs FOXD3, HIST1H1D, HIST1H2BH and HMGA2. We can see that the experiment group (LGG samples) has a lower average methylation level of sites with the FOX motif plotted and a higher average expression of the TFs. And for each relevant TF we will use the clinical data to analyze the survival curves for the 30% patients with higher expression of that TF versus the 30% with lower expression. The code below shows that analysis. TCGAanalyze_survival(clinical, "tf_groups", 25 legend=paste0(TF," Exp level"), 26 filename = paste0(TF,".pdf")) 27 } 28 29 # get clinical patient data for GBM samples 30 gbm_clin <-GDCquery_clinic("TCGA-GBM","Clinical") 31 32 # get clinical patient data for LGG samples 33 lgg_clin <-GDCquery_clinic("TCGA-LGG","Clinical") 34 35 # Bind the results, as the columns might not be the same, 36 # we will will plyr rbind.fill, to have all columns from both files 37 clinical <-plyr::rbind.fill(gbm_clin,lgg_clin) 38 # Call the function we created 39 TCGAsurvival_TFplot("ELF4",mee,clinical) 40 TCGAsurvival_TFplot("HMGA2",mee,clinical) 41 TCGAsurvival_TFplot("FOXD3",mee,clinical) The Figure 20, shows that the samples with lower expression of some of these TFs have a better survival than those with higher expression.

Conclusion
This workflow outlines how one can use specific Bioconductor packages for the analysis of cancer genomics and epigenomics data derived from the TCGA. In addition, we highlight the importance of using ENCODE and Roadmap data to inform on the biology of the non-coding elements defined by functional roles in gene regulation. We introduced TCGAbiolinks and RTCGAToolbox Bioconductor packages in order to illustrate how one can acquire TCGA specific data, followed by key steps for genomics analysis using GAIA package, for transcriptomic analysis using TCGAbiolinks, dnet, pathview packages and for DNA methylation analysis using TCGAbiolinks package. An inference of gene regulatory networks was also introduced by MINET package. Finally, we introduced Bioconductor packages AnnotationHub, ChIPSeeker, ComplexHeatmap, and ELMER to illustrate how one can acquire ENCODE/Roadmap data and integrate these data with the results obtained from analyzing TCGA data to identify and characterize candidate regulatory enhancers associated with cancer.

Data and software availability
This workflow depends on various packages from version 3.2 of the Bioconductor project, running on R version 3.2.2 or higher. It requires a number of software packages, including AnnotationHub, ChIPSeeker, ELMER, ComplexHeatmap, GAIA, rGADEM, MotIV, MINET, RTCGAToolbox and TCGAbiolinks. Version numbers for all packages used are in section "Session Information". Listing 32 shows how to install all the required packages.
The "broadinstitute website" should be changed to "Broad Institute website" and include a link.
The "nce" is not defined as far as I could see, and the sentence just before the "Retrieving known interactions" seems unfinished. It is not clear to me how the information that "There are 3,941 unique interactions between the 2,901 differentially expressed genes" was obtained, since the function in the code chunk above is not run until later.
It would be easier to get an overview of the text of the headers were more consistent. For example, the gene regulatory network section has subheaders that are larger than the main section header.
Finally, as mentioned by one of the other reviewers, the paper would benefit from modifications to the text in order to make it more readable. There are many grammatical errors, long and complicated sentences, and sentences that are unfinished and thus difficult to understand. I believe that if the authors can make sure that all the code in the workflow runs as expected, and the text and code is clarified and harmonized throughout the manuscript, this can be a very valuable addition to the literature.

I have read this submission. I believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard, however I have significant reservations, as outlined above.
No competing interests were disclosed. As it has been pointed out already by the first reviewer, it is important to verify that the pipeline is updated according to the data migration to the GCD server.
Apart from this, I find this manuscript very nice and the workflow an important contribution in the field, allowing to make accessible large datasets with genomics profiling of cancer patients to the community and even for not advanced R-users. I would like to compliment with the authors to make such a comprehensive workflow available.
Some minor points before final publication is indexed: The main issue with this manuscript is that it will need an extensive revision by a native English speaker since some sentences are hard to read and very often they risk to convey the wrong message or to be misinterpreted. One example from the abstract, "we provide a series of F1000Research message or to be misinterpreted. One example from the abstract, "we provide a series of biologically focused integrative downstream analyses of different molecular data" includes to many adjectives. Few lines below "we provide a workplan to identify candidate biologically relevant functional epigenomic elements associated with cancer", this sentence does not sound that right in English.
Be also careful in the style use to call the software or Bioconductor packages in the text since sometimes they are italics and other times they are not.
There is also a lot of redundancy on some of the statements in the introduction that if removed will improve the readability. For example the sentence at the end of the first paragraph of the introduction starting with " Here we describe a comprehensive workflow that integrates many bioconductor packages..." is probably not needed since the same concept is provided at the end of the introduction in the same page.
In the Introduction when TCGA data are illustrated I would also mention metastatic samples (they are there, isn't it? or may be I don't recall correctly...) I suggest to change the title of the first section of the methods from 'Experimental data' to 'Access to the data' or something similar so that it will better reflect the content of the paragraph.
I felt that sometimes the manuscript sounds more a user guide than an article, so I would suggest taking as an example the section on methylation to improve the discussion of the analyses outcome and aims in the other paragraphs.
Also please be careful to make things accessible to any readers for example at page 4, the authors introduce the GISTIC data without explaining them in details.
The sentence at page 4 "the data used in this workflow are published data and freely available" might be redundant so the authors could remove it.
Page 4 summarizedExperiment object is missing the reference to the original paper about summarizedExperiment format.
Page 5: in the text the authors illustrate assay, rowRanges and colData but in the example the order is different, i.e. assay, colData, rowRanges -I would suggest to keep consistency in the order between main text and examples (it also holds for other places in the paper).
I am not sure I can understand the sentence at page 5, "the users should use all the follow up files in your analyses, not just the latest version". Not sure that Listing 5 also fits there... why not mention it before together with the others summ.exp. accessors.
In writing it is important that the authors really pay attention to convey the message that it is an F1000Research "Page 19, I might have missed something but I cannot find in maintext any explanation of the DEA pipeline and more comments would be needed there for a new user to really appreciate it and its importance." For DEA pipeline we used the TCGAanalyze_DEA function from our package TCGAbiolinks, that allows user to perform Differentially expression analysis (DEA), using edgeR package to identify differentially expressed genes (DEGs). It is possible to do a two-class analysis. TCGAanalyze_DEA performs DEA using following functions from edgeR: edgeR::DGEList converts the count matrix into an edgeR object. edgeR::estimateCommonDisp each gene gets assigned the same dispersion estimate. edgeR::exactTest performs pair-wise tests for differential expression between two groups. edgeR::topTags takes the output from exactTest(), adjusts the raw p-values using the False Discovery Rate (FDR) correction, and returns the top differentially expressed genes. "On PEA, I was wondering why for example the R package for Reactome was not integrated here and if the authors could comment on this." Thank you for this suggestion we are working to integrate Reactome as source of genes annotated within pathways.
Our PEA provided only one plot to show Gene ontologies and pathways enriched by a list of genes to have an overview of top biological functions and pathways altered by molecules inside the gene signature. We focused mainly on Bioconductor packages, but in particular we are interested to add wrapping to some functionalities of the ReactomePA package.
"Page 23, why for protein-protein interactions only BioGrid is used and not for example a more comprehensive resource such as I2D (Interologous Interaction Database)." Different sources for protein-protein interactions are available. We used BioGrid as an example, but the users can chose their preferred database such as I2D

F1000Research
This review comes at a very inopportune moment. The entire software pipeline is based on the TCGAbiolinks tool kit, which downloads files from the TCGA DCC service. Unfortunately, just as this paper was being sent for review, the NCI began their migration to the GDC service. This means that all of the data services at DCC TCGA data portal (https://tcga-data.nci.nih.gov) are no longer active, and users are being directed to the GDC NCI data portal (https://gdc.nci.nih.gov/). This means that the TCGAbiolinks tool kit is broken. My attempts to run some of the examples listed in the paper where stopped by this issue. Many parts of the TCGAbiolinks will have to be fixed and/or re-written to deal with this change. This isn't the fault of the authors, but it does render the example and pipelines described in the paper inoperable.
There are a few technical issues that need to be addressed. The package 'RTCGAtoolbox' is mentioned, but actually its name is 'RTCGAToolbox'. (lower-case T to upper-case), and because the bioconductor website is case-sensitive the provided URL ( http://bioconductor.org/packages/RTCGAtoolbox/ ) doesn't work. On page 3, 'TCGA Wikipedia' should be 'TCGA Wiki'. In the methods section on page 3, where the various data levels of TCGA are explained, on the most important aspects about the different levels is of access requirements. Data levels 1 and 2 cannot be accessed without dbGap, while level 3 and 4 are generally available for public access.
There is also the question of whether or not this paper fits the criteria of the Software Tools Article guidelines. The authors have presented a loose set of examples that utilize various existing, and previously published, tools. In the guidelines of F1000Research's Software Tool Articles, the criteria for a paper is "We welcome the description of new software tools. A Software Tool Article should include the rationale for the development of the tool and details of the code used for its construction." The authors refer to the code included in the paper as a workflow, but reads more like a series of point example that the reader can copy and alter for their own research. In introducing a new piece of software, one expects that the authors of the article are responsible for the software being presented. And while this article mentioned an extensive number of R packages, I believe only one of them, TCGAbiolinks, was written by the authors, and was published last year. While the analysis they present is very detailed and covers an expansive number of topics, I'm not sure if this should be classified as 'Software Tool Article'. On my first read through, my assumption was that the authors were responsible for all of the tools mentioned in the abstract. And this wasn't clarified in the text. For example, in the conclusion section on page 43, they include the phrase 'We introduced TCGAbiolinks and RTCGAtoolbox bioconductor packages in order to illustrate how one can acquire TCGA specific data', it wasn't until later, that I realized that the author of the RTCGAToolbox wasn't in the author list of the paper.
Of course, this is at the discretion of the editors. I would feel that this material would better be presented as a review article covering the different methods of TCGA data analysis, with better notation and attributions about the source of different pieces of software.
I have read this submission. I believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard, however I have significant reservations, as outlined above.
No competing interests were disclosed. Thank you for your comments and suggestions. We made several changes in the version 2 of the workflow, some of the changes and answers to your points are below: TCGAbiolinks package was entirely redesigned to query, download and prepare data from the GDC NCI data portal ( ). All the codes are working and we will https://gdc.nci.nih.gov/ submit very soon the workflow to Bioconductor. For the moment, the RMarkdown can be found in our Github repository ( ). https://github.com/BioinformaticsFMRP/f1000_TCGA_Workflow/blob/master/f1000.Rmd We also highlighted the difference between open (TCGA level 3 and 4 data) and controlled data (TCGA level 1 and 2 data) and we added some useful sources that can help the user request access to the controlled data.
Despite the description of a software article, we decided by this type of article simply because there were already other workflows with that choice and we did not find any other possibility that best described this type of article. In addition, our main focus was on using the tools (in a reasonable time) rather than analyzing the results. The analysis itself can be verified with some articles of our group, for example, the DNA methylation analysis was We corrected links and references such as . RTCGAToolbox No competing interests were disclosed. Competing Interests: