---
title: 'TCGA Workflow: Analyze cancer genomics and epigenomics data using Bioconductor
  packages'
date: "`r Sys.Date()`"
author:
- name: Tiago Chedraoui Silva
  affiliation: 
  - &id Department of Genetics, Ribeirao Preto Medical School, University of Sao Paulo, Ribeirao Preto, Brazil
  email: tiagochst@gmail.com
- name: Antonio Colaprico
  affiliation:  
  - &id2 Interuniversity Institute of Bioinformatics in Brussels , Brussels, Belgium
  - &id3 Machine Learning Group (MLG), University Libre de Bruxelles, Brussels, Belgium
  email: antonio.colaprico@gmail.com
- name: Catharina Olsen
  affiliation: 
  - *id2
  - *id3
  email: colsen@ulb.ac.be
- name: Fulvio D’Angelo
  affiliation: 
  - Institute for Cancer Genetics, Columbia University Medical Center, New York, New York 10032, USA.
  - BIOGEM Istituto di Ricerche Genetiche ‘G. Salvatore’, Campo Reale, 83031 Ariano Irpino, Italy. 
  email: fulvio.dangelo@biogem.it
- name: Gianluca Bontempi
  affiliation: *id3
  email: gbonte@ulb.ac.be
- name: Michele Ceccarelli
  affiliation: 
  - Department of Science and Technology, University of Sannio, Benevento, Italy.
  - BIOGEM Istituto di Ricerche Genetiche "G. Salvatore", Ariano Irpino, Italy.
  email: ceccarelli@unisannio.it
- name: Houtan Noushmehr
  affiliation: 
  - Department of Neurosurgery, Henry Ford Hospital, Detroit, Detroit, MI, USA
  - *id
  email: houtana@gmail.com
output: 
    BiocStyle::html_document:
        toc: true
        number_sections: false
        toc_depth: 2
        toc_float: true
        highlight: breezedark
        theme: lumen
        df_print: kable
bibliography: bibliography.bib

vignette: >
  %\VignetteIndexEntry{'TCGA Workflow: Analyze cancer genomics and epigenomics data using Bioconductor packages'}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---

<style> body { width: 100%; text-align: justify} </style>

```{r setup, include=FALSE, echo=F, warning= F, message=F}
knitr::opts_chunk$set(
  message = FALSE, 
  warning = FALSE, 
  error = FALSE, 
  tidy = FALSE,
  fig.align = "center", 
  dpi = 150, 
  cache = FALSE,
  progress = FALSE, 
  quite = TRUE
)
```

```{r, echo=FALSE, results="hide", warning=FALSE,message=FALSE}
library(TCGAWorkflow)
```


# Environment

**R version**: `r R.version.string`

**Bioconductor version**: `r BiocManager::version()`

**Package**: `r packageVersion("TCGAWorkflow")`


# About

This workflow is based on the article: [TCGA Workflow: Analyze cancer genomics and epigenomics data using Bioconductor packages](https://f1000research.com/articles/5-1542/v2) [@10.12688/f1000research.8923.2].
Due to time and space limitations,  we downloaded only a subset of the data, for a real analysis please use all data available. The data used in the examples are available in the package `r Biocpkg("TCGAWorkflowData")`. 


## Installation

To be able to execute all the steps of this workflow please install it with the following code:

```{R, eval=FALSE, include=TRUE}
if (!"BiocManager" %in% rownames(installed.packages()))
  install.packages("BiocManager")
BiocManager::install("TCGAWorkflow")
BiocManager::install("TCGAWorkflowData")
```

## Loading packages

At the beginning of each section, the packages required to execute the code will be loaded. However, 
the following packages are required for all sections.

- TCGAWorkflowData: this package contains the data necessary to execute each of the analysis
steps. This is a subset of the downloaded to make the example faster. For a real analysis, please
use all the data available.
- DT: we will use it to visualize the results

```{R, eval=TRUE, message=FALSE, warning=FALSE, include=TRUE}
library(TCGAWorkflowData)
library(DT)
```

```{R, eval=TRUE, message=FALSE, warning=FALSE, include=FALSE}
library(TCGAbiolinks)
```

# Abstract

Biotechnological advances in sequencing have led to an explosion of
publicly available data via large international consortia such as [The
Cancer Genome Atlas (TCGA)](http://cancergenome.nih.gov/), [The
Encyclopedia of DNA Elements (ENCODE)](http://www.encodeproject.org/),
and [The NIH Roadmap Epigenomics Mapping Consortium
(Roadmap)](http://www.roadmapepigenomics.org/). 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](http://www.bioconductor.org/)
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). 

All the  package landing pages used in this workflow can be found through 
the [biocViews interface](http://www.bioconductor.org/packages/release/BiocViews.html#___Software).

**Keywords:** Epigenomics, Genomics, Cancer, non-coding, TCGA, ENCODE, Roadmap, Bioinformatics.

# 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. Here we describe a comprehensive workflow that
integrates many Bioconductor packages in order to analyze and integrate
the multiplicity 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](http://cancergenome.nih.gov/)), The Encyclopedia of DNA
Elements ([ENCODE](http://www.encodeproject.org/)), and The NIH Roadmap
Epigenomics Mapping Consortium
([Roadmap](http://www.roadmapepigenomics.org/)) [@Hawkins]. 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
[@weinstein2013cancer].

-   **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 include immortalized cell lines,
tissues, primary cells and stem cells [@encode2011user].

-   **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 [@Fingerman; @Bernstein].

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
a 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 [@kannan2015public] in
accessing cancer associated data,
[Bioconductor](http://www.bioconductor.org/) represents 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
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](http://www.bioconductor.org/) packages
to work with genomic and epigenomics data.


# Methods
## Access to the data

TCGA data is accessible via the NCI Genomic Data Commons (GDC) [data
portal](https://gdc-portal.nci.nih.gov/), and the
[Broad Institute’s GDAC Firehose](gdac.broadinstitute.org). The GDC Data
Portal provides access to the subset of TCGA data that has been
harmonized against GRCh38 (hg38) using GDC Bioinformatics Pipelines
which provides methods to the standardization of biospecimen and
clinical data, the re-alignment of DNA and RNA sequence data against a
common reference genome build GRCh38, and the generation of derived
data. 

The previously stored data in CGHub, TCGA Data Portal and Broad
Institute’s GDAC Firehose, were provided as different levels or tiers
that were defined in terms of a specific combination of both processing
level (raw, normalized, integrated) and access level (controlled or open
access). Level 1 indicated raw and controlled data, level 2 indicated
processed and controlled data, level 3 indicated Segmented or
Interpreted Data and open access and level 4 indicated region of
interest and open access data. While the TCGA data portal provided level
1 to 3 data, Firehose only provides level 3 and 4. An explanation of the
different levels can be found at [TCGA
Wiki](https://wiki.nci.nih.gov/display/TCGA/Data+level). However, the
GDC data portal no longer uses this based classification model in
levels. Instead, a new data model was created, its documentation can be
found in [GDC
documentation](https://gdc.nci.nih.gov/developers/gdc-data-model/gdc-data-model-components).
In this new model, data can be open or controlled access. While the GDC
open access data does not require authentication or authorization to
access it and generally includes high-level genomic data that is not
individually identifiable, as well as most clinical and all biospecimen
data elements, the GDC controlled access data requires dbGaP
authorization and eRA Commons authentication and generally includes
individually identifiable data such as low-level genomic sequencing
data, germline variants, SNP6 genotype data, and certain clinical data
elements. The process to obtain access to controlled data is found in
[GDC web
site](https://gdc.nci.nih.gov/access-data/obtaining-access-controlled-data).

Finally, the data provided by [GDC data
portal](https://gdc-portal.nci.nih.gov/) can be
accessed using Bioconductor package
[TCGAbiolinks](http://bioconductor.org/packages/TCGAbiolinks/), while
the data provided by Firehose can be accessed by Bioconductor package
[RTCGAToolbox](http://bioconductor.org/packages/RTCGAToolbox/).

The next steps describe how one could use
[TCGAbiolinks](http://bioconductor.org/packages/TCGAbiolinks/) &
[RTCGAToolbox](http://bioconductor.org/packages/RTCGAToolbox/) to
download clinical, genomics, transcriptomics, epigenomics data, as well
as subtype information and GISTIC results (i.e., identified genes
targeted by somatic copy-number alterations (SCNAs) that drive cancer
growth). All the data used in this workflow has as reference the Genome
Reference Consortium human genome (build 37 - hg19).

### Downloading data from TCGA data portal

The Bioconductor package
[TCGAbiolinks](http://bioconductor.org/packages/TCGAbiolinks/)
[@TCGAbiolinks] has three main functions *GDCquery*, *GDCdownload* and
*GDCprepare* that should sequentially be used to respectively search,
download and load the data as an R object.

*GDCquery* uses [GDC
API](https://gdc.cancer.gov/developers/gdc-application-programming-interface-api)
to search the data for a given project and data category and filters the
results by samples, sample type, file type and others features if
requested by the user. This function returns an object with a summary
table with the results found (samples, files and other useful
information) and the arguments used in the query. The most important
*GDCquery* arguments are *project* which receives a GDC project
(TCGA-USC, TCGA-LGG, TARGET-AML, etc), *data.category* which receives a
data category (Transcriptome Profiling, Copy Number Variation, DNA
methylation, Gene expression, etc), *data.type* which receives a data
type (Gene expression quantification, Isoform Expression Quantification,
miRNA Expression Quantification, Copy Number Segment, Masked Copy Number
Segment, etc), *workflow.type*, which receives a GDC workflow type
(STAR - Counts), and *platform*,
which receives a platform for the searches in the legacy database
(HumanMethylation27, Genome\_Wide\_SNP\_6, IlluminaHiSeq\_RNASeqV2,
etc). A complete list of possible entries for arguments can be found in
the [TCGAbiolinks
vignette](https://bioconductor.org/packages/release/bioc/vignettes/TCGAbiolinks/inst/doc/query.html#Harmonized_data_options).
Listing 1 shows an example of this function.

After the search step, the user will be able to download the data using
the *GDCdownload* function which can use either the GDC API to download
the samples, or the [gdc client
tools](https://gdc.cancer.gov/access-data/gdc-data-transfer-tool). The
downloaded data will be saved in a directory with the project name and a
sub-folder with the data.category, for example
"TCGA-GBM/DNA\_methylation".

Finally, *GDCprepare* transforms the downloaded data into a
[summarizedExperiment](http://bioconductor.org/packages/SummarizedExperiment/)
object [@huber2015orchestrating] or a data frame. If
*SummarizedExperiment* is set to TRUE, TCGAbiolinks will add to the
object sub-type information, which was defined by The Cancer Genome
Atlas (TCGA) Research Network reports (the full list of papers can be
seen in [TCGAquery\_subtype
section](http://bioconductor.org/packages/devel/bioc/vignettes/TCGAbiolinks/inst/doc/tcgaBiolinks.html#tcgaquery_subtype-working-with-molecular-subtypes-data.)
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](https://bioconductor.org/packages/release/bioc/vignettes/TCGAbiolinks/inst/doc/index.html).


```{R, eval=FALSE, include=TRUE}
library(TCGAbiolinks)
query_met_gbm <- GDCquery(
  project = "TCGA-GBM", 
  data.category = "DNA Methylation",
  data.type = "Methylation Beta Value",
  platform = "Illumina Human Methylation 450", 
  barcode = c("TCGA-76-4926-01B-01D-1481-05", "TCGA-28-5211-01C-11D-1844-05")
)
GDCdownload(query_met_gbm)

met_gbm_450k <- GDCprepare(
  query = query_met_gbm,
  summarizedExperiment = TRUE
)

query_met_lgg <- GDCquery(
  project = "TCGA-LGG", 
  data.category = "DNA Methylation",
  data.type = "Methylation Beta Value",
  platform = "Illumina Human Methylation 450",
  barcode = c("TCGA-HT-7879-01A-11D-2399-05", "TCGA-HT-8113-01A-11D-2399-05")
)
GDCdownload(query_met_lgg)
met_lgg_450k <- GDCprepare(
  query = query_met_lgg,
  summarizedExperiment = TRUE
)

met_lgg_450k$days_to_death <- NA
met_lgg_450k$year_of_death <- NA
met_gbm_lgg <- SummarizedExperiment::cbind(
  met_lgg_450k, 
  met_gbm_450k
)


# A total of 2.27 GB
query_exp_lgg <- GDCquery(
  project = "TCGA-LGG",
  data.category = "Transcriptome Profiling",
  data.type = "Gene Expression Quantification", 
  workflow.type = "STAR - Counts"
)

GDCdownload(query_exp_lgg)
exp_lgg <- GDCprepare(
  query = query_exp_lgg
)

query_exp_gbm <- GDCquery(
  project = "TCGA-GBM",
  data.category = "Transcriptome Profiling",
  data.type = "Gene Expression Quantification", 
  workflow.type = "STAR - Counts"
)
GDCdownload(query_exp_gbm)
exp_gbm <- GDCprepare(
  query = query_exp_gbm
)

# The following clinical data is not available in GBM
missing_cols <- setdiff(colnames(colData(exp_lgg)),colnames(colData(exp_gbm)))
for(i in missing_cols){
  exp_lgg[[i]] <- NULL
}

exp_gbm_lgg <- SummarizedExperiment::cbind(
  exp_lgg, 
  exp_gbm
)
```

```{R, eval=FALSE, include=TRUE}
#-----------------------------------------------------------------------------
#                   Data.category: Copy number variation aligned to hg38
#-----------------------------------------------------------------------------
query <- GDCquery(
  project = "TCGA-ACC",
  data.category = "Copy Number Variation",
  data.type = "Copy Number Segment",
  barcode = c( "TCGA-OR-A5KU-01A-11D-A29H-01", "TCGA-OR-A5JK-01A-11D-A29H-01")
)
GDCdownload(query)
data <- GDCprepare(query)

query <- GDCquery(
  project = "TCGA-ACC",
  data.category = "Copy Number Variation",
  data.type = "Masked Copy Number Segment",
  sample.type = c("Primary Tumor")
) # see the barcodes with getResults(query)$cases
GDCdownload(query)
data <- GDCprepare(query)
```

If a [SummarizedExperiment object](http://www.nature.com/nmeth/journal/v12/n2/fig_tab/nmeth.3252_F2.html) 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)
[@huber2015orchestrating; @SummarizedExperiment]. An example is shown in
listing below.


```{R, eval=TRUE, include=TRUE}
library(SummarizedExperiment)

# Load object from TCGAWorkflowData package
# This object will be created in subsequent sections for enhanced clarity and understanding.
data(TCGA_GBM_Transcriptome_20_samples) 

# get gene expression matrix
data <- assay(exp_gbm)
datatable(
  data = data[1:10,], 
  options = list(scrollX = TRUE, keys = TRUE, pageLength = 5), 
  rownames = TRUE
)

# get genes information
genes.info <- rowRanges(exp_gbm)
genes.info

# get sample information
sample.info <- colData(exp_gbm)
datatable(
  data = as.data.frame(sample.info), 
  options = list(scrollX = TRUE, keys = TRUE, pageLength = 5), 
  rownames = FALSE
)

```

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 below. This function has two
arguments *project* ("TCGA-GBM","TARGET-AML",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),
admin (batch number, project code, disease code, Biospecimen Core
Resource).

```{R, eval=TRUE, include=TRUE}
# get indexed clinical patient data for GBM samples
gbm_clin <- GDCquery_clinic(
  project = "TCGA-GBM", 
  type = "Clinical"
)

# get indexed clinical patient data for LGG samples
lgg_clin <- GDCquery_clinic(
  project = "TCGA-LGG", 
  type = "Clinical"
)

# Bind the results, as the columns might not be the same,
# we will will plyr rbind.fill, to have all columns from both files
clinical <- plyr::rbind.fill(
  gbm_clin,
  lgg_clin
)
```

```{r  echo = TRUE, message = FALSE, warning = FALSE}
datatable(
  clinical[1:10,], 
  options = list(scrollX = TRUE, keys = TRUE), 
  rownames = FALSE
)
```


```{r results = 'hide', echo=TRUE, message=FALSE, warning=FALSE}
# Fetch clinical data directly from the clinical XML files.
# if barcode is not set, it will consider all samples.
# We only set it to make the example faster
query <- GDCquery(
  project = "TCGA-GBM",
  data.format = "bcr xml",
  data.category = "Clinical",
  barcode = c("TCGA-08-0516","TCGA-02-0317")
) 
GDCdownload(query)
clinical <- GDCprepare_clinic(
  query = query, 
  clinical.info = "patient"
)
```

```{r  echo = TRUE, message = FALSE, warning = FALSE}
datatable(
  data = clinical, 
  options = list(scrollX = TRUE, keys = TRUE), 
  rownames = FALSE
)
```

```{r results = 'hide', echo=TRUE, message=FALSE, warning=FALSE}
clinical_drug <- GDCprepare_clinic(
  query = query, 
  clinical.info = "drug"
)
```
```{r  echo = TRUE, message = FALSE, warning = FALSE}
clinical_drug |>
  datatable(
    options = list(scrollX = TRUE, keys = TRUE), 
    rownames = FALSE
  )
```

```{r results = 'hide', echo=TRUE, message=FALSE, warning=FALSE}
clinical_radiation <- GDCprepare_clinic(
  query = query, 
  clinical.info = "radiation"
)
```

```{r  echo = TRUE, message = FALSE, warning = FALSE}
clinical_radiation |> 
  datatable(
    options = list(scrollX = TRUE, keys = TRUE), 
    rownames = FALSE
  )
```

```{r results = 'hide', echo=TRUE, message=FALSE, warning=FALSE}
clinical_admin <- GDCprepare_clinic(
  query = query, 
  clinical.info = "admin"
)
```

```{r  echo = TRUE, message = FALSE, warning = FALSE}
clinical_admin |>
  datatable(
    options = list(scrollX = TRUE, keys = TRUE), 
    rownames = FALSE
  )
```

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 code below shows how to download Somatic MAFs data using
[TCGAbiolinks](http://bioconductor.org/packages/TCGAbiolinks/).

```{R, eval=FALSE, include=TRUE}
query <- GDCquery(
  project = c("TCGA-LGG","TCGA-GBM"), 
  data.category = "Simple Nucleotide Variation", 
  access = "open",
  data.type = "Masked Somatic Mutation", 
  workflow.type = "Aliquot Ensemble Somatic Variant Merging and Masking"
)
GDCdownload(query)
maf <- GDCprepare(query)
```

```{r  echo = TRUE, message = FALSE, warning = FALSE}
data(maf_lgg_gbm)
maf[1:10,] |>
  datatable(
    options = list(scrollX = TRUE, keys = TRUE), 
    rownames = FALSE
  )
```

Finally, the Cancer Genome Atlas (TCGA) Research Network has reported
integrated genome-wide studies of various diseases (ACC
[@zheng2016comprehensive], BRCA [@cancer2012comprehensive_brca], COAD
[@cancer2012comprehensive_colon], GBM [@Cell], HNSC
[@cancer2015comprehensive], KICH [@davis2014somatic], KIRC
[@cancer2013comprehensive], KIRP [@cancer2016comprehensive], LGG
[@Cell], LUAD [@cancer2014comprehensive_lung], LUSC
[@cancer2012comprehensive_lusc], PRAD [@cancer2015molecular_prad], READ
[@cancer2012comprehensive_colon], SKCM [@cancer2015genomic_skcm], STAD
[@cancer2014comprehensive_stad], THCA [@cancer2014integrated_thca] and
UCEC [@cancer2014comprehensive_gastric]) which classified them in
different subtypes. This classification can be retrieved using the
*TCGAquery\_subtype* function or by accessing the samples information in
the SummarizedExperiment object that created by the *GDCprepare*
function, which automatically incorporates it into the object.

```{R, eval=TRUE, include=TRUE}
gbm_subtypes <- TCGAquery_subtype(
  tumor = "gbm"
)
```

```{r  echo = TRUE, message = FALSE, warning = FALSE}
datatable(
  gbm_subtypes[1:10,], 
  options = list(scrollX = TRUE, keys = TRUE), 
  rownames = FALSE
)
```


### Downloading data from Broad TCGA GDAC

The Bioconductor package
[RTCGAToolbox](http://bioconductor.org/packages/RTCGAToolbox/)
[@samur2014rtcgatoolbox] provides access to Firehose Level 3 and 4 data
through the function *getFirehoseData*. The following arguments allow
users to select the version and tumor type of interest:

-   dataset - Tumor to download. A complete list of possibilities can be
view 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.

These arguments can be used to select the data type to download:
RNAseq\_Gene, Clinic, miRNASeq\_Gene, ccRNAseq2\_Gene\_Norm, CNA\_SNP,
CNV\_SNP, CNA\_Seq, CNA\_CGH, Methylation, Mutation, mRNA\_Array ,
miRNA\_Array, and RPPA.

By default,
[RTCGAToolbox](http://bioconductor.org/packages/RTCGAToolbox/) allows
users to download up to 500 MB worth of data. To increase the size of
the download, users are encouraged to use *fileSizeLimit* argument. An
example is found in listing below. The *getData* function
allows users to access the downloaded data.

```{R, eval=FALSE, include=TRUE}
library(RTCGAToolbox)
```

```{R, eval=FALSE, include=TRUE}
# Get the last run dates
lastRunDate <- getFirehoseRunningDates()[1]

# get DNA methylation data, RNAseq2 and clinical data for GBM
gbm_data <- getFirehoseData(
  dataset = "GBM",
  runDate = lastRunDate, 
  gistic2Date = getFirehoseAnalyzeDates(1),
  Methylation = FALSE,  
  clinical = TRUE,
  RNASeq2GeneNorm  = FALSE, 
  Mutation = TRUE,
  fileSizeLimit = 10000
)

gbm_mut <- getData(gbm_data,"Mutation")
gbm_clin <- getData(gbm_data,"clinical")
```

Finally, [RTCGAToolbox](http://bioconductor.org/packages/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) [@mermel2011gistic2].

```{R, eval=FALSE, message=FALSE, warning=FALSE, include=TRUE}
# Download GISTIC results
lastanalyzedate <- getFirehoseAnalyzeDates(1)
gistic <- getFirehoseData(
  dataset = "GBM",
  GISTIC = TRUE, 
  gistic2Date = lastanalyzedate
)

# get GISTIC results
gistic_allbygene <- getData(
  object = gistic, 
  type = "GISTIC", 
  platform = "AllByGene"
)
gistic_thresholedbygene <- getData(
  object = gistic, 
  type = "GISTIC", 
  platform = "ThresholdedByGene"
)
```

```{R, eval=TRUE, message=FALSE, warning=FALSE, include=TRUE}
data(gbm_gistic)
gistic_allbygene %>% head() %>% gt::gt()
gistic_thresholedbygene %>% head() %>% gt::gt()
```

## 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 of 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 the cancer genome. We analyzed GBM 
segmented CNV from SNP array (Affymetrix Genome-Wide Human SNP Array
6.0).

### Visualizing multiple genomic alteration events

In order to visualize multiple genomic alteration events, we recommend
using `maftools` plot which is provided by Bioconductor package
[maftools](http://bioconductor.org/packages/maftools/)
[@mayakonda2016maftools]. The listing below shows how to download
mutation data using *GDCquery\_maf* (line 4) and prepare it to use with maftools.

The function `read.maf` is used to prepare the MAF data to be used with maftools.
We also added clinical information that will be used in survival plots.

```{r results='hide', echo=TRUE, message=FALSE, warning=FALSE}
library(maftools)
# recovering data from TCGAWorkflowData package.
data(maf_lgg_gbm)

# To prepare for maftools we will also include clinical data
# For a mutant vs WT survival analysis 
# get indexed clinical patient data for GBM samples
gbm_clin <- GDCquery_clinic(project = "TCGA-GBM", type = "Clinical")
# get indexed clinical patient data for LGG samples
lgg_clin <- GDCquery_clinic(project = "TCGA-LGG", type = "Clinical")
# Bind the results, as the columns might not be the same,
# we will will plyr rbind.fill, to have all columns from both files
clinical <- plyr::rbind.fill(gbm_clin,lgg_clin)
colnames(clinical)[grep("submitter_id",colnames(clinical))] <- "Tumor_Sample_Barcode"

# we need to create a binary variable 1 is dead 0 is not dead
plyr::count(clinical$vital_status)
clinical$Overall_Survival_Status <- 1 # dead
clinical$Overall_Survival_Status[which(clinical$vital_status != "Dead")] <- 0

# If patient is not dead we don't have days_to_death (NA)
# we will set it as the last day we know the patient is still alive
clinical$time <- clinical$days_to_death
clinical$time[is.na(clinical$days_to_death)] <- clinical$days_to_last_follow_up[is.na(clinical$days_to_death)]

# Create object to use in maftools
maf <- read.maf(
  maf = maf, 
  clinicalData = clinical, 
  isTCGA = TRUE
)
```

We can plot a MAF summary.
```{r echo=TRUE, message=FALSE, warning=FALSE,fig.width=10}
plotmafSummary(
  maf = maf,
  rmOutlier = TRUE,
  addStat = 'median',
  dashboard = TRUE
)
```

We can draw oncoplot with the top 20 most mutated genes and add metadata information such as molecular subtypes information.

```{r echo=TRUE, message=FALSE, warning=FALSE,fig.height=10,fig.width=15,eval=FALSE}
oncoplot(
  maf = maf,
  top = 20,
  legendFontSize = 8,
  clinicalFeatures = c("tissue_or_organ_of_origin")
)
```

We can also perform survival analysis by grouping samples from MAF based on mutation status of given gene(s).

```{r echo=TRUE, message=FALSE, warning=FALSE}
plot <- mafSurvival(
  maf = maf,
  genes = "TP53",
  time = 'time',
  Status = 'Overall_Survival_Status',
  isTCGA = TRUE
)
```


## Transcriptomic analysis

### Pre-Processing Data

The
[LGG](https://gdc-portal.nci.nih.gov/legacy-archive/search/f?filters=%7B%22op%22:%22and%22,%22content%22:%5B%7B%22op%22:%22in%22,%22content%22:%7B%22field%22:%22cases.project.project_id%22,%22value%22:%5B%22TCGA-LGG%22%5D%7D%7D,%7B%22op%22:%22in%22,%22content%22:%7B%22field%22:%22files.data_category%22,%22value%22:%5B%22Gene%20expression%22%5D%7D%7D,%7B%22op%22:%22in%22,%22content%22:%7B%22field%22:%22files.data_type%22,%22value%22:%5B%22Gene%20expression%20quantification%22%5D%7D%7D,%7B%22op%22:%22in%22,%22content%22:%7B%22field%22:%22files.platform%22,%22value%22:%5B%22Illumina%20HiSeq%22%5D%7D%7D%5D%7D)
and
[GBM](https://gdc-portal.nci.nih.gov/legacy-archive/search/f?filters=%7B%22op%22:%22and%22,%22content%22:%5B%7B%22op%22:%22in%22,%22content%22:%7B%22field%22:%22cases.project.project_id%22,%22value%22:%5B%22TCGA-GBM%22%5D%7D%7D,%7B%22op%22:%22in%22,%22content%22:%7B%22field%22:%22files.data_category%22,%22value%22:%5B%22Gene%20expression%22%5D%7D%7D,%7B%22op%22:%22in%22,%22content%22:%7B%22field%22:%22files.data_type%22,%22value%22:%5B%22Gene%20expression%20quantification%22%5D%7D%7D,%7B%22op%22:%22in%22,%22content%22:%7B%22field%22:%22files.platform%22,%22value%22:%5B%22Illumina%20HiSeq%22%5D%7D%7D%5D%7D)
data used for following transcriptomic analysis were downloaded using
TCGAbiolinks. We downloaded only primary 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 filename including both the cancer name and the
name of the platform used for gene expression data.



```{R, eval=FALSE, include=TRUE}
query_exp_lgg <- GDCquery(
  project = "TCGA-LGG",
  data.category = "Transcriptome Profiling",
  data.type = "Gene Expression Quantification", 
  workflow.type = "STAR - Counts"
)
# Get only first 20 samples to make example faster
query_exp_lgg$results[[1]] <- query_exp_lgg$results[[1]][1:20,]
GDCdownload(query_exp_lgg)
exp_lgg <- GDCprepare(
  query = query_exp_lgg
)

query_exp_gbm <- GDCquery(
  project = "TCGA-GBM",
  data.category = "Transcriptome Profiling",
  data.type = "Gene Expression Quantification", 
  workflow.type = "STAR - Counts"
)
# Get only first 20 samples to make example faster
query_exp_gbm$results[[1]] <- query_exp_gbm$results[[1]][1:20,]
GDCdownload(query_exp_gbm)
exp_gbm <- GDCprepare(
  query = query_exp_gbm
)
```

To pre-process the data, first, we searched for possible outliers using
the *TCGAanalyze\_Preprocessing* function, which performs an Array Array
Intensity correlation AAIC.
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](http://bioconductor.org/packages/EDASeq/) package, we
normalized mRNA transcripts.

This function implements Within-lane normalization procedures to adjust
for GC-content effect (or other gene-level effects) on read counts:
loess robust local regression, global-scaling, and full-quantile
normalization [@risso2011gc] and between-lane normalization procedures
to adjust for distributional differences between lanes (e.g., sequencing
depth): global-scaling and full-quantile normalization
[@bullard2010evaluation].

```{r results='asis', echo=TRUE, message=FALSE,warning=FALSE}
data("TCGA_LGG_Transcriptome_20_samples")
data("TCGA_GBM_Transcriptome_20_samples")

exp_lgg_preprocessed <- TCGAanalyze_Preprocessing(
  object = exp_lgg,
  cor.cut = 0.6,    
  datatype = "unstranded",
  filename = "LGG_IlluminaHiSeq_RNASeqV2.png"
)

exp_gbm_preprocessed <- TCGAanalyze_Preprocessing(
  object = exp_gbm,
  cor.cut = 0.6, 
  datatype = "unstranded",
  filename = "GBM_IlluminaHiSeq_RNASeqV2.png"
)
exp_preprocessed <- cbind(
  exp_lgg_preprocessed, 
  exp_gbm_preprocessed
)

exp_normalized <- TCGAanalyze_Normalization(
  tabDF = cbind(exp_lgg_preprocessed, exp_gbm_preprocessed),
  geneInfo = TCGAbiolinks::geneInfoHT,
  method = "gcContent"
) # 60513   40

exp_filtered <- TCGAanalyze_Filtering(
  tabDF = exp_normalized,
  method = "quantile",
  qnt.cut =  0.25
)  # 44630   40

exp_filtered_lgg <- exp_filtered[
  ,substr(colnames(exp_filtered),1,12) %in% lgg_clin$bcr_patient_barcode
]

exp_filtered_gbm <-   exp_filtered[
  ,substr(colnames(exp_filtered),1,12) %in% gbm_clin$bcr_patient_barcode
]

diff_expressed_genes <- TCGAanalyze_DEA(
  mat1 = exp_filtered_lgg,
  mat2 = exp_filtered_gbm,
  Cond1type = "LGG",
  Cond2type = "GBM",
  fdr.cut = 0.01 ,
  logFC.cut = 1,
  method = "glmLRT"
)
```

```{r results='asis', echo=TRUE, message=FALSE,warning=FALSE}
# Number of differentially expressed genes (DEG)
nrow(diff_expressed_genes)
```


### EA: enrichment analysis
In order to understand the underlying biological process of DEGs we
performed an enrichment analysis using *TCGAanalyze\_EA\_complete*
function.

```{r results='hide', echo=TRUE, message=FALSE,warning=FALSE,fig.height=10, fig.width=10}
#-------------------  4.2 EA: enrichment analysis             --------------------
ansEA <- TCGAanalyze_EAcomplete(
  TFname = "DEA genes LGG Vs GBM", 
  RegulonList = diff_expressed_genes$gene_name
)

TCGAvisualize_EAbarplot(
  tf = rownames(ansEA$ResBP),
  filename = NULL,
  GOBPTab = ansEA$ResBP,
  nRGTab = diff_expressed_genes$gene_name,
  nBar = 20
)
```

```{r results='hide', echo=TRUE, message=FALSE,warning=FALSE,fig.height=10, fig.width=10}
TCGAvisualize_EAbarplot(
  tf = rownames(ansEA$ResBP),
  filename = NULL,
  GOCCTab = ansEA$ResCC,
  nRGTab = diff_expressed_genes$gene_name,
  nBar = 20
)
```

```{r results='hide', echo=TRUE, message=FALSE,warning=FALSE,fig.height=10, fig.width=10}
TCGAvisualize_EAbarplot(
  tf = rownames(ansEA$ResBP),
  filename = NULL,
  GOMFTab = ansEA$ResMF,
  nRGTab = diff_expressed_genes$gene_name,
  nBar = 20
)
```

```{r results='hide', echo=TRUE, message=FALSE,warning=FALSE,fig.height=12, fig.width=15}
TCGAvisualize_EAbarplot(
  tf = rownames(ansEA$ResBP),
  filename = NULL,
  PathTab = ansEA$ResPat,
  nRGTab = rownames(diff_expressed_genes),
  nBar = 20
)
```

The plot shows canonical pathways significantly overrepresented (enriched)
by the DEGs (differentially expressed genes)  with the number of genes for the main categories of
three ontologies (GO:biological process, GO:cellular component, and GO:molecular function, respectively).
The most statistically significant canonical pathways identified in DEGs list are listed according
to their p-value corrected FDR (-Log) (colored bars) and the ratio of list genes found in each
pathway over the total number of genes in that pathway (ratio, red line).]

*TCGAanalyze\_EAbarplot* outputs a bar chart as shown in figure
\ref{eabarplot} 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  shows canonical pathways significantly
overrepresented (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](http://bioconductor.org/packages/pathview/)
[@luo2013pathview] can be used. Listing below shows an
example how to use it. It can receive, for example, a named vector of
genes with their expression level, the pathway.id which can be found in
[KEGG database](http://www.genome.jp/kegg/pathway.html), the species
(’hsa’ for Homo sapiens) and the limits for the gene expression.

```{r results='asis', echo=TRUE, message=FALSE,warning=FALSE}
library(SummarizedExperiment)

# DEGs TopTable
dataDEGsFiltLevel <- TCGAanalyze_LevelTab(
  FC_FDR_table_mRNA = diff_expressed_genes,
  typeCond1 = "LGG",
  typeCond2 = "GBM",
  TableCond1 = exp_filtered[,colnames(exp_filtered_lgg)],
  TableCond2 = exp_filtered[,colnames(exp_filtered_gbm)]
)

dataDEGsFiltLevel$GeneID <- 0

library(clusterProfiler)
# Converting Gene symbol to geneID
eg = as.data.frame(
  bitr(
    dataDEGsFiltLevel$mRNA,
    fromType = "ENSEMBL",
    toType = c("ENTREZID","SYMBOL"),
    OrgDb = "org.Hs.eg.db"
  )
)
eg <- eg[!duplicated(eg$SYMBOL),]
eg <- eg[order(eg$SYMBOL,decreasing=FALSE),]

dataDEGsFiltLevel <- dataDEGsFiltLevel[dataDEGsFiltLevel$mRNA %in% eg$ENSEMBL,]
dataDEGsFiltLevel <- dataDEGsFiltLevel[eg$ENSEMBL,]
rownames(dataDEGsFiltLevel) <- eg$SYMBOL

all(eg$SYMBOL == rownames(dataDEGsFiltLevel))
dataDEGsFiltLevel$GeneID <- eg$ENTREZID

dataDEGsFiltLevel_sub <- subset(dataDEGsFiltLevel, select = c("GeneID", "logFC"))
genelistDEGs <- as.numeric(dataDEGsFiltLevel_sub$logFC)
names(genelistDEGs) <- dataDEGsFiltLevel_sub$GeneID
library(pathview)
# pathway.id: hsa05214 is the glioma pathway
# limit: sets the limit for gene expression legend and color
hsa05214 <- pathview::pathview(
  gene.data  = genelistDEGs,
  pathway.id = "hsa05214",
  species    = "hsa",
  limit = list(gene = as.integer(max(abs(genelistDEGs))))
)

```

The red genes are up-regulated and the green genes are down-regulated in
the LGG samples compared to GBM.

![\label{gliomapathview} Pathways enrichment analysis: glioma pathway. Red defines genes that are up-regulated and green defines genes that are down-regulated.](hsa05214.pathview.png)

### 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 [@margolin2006aracne], CLR [@faith2007large], MRNET
[@meyer2007information] and C3NET [@altay2010inferring]. These methods
are based on mutual inference and use different heuristics to infer the
edges of the network. These methods have been made available via
Bioconductor/CRAN packages
([MINET](http://bioconductor.org/packages/minet/) [@Meyer2008] and
[c3net](https://cran.r-project.org/web/packages/c3net/index.html),
[@altay2010inferring] 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 [@Stark01012006] or GeneMANIA [@montojo2010genemania]. 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 below.

validated   not validated/non-existing
-------------- ----------- ----------------------------
inferred           TP                   FP
not inferred       FN                   TN

: Confusion matrix, comparing inferred network to network of validated
interactions.<span data-label="table::confusion_matrix"></span>

Different quality measures can then be computed such as the false
positive rate $$fpr=\frac{FP}{FP+TN},$$ the true positive rate (also
called recall) $$tpr=\frac{TP}{TP+FN}$$ and the precision
$$p=\frac{TP}{TP+FP}.$$ 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](http://thebiogrid.org/).

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 one gene set for the LGG data.

```{R, eval=FALSE, include=TRUE}
### read biogrid info (available in TCGAWorkflowData as "biogrid")
### Check last version in https://thebiogrid.org/download.php 
file <- "https://downloads.thebiogrid.org/Download/BioGRID/Latest-Release/BIOGRID-ALL-LATEST.tab2.zip"
if(!file.exists(gsub("zip","txt",basename(file)))){
  downloader::download(file,basename(file))
  unzip(basename(file),junkpaths =TRUE)
}

tmp.biogrid <- vroom::vroom(
  dir(pattern = "BIOGRID-ALL.*\\.txt")
)
```

```{R, eval=TRUE, include=TRUE}
### plot details (colors & symbols)
mycols <- c('#e41a1c','#377eb8','#4daf4a','#984ea3','#ff7f00','#ffff33','#a65628')

### load network inference libraries
library(minet)
library(c3net)

### deferentially identified genes using TCGAbiolinks
# we will use only a subset (first 50 genes) of it to make the example faster
names.genes.de <- rownames(diff_expressed_genes)[1:30]

data(biogrid)
net.biogrid.de <- getAdjacencyBiogrid(tmp.biogrid, names.genes.de)

mydata <- exp_filtered_lgg[names.genes.de, ]

### infer networks
t.mydata <- t(mydata)
net.aracne <- minet(t.mydata, method = "aracne")
net.mrnet <- minet(t.mydata)
net.clr <- minet(t.mydata, method = "clr")
net.c3net <- c3net(mydata)

### validate compared to biogrid network
tmp.val <- list(
  validate(net.aracne, net.biogrid.de), 
  validate(net.mrnet, net.biogrid.de),
  validate(net.clr, net.biogrid.de), 
  validate(net.c3net, net.biogrid.de)
)

### plot roc and compute auc for the different networks
dev1 <- show.roc(tmp.val[[1]],cex=0.3,col=mycols[1],type="l")
res.auc <- auc.roc(tmp.val[[1]])
for(count in 2:length(tmp.val)){
  show.roc(tmp.val[[count]],device=dev1,cex=0.3,col=mycols[count],type="l")
  res.auc <- c(res.auc, auc.roc(tmp.val[[count]]))
}

legend(
  "bottomright", 
  legend = paste(c("aracne","mrnet","clr","c3net"), signif(res.auc,4), sep=": "),
  col = mycols[1:length(tmp.val)],
  lty = 1, 
  bty = "n" 
)
# Please, uncomment this line to produce the pdf files.
# dev.copy2pdf(width=8,height=8,device = dev1, file = paste0("roc_biogrid_",cancertype,".pdf"))
```


![ROC with corresponding AUC for inferred GBM networks compared to
BioGrid interactions\label{fig::roc_GBM}](roc_biogrid_GBM.png)

In Figure above, the obtained ROC curve and the
corresponding area under the curve (AUC) are presented. It can be observed
that CLR and MRNET perform best when comparing the inferred network with
known interactions from the BioGrid database.


## 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
[@phillips2008role].

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 [@deaton2011cpg].

Thus, the investigation of the DNA methylation is crucial to
understanding regulatory gene networks in cancer as the DNA methylation
represses transcription [@robertson2005dna]. 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](http://bioconductor.org/packages/TCGAbiolinks/)
[@TCGAbiolinks]. 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.
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](https://support.illumina.com/content/dam/illumina-support/documents/myillumina/90666eaa-0c66-48b4-8199-3be99b2b3ef9/genomestudio_methylation_v1.8_user_guide_11319130_b.pdf).
Therefore, the rs probes were removed. Also, probes in chromosomes X, Y were removed to eliminate
potential artifacts originating from the presence of a different
proportion of males and females [@marabita2013evaluation]. The last
pre-processing steps were to remove probes with at least one NA value.

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

```{R, eval=FALSE, include=TRUE}
#----------------------------
# Obtaining DNA methylation
#----------------------------
# Samples
lgg.samples <- matchedMetExp("TCGA-LGG", n = 10)
gbm.samples <- matchedMetExp("TCGA-GBM", n = 10)
samples <- c(lgg.samples,gbm.samples)

#-----------------------------------
# 1 - Methylation
# ----------------------------------
# For DNA methylation it is quicker in this case to download the tar.gz file
# and get the samples we want instead of downloading files by files
query <- GDCquery(
  project = c("TCGA-LGG","TCGA-GBM"),
  data.category = "DNA Methylation",
  platform = "Illumina Human Methylation 450",
  data.type = "Methylation Beta Value",
  barcode = samples
)
GDCdownload(query)
met <- GDCprepare(
  query = query, 
  save = FALSE
)

# We will use only chr9 to make the example faster
met <- subset(met,subset = as.character(seqnames(met)) %in% c("chr9"))
# This data is avaliable in the package (object elmerExample)
```

```{r echo=TRUE, message=FALSE,warning=FALSE, fig.height=8,fig.width=8}
data(elmerExample)
#----------------------------
# Mean methylation
#----------------------------
# Plot a barplot for the groups in the disease column in the
# summarizedExperiment object

# remove probes with NA (similar to na.omit)
met <- met[rowSums(is.na(assay(met))) == 0,]

df <- data.frame(
  "Sample.mean" = colMeans(assay(met), na.rm = TRUE),
  "groups" = met$project_id
)

library(ggpubr)
ggpubr::ggboxplot(
  data = df,
  y = "Sample.mean",
  x = "groups",
  color = "groups",
  add = "jitter",
  ylab = "Mean DNA methylation (beta-values)",
  xlab = ""
) + stat_compare_means() 
```

The figure above 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.

### 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\_DMC* function
(see listing below). 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 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.15 and an adjusted p-value of less than $0.05$.

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. 

```{r results='hide', echo=TRUE, message=FALSE,warning=FALSE}
#------- Searching for differentially methylated CpG sites     ----------
dmc <- TCGAanalyze_DMC(
  data = met,
  groupCol = "project_id", # a column in the colData matrix
  group1 = "TCGA-GBM", # a type of the disease type column
  group2 = "TCGA-LGG", # a type of the disease column
  p.cut = 0.05,
  diffmean.cut = 0.15,
  save = FALSE,
  legend = "State",
  plot.filename = "LGG_GBM_metvolcano.png",
  cores = 1 # if set to 1 there will be a progress bar
)
```

The figure below shows the
volcano plot produced by listing below. This plot aids the
user in selecting relevant thresholds, as we search for candidate
biological DMRs.

![Volcano plot: searching for differentially methylated
CpG sites (x-axis:difference of mean DNA methylation, y-axis:
statistical
significance)](LGG_GBM_metvolcano.png)

To visualize the level of DNA methylation of these probes across all
samples, we use heatmaps that can be generated by the Bioconductor
package
[complexHeatmap](http://bioconductor.org/packages/ComplexHeatmap/)
[@ComplexHeatmap]. To create a heatmap using the
[complexHeatmap](http://bioconductor.org/packages/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 below shows
the code to produce the heatmap of a DNA methylation datum.


```{r results='hide', echo=TRUE, message=FALSE,warning=FALSE}
#--------------------------
# DNA Methylation heatmap
#-------------------------
library(ComplexHeatmap)
clinical <- plyr::rbind.fill(
  gbm_clin,
  lgg_clin
)

# get the probes that are Hypermethylated or Hypomethylated
# met is the same object of the section 'DNA methylation analysis'
status.col <- "status"
probes <- rownames(dmc)[grep("hypo|hyper",dmc$status,ignore.case = TRUE)]
sig.met <- met[probes,]


# top annotation, which samples are LGG and GBM
# We will add clinical data as annotation of the samples
# we will sort the clinical data to have the same order of the DNA methylation matrix
clinical.ordered <- clinical[match(substr(colnames(sig.met),1,12),clinical$bcr_patient_barcode),]

ta <- HeatmapAnnotation(
  df = clinical.ordered[, c("primary_diagnosis", "gender", "vital_status", "race")],
  col = list(
    disease = c("LGG" = "grey", "GBM" = "black"),
    gender = c("male" = "blue", "female" = "pink")
  )
)

# row annotation: add the status for LGG in relation to GBM
# For exmaple: status.gbm.lgg Hypomethyated means that the
# mean DNA methylation of probes for lgg are hypomethylated
# compared to GBM ones.
ra = rowAnnotation(
  df = dmc[probes, status.col],
  col = list(
    "status.TCGA.GBM.TCGA.LGG" =
      c(
        "Hypomethylated" = "orange",
        "Hypermethylated" = "darkgreen"
      )
  ),
  width = unit(1, "cm")
)

heatmap  <- Heatmap(
  matrix = assay(sig.met),
  name = "DNA methylation",
  col = matlab::jet.colors(200),
  show_row_names = FALSE,
  cluster_rows = TRUE,
  cluster_columns = FALSE,
  show_column_names = FALSE,
  bottom_annotation = ta,
  column_title = "DNA Methylation"
) 
# Save to pdf
png("heatmap.png",width = 600, height = 400)
draw(heatmap, annotation_legend_side =  "bottom")
dev.off()
```


![](heatmap.png)


### 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](http://bioconductor.org/packages/rGADEM/)
[@droit2015rgadem; @li2009gadem] 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.
An object will be returned which
contains all relevant information about your motif analysis (i.e.,
sequence consensus, PWM, chromosome, p-value, etc).

Using Bioconductor package
[motifStack](http://bioconductor.org/packages/motifStack/)
[@ou2013motifstack] it is possible to generate a graphic representation
of multiple motifs with different similarity scores.
```{r results='asis', echo=TRUE, message=FALSE,warning=FALSE}
library(rGADEM)
library(BSgenome.Hsapiens.UCSC.hg19)
library(motifStack)
library(SummarizedExperiment)
library(dplyr)

probes <- rowRanges(met)[rownames(dmc)[grep("hypo|hyper",dmc$status,ignore.case = TRUE)],]

# Get hypo/hyper methylated probes and make a 200bp window 
# surrounding each probe.
sequence <- GRanges(
  seqnames = as.character(seqnames(probes)),
  IRanges(
    start = ranges(probes) %>% as.data.frame() %>% dplyr::pull("start") - 100,
    end = ranges(probes) %>% as.data.frame() %>% dplyr::pull("end") + 100), 
  strand = "*"
)
#look for motifs
gadem <- GADEM(sequence, verbose = FALSE, genome = Hsapiens)

# How many motifs were found?
nMotifs(gadem)

# get the number of occurrences
nOccurrences(gadem)

# view all sequences consensus
consensus(gadem)

# Print motif
pwm <- getPWM(gadem)
pfm  <- new("pfm",mat = pwm[[1]],name = "Novel Site 1")
plotMotifLogo(pfm)

# Number of instances of motif 1?
length(gadem@motifList[[1]]@alignList)

```

## 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 throughput data
[@phillips2008role; @shi2014integrative; @rhodes2005integrative]. In
this section, we will introduce a Bioconductor package called
[ELMER](http://bioconductor.org/packages/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](http://www.encodeproject.org/) and
[Roadmap](http://www.roadmapepigenomics.org/).

### 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 the table below.


```{r table4, echo=FALSE, message=FALSE, warnings=FALSE, results='asis'}
tabl <- "
|                  Histone marks                 |                                                   Role                                                  |
|:----------------------------------------------:|:-------------------------------------------------------------------------------------------------------:|
| Histone H3 lysine 4 trimethylation (H3K4me3)   | Promoter regions [@heintzman2007distinct,@bernstein2005genomic]                                      |
| Histone H3 lysine 4 monomethylation (H3K4me1)  | Enhancer regions [@heintzman2007distinct]                                                           |
| Histone H3 lysine 36 trimethylation (H3K36me3) | Transcribed regions                                                                                     |
| Histone H3 lysine 27 trimethylation (H3K27me3) | Polycomb repression [@bonasio2010molecular]                                                         |
| Histone H3 lysine 9 trimethylation (H3K9me3)   | Heterochromatin regions  [@peters2003partitioning]                                                  |
| Histone H3 acetylated at lysine 27 (H3K27ac)   | Increase activation of genomic elements [@heintzman2009histone,@rada2011unique,@creyghton2010histone] |
| Histone H3 lysine 9 acetylation  (H3K9ac)      | Transcriptional activation [@nishida2006histone]                                                    |
"
cat(tabl) 
```

Besides, ChIP-seq data exists in the ROADMAP database and can be obtained through the
[AnnotationHub](http://bioconductor.org/packages/AnnotationHub/) package [@annotationHub]
or from [Roadmap web portal](http://egg2.wustl.edu/roadmap/webportal/processed_data.html).
The table below shows the description of all the roadmap files that are available through AnnotationHub.

```{r table5, echo=FALSE, message=FALSE, warnings=FALSE, results='asis'}
tabl <- "  
|                File                |                               Description                              |
|:----------------------------------:|:----------------------------------------------------------------------:|
| fc.signal.bigwig                   | Bigwig File containing  fold enrichment signal tracks                  |
| pval.signal.bigwig                 | Bigwig File containing -log10(p-value) signal tracks                   |
| hotspot.fdr0.01.broad.bed.gz       | Broad domains on enrichment for  DNase-seq for consolidated epigenomes |
| hotspot.broad.bed.gz               | Broad domains on enrichment for DNase-seq for consolidated epigenomes  |
| broadPeak.gz                       | Broad ChIP-seq peaks for consolidated  epigenomes                      |
| gappedPeak.gz                      | Gapped ChIP-seq peaks for consolidated   epigenomes                    |
| narrowPeak.gz                      | Narrow ChIP-seq peaks for consolidated epigenomes                      |
| hotspot.fdr0.01.peaks.bed.gz       | Narrow DNasePeaks for   consolidated epigenomes                        |
| hotspot.all.peaks.bed.gz           | Narrow DNasePeaks for  consolidated epigenomes                         |
| .macs2.narrowPeak.gz               | Narrow DNasePeaks for consolidated epigenomes                          |
| coreMarks_mnemonics.bed.gz        | 15 state chromatin segmentations                                       |
| mCRF_FractionalMethylation.bigwig | MeDIP/MRE(mCRF) fractional methylation calls                           |
| RRBS_FractionalMethylation.bigwig | RRBS fractional methylation calls                                      |
| WGBS_FractionalMethylation.bigwig | Whole genome bisulphite fractional methylation calls                   |
"
cat(tabl) 
```

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 [@ChIP-seqbioc]. The first step is shown in listing
below 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 the brain (keywords can be seen in the [summary
table](https://docs.google.com/spreadsheets/d/1yikGx4MsO9Ei36b64yOy9Vb6oPC5IBGlFbYEt-N6gOM/edit#gid=15))
and narrowPeak is selecting the type of file. The data downloaded is a
processed data from an integrative Analysis of 111 reference human
epigenomes [@kundaje2015integrative].

```{r results='hide', eval=TRUE, echo=TRUE, message=FALSE,warning=FALSE}
library(ChIPseeker)
library(pbapply)
library(ggplot2)
```

```{r results='hide', eval=FALSE, echo=TRUE, message=FALSE,warning=FALSE}
#------------------ Working with ChipSeq data ---------------
# Step 1: download histone marks for a brain and non-brain samples.
#------------------------------------------------------------
# loading annotation hub database
library(AnnotationHub)
ah = AnnotationHub()

# Searching for brain consolidated epigenomes in the roadmap database
bpChipEpi_brain <- query(ah , c("EpigenomeRoadMap", "narrowPeak", "chip", "consolidated","brain","E068"))
# Get chip-seq data
histone.marks <- pblapply(names(bpChipEpi_brain), function(x) {ah[[x]]})
names(histone.marks) <- names(bpChipEpi_brain) 
# OBS: histone.marks is available in TCGAWorkflowData package
```


The [Chipseeker](http://bioconductor.org/packages/ChIPseeker/) package
[@yu2015chipseeker] implements functions that use 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,
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 below
shows an example of code to plot the average profile. 

To help the user better understand the regions found in the DMR
analysis, we downloaded histone marks specific to brain tissue using
the AnnotationHub package that can access the Roadmap database. 
Next, the Chipseeker was used to visualize how
histone modifications are enriched onto hypomethylated and
hypermethylated regions, (listing below). The
enrichment heatmap and the average profile of peaks binding to those
regions.

```{r getTagMatrix, results='hide', echo=TRUE, message=FALSE,warning=FALSE}
data(histoneMarks)
# Create a GR object based on the hypo/hypermethylated probes.
probes <- keepStandardChromosomes(
  rowRanges(met)[rownames(dmc)[dmc$status %in% c("Hypermethylated in TCGA-GBM", "Hypomethylated in TCGA-GBM")],]
)
# Defining a window of 3kbp - 3kbp_probe_3kbp
# to make it work with ChIPseeker package version "1.31.3.900"
attributes(probes)$type <- "start_site"
attributes(probes)$downstream <- 3000
attributes(probes)$upstream <- 3000
probes <- GenomicRanges::resize(probes,6001,fix = "center") 

### Profile of ChIP peaks binding to TSS regions
# First of all, to calculate the profile of ChIP peaks binding to TSS regions, we should
# prepare the TSS regions, which are defined as the flanking sequence of the TSS sites.
# Then align the peaks that are mapping to these regions and generate the tagMatrix.
tagMatrixList <- pbapply::pblapply(histone.marks, function(x) {
  getTagMatrix(keepStandardChromosomes(x), windows = probes, weightCol = "score")
})
# change names retrieved with the following command: basename(bpChipEpi_brain$title)
names(tagMatrixList) <- c("H3K4me1","H3K4me3", "H3K9ac", "H3K9me3", "H3K27ac",  "H3K27me3", "H3K36me3")
```

To plot the enrichment heatmap use the function `tagHeatmap`

```{R tagHeatmap, eval=FALSE, include=TRUE}
tagHeatmap(tagMatrixList)
```

```{r results='asis', echo=FALSE, fig.width = 20, message=FALSE,warning=FALSE}
tagHeatmap(tagMatrixList)
```

To plot the  average profile of peaks binding to those region use `plotAvgProf`:
```{R plotAvgProf, eval=FALSE, include=TRUE}
p <- plotAvgProf(tagMatrixList, xlim = c(-3000,3000), xlab = "Genomic Region (5'->3', centered on CpG)")
# We are centreing in the CpG instead of the TSS. So we'll change the labels manually
p <- p + scale_x_continuous(
  breaks = c(-3000,-1500,0,1500,3000),
  labels = c(-3000,-1500,"CpG",1500,3000)
)

library(ggthemes)
p + theme_few() + scale_colour_few(name = "Histone marks") +  guides(colour = guide_legend(override.aes = list(size=4)))
```

```{r results='asis', echo=FALSE, message=FALSE,warning=FALSE}
p <- plotAvgProf(tagMatrixList, xlim = c(-3000,3000), xlab = "Genomic Region (5'->3', centered on CpG)")
# We are centreing in the CpG instead of the TSS. So we'll change the labels manually

if (requireNamespace("ggplot2", quietly = TRUE)) {
  p <- p + ggplot2::scale_x_continuous(
    breaks = c(-3000,-1500,0,1500,3000),
    labels = c(-3000,-1500,"CpG",1500,3000)
  )
}

if (requireNamespace("ggthemes", quietly = TRUE)){ 
  p + ggthemes::theme_few() + 
    ggthemes::scale_colour_few(name="Histone marks") + 
    ggplot2::guides(colour = guide_legend(override.aes = list(size=4)))
}
```

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 do not show an
enrichment nearby the position 0, and the average profile also does not
show a peak at position 0.

## Enhancer Linking by Methylation/Expression Relationship

Recently, many studies suggest that enhancers play a major role as
regulators of cell-specific phenotypes leading to alteration in
transcriptomes related to diseases
[@giorgio2015large; @groschel2014single; @sur2012mice; @yao2015demystifying].
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)](http://bioconductor.org/packages/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
the 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 datasets and custom DNA methylation and gene expression data
sets [@yao2015inferring,@ChedraouiSilva148726].

[ELMER](http://bioconductor.org/packages/ELMER/) analysis have 5 main
steps:

1. Identify distal probes on HM450K or EPIC.

2. Identification of distal probes with significant differential DNA methylation (i.e. DMCs) in tumor vs. normal samples.

3. Identification of putative target gene(s) for differentially methylated distal probes.

4. Identification of enriched motifs within a set of probes in significant probe-gene pairs.

5. Identification of master regulator Transcription Factors (TF) for each enriched motif.

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

### Preparing the data for ELMER package {#preparing-the-data-for-elmer-package .unnumbered}


The listing below shows how to use
[TCGAbiolinks](http://bioconductor.org/packages/TCGAbiolinks/)
[@TCGAbiolinks] to search, download and prepare the data for the
[ELMER](http://bioconductor.org/packages/ELMER/) package. Due to time
and memory constraints, we will use in this example only data from 10
LGG patients and 10 GBM patients that have both DNA methylation and gene
expression data. These samples are the same used in the previous steps.


To perform ELMER analyses, we need to populate a **MultiAssayExperiment** with a DNA methylation matrix or **SummarizedExperiment** object from HM450K or EPIC platform; a gene expression matrix or SummarizedExperiment object for the same samples; a matrix mapping DNA methylation samples to gene expression samples; and a matrix with sample metadata (i.e. clinical data, molecular subtype, etc.). If TCGA data are used, the last two matrices will be automatically generated.
If using non-TCGA data,  the matrix with sample metadata should be provided with at least a column with a patient identifier and another one identifying its group which will be used for analysis, if samples in the methylation and expression matrices are not ordered and with same names, a matrix mapping for each patient identifier their DNA methylation samples and their gene expression samples should be provided to the *createMAE* function.
Based on the genome of reference selected, metadata for the DNA methylation probes, such as genomic coordinates, will be added from http://zwdzwd.github.io/InfiniumAnnotation [@zhou2016comprehensive]; 
and metadata for gene expression and annotation is added from Ensembl database [@yates2015ensembl] using [biomaRt](http://bioconductor.org/packages/biomaRt/) [@durinck2009mapping]. 

```{R, eval=FALSE, include=TRUE}
#----------- 8.3 Identification of Regulatory Enhancers   -------
library(TCGAbiolinks)
# Samples: primary solid tumor w/ DNA methylation and gene expression
lgg.samples <- matchedMetExp("TCGA-LGG", n = 10)
gbm.samples <- matchedMetExp("TCGA-GBM", n = 10)
samples <- c(lgg.samples,gbm.samples)

#-----------------------------------
# 1 - Methylation
# ----------------------------------
query_met <- GDCquery(
  project = c("TCGA-LGG","TCGA-GBM"),
  data.category = "DNA Methylation",
  platform = "Illumina Human Methylation 450",
  data.type = "Methylation Beta Value",
  barcode = samples
)
GDCdownload(query_met)
met <- GDCprepare(query_met, save = FALSE)
met <- subset(met,subset = as.character(GenomicRanges::seqnames(met)) %in% c("chr9"))

#-----------------------------------
# 2 - Expression
# ----------------------------------
query_exp <- GDCquery(
  project = c("TCGA-LGG","TCGA-GBM"),
  data.category = "Gene expression",
  data.type = "Gene expression quantification",
  platform = "Illumina HiSeq", 
  file.type  = "results", 
  legacy = TRUE, 
  barcode =  samples
)
GDCdownload(query_exp)
exp <- GDCprepare(query_exp, save = FALSE)
save(exp, met, gbm.samples, lgg.samples, file = "elmer.example.rda", compress = "xz")
```

Probes from HumanMethylationEPIC (EPIC) array and Infinium HumanMethylation450 (HM450) array are removed from the analysis if they have either internal SNPs close to the $3'$ end of the probe; non-unique mapping to the bisulfite-converted genome; or off-target hybridization due to partial overlap with non-unique elements [@doi:10.1093/nar/gkw967]. This probe metadata information is
included in [ELMER.data](https://github.com/tiagochst/ELMER.data) package, populated from the source file at http://zwdzwd.github.io/InfiniumAnnotation [@doi:10.1093/nar/gkw967].
To limit ELMER to the analysis of distal elements, probes located in regions of $\pm2 kb$ around transcription start sites (TSSs) were removed. The distal elements are retrieved using the *get.feature.probe* function.


```{r elmer, message = FALSE,warning = FALSE}
library(ELMER)
library(MultiAssayExperiment)
library(GenomicRanges)
distal.probes <- get.feature.probe(
  genome = "hg19", 
  met.platform = "450K"
)

# Recover the data created in the last step
data(elmerExample)
rownames(exp) <- values(exp)$ensembl_gene_id
mae <- createMAE(
  exp = assay(exp), 
  met = met,
  save = TRUE,
  linearize.exp = TRUE,
  save.filename = "mae.rda",
  filter.probes = distal.probes,
  met.platform = "450K",
  genome = "hg19",
  TCGA = TRUE
)
mae
```

## ELMER analysis

After preparing the data into an MAE object, we executed the five
[ELMER](http://bioconductor.org/packages/ELMER/) steps for the hypo
direction, which means where are searching for distal enhancer probes hypomethylated in the GBM group compared to the LGG group. The
code is shown below. A description of how these distal enhancer probes
are identified is found in the [ELMER.data](https://bioconductor.org/packages/release/data/experiment/vignettes/ELMER.data/inst/doc/vignettes.pdf) vignette. As the number of samples is small we loosened some of the cut-off values.


```{r, warning=FALSE, results="hide"} 
cores <- 1 # you can use more cores if you want
group.col <- "project_id"
group1 <- "TCGA-GBM"
group2 <- "TCGA-LGG"

# Available directions are hypo and hyper, we will use only hypo
# due to speed constraint
direction <- "hypo"
dir.out <- paste0("elmer/",direction)
dir.create(dir.out, showWarnings = FALSE, recursive = TRUE)
```

### Supervised and unsupervised analysis mode


ELMER is designed to identify differences between two sets of samples within a given dataset. But, before we start our analysis, it is better to understand the two analysis mode available in ELMER 2.0 (supervised and unsupervised mode) which affects how the samples are selected for each comparison.

In original version  [@yao2015inferring], the first step - identification of differentially methylated CpG probes (DMCs) - was hard-coded to identify DMCs between non-cancer vs. cancer samples, and the subsequent step was unsupervised, identifying changes within any subset of tumors. In ELMER 2.0, we generalize these strategies so that they are applicable to any paired dataset, including disease vs. healthy tissue for any disease type. 

In the Unsupervised mode, as in ELMER 1.0, it is assumed that at least one of our groups is heterogenous (i.e Breast Invasive Carcinoma primary solid tumors samples might belong to basal,Luminal A, Her2 molecular subtypes), for that reason  for each comparison samples from two most extreme quintiles based on DNA methylation levels were used, generating both the $M$ (methylated) and the $U$ (unmethylated) groups. In the Supervised mode, it is assumed that our groups are homogenous, for that reason the $U$ and $M$ groups are defined strictly by sample group labels, and all samples in each group are used [@ChedraouiSilva148726].

The next steps we will perform ELMER analysis using the supervised mode, which will use all samples from GBM and compare to all LGG samples in all analysis.

### Identification of differentially methylated CpGs (DMCs)

For each distal probe, samples of each group (group 1 and group 2) are ranked by their DNA methylation beta values, those samples in the lower quintile (20\% samples with the lowest methylation levels) of each group are used to identify if the probe is hypomethylated in group 1 compared to group 2, using an unpaired one-tailed t-test. The 20\% is a parameter to the *diff.meth* function called *minSubgroupFrac*. For the (ungrouped) cancer case, this is set to 20\% as in [@yao2015inferring], because we typically wanted to be able to detect a specific molecular subtype among the tumor samples; these subtypes often make up only a minority of samples, and 20\% was chosen as a lower bound for the purposes of statistical power (high enough sample numbers to yield t-test p-values that could overcome multiple hypothesis corrections, yet low enough to be able to capture changes in individual molecular subtypes occurring in 20\% or more of the cases.) This number can be set arbitrarily as an input to the *diff.meth* function and should be tuned based on sample sizes in individual studies. 

As we are dealing with a low number of samples, and we are comparing two cancer types the *minSubgroupFrac* parameter is set to 100\% to use all samples in this analysis. 


```{r, warning=FALSE, results="hide"} 
#--------------------------------------
# STEP 3: Analysis                     |
#--------------------------------------
# Step 3.1: Get diff methylated probes |
#--------------------------------------
message("Get diff methylated probes")
Sig.probes <- get.diff.meth(
  data = mae, 
  group.col = group.col,
  group1 = group1,
  group2 =  group2,
  minSubgroupFrac = 1.0, # Use all samples
  sig.dif = 0.2, # defualt is 0.3
  diff.dir = direction, # Search for hypomethylated probes in group 1
  cores = cores, 
  dir.out = dir.out, 
  pvalue = 0.1
)
```

```{r, warning=FALSE} 
datatable(
  data = Sig.probes[1:10,], 
  options = list(scrollX = TRUE, keys = TRUE, pageLength = 5), 
  rownames = TRUE
)
```

### Identification of putative target gene(s) for differentially methylated distal probes

For each differentially methylated distal probe (DMC), the closest 10 upstream 
genes and the closest 10 downstream genes are tested for inverse correlation between 
methylation of the probe and expression of the gene (the number 10 can be changed using the *numFlankingGenes* parameter). To select these genes, the probe-gene distance is defined as the distance from the probe to the transcription 
start site specified by the ENSEMBL gene level annotations [@yates2015ensembl] accessed via
the R/Bioconductor package [biomaRt](http://bioconductor.org/packages/biomaRt/) [@durinck2009mapping,@durinck2005biomart]. By choosing a constant number of genes to test for each probe, our goal is to avoid systematic false positives for probes in gene-rich regions. This is especially important given the highly non-uniform gene density of mammalian genomes.
Thus, exactly 20 statistical tests were performed for each probe, as follows. 

For each probe-gene pair, the samples (all samples from both groups) are divided into two 
groups: the M group, which consisted of the upper methylation quintile (the 20\%
of samples with the highest methylation at the enhancer probe), and the U group, 
which consists of the lowest methylation quintile (the 20\% of samples with the 
lowest methylation.) The 20\% ile cutoff is a configurable parameter *minSubgroupFrac* in the *get.pair* function.
As with its usage in the *diff.meth* function, the default value of 20\% is a balance, allowing for the identification of changes in a  molecular subtype making up a minority (i.e. 20\%) of cases, while also yielding 
enough statistical power to make strong predictions. For larger sample sizes or other experimental designs, this could be set even lower.

For each candidate probe-gene pair, 
the Mann-Whitney U test is used to test the null hypothesis that overall gene 
expression in group M is greater than or equal than that in group U. 
This non-parametric test was used in order to minimize the effects 
of expression outliers, which can  occur across a very wide dynamic range. 
For each probe-gene pair tested, the raw p-value $P_r$ is corrected for multiple 
hypothesis using a permutation approach as follows.
The gene in the pair is held constant, and *x* random methylation probes are 
chosen to perform the same one-tailed U test, generating a set of *x* permutation
p-values $P_p$. We chose the x random probes only from among those that were 
"distal" (farther than $2kb$ from an annotated transcription start site), in order 
to draw these null-model probes from the same set as the probe being tested [@sham2014statistical]. 
An empirical p-value $P_e$ value was calculated using the following formula 
(which introduces a pseudo-count of 1):

\begin{equation}
P_e = \frac{num(P_p \leq P_r)+ 1}{x+1}
\end{equation}

Notice that in the *Supervised* mode, no additional filtering is necessary to ensure that the M and U group segregate by sample group labels.  The two sample groups are segregated by definition, since these probes were selected for their differential methylation, with the same directionality, between the two groups. 

In our example, we will reduce the number of closest genes evaluated to 5 (the 2 closest upstream 
genes and the closest 2 downstream genes). The *mode* is set to "supervised" which will set group 1 (GBM samples) as the $U$ group, and the group 2 (LGG samples) as $M$ group.
Also, to make the example faster, the number of permutation is reduced to 5 (*permu.size* = 5), the raw p-value cut-off is loosened to 0.1 and the empirical p-value to 0.5, but it is recommended in a real analysis to use the default values.

```{r, warning=FALSE, results="hide"} 
#-------------------------------------------------------------
# Step 3.2: Identify significant probe-gene pairs            |
#-------------------------------------------------------------
# Collect nearby 20 genes for Sig.probes
message("Get nearby genes")
nearGenes <- GetNearGenes(
  data = mae,
  numFlankingGenes = 4, # default is 20 genes
  probes = Sig.probes$probe
)
```

```{R}
length(Sig.probes$probe)
dim(nearGenes)
head(nearGenes)
```

```{r, warning=FALSE, results="hide"} 
message("Get anti correlated probes-genes")
pair <- get.pair(
  data = mae,
  group.col = group.col,
  group1 = group1,
  group2 =  group2,
  nearGenes = nearGenes,
  mode = "supervised",
  minSubgroupFrac = 1, # % of samples to use in to create groups U/M
  raw.pvalue = 0.1,   # defualt is 0.001
  Pe = 0.5, # Please set to 0.001 to get significant results
  filter.probes = TRUE, # See preAssociationProbeFiltering function
  filter.percentage = 0.05,
  save = FALSE, # Create CVS file
  filter.portion = 0.3,
  dir.out = dir.out,
  diff.dir = direction,
  cores = cores,
  label = direction
)
```

```{r, warning=FALSE} 
datatable(
  pair[1:10,], 
  options = list(scrollX = TRUE, keys = TRUE, pageLength = 5), 
  rownames = TRUE
)
```

### Motif enrichment analysis

To identify enriched motifs and potential upstream regulatory TFs,
HOCOMOCO v11 [@kulakovskiy2016hocomoco] TF binding models were used as input for HOMER [@heinz2010simple] to find motif occurrences in a $\pm 250bp$ region around each probe from EPIC and HM450 arrays. Transcription factor (TF) binding models are available at http://hocomoco.autosome.ru/downloads (using the HOMER specific format with threshold score levels corresponding to p-value 
$>1^{-4}$). 

A motif enrichment analysis using Fisher's exact test and Benjamini-Hochberg multiple hypothesis testing correction is performed using as correction background all distal probes.

A probe set was considered significantly enriched 
for a particular motif if the 95\% confidence interval of the Odds Ratio was 
greater than $1.1$ (specified by option *lower.OR*, $1.1$ is default), and the motif 
occurred at least 10 times (specified by option *min.incidence*, $10$ is default) in 
the probe set. 

```{r, warning=FALSE, results="hide"} 
#-------------------------------------------------------------
# Step 3.3: Motif enrichment analysis on the selected probes |
#-------------------------------------------------------------
enriched.motif <- get.enriched.motif(
  data = mae,
  probes = pair$Probe, 
  dir.out = dir.out, 
  label = direction,
  pvalue = 1, # default is FDR < 0.05
  min.incidence = 10,
  lower.OR = 1.1
)
# One of the output from the  previous function is a file with the motif, OR and Number of probes
# It will be used for plotting purposes
motif.enrichment <- read.csv(paste0(dir.out,"/getMotif.",direction, ".motif.enrichment.csv"))
```

```{r, warning=FALSE} 
head(enriched.motif[names(enriched.motif)[1]]) ## probes in the given set that have the first motif.
motif.enrichment %>% head %>% gt::gt()
```

### Identification of master regulator TFs

When a group of enhancers is coordinately altered in a specific sample subset, this is often the result of an altered upstream *master regulator* transcription factor in the gene regulatory network. *ELMER* tries to identify such transcription factors corresponding to each of the TF binding motifs enriched from the previous analysis step.
For each enriched motif, *ELMER* takes the average DNA methylation of all distal probes (in significant probe-gene pairs) that contain that motif occurrence (within a $\pm 250bp$ region) and compares this average DNA methylation to the expression of each gene annotated as a human TF.


In the Unsupervised mode, a statistical test is performed for each motif-TF pair, as follows. All samples 
are divided into two groups: the $M$ group, which consists 
of the 20\% of samples with the highest average methylation at all motif-adjacent
probes, and the $U$ group, which consisted of the 20\%  of samples with the lowest 
methylation. This step is performed by the get.TFs function, which takes *minSubgroupFrac* as an input parameter, again with a default of 20\%.
For each candidate motif-TF pair, the Mann-Whitney U test is used to test 
the null hypothesis that overall gene expression in group $M$ is greater or equal 
than that in group $U$. This non-parametric test was used in order to minimize the 
effects of expression outliers, which can occur across a very wide dynamic range. 
For each motif tested, this results in a raw p-value ($P_r$) for each of the human TFs.

All TFs are ranked by their $-log_{10}(Pr)$ values, and those falling within the top 5\% of
this ranking were considered candidate upstream regulators. The best upstream 
TFs which are known to recognize to specific binding motif are automatically extracted as putative 
regulatory TFs, and rank ordered plots are created to visually inspect these relationships, as shown in the example below. Because the same motif can be recognized by many transcription factors of the same binding domain family, we define these relationships at both the family and subfamily classification level using the 
classifications from TFClass database [@wingender2013tfclass]. 

The Supervised  mode uses the same approach as described for the identification of putative target gene(s) step. The $U$ and $M$ groups are one of the the label group of samples and the *minSubgroupFrac* parameter is set to 100\% to use all samples from both groups in the statistical test. 



```{r, warning=FALSE, results="hide"} 
#-------------------------------------------------------------
# Step 3.4: Identifying regulatory TFs                        |
#-------------------------------------------------------------
TF <- get.TFs(
  data = mae, 
  group.col = group.col,
  group1 = group1,
  group2 =  group2,
  mode = "supervised",
  enriched.motif = enriched.motif,
  dir.out = dir.out, 
  cores = cores,
  save.plots = FALSE,
  diff.dir = direction,
  label = direction
)

# One of the output from the previous function is a file with the raking of TF,
# for each motif. It will be used for plotting purposes
TF.meth.cor <- get(load(paste0(dir.out,"/getTF.",direction,".TFs.with.motif.pvalue.rda")))
```


The output of this step is a data frame with the following columns:

- *motif*:  enriched motif name.
- *top_5percent_TFs*: the top 5\% TFs ranked
- *potential.TFs.family*: TF from the "top_5percent_TFs"  that belongs to the same family as the TF of the motif,  
- *top.potential.TFs.family* is the highest ranked TF belonging to the same family as the TF of the motif (same as the first TF from *potential.TFs.family* column). 
- *potential.TFs.subfamily* and *top.potential.TFs.subfamily* are the same as *potential.TFs.family* and *top.potential.TFs.family*  but considering the subfamily classification instead (you can check the classification at http://hocomoco11.autosome.ru/human/mono?full=true). 


```{r, warning=FALSE} 
datatable(
  TF, 
  options = list(scrollX = TRUE, keys = TRUE, pageLength = 5), 
  rownames = FALSE
)
datatable(
  TF.meth.cor[1:10,1:6], 
  options = list(scrollX = TRUE, keys = TRUE, pageLength = 5), 
  rownames = TRUE
)
```

### ELMER visualization functions

The pairs identified by ELMER can be visualized through a heatmap using the function `heatmapPairs`.
The first heatmap shows the DNA methylation levels, while the second heatmap shows the gene expression
levels of the target genes. The last heatmap shows the distance between the probe and the gene TSS.

```{r, echo=TRUE, message=FALSE,eval = FALSE, warnings=FALSE, results='asis',fig.height=6, fig.width=10}
heatmapPairs(
  data = mae, 
  group.col = group.col,
  group1 = group1, 
  group2 = group2, 
  annotation.col = c("gender"),
  pairs = pair,
  filename =  NULL
)
```


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 motif found.

The list of enriched motifs for the hypo direction (probes hypomethylated in
GBM group compared to the LGG group) with lower OR >= 1.1 is found in the Figure
below. 


```{r, echo=TRUE, message=FALSE, warnings=FALSE, results='asis',fig.height=7,fig.width=15}
motif.enrichment.plot(
  motif.enrichment = motif.enrichment, 
  save = FALSE,
  significant = list(lowerOR = 1.1)
) # Filter motifs in the plot lowerOR > 1.3
```

After finding the enriched motifs,
[ELMER](http://bioconductor.org/packages/ELMER/) identifies regulatory
transcription factors (TFs) whose expression is associated with DNA
methylation at motifs. [ELMER](http://bioconductor.org/packages/ELMER/)
automatically creates a TF ranking plot for each enriched motif. This
plot shows the TF ranking plots based on the association score
$(-log(p-value))$ between TF expression and DNA methylation of the
motif. This plot represents the rank of p-values of all TF, the top 3, 
the TF classified in the same family and sub-family according to TFClass database
are highlighted. The dashed line represents the top 5% of all TF.

```{R, eval=FALSE, include=TRUE}
grid:TF.rank.plot(motif.pvalue=TF.meth.cor, motif=TF$motif[1], save=FALSE)
```	

```{r results='asis', echo=FALSE, message=FALSE,warning=FALSE,  fig.align='center',fig.height=8,fig.width=8}
gg <- TF.rank.plot(motif.pvalue = TF.meth.cor, motif=TF$motif[1], save=FALSE)
grid::grid.draw(gg[[1]])
```	

Also, for each motif, we can take a look at the three most relevant TFs.
For example, Figure below shows the average DNA methylation level
of sites with the first motif plotted against the expression of some top TFs
associated with it. We can see that the GBM samples have a lower average methylation level of sites with
the motif plotted and a higher average expression of the TFs.

```{r results='asis', echo=TRUE, message=FALSE, warning=FALSE}
png("TF.png",width = 800, height = 400)
scatter.plot(
  data = mae, 
  category = group.col, 
  save = FALSE, 
  lm_line = TRUE,
  byTF = list(
    TF = unlist(stringr::str_split(TF[1,"top_5percent_TFs"],";"))[1:4], 
    probe = enriched.motif[[TF$motif[1]]]
  )
)
dev.off()
```

![](TF.png){ width=100%}



# 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](http://bioconductor.org/packages/TCGAbiolinks/) and
[RTCGAToolbox](http://bioconductor.org/packages/RTCGAToolbox/)
Bioconductor packages in order to illustrate how one can acquire TCGA
specific data, followed by key steps for genomics analysis using
[GAIA](http://bioconductor.org/packages/gaia/) package and
[maftools](http://bioconductor.org/packages/maftools/), for
transcriptomic analysis using
[TCGAbiolinks](http://bioconductor.org/packages/TCGAbiolinks/), 
[pathview](http://bioconductor.org/packages/pathview/) packages and for
DNA methylation analysis using
[TCGAbiolinks](http://bioconductor.org/packages/TCGAbiolinks/) package.
An inference of gene regulatory networks was also introduced by
[MINET](http://bioconductor.org/packages/minet/) package. Finally, we
introduced Bioconductor packages
[AnnotationHub](http://bioconductor.org/packages/AnnotationHub/),
[ChIPSeeker](http://bioconductor.org/packages/ChIPseeker/),
[ComplexHeatmap](http://bioconductor.org/packages/ComplexHeatmap/), and
[ELMER](http://bioconductor.org/packages/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.

# Session Information

```{r}
pander::pander(sessionInfo(), compact = FALSE)
```

# Author contributions {#author-contributions .unnumbered}

HN conceived the study. HN, MC and GB provided direction on the design
of the Transcriptomics, Genomics, master regulatory networks and DNA
methylation workflows. TCS developed and tested sections "Experimental
data", "DNA methylation analysis", "Motif analysis" and "Integrative
analysis". AC developed and tested section "Transcriptomic analysis". CO
developed and tested the section "Inference of gene regulatory
networks". FDA developed and tested section "Genomic analysis". TCS, AC,
CO, and FDA prepared the first draft of the manuscript. All authors were
involved in the revision of the draft manuscript and have agreed to the
final content. Also, AC, TS, CO, MC, GB, and HN are authors of the
TCGAbiolinks package and MC is the author of the GAIA package.

# Competing interests {#competing-interests .unnumbered}

No competing interests were disclosed.

# Grant information {#grant-information .unnumbered}

The project was supported by the São Paulo Research Foundation
([FAPESP](http://www.fapesp.br/)) (2015/02844-7 and 2016/01389-7 to
T.C.S. & H.N. and 2015/07925-5 to H.N.), the [BridgeIRIS
project](http://mlg.ulb.ac.be/BridgeIRIS), funded by INNOVIRIS, Region
de Bruxelles Capitale, Brussels, Belgium, and by GENomic profiling of
Gastrointestinal Inflammatory-Sensitive CANcers
([GENGISCAN](http://mlg.ulb.ac.be/GENGISCAN)), Belgian FNRS PDR
(T100914F to G.B.). Funding for open access charge: São Paulo Research
Foundation (FAPESP) (2015/07925-5).

# Acknowledgements {#acknowledgements .unnumbered}

We are grateful to all the authors of the packages used in this article.
Also, we would like to thank The GDC Support Team for the provided help,
which was necessary in order to update the TCGAbiolinks package to use
GDC API.

# References