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