# 11 - ELMER: Use case

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

    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 Supervised analysis

    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 Unsupervised analysis

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