TCGA Workflow: Analyze cancer genomics and epigenomics data using Bioconductor packages [version 1; peer review: 1 approved, 1 approved with reservations]

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 opensource 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 comprehensive tool that provides a complete integrative analysis harnessing 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 downstream 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 and by using Roadmap and ENCODE data, we provide a workplan to identify candidate 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, Open Peer Review


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. In order to get insight in the 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, the understanding of 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. Here we describe a comprehensive workflow that integrates many Bioconductor packages in order to analyze and integrate the molteplicity of molecular observation layers in large scale cancer dataset.
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 that include: exome (variant analysis), single nucleotide polymorphism (SNP), DNA methylation, transcriptome (mRNA), microRNA (miRNA), proteome and clinical information. Sample types available at TCGA are: primary solid tumors, recurrent solid tumors, blood derived normal and tumor, and solid tissue normal 2 .
The next steps describes how one could use TCGAbiolinks & RTCGAtoolbox to download clinical, genomics, transcriptomics, epigenomics data, as well as subtype information and GISTIC results (identified genes targeted by somatic copy-number alterations (SCNAs) that drive cancer growth). Just to reiterate, the data used in this workflow are published data and freely available.
Downloading data from TCGA data portal. The Bioconductor package TCGAbiolinks 7 has three main functions TGCAquery, TCGAdownload and TCGAprepare that should sequentially be used to respectively search, download and load the data as an R object.
TGCAquery searches in a pre-processed TCGA database and returns a summary table with the found files, samples, version and other useful information. The most important TGCAquery arguments are tumor which receives one or multiple tumor types (USC, LGG, SKCM, KICH, CHO, etc), platform which receives the platform (HumanMethylation27, Genome_Wide_SNP_6, IlluminaHiSeq_RNASeqV2, etc), version which receives the version of the data to be downloaded if the user wants an older version and samples which receives a list of TCGA barcodes (ex. "TCGA-CS-4938") to filter the search results. A complete list of possible entries for arguments can be found in the TCGAbiolinks vignette. Lines 6 and 13 of Listing 1 show an example of this function.
After searching, the user will be able to download the data with TCGAdownload. An important feature of this function is the ability to filter the data using the arguments type if the user wants to specify file tumor type and samples if user wants to specify samples (list of TCGA barcodes). For example, lines 15 and 18 of Listing 1 are used to select a specific tumor type to download and prepare the data respectively. The platforms and their possible inputs for the type argument is shown below: Finally, TCGAprepare transforms the downloaded data into a summarizedExperiment object or a data frame. If summarizedExperiment is set to TRUE, TCGAbiolinks will add metadata to the object in order to help the user when working with the data. Also, if the user sets the argument add.subtype to TRUE the summarizedExperiment will receive subtype information defined by The Cancer Genome Atlas (TCGA) Research Network reports (the full list of papers can be seen in TCGAquery_subtype section in TCGAbiolinks vignette), Likewise, if the user sets the argument add.clinical to TRUE the summarizedExperiment will receive clinical information. Lines 8-11 and 18-22 of Listing 1 illustrates this function.

Listing 1. Downloading DNA methylation and gene expression data from TCGA 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) 8,9 . An example is shown in Listing 2.

Listing 2. summarizedExperiment accessors
Clinical data can be obtained using the function TCGAquery_clinical which can be used as described in Listing 3. This function has three arguments tumor, clinical_data_type and samples. The clinical_data_type argument is always required and should be accompanied by at least one of the other two parameters. Examples for the argument clinical_data_type are: "clinical_drug", "clinical_patient", and "clinical_radiation" (a complete list and description can be found in the section 'Working with clinical data.' of the TCGAbiolinks vignette).
An important note about the clinical data is that follow-up data for TCGA patients are contained in the 'clinical_ follow_up' files for each cancer type and to obtain all available disease progression information, the users should use all the follow_up files in your analyses, not just the latest version.

Listing 3. Downloading clinical data with TCGAbiolinks
Mutation information is stored in Mutation Annotation Format (MAF) files which contain different mutation types (somatic or germline) and states (validated or putative). A summary of all the Mutation Annotation Format (MAF) can be accessed at TCGA wiki. To download these data using TCGAbiolinks, TCGAquery_maf function is provided. It will download the non-obsolete tables from TCGA wiki, remove the protected entries and ask the user which file s/he wants to download (see Listing 4). It will then download and return a data frame with the data. 20 # To access the data you should use the getData function 21 # or simply access with @ (for example gbm.data@Clinical) 22 gbm.mut <− getData(gbm.data,"Mutations") 23 gbm.clin <− getData(gbm.data,"Clinical") 24 gbm.gistic <− getData(gbm.data,"GISTIC") Listing 6. Downloading TCGA data files with RTCGAtoolbox Finnaly, RTCGAtoolbox can access level 4 data, which can be handy when the user requires GISTIC results. GISTIC is used to identify genes targeted by somatic copy-number alterations (SCNAs) 16 (see Listing 7).

Genomic analysis
Copy number variations (CNV) has 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. CNV 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. Data for selected samples were downloaded and prepared in two separate rse objects (RangedSummarizedExperiment).

Identification of recurrent CNV in cancer.
Cancer related CNV have to be present in many of the analyzed genomes. The most significant recurrent CNV were identified using GAIA 17 , an iterative procedure where a statistical hypothesis framework is extended to take into account within-sample homogeneity. GAIA is based on a conservative permutation test allowing the estimation of the probability distribution of the contemporary mutations expected for non-driver markers. Segmented data retrieved from TCGA were used to generate a matrix including all needed information about the observed aberrant regions. Furthermore, GAIA requires genomic probes metadata (specific for each CNV technology), that can be downloaded from broadinstitute website. colnames(cnvMatrix) <-c("Sample.Name", "Chromosome", "Start", "End", "Num.of.Markers", "Aberration") 19 # Substitute Chromosomes "X" and "Y" with "23" and "24" 20 xidx <-which(cnvMatrix$Chromosome=="X") 21 yidx <-which(cnvMatrix$Chromosome=="Y") 22 cnvMatrix save(results, RecCNV, threshold, file = paste0(cancer,"_CNV_results.rda")) 85 } Listing 9. 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.

Listing 10. Gene annotation of recurrent CNV
Visualizing multiple genomic alteration events. In order to visualize multiple genomic alteration events we recommend using OncoPrint plot which is provided by bioconductor package complexHeatmap 18 . The Listing 11 shows how to download mutation data using TCGAquery_maf (line 4), then we filtered the genes to obtain genes with mutations found among glioma specific pathways (lines 6 -12). The following steps prepared the data into a matrix to fit oncoPrint function. We defined SNPs as blue, insertions as green and deletions as red. The upper barplot indicates the number of genetic mutation per patient, while the right barplot shows the number of genetic mutations per gene. Also, it is possible to add annotations to rows or columns. In the columns case, if an insertion is made at the top, will remove the barplot. The final result for adding the annotation to the bottom is highlighted in Figure 2.

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 mutations (selecting curated genetic variations retrieved from TCGA that are identified in at least two tumor samples) in LGG (see Listing 13). 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.

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 object (RangedSummarizedExperiment) saving them as an R object with a filename including both the name of the cancer and the name of the plaftorm used for gene expression data (see Listing 13).

Listing 13. Searching, downloading and preparing RNA-seq data with TCGAbiolinks
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 14). In this way 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.
This function implements Within-lane normalization procedures to adjust for GC-content effect (or other genelevel effects) on read counts: loess robust local regression, global-scaling, and full-quantile normalization 19 and between-lane normalization procedures to adjust for distributional differences between lanes (e.g., sequencing depth): global-scaling and full-quantile normalization 20 .
EA: enrichment analysis. In order to understand the underlying biological process from DEGs we performed an enrichment analysis using TCGAanalyze_EA_complete function (see Listing 15).

Listing 15. 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 (GO:biological process, GO:cellular component, and GO:molecular function. The Figure 5 shows canonical pathways significantly overrepresented (enriched) by the DEGs. The most statistically significant canonical pathways identified in DEGs list are listed 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).

Listing 16. Pathways enrichment analysis with pathview package
The red genes are up-regulated and the green genes are down-regulated in the LGG samples compared to GBM.

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 22 , CLR 23 , MRNET 24 and C3NET 25 . 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 26 and c3net 25 , 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 27 or GeneMANIA 28 . 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. Different quality measures can then be computed such as the false positive rate  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
We obtained a set of known interactions from the BioGrid database. There are 3,941 unique interactions between the 2,901 differentially expressed genes.

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 29 . In mammals, DNA methylation is found sparsely but globally, distributed in definite CpG sequences throughout the entire genome. There is however an exception, the CpG islands (CGIs) which are short interspersed DNA sequences with are enriched for GC. These CpG islands are normally found in sites of transcription initiation and their methylation can lead to gene silencing 30 Thus, the investigation of the DNA methylation is crucial to understanding regulatory gene networks in cancer as the DNA methylation represses transcription 31 . 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 7 . 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-7 of the Listing 17 describes how to make the data acquisition). We started by checking the mean DNA methylation of different groups of samples, then a DMR (Differentially methylated region) analysis is performed in which we search for regions that have possible biological significance, for example, 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 have 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 17 lines 43). Also, probes in chromosomes X, Y were removed to eliminate potential artifacts originating from the presence of a different proportion of males and females 32 . The last pre-processing steps were to remove probes with at least one NA value (see Listing 17 lines 40).
After this pre-processing step, using the function TCGAvisualize_meanMethylation provided we can take a 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 17 lines 46-50). 39 # remove probes with NA (similar to na.omit) 40 met <-subset(met,subset = (rowSums(is.na(assay(met))) == 0)) 41 42 # remove probes in chromossomes X, Y and NA 43 met <-subset(met,subset = !as.character(seqnames(met)) %in% c("chrNA","chrX","chrY")) 44 45 Listing 17. Visualizing the DNA mean methylation of groups Figure 8 illustrates a mean DNA methylation plot for each sample in the GBM group (140 samples) and a mean DNA methylation for each sample in the LGG group. Genome-wide view of the data highlights a difference between the groups of tumors (p-value = 6.1 × 10 −06 ).
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 18). 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, it calculates the difference between the mean DNA methylation (mean of the beta-values) of each group for each probe. Second, it tests 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 18. This plot aids the user in selecting relevant thresholds, as we search for candidate biological DMRs.

Listing 18. 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 18 . 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 19 shows the code to produce the heatmap of a DNA methylation data ( Figure 10).

Listing 19. 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 sequence short nucleotide sequences (6-15 bp) might have a biological significance as it can be used to control the expression of genes. These sequences are called Regulatory motifs. The bioconductor package rGADEM 33,34 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 20). An object will be returned which contains all relevant information about your motif analysis (sequence consensus, pwm, chromosome, pvalue…).
Using bioconductor package motifStack 35 it is possible to generate a graphic representation of multiple motifs with different similarity scores (see Figure 11).

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 29,40,41 . 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 possible question one might ask is whether nearby genes also undergo a change in its expression, either an increase or a decrease. DNA methylation at promoters of genes have 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 applied for a study of DNA methylation and gene expression 42 . 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 43 .
The function TCGAvisualize_starburst creates a Starburst plot for comparison of DNA methylation and gene expression. 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 22 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.

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 exist in the ROAD-MAP database and can be obtained through the AnnotationHub package 44 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 53 .
The first step shown in Listing 23 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 (a table for the list is found in this summary table) 54 .and narrowPeak is select ing the type of file. The data downloaded are processed data from an integrative Analysis of 111 reference human epigenomes 54 .

Listing 24. 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.  H3K27me3, H3K36me3, H3K4me1, H3K4me3, H3K9ac,  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 related to diseases [56][57][58][59] . In order to investigate regulatory enhancers that can be located at long distances upstream or downstream of target genes bioconductor offers 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 60 .
ELMER analysis has five main steps: 1. Identify distal enhancer probes on HM450K.
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.
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 26. Preparing TCGA data for ELMER's mee object
Finally, the ELMER input is a mee object that contains a DNA methylation matrix, a 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 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 two groups the experiment group (LGG) and control samples (GBM).

Listing 28. 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 found for the hyper direction (probes hypomethylated in LGG group compared to the GBM group) is found in Figure 17. To select the motifs we select the 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 values are the default from the ELMER package.
The analysis found 14 enriched motifs for the hyper direction and no enriched motifs for the hypo direction.
After finding these list of enriched motifs, ELMER identifies regulatory TFs whose expression associate with DNA methylation at motifs and for each enriched motif a TF ranking plot is created automatically by ELMER. This plot shows the TF ranking plots based on the score (−log(Pvalue)) of association between TF expression and DNA methylation of the motif. We can see in Figure 18 that the top three associated TFs that are associated with that AP1 motif are POLR3K, DLX3 and NEUROD2.
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 on the three most relevant transcription factors. For example, for the AP1 motif the average DNA methylation level of sites with the AP1 motif plotted against the expression of the transcription factors WT1, ZNF208, ATF4 and DDX5 is show in Figure 19. We can see that the experiment group (GBM samples) has a lower average methylation level of sites with the AP1 motif plotted and a higher expression of the transcription factors.

Listing 29. Visualizing the average DNA methylation level of sites with a chosen motif vs TF expression
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 transcription factor versus the 30% with lower expression. The code below shows that analysis.

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 with the results obtained from analyzing TCGA data in order 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 31 shows how to install all the required packages.

Listing 31. Installing packages
All data used in this workflow is freely available and can be accessed using a R/Bioconductor package. There are two main sources of data: The Cancer Genome Atlas (TCGA) and a supplementary data repository with processed datasets from the Roadmap Epigenomics Project and from The Encyclopedia of DNA Elements (ENCODE) project 54 . For the first, a summary of the data available can be seen in https://tcga-data.nci.nih.gov/tcga/ and its data can be accessed using the R/Bioconductor TCGAbiolinks package. For the second, a summary of the data available can be found in this spread sheet and the data can be accessed using the R/Bioconductor AnnotationHub package. 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 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. First, we would like to thank you for your review and for providing a detailed feedback for our workflow. We have made several changes in this new version (2) of the workflow. The changes and responses of all authors to these points are described below: "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." TCGAbiolinks package was entirely redesigned to query, download and prepare data from the GDC NCI data portal (https://gdc.nci.nih.gov/). All the codes are working and we will 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 ).
○ ○ "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." We applied the following pattern: name of the packages normal with a link to the package, name of functions and objects in italics.
○ ○ "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." 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 described in  In writing it is important that the authors really pay attention to convey the message that it is an integrative approach ... at the very first reading of the manuscript I felt that it was more a list of tools. and the integrative part of it was missing" We agree with you. The article starts with specific analysis, which might cover 65% of the content and unfortunately it might lead to the feeling of the missing integrative approach. Also, we decided not to go much deeper into the analyzes because they already exist in the referenced articles and we focus more on the execution of the analyzes with the tools of the bioconductor project, something that was not covered by the cited articles, which might also give the feeling of a "list of tools". But, the integrative analysis is shown strongly in the last two sections. We may highlight that ELMER, for example, uses DNA methylation, gene expression, histone marks, motif information and clinical data for its analysis. Also, after the DNA methylation analysis available in TCGAbiolinks package, we do a motif analysis on those regions and integrated the roadmap histone marks data to evaluate them. Finally, another integrative analysis presented in the TCGAbiolinks package is the starburst plot which integrates the differential expression analysis results with the DNA methylation results and we try to identify if the nearest gene might have been affected by the change of the DNA methylation.
○ ○ "Page 18: a reference is missing for the AAIC method and an explanation to the reader about the choice of the 0.6 correlation cutoff." Array-array intensity correlation plot (AAIC) is a re-adaptation of a affyQCReport::correlationPlot working with data from gene expression of class RangedSummarizedExperiment as output of GDCprepare.Reference from https://www.bioconductor.org/packages/devel/bioc/vignettes/affyQCReport/inst/doc/affyQCRep . AAIC 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.

○
The 0.6 threshold came out from unsupervised hierarchical clustering with ward methodology that obtained distinct groups of samples and first one had pearson correlation lower than 0.6 considering them as possible outliers.
○ ○ "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. 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. 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 (https://gdc.nci.nih.gov/). All the codes are working and we will 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 described in