```{r setup, include=FALSE}
knitr::opts_chunk$set(dpi = 300)
```

# Instructor names and contact information

1. Benjamin P. Berman
2. Tiago C. Silva

# Workshop Description

This workshop demonstrates how to perform ELMER analysis using matched RNA-seq
and DNA methylation data.

You can find a detailed information about ELMER at http://bioconductor.org/packages/3.10/bioc/vignettes/ELMER/inst/doc/index.html.

Articles about ELMER:

* Tiago C Silva, Simon G Coetzee, Nicole Gull, Lijing Yao, Dennis J Hazelett, Houtan Noushmehr, De-Chen Lin, Benjamin P Berman, ELMER v.2: an R/Bioconductor package to reconstruct gene regulatory networks from DNA methylation and transcriptome profiles, Bioinformatics, Volume 35, Issue 11, 1 June 2019, Pages 1974–1977, https://doi.org/10.1093/bioinformatics/bty902
* Yao, Lijing, et al. "Inferring regulatory element landscapes and transcription factor networks from cancer methylomes." Genome biology 16.1 (2015): 105. https://doi.org/10.1186/s13059-015-0668-3
* Yao, Lijing, Benjamin P. Berman, and Peggy J. Farnham. "Demystifying the secret mission of enhancers: linking distal regulatory elements to target genes." Critical reviews in biochemistry and molecular biology 50.6 (2015): 550-573. https://doi.org/10.3109/10409238.2015.1087961

## Pre-requisites

* Basic knowledge of R syntax
* Familiarity with the SummarizedExperiment classes
* Familiarity with ’omics data types including DNA methylation and gene expression

## Workshop Participation

Students will have a chance to run ELMER analysis on a provided `MultiAssayExperiment` 
object created from TCGA data from GDC data portal.

## Goals and objectives

* gain familiarity ELMER input data, a MultiAssayExperiment object
* Execute ELMER analysis on real data and understand its meaning

#  R/Bioconductor packages used

* [ELMER](http://bioconductor.org/packages/ELMER/)
* [MultiAssayExperiment](http://bioconductor.org/packages/MultiAssayExperiment/)


```{r, message=FALSE}
library("ELMER")
library("MultiAssayExperiment")
```

# Retrieving the data

The RNA-seq data, DNA methylation data and patients metadata (i.e age, gender) used in this workshop is 
structured as a MultiAssayExperiment. You can read more about a MultiAssayExperiment object at 	http://bioconductor.org/packages/MultiAssayExperiment/ 
and https://bioconductor.github.io/BiocWorkshops/workflow-for-multi-omics-analysis-with-multiassayexperiment.html.

Through the next section we will:

1. load the data
2. Verify the RNA-seq data
3. Verify the DNA methylation data
4. Verify the samples metadata

## Download the data

The data is available in this [google drive](https://drive.google.com/drive/folders/1N9GKy8rEdmNsPa3W1vZ7WnBaUlXky2L4?usp=sharing).

## Loading the data

```{r}
mae <- readRDS("Data/TCGA_ESCA_MAE_distal_regions.rds")
mae
```

## Verify the RNA-seq data

```{r}
# check experiments
experiments(mae)

# Get the Gene expression object
rna.seq <- mae[[2]]

# nb of genes
nrow(rna.seq)

# nb of samples
ncol(rna.seq)

# Check genes metadata
rowRanges(rna.seq)

# Check genes expression
assay(rna.seq)[1:4,1:4]
```

## Verify the DNA methylation data

```{r}
# check experiments
experiments(mae)

# Get the DNA methylation object
dna.met <- mae[[1]]

# nb of DNA methylation probes
nrow(dna.met)

# nb of samples
ncol(dna.met)

# Check DNA methylation probes metadata
rowRanges(dna.met)[,1:4]

# Check DNA methylation beta values
assay(dna.met)[1:4,1:4]
```


## Verify the samples metadata

```{r}
# metadata can be accessed using the function colData
# we will check the first 4 samples and their 4 columns
colData(mae)[1:4,1:4]

# or you can also access the metadata directly with the $
mae$primary_diagnosis[1:4]

# you can summarize the number of tissue type with the function table
table(mae$definition)

# or summarize the numbers of samples using two features
table(mae$primary_diagnosis,mae$definition)
```

# Selecting the samples

Since we will use only Primary solid Tumor samples we will remove the Metastatic and 
Solid Tissue Normal samples from our object
```{r}
mae <- mae[, mae$definition == "Primary solid Tumor"]
```


# Performing ELMER analysis

We will first define some parameters used in the majority of ELMER function
```{r}
# which is the column of the samples metadata defining our groups
group.col <- "primary_diagnosis"

# which are the grops withing the group.col that we want to compare 
group1 <- "Adenocarcinoma, NOS"
group2 <- "Squamous cell carcinoma, NOS"

# Identify probes hypo in group1 or hypo in group2 ? 
direction <- "hypo" # hypo in group 1

cores <- 2
mode <- "supervised"

# Where do we want to save the results ?
dir.out <- paste0("analysis/",group.col,"-",
                  gsub("[[:punct:]]| ", "_", group1),
                  "_vs_",
                  gsub("[[:punct:]]| ", "_", group2),
                  "_", mode,
                  "/", direction)
dir.create(dir.out, showWarnings = FALSE,recursive = T)
```
## Identifying hypo methylated distal probes between ESAD samples compared to ESCC samples

```{r}
# time expected with one core: 4 minuntes
diff.probes <- get.diff.meth(data = mae, 
                             group.col = group.col,
                             group1 = group1,
                             group2 = group2,
                             diff.dir = direction, # Get probes hypometh. in group 1 
                             cores = cores, 
                             mode = mode,
                             pvalue = 0.01, 
                             sig.dif = 0.3,
                             dir.out =  dir.out, 
                             save = TRUE)
# check the results
head(diff.probes)
```

## Correlate RNA-seq expression and DNA methylation

Now that we have our hypomethylated probes we will look for the 20 nearest genes  
to each probes to identify  if any of these genes could being affected by the 
hypomethylation in the distal regulatory region.


```{r}
# time expected less than one minute
nearGenes <- GetNearGenes(data = mae, 
                          probes =  diff.probes$probe, 
                          numFlankingGenes = 20)
head(nearGenes)
```



Now that we have our hypomethylated probes and the 20 nearest genes we want to identify 
if any of these genes could being affected by the hypomethylatio in the distal regulatory region.


```{r}
# time expected  one minute
pair <- get.pair(data = mae,
                 nearGenes = nearGenes,
                 group.col = group.col,
                 group1 = group1,
                 group2 = group2,
                 mode = mode,
                 diff.dir = direction, 
                 raw.pvalue = 0.001,   
                 Pe = 0.001, 
                 filter.probes = TRUE,
                 filter.percentage = 0.05,
                 filter.portion = 0.3,
                 dir.out =  dir.out,
                 cores = cores, 
                 label = direction)

head(pair)
```

We can plot the correlation heatmap of DNA methylation and gene expression
with the `heatmapPairs` function. We will plot only the first 100 pairs.

```{r, fig.width=10, fig.height=4}
heatmapPairs(data = mae, 
             group.col = group.col,
             group1 = group1,
             group2 = group2,
             subset = TRUE,
             pairs = pair[1:100,],
             filename =  NULL)
```

## Identifying TF enriched motifs on the hypomethylated regions

We have hypomethylated regions that has a distal gene upregulated.
We want to identify which is the MR TF binding to those regions and regulating the distal gene.
First we identify which are the enriched motifs for those hypomethylated regions using 
the function `get.enriched.motif`. This will output the enriched motifs and the respective
hypomethylated probes with the motif signature around it 
(since a DNA methylation probes is 1bp, we extended the search window to $\pm$ 250 bp).

```{r,echo=TRUE, results='hide', fig.keep='all', message = FALSE}
# time expected  less than one minute
enriched.motif <- get.enriched.motif(data = mae,
                                     probes = unique(pair$Probe), 
                                     dir.out = dir.out,
                                     label = direction,
                                     min.incidence = 10,
                                     lower.OR = 1.5)
```

```{r,echo=TRUE}
# check enriched motifs
names(enriched.motif)[1:4]

# and the paired hypomethylated probes with the motif signature
enriched.motif$GATA4_HUMAN.H11MO.0.A
```


```{r}
motif.result.file <- file.path(dir.out,"getMotif.hypo.motif.enrichment.csv")

motif.enrichment.plot(motif.enrichment = motif.result.file,
                      significant = list(lowerOR = 2.0), 
                      label = "hypo", 
                      summary = TRUE,
                      save = FALSE)  
```

## Inferring MR TF binding to hypomethylated regions

Now that we have our enriched motifs, we have the following question: 
- Which is the candidate master regulator Transcription factor binding to those regions ?

It is important to highlight that since TF within the same family and subfamily have 
very close motifs, they will likely be indentified to bind the same regions.
You can check that by verifying the intersection of probes from our enrichment analysis results.

```{r}
length(enriched.motif$GATA4_HUMAN.H11MO.0.A)
length(enriched.motif$GATA6_HUMAN.H11MO.0.A)
length(intersect(enriched.motif$GATA4_HUMAN.H11MO.0.A,enriched.motif$GATA6_HUMAN.H11MO.0.A))
```

So more than 70% of probes are enriched for both GATA4 and GATA6.
How could we infer which of those MR might be the one with higher impact ?

For that ELMER correlates the TF expression vs the average mean methylation of the hypomethylated probes
with a motif signature and rank them. We expect that a TF with higher expression and lower methylation in the 
binding regions will be a better MR TF candidate than TF that are for example equally expressed on both groups.

```{r}
# time expected  ~5 minutes
TF <- get.TFs(data = mae, 
              group.col = group.col,
              mode = mode,
              group1 = group1,
              group2 = group2,
              diff.dir = direction,
              enriched.motif = enriched.motif,
              dir.out = dir.out,
              label = direction, 
              cores = cores)
# Check the results for the GATA4_HUMAN.H11MO.0.A motif
TF["GATA4_HUMAN.H11MO.0.A",1:5]
```

### Visualizing the results with TF ranking plot

ELMER provides an easy way to visualize the rank of all human TF evaluated with the
function `TF.rank.plot`. 

```{r,echo=FALSE, results='hide', fig.keep='all', message = FALSE}
load(file.path(dir.out,"getTF.hypo.TFs.with.motif.pvalue.rda"))
motif <-  "GATA4_HUMAN.H11MO.0.A"
TF.rank.plot(motif.pvalue = TF.meth.cor, 
             motif = motif,
             save = FALSE) 
```

# Summarizing several ELMER analysis

## MR TF heatmap 


When you perform several ELMER analysis with different arguments, you might want to summarize and visualize those several results. For that, we have the `ELMER:::get.tabs` function , which accepts as input 
the directories of the results and if you want to select the MR TF within the 
family of subfamily classification.

In the `Data/analysis_april_2019` folder we provided 4 ELMER analysis. The first step is to retrieve the directory names with the results. Also, it is a golden rule when running several ELMER runs to add a understable name to the `dir.out` directories, such that you easily know which analysis were ran, with which arguments.
In the example below, you can identify the tumor type, type of ELMER run (supervised/unsupervised), groups compared and direction of the analysis (hypo or hyper in group 1).

In the directory `analysis_april_2019` we provide 4 results from different ELMER runs.
```{r}
analsysis.dir <- unique(dirname(dir("Data/analysis_april_2019/",recursive = T,full.names = T)))
analsysis.dir
```

With those directory names we will use the function `ELMER:::get.tabs` to summarize the results.
```{r,echo=TRUE,warning=FALSE}
classification <- "subfamily"
tabs <- ELMER:::get.tabs(analsysis.dir,classification)
names(tabs)

# To be more readable we will change the path names to the analysis name
analysis.names <- c( "COAD-READ (vs Normal)",
                     "EAC (vs Normal)",
                     "EAC (vs ESCC)",
                     "PAAD (vs Normal)")
for (i in 1:3) colnames(tabs[[i]]) <- analysis.names
```

This function creates 4 matrices using some of the files in each analysis directory:
1.  a binary matrix saying if the TF was found or not in the analysis: this function parses each `getTF.hypo.significant.TFs.with.motif.summary.csv` (create after you run the `get.TFs` function) file that has the MR TF for each enriched motif.

```{r,echo=TRUE,warning=FALSE}
# 1.  a binary matrix saying if the TF was found or not in the analysis
tab.binary <- tabs$tab
head(tab.binary)
```

2. a matrix  with the MR TF p-value found in the analysis: this uses `getTF.hypo.significant.TFs.with.motif.summary.csv` to identify the MR TF and `getTF.hypo.TFs.with.motif.pvalue.rda`, which is a matrix with TF as columns and enriched motifs as rows and the contains the p-value of the wilcoxon test of the TF expression for the Unmethylated vs Methylated groups.

```{r,echo=TRUE,warning=FALSE}

# 2. a matrix  with the MR TF p-value found in the analysis
tab.pval <- tabs$tab.pval
head(tab.pval)
```

3. a matrix  with the highest motif OR for the MR TF found in the analysis: this function uses the file `getTF.hypo.significant.TFs.with.motif.summary.csv`  to get the MR TF and then uses `getMotif.hypo.motif.enrichment.csv` to get the motifs OR to which the MR TF is known to bind and returns the highest OR.

```{r,echo=TRUE,warning=FALSE}

# 3. a matrix  with the highest motif OR for the MR TF found in the analysis
tab.or <- tabs$tab.or
head(tab.or)
```

4. a matrix with 4 collums: MR TF, motif, FDR and analysis it was identified: This function uses the same files as the previous matrices, but returns the results in a different format.

```{r,echo=TRUE,warning=FALSE}
# 4. a matrix  with 4 collums: MR TF, motif, FDR and analysis it was identified.
tf.or.table <- tabs$tf.or.table

# Rename analysis
tf.or.table$analysis <- gsub("ESCA","EAC",paste0(basename(dirname(dirname(dirname(tf.or.table$analysis)))), 
                               " (vs ",
                               ifelse(grepl("Normal",tf.or.table$analysis),"Normal","ESCC"),")"))
head(tf.or.table)
```

## Creating Summary plots

With those summarized results we can create two summary plots: 

1. A heatmap with the MR TFs identified in each analysis. 
2. A scatter plot of the TF expression vs the avg DNA methylation of the TFBS. 

### MR TF heatmap
The code below creates the heatmap using complexHeatmap package for the binary MR TF matrix and the p-value one.

```{r, echo = TRUE, warning = FALSE, results = 'hide', fig.keep = 'all', message = FALSE}
library(ComplexHeatmap)
library(vegan)

col <- ifelse(classification == "family","top.potential.TF.family", "top.potential.TF.subfamily")

analysis.colors <- c(
  "COAD-READ (vs Normal)" = '#ff964f',
  "EAC (vs Normal)" = '#94e894',
  "EAC (vs ESCC)" = '#eaadb9',
  "PAAD (vs Normal)" = '#9ab9f9'                   
)

# get top TFs in each anaylsis
labels <- ELMER:::get.top.tf.by.pval(tab.pval,top = nrow(tab.pval))

hb = HeatmapAnnotation(df = data.frame("Mode" = rep("unsupervised",4),
                                       "Analysis" = analysis.names),
                       col = list("Analysis" = analysis.colors,
                                  "Mode" = c("unsupervised" = "#ebd540",
                                             "supervised" = "#718da5")),
                       show_annotation_name = FALSE,
                       show_legend = T,
                       annotation_name_side = "left",
                       annotation_name_gp = gpar(fontsize = 12))

# Binary heatmap
ht_list_binary <- Heatmap(as.matrix(tab.binary),
                   name = "Binary heatmap" ,
                   col = c("1" = "black", "0" = "white"),
                   column_names_gp = gpar(fontsize = 8),
                   show_column_names = FALSE,
                   use_raster = TRUE,
                   na_col = "white",
                   raster_device = c("png"),
                   raster_quality = 2,
                   show_row_names = F,
                   show_heatmap_legend = TRUE,
                   cluster_columns = FALSE,
                   cluster_rows = TRUE,
                   show_column_dend = FALSE,
                   show_row_dend = FALSE,
                   clustering_distance_rows = function(x) {
                     vegan::vegdist(tab.binary,method = "jaccard",binary = TRUE)
                   },
                   clustering_method_rows = "average",
                   top_annotation = hb,
                   width = unit(5, "cm"),
                   row_names_gp = gpar(fontsize = 1),
                   column_title = paste0("Binrary MR TF Heatmap"), 
                   column_title_gp = gpar(fontsize = 10), 
                   row_title_gp = gpar(fontsize = 16)) 

ht_list_binary <-  ht_list_binary + 
  rowAnnotation(link = anno_mark(at = match(labels,rownames(tab.pval)), 
                                 labels = labels,labels_gp = gpar(fontsize = 10)),
                width = unit(1, "cm") + max_text_width(labels)
  ) 
draw(ht_list_binary,
     newpage = FALSE, 
     column_title_gp = gpar(fontsize = 12, fontface = "bold"),
     heatmap_legend_side = "left",
     show_heatmap_legend = TRUE,
     annotation_legend_side = "right")

# P-value heatmap
ht_list <- Heatmap(scale(as.matrix(-log10(tab.pval))),
                   name = "column z-score FDR",
                   col = colorRampPalette(c('#3361A5', '#248AF3', '#14B3FF', 
                                            '#88CEEF', '#C1D5DC', '#EAD397', 
                                            '#FDB31A','#E42A2A', '#A31D1D'))(100),
                   column_names_gp = gpar(fontsize = 8),
                   show_column_names = FALSE,
                   use_raster = TRUE,
                   na_col = "white",
                   raster_device = c("png"),
                   raster_quality = 2,
                   show_row_names = F,
                   show_heatmap_legend = TRUE,
                   cluster_columns = FALSE,
                   cluster_rows = FALSE,
                   row_order = row_order(ht_list_binary),
                   top_annotation = hb,
                   width = unit(5, "cm"),
                   row_names_gp = gpar(fontsize = 1),
                   column_title = paste0("MR TF Heatmap"), 
                   column_title_gp = gpar(fontsize = 10), 
                   row_title_gp = gpar(fontsize = 16)) 

ht_list <-  ht_list + 
  rowAnnotation(link = anno_mark(at = match(labels,rownames(tab.pval)), 
                                 labels = labels,labels_gp = gpar(fontsize = 10)),
                width = unit(1, "cm") + max_text_width(labels)
  ) 

draw(ht_list,
     newpage = TRUE,     
     column_title_gp = gpar(fontsize = 12, fontface = "bold"),
     heatmap_legend_side = "left",
     show_heatmap_legend = TRUE,
     annotation_legend_side = "right")
```

## MR TF expression vs DNA methylation at binding motifs

The main difference from the last plot is the requirement of the MAE object with the DNA methylation and gene expression data. 
For a given MR TF, we need to indentify the top enriched motif it binds and then for each enriched motifs its distal paired differently methylated probes using the `getMotif.hypo.enriched.motifs.rda` object within the results folder.

Also, since we only provide one object used for two analysis ("EAC (vs Normal)"  and "EAC (vs ESCC)" )
we will only plot the  MR TF expression vs DNA methylation at binding motifs
for them.


```{r,echo=TRUE, results='hide', fig.keep='all', message = FALSE,fig.height=8,fig.width=15,warning=FALSE}
analsysis.dir.esca <- unique(dirname(dir("Data/analysis_april_2019/ESCA/",recursive = T,full.names = T)))
mae <- readRDS("Data/TCGA_ESCA_MAE_distal_regions.rds")

tab.pval.esca <- tab.pval[,2:3]

suppressMessages({
library(ELMER)
library(grid)
library(ggplot2)
library(ggpubr)
library(dplyr)
})


analysis.colors <- c(
  "COAD-READ (vs Normal)" = '#ff964f',
  "EAC (vs Normal)" = '#94e894',
  "EAC (vs ESCC)" = '#eaadb9',
  "PAAD (vs Normal)" = '#9ab9f9'                   
)

# Top has the number of MRTF to plot
for (top in c(5)) {
  
  scatter.list <- list()
  scatter.grid <- list()
  scatter.labels <- c()
  
  for (path in analsysis.dir.esca) {
    # Load of object with enriched motifs and the probes mapped to it
    load(dir(path = path, 
             pattern = "getMotif.hypo.enriched.motifs.rda",
             recursive = T,full.names = T)    
    )
    i <- which(path == analsysis.dir.esca)
    
    tumor <- basename(dirname(dirname(dirname(path))))
    mae.plot <- mae
      
     # Change the data to plot based on the analysis performed  
    if (grepl("Squamous",path)) {
      category <- "primary_diagnosis"
      # keep only samples used in the analysis
      mae.plot <- mae.plot[,mae.plot$primary_diagnosis %in% c("Adenocarcinoma, NOS", 
                                                              "Squamous cell carcinoma, NOS") & 
                             mae.plot$definition %in% c("Primary solid Tumor")]          
    } else {
      category <- "definition"
      # keep only samples used in the analysis
      mae.plot <- mae.plot[,mae.plot$primary_diagnosis %in% c("Adenocarcinoma, NOS") & 
                             mae.plot$definition %in% c("Primary solid Tumor","Solid Tissue Normal")]
    }
    # get top Tfs
    topTFs <- rownames(tab.pval)[head(sort(DelayedArray::rowMins(tab.pval.esca[,i,drop = F],na.rm = T), 
                                           index.return = TRUE, decreasing = F)$ix, n = top)]
    topTFs <- topTFs[1:top]
    for (j in topTFs) {
      TF <- j
      scatter.labels <- c(scatter.labels,colnames(tab.pval.esca)[i])
      
      motif <- tf.or.table[tf.or.table$TF == j & 
                             tf.or.table$analysis == colnames(tab.pval.esca)[i],"motif"][1]
      
      # Change the dots color based on the analysis performed  
      if (grepl("Squamous",path))   color.value <- c("#FF0000","#0000FF")
      if (grepl("Adenocarcinoma.Solid.Tissue.Normal",path)) color.value <- c("#FF0000","#176d55")
      
      if (is.na(motif)) next
        
        # Plot TF expression vs avg DNA methylation of the TFBS
        suppressMessages({
          scatter <- scatter.plot(data = mae.plot,
                              correlation = TRUE,
                              byTF = list(TF = TF,
                                          probe = enriched.motif[[motif]]), 
                              category = category,
                              ylim = c(0,25), 
                              dots.size = 0.8,
                              color.value = color.value,
                              save = FALSE, 
                              lm_line = TRUE) +
            ylab(bquote(.(TF) ~ log[2](RSEM + 1))) + 
            xlab(paste0("Avg. DNA met. at \n", length(enriched.motif[[motif]]),
                    " CpGs w/ \n enriched motif ",
                    sub("_HUMAN.H11MO.[0-1].[A-D]", "",  motif))) + 
            theme(legend.position = "none", 
                  plot.title = element_blank(),
                  axis.title.y = element_text(size = 8, face = "bold"),
                  axis.title.x = element_text(size = 8, face = "bold"),
                  strip.background = element_rect(fill = analysis.colors[i + 1]),
                  panel.grid.major = element_blank(),
                  panel.grid.minor = element_blank(),
                  # Remove panel background
                  panel.background = element_blank()
        )
      })
      
      # If this is the first plot of the analysis we will add a labels (rectangle + text)  
      if(j == topTFs[1]){
        gt <- gtable::gtable(widths = unit(0.5,"cm"), heights = unit(3, "in")) %>%
          gtable::gtable_add_grob(grid::rectGrob(gp=gpar(fill= analysis.colors[i + 1])), 1, 1, name = 1) %>% 
          gtable::gtable_add_grob(grid::textGrob(colnames(tab.pval.esca)[i], 
                                                 rot = 90, 
                                                 gp = grid::gpar(fontsize = 8)), 1, 1, name = 2)
        scatter <-  scatter %>% ggpubr::annotate_figure(left = gt)
      }
      scatter.list[[paste0(i,j)]] <-  scatter
      
    }
    scatter.grid[[i]] <-  ggpubr::ggarrange(plotlist = scatter.list[grep(paste0("^",i),names(scatter.list))], 
                                            ncol = top,
                                            # the first column should be wider since it has the label
                                            widths = c(1.15,rep(1,top - 1)), 
                                            nrow = 1, 
                                            common.legend = T,
                                            legend = "bottom")
  }     
}
options(repr.plot.width=20, repr.plot.height=8)

# we use ggarrange to arrange the plots together
ggpubr::ggarrange(plotlist = scatter.grid, 
                  ncol = 1,
                  # the first column should be wider since it has the label
                  widths = c(1.15,rep(1,top - 1)), 
                  nrow = 2, # Number of analysis - rows of the plot
                  heights = c(rep(1.5,length(scatter.list) - 1), 2), 
                  common.legend = F,
                  legend = "bottom")
```

# Code to create the data used in the workshop

The function `getTCGA` uses TCGAbiolinks to retrieve TCGA data in ELMER accept format.


```{r, eval=FALSE}
genome <- "hg38"

# Get RNA-seq and DNA methylation from GDC using TCGAbiolinks
# Data will be saved in Data folder by default
getTCGA(disease = "ESCA", genome = genome, Meth = TRUE,RNA = TRUE)

# We will keep only the distal probes
distal.probes <- get.feature.probe(feature = NULL,
                                   genome = genome,
                                   met.platform = "450K") 

mae <- createMAE(exp = paste0("Data/",tumor,"/",tumor,"_RNA_hg38_no_filter.rda"), 
                   met = paste0("Data/",tumor,"/",tumor,"_Meth_hg38_no_filter.rda"), 
                   met.platform = "450K",
                   genome = "hg38",
                   linearize.exp = TRUE,
                   filter.probes = distal.probes,
                   save = FALSE,
                   TCGA = TRUE) 
mae$group <- paste0(mae$name,"-",mae$definition)
mae <- mae[,!mae$is_ffpe] # remove FFPE samples from analysis
saveRDS(object = mae, file = "TCGA_ESCA_MAE_distal_regions.rds")
```

# Session information
```{r}
sessionInfo()
```

# Workshop materials

## Workshops HTMLs

* ELMER data Workshop HTML: http://rpubs.com/tiagochst/elmer-data-workshop-2019
* ELMER analysis Workshop HTML: http://rpubs.com/tiagochst/ELMER_workshop
* ATAC-seq Workshop HTML: http://rpubs.com/tiagochst/atac_seq_workshop


## Workshop videos

We have a set of recorded videos, explaining some of the workshops.

* All videos playlist: https://www.youtube.com/playlist?list=PLoDzAKMJh15kNpCSIxpSuZgksZbJNfmMt
* ELMER algorithm: https://youtu.be/PzC31K9vfu0
* ELMER data: https://youtu.be/R00wG--tGo8
* ELMER analysis part1 : https://youtu.be/bcd4uyxrZCw
* ELMER analysis part2: https://youtu.be/vcJ_DSCt4Mo
* ELMER summarizing several analysis: https://youtu.be/moLeik7JjLk
* ATAC-Seq workshop: https://youtu.be/3ftZecz0lU4

