inst/doc/usecase.R

In ELMER: Inferring Regulatory Element Landscapes and Transcription Factor Networks Using Cancer Methylomes

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