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

## ----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"
# )
# 

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

