---
title: "11 - ELMER: Use case"
output: 
  html_document:
    self_contained: true
    number_sections: no
    theme: flatly
    highlight: tango
    mathjax: null
    toc: true
    toc_float: true
    toc_depth: 2
    css: style.css
fontsize: 11pt
bibliography: bibliography.bib    
vignette: >
  %\VignetteIndexEntry{"11 - ELMER: Use case"}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---

<br> 

# Introduction

This section contains the complete ELMER code for the following analysis: 

* Vignette example
* BRCA Supervised analysis
* BRCA Unsupervised analysis

# Vignette example

Below is the complete code that was explained in the other sections.

```{r, eval=FALSE}
library(MultiAssayExperiment)
library(ELMER.data)
library(ELMER)
# get distal probes that are 2kb away from TSS on chromosome 1
distal.probes <- get.feature.probe(
  genome = "hg19", 
  met.platform = "450K", 
  rm.chr = paste0("chr",c(2:22,"X","Y"))
)
data(LUSC_RNA_refined,package = "ELMER.data") # GeneExp
data(LUSC_meth_refined,package = "ELMER.data") # Meth

mae <- createMAE(
  exp = GeneExp, 
  met = Meth,
  save = TRUE,
  linearize.exp = TRUE,
  save.filename = "mae.rda",
  filter.probes = distal.probes,
  met.platform = "450K",
  genome = "hg19",
  TCGA = TRUE
)

group.col <- "definition"
group1 <-  "Primary solid Tumor"
group2 <- "Solid Tissue Normal"
dir.out <- "result"
diff.dir <-  "hypo" # Search for hypomethylated probes in group 1

sig.diff <- get.diff.meth(
  data = mae, 
  group.col = group.col,
  group1 = group1,
  group2 = group2,
  minSubgroupFrac = 0.2,
  sig.dif = 0.3,
  diff.dir = diff.dir,
  cores = 1, 
  dir.out = dir.out,
  pvalue = 0.01
)


nearGenes <- GetNearGenes(
  data = mae, 
  probes = sig.diff$probe, 
  numFlankingGenes = 20
) # 10 upstream and 10 dowstream genes

pair <- get.pair(
  data = mae,
  group.col = group.col,
  group1 = group1,
  mode = "unsupervised",
  group2 = group2,
  nearGenes = nearGenes,
  diff.dir = diff.dir,
  minSubgroupFrac = 0.4, # % of samples to use in to create groups U/M
  permu.dir = file.path(dir.out,"permu"),
  permu.size = 100, # Please set to 100000 to get significant results
  raw.pvalue = 0.05,   
  Pe = 0.01, # Please set to 0.001 to get significant results
  filter.probes = TRUE, # See preAssociationProbeFiltering function
  filter.percentage = 0.05,
  filter.portion = 0.3,
  dir.out = dir.out,
  cores = 1,
  label = diff.dir
)

# Identify enriched motif for significantly hypomethylated probes which 
# have putative target genes.
enriched.motif <- get.enriched.motif(
  data = mae,
  probes = pair$Probe, 
  dir.out = dir.out, 
  label = diff.dir,
  min.incidence = 10,
  lower.OR = 1.1
)

TF <- get.TFs(
  data = mae, 
  mode = "unsupervised",
  group.col = group.col,
  group1 = group1,
  group2 = group2,
  enriched.motif = enriched.motif,
  dir.out = dir.out,
  cores = 1, 
  label = diff.dir
)


```

# BRCA Unsupervised analysis

```{r, eval=FALSE}
library(stringr)
library(TCGAbiolinks)
library(dplyr)
library(ELMER)
library(MultiAssayExperiment)
library(parallel)
library(readr)
dir.create("~/paper_elmer/",showWarnings = FALSE)
setwd("~/paper_elmer/")

file <- "mae_BRCA_hg38_450K_no_ffpe.rda"
if(file.exists(file)) {
  mae <- get(load(file))
} else {
  getTCGA(
    disease = "BRCA", # TCGA disease abbreviation (BRCA,BLCA,GBM, LGG, etc)
    basedir = "DATA", # Where data will be downloaded
    genome  = "hg38"
  ) # Genome of refenrece "hg38" or "hg19"
  
  distal.probes <- get.feature.probe(
    feature = NULL,
    genome = "hg38", 
    met.platform = "450K"
  )
  
  
  mae <- createMAE(
    exp = "~/paper_elmer/Data/BRCA/BRCA_RNA_hg38.rda", 
    met = "~/paper_elmer/Data/BRCA/BRCA_meth_hg38.rda", 
    met.platform = "450K",
    genome = "hg38",
    linearize.exp = TRUE,
    filter.probes = distal.probes,
    met.na.cut = 0.2,
    save = FALSE,
    TCGA = TRUE
  ) 
  # Remove FFPE samples from the analysis
  mae <- mae[,!mae$is_ffpe]
  
  # Get molecular subytpe information from cell paper and more metadata (purity etc...)
  # https://doi.org/10.1016/j.cell.2015.09.033
  file <- "http://ars.els-cdn.com/content/image/1-s2.0-S0092867415011952-mmc2.xlsx"
  downloader::download(file, basename(file))
  subtypes <- readxl::read_excel(basename(file), skip = 2)
  
  subtypes$sample <- substr(subtypes$Methylation,1,16)
  meta.data <- merge(colData(mae),subtypes,by = "sample",all.x = T)
  meta.data <- meta.data[match(colData(mae)$sample,meta.data$sample),]
  meta.data <- S4Vectors::DataFrame(meta.data)
  rownames(meta.data) <- meta.data$sample
  stopifnot(all(meta.data$patient == colData(mae)$patient))
  colData(mae) <- meta.data
  save(mae, file = "mae_BRCA_hg38_450K_no_ffpe.rda")
}
dir.out <- "BRCA_unsupervised_hg38/hypo"
cores <- 10
diff.probes <- get.diff.meth(
  data = mae, 
  group.col = "definition",
  group1 = "Primary solid Tumor",
  group2 = "Solid Tissue Normal",
  diff.dir = "hypo", # Get probes hypometh. in group 1 
  cores = cores,
  minSubgroupFrac = 0.2, # % group samples  used. 
  pvalue = 0.01, 
  sig.dif = 0.3,
  dir.out = dir.out,
  save = TRUE
)

# For each differently methylated probes we will get the 
# 20 nearby genes (10 downstream and 10 upstream)
nearGenes <- GetNearGenes(
  data = mae, 
  probes =  diff.probes$probe, 
  numFlankingGenes = 20
)

# This step is the most time consuming. Depending on the size of the groups
# and the number of probes found previously it migh take hours
Hypo.pair <- get.pair(
  data = mae,
  nearGenes = nearGenes,
  group.col = "definition",
  group1 = "Primary solid Tumor",
  group2 = "Solid Tissue Normal",
  permu.dir = paste0(dir.out,"/permu"),
  permu.size = 10000, 
  mode = "unsupervised",
  minSubgroupFrac = 0.4, # 40% of samples to create U and M
  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 = "hypo"
)
# Number of pairs: 2950 


enriched.motif <- get.enriched.motif(
  data = mae,
  min.motif.quality = "DS",
  probes = unique(Hypo.pair$Probe), 
  dir.out = dir.out, 
  label = "hypo",
  min.incidence = 10,
  lower.OR = 1.1
)
TF <- get.TFs(
  data = mae, 
  group.col = "definition",
  group1 = "Primary solid Tumor",
  group2 = "Solid Tissue Normal",
  minSubgroupFrac = 0.4, # Set to 1 if supervised mode
  enriched.motif = enriched.motif,
  dir.out = dir.out, 
  cores = cores, 
  label = "hypo"
)

```

# BRCA Supervised analysis

```{r, eval=FALSE}
library(stringr)
library(TCGAbiolinks)
library(dplyr)
library(ELMER)
library(MultiAssayExperiment)
library(parallel)
library(readr)
#-----------------------------------
# 1 - Samples
# ----------------------------------
dir.create("~/paper_elmer/",showWarnings = FALSE)
setwd("~/paper_elmer/")

file <- "mae_BRCA_hg38_450K_no_ffpe.rda"
if(file.exists(file)) {
  mae <- get(load(file))
} else {
  getTCGA(
    disease = "BRCA", # TCGA disease abbreviation (BRCA,BLCA,GBM, LGG, etc)
    basedir = "DATA", # Where data will be downloaded
    genome  = "hg38"
  ) # Genome of refenrece "hg38" or "hg19"
  
  distal.probes <- get.feature.probe(
    feature = NULL,
    genome = "hg38", 
    met.platform = "450K"
  ) 
  
  mae <- createMAE(
    exp = "DATA/BRCA/BRCA_RNA_hg38.rda", 
    met = "DATA/BRCA/BRCA_meth_hg38.rda", 
    met.platform = "450K",
    genome = "hg38",
    linearize.exp = TRUE,
    filter.probes = distal.probes,
    met.na.cut = 0.2,
    save = FALSE,
    TCGA = TRUE
  ) 
  # Remove FFPE samples from the analysis
  mae <- mae[,!mae$is_ffpe]
  
  # Get molecular subytpe information from cell paper and more metadata (purity etc...)
  # https://doi.org/10.1016/j.cell.2015.09.033
  file <- "http://ars.els-cdn.com/content/image/1-s2.0-S0092867415011952-mmc2.xlsx"
  downloader::download(file, basename(file))
  subtypes <- readxl::read_excel(basename(file), skip = 2)
  
  subtypes$sample <- substr(subtypes$Methylation,1,16)
  meta.data <- merge(colData(mae),subtypes,by = "sample",all.x = T)
  meta.data <- meta.data[match(colData(mae)$sample,meta.data$sample),]
  meta.data <- S4Vectors::DataFrame(meta.data)
  rownames(meta.data) <- meta.data$sample
  stopifnot(all(meta.data$patient == colData(mae)$patient))
  colData(mae) <- meta.data
  save(mae, file = "mae_BRCA_hg38_450K_no_ffpe.rda")
}

cores <- 6
direction <- c( "hypo","hyper")
genome <- "hg38"
group.col  <- "PAM50"
groups <- t(combn(na.omit(unique(colData(mae)[,group.col])),2))
for(g in 1:nrow(groups)) {
  group1 <- groups[g,1]
  group2 <- groups[g,2]
  for (j in direction){
    tryCatch({
      message("Analysing probes ",j, "methylated in ", group1, " vs ", group2)
      dir.out <- paste0("BRCA_supervised_",genome,"/",group1,"_",group2,"/",j)
      dir.create(dir.out, recursive = TRUE)
      #--------------------------------------
      # STEP 3: Analysis                     |
      #--------------------------------------
      # Step 3.1: Get diff methylated probes |
      #--------------------------------------
      Sig.probes <- get.diff.meth(
        data       = mae,
        group.col  = group.col,
        group1     = group1,
        group2     = group2,
        sig.dif    = 0.3,
        minSubgroupFrac = 1,
        cores      = cores,
        dir.out    = dir.out,
        diff.dir   = j,
        pvalue     = 0.01
      )
      if(nrow(Sig.probes) == 0) next
      #-------------------------------------------------------------
      # Step 3.2: Identify significant probe-gene pairs            |
      #-------------------------------------------------------------
      # Collect nearby 20 genes for Sig.probes
      nearGenes <- GetNearGenes(
        data  = mae,
        probe = Sig.probes$probe
      )
      
      pair <- get.pair(
        data       = mae,
        nearGenes  = nearGenes,
        group.col  = group.col,
        group1     = group1,
        group2     = group2,
        permu.dir  = paste0(dir.out,"/permu"),
        dir.out    = dir.out,
        mode       = "supervised", 
        diff.dir   = j,
        cores      = cores,
        label      = j,
        permu.size = 10000,
        raw.pvalue = 0.001
      )
      
      Sig.probes.paired <- readr::read_csv(
        paste0(dir.out,
               "/getPair.",j,
               ".pairs.significant.csv")
      )[,1, drop = TRUE]
      
      
      #-------------------------------------------------------------
      # Step 3.3: Motif enrichment analysis on the selected probes |
      #-------------------------------------------------------------
      if(length(Sig.probes.paired) > 0 ){
        #-------------------------------------------------------------
        # Step 3.3: Motif enrichment analysis on the selected probes |
        #-------------------------------------------------------------
        enriched.motif <- get.enriched.motif(
          probes  = Sig.probes.paired,
          dir.out = dir.out,
          data    = mae,
          label   = j,
          plot.title =  paste0("BRCA: OR for paired probes ",
                               j, "methylated in ",
                               group1, " vs ",group2)
        )
        motif.enrichment <- readr::read_csv(
          paste0(dir.out,
                 "/getMotif.",j,
                 ".motif.enrichment.csv")
        )
        if(length(enriched.motif) > 0){
          #-------------------------------------------------------------
          # Step 3.4: Identifying regulatory TFs                        |
          #-------------------------------------------------------------
          print("get.TFs")
          
          TF <- get.TFs(
            data           = mae,
            enriched.motif = enriched.motif,
            dir.out        = dir.out,
            mode           = "supervised",
            group.col      = group.col,
            group1         = group1,
            diff.dir       = j,
            group2         = group2,
            cores          = cores,
            label          = j
          )
          TF.meth.cor <- get(
            load(paste0(dir.out, "/getTF.",j, ".TFs.with.motif.pvalue.rda"))
          )
          save(
            mae, TF, enriched.motif, Sig.probes.paired,
            pair, nearGenes, Sig.probes, motif.enrichment,
            TF.meth.cor,
            file = paste0(dir.out,"/ELMER_results_",j,".rda")
          )
        }
      }
    }, error = function(e){
      message(e)
    })
  }
}
```