inst/doc/MethReg.R

In MethReg: Assessing the regulatory potential of DNA methylation regions or sites on gene transcription

## ----settings, include = FALSE--------------------------------------------------------------------
options(width = 100)
knitr::opts_chunk$set(collapse = TRUE, comment = "#>",class.source = "whiteCode")
library(dplyr)

## ----sesame, include = FALSE----------------------------------------------------------------------
library(sesameData)

## ---- eval = FALSE--------------------------------------------------------------------------------
#  if (!"BiocManager" %in% rownames(installed.packages()))
#       install.packages("BiocManager")
#  BiocManager::install("MethReg", dependencies = TRUE)

## ----setup, eval = TRUE---------------------------------------------------------------------------
library(MethReg)

## ----workflow, fig.cap = "MethReg workflow", echo = FALSE, fig.width = 13-------------------------
jpeg::readJPEG("workflow_methReg.jpg") %>% grid::grid.raster()

## ----warning=FALSE--------------------------------------------------------------------------------
data("dna.met.chr21")

## -------------------------------------------------------------------------------------------------
dna.met.chr21[1:5,1:5]

## -------------------------------------------------------------------------------------------------
dna.met.chr21.se <- make_dnam_se(
  dnam = dna.met.chr21,
  genome = "hg38",
  arrayType = "450k",
  betaToM = FALSE, # transform beta to m-values 
  verbose = FALSE # hide informative messages
)

## -------------------------------------------------------------------------------------------------
dna.met.chr21.se
SummarizedExperiment::rowRanges(dna.met.chr21.se)[1:4,1:4]

## -------------------------------------------------------------------------------------------------
data("gene.exp.chr21.log2")
gene.exp.chr21.log2[1:5,1:5]

## -------------------------------------------------------------------------------------------------
gene.exp.chr21.se <- make_exp_se(
  exp = gene.exp.chr21.log2,
  genome = "hg38",
  verbose = FALSE
)
gene.exp.chr21.se
SummarizedExperiment::rowRanges(gene.exp.chr21.se)[1:5,]

## ----plot, fig.cap = "Different target linking strategies", echo = FALSE--------------------------
png::readPNG("mapping_target_strategies.png") %>% grid::grid.raster()

## ---- message = FALSE, results = "hide"-----------------------------------------------------------
triplet.promoter <- create_triplet_distance_based(
  region = dna.met.chr21.se,
  target.method = "genes.promoter.overlap",
  genome = "hg38",
  target.promoter.upstream.dist.tss = 2000,
  target.promoter.downstream.dist.tss = 2000,
  motif.search.window.size = 500,
  motif.search.p.cutoff  = 1e-08,
  cores = 1  
)

## ---- message = FALSE, results = "hide"-----------------------------------------------------------
# Map probes to genes within 500kb window
triplet.distal.window <- create_triplet_distance_based(
  region = dna.met.chr21.se,
    genome = "hg38", 
    target.method = "window",
    target.window.size = 500 * 10^3,
    target.rm.promoter.regions.from.distal.linking = TRUE,
    motif.search.window.size = 500,
    motif.search.p.cutoff  = 1e-08,
    cores = 1
)

## ---- message = FALSE, results = "hide"-----------------------------------------------------------
# Map probes to 5 genes upstream and 5 downstream
triplet.distal.nearby.genes <- create_triplet_distance_based(
  region = dna.met.chr21.se,
    genome = "hg38", 
    target.method = "nearby.genes",
    target.num.flanking.genes = 5,
    target.window.size = 500 * 10^3,
    target.rm.promoter.regions.from.distal.linking = TRUE,
    motif.search.window.size = 500,
    motif.search.p.cutoff  = 1e-08,
    cores = 1  
)

## ---- eval = FALSE--------------------------------------------------------------------------------
#  if (!"BiocManager" %in% rownames(installed.packages()))
#       install.packages("BiocManager")
#  BiocManager::install("remap-cisreg/ReMapEnrich", dependencies = TRUE)

## ---- eval = FALSE--------------------------------------------------------------------------------
#  library(ReMapEnrich)
#  remapCatalog2018hg38 <- downloadRemapCatalog("/tmp/", assembly = "hg38")
#  remapCatalog <- bedToGranges(remapCatalog2018hg38)

## ---- eval = FALSE--------------------------------------------------------------------------------
#  #-------------------------------------------------------------------------------
#  # Triplets promoter using remap
#  #-------------------------------------------------------------------------------
#  triplet.promoter.remap <- create_triplet_distance_based(
#    region = dna.met.chr21.se,
#    genome = "hg19",
#    target.method =  "genes.promoter.overlap",
#    TF.peaks.gr = remapCatalog,
#    motif.search.window.size = 500,
#    max.distance.region.target = 10^6,
#  )

## ---- eval = FALSE--------------------------------------------------------------------------------
#  if (!"BiocManager" %in% rownames(installed.packages()))
#       install.packages("BiocManager")
#  BiocManager::install("dorothea", dependencies = TRUE)

## -------------------------------------------------------------------------------------------------
regulons.dorothea <- dorothea::dorothea_hs
regulons.dorothea %>% head

## -------------------------------------------------------------------------------------------------
rnaseq.tf.es <- get_tf_ES(
  exp = gene.exp.chr21.se %>% SummarizedExperiment::assay(),
  regulons = regulons.dorothea
)

## ---- message = FALSE, results = "hide"-----------------------------------------------------------
  triplet.regulon <- create_triplet_regulon_based(
    region = dna.met.chr21.se,
    genome = "hg38",  
    motif.search.window.size = 500,
    tf.target = regulons.dorothea,
    max.distance.region.target = 10^6 # 1Mbp
  ) 

## -------------------------------------------------------------------------------------------------
triplet.regulon %>% head

## -------------------------------------------------------------------------------------------------
str(triplet.promoter)
triplet.promoter$distance_region_target_tss %>% range
triplet.promoter %>% head

## ----interaction_model, message = FALSE, results = "hide", eval = TRUE----------------------------
results.interaction.model <- interaction_model(
    triplet = triplet.promoter, 
    dnam = dna.met.chr21.se,
    exp = gene.exp.chr21.se,
    sig.threshold = 0.05,
    fdr = TRUE,
    filter.correlated.tf.exp.dnam = TRUE,
    filter.triplet.by.sig.term = TRUE
)

## -------------------------------------------------------------------------------------------------
# Results for quartile model
results.interaction.model %>% dplyr::select(
  c(1,4,5,grep("quant",colnames(results.interaction.model)))
  ) %>% head

## ----stratified_model, message = FALSE, warning = FALSE, results = "hide", eval = TRUE------------
results.stratified.model <- stratified_model(
    triplet = results.interaction.model,
    dnam = dna.met.chr21.se,
    exp = gene.exp.chr21.se
)

## -------------------------------------------------------------------------------------------------
results.stratified.model %>% head

## ----plot_interaction_model, eval = TRUE, message = FALSE, results = "hide", warning = FALSE------
plots <- plot_interaction_model(
    triplet.results = results.interaction.model[1,], 
    dnam = dna.met.chr21.se, 
    exp = gene.exp.chr21.se
)

## ---- fig.height = 8, fig.width = 13, eval = TRUE, fig.cap = "An example output from MethReg."----
plots

## ----scenarios, fig.cap =  "Scenarios modeled by MethReg.", echo = FALSE, fig.width=13------------
png::readPNG("scenarios.png")  %>% grid::grid.raster()

## ----residuals, results = "hide", eval = FALSE----------------------------------------------------
#  data("gene.exp.chr21.log2")
#  data("clinical")
#  metadata <- clinical[,c("sample_type","gender")]
#  
#  gene.exp.chr21.residuals <- get_residuals(gene.exp.chr21, metadata) %>% as.matrix()

## ---- eval = FALSE--------------------------------------------------------------------------------
#  gene.exp.chr21.residuals[1:5,1:5]

## ---- results = "hide", eval = FALSE--------------------------------------------------------------
#  data("dna.met.chr21")
#  dna.met.chr21 <- make_se_from_dnam_probes(
#    dnam = dna.met.chr21,
#    genome = "hg38",
#    arrayType = "450k",
#    betaToM = TRUE
#  )
#  dna.met.chr21.residuals <- get_residuals(dna.met.chr21, metadata) %>% as.matrix()

## ---- eval = FALSE--------------------------------------------------------------------------------
#  dna.met.chr21.residuals[1:5,1:5]

## ---- message = FALSE, results = "hide", eval = FALSE---------------------------------------------
#  results <- interaction_model(
#      triplet = triplet,
#      dnam = dna.met.chr21.residuals,
#      exp = gene.exp.chr21.residuals
#  )

## -------------------------------------------------------------------------------------------------
regulons.dorothea <- dorothea::dorothea_hs
regulons.dorothea %>% head

## ---- message = FALSE, results = "hide"-----------------------------------------------------------
rnaseq.tf.es <- get_tf_ES(
  exp = gene.exp.chr21.se %>% SummarizedExperiment::assay(),
  regulons = regulons.dorothea
)

## -------------------------------------------------------------------------------------------------
rnaseq.tf.es[1:4,1:4]

## ---- message = FALSE, results = "hide"-----------------------------------------------------------
regulons.dorothea <- dorothea::dorothea_hs
regulons.dorothea$tf <- MethReg:::map_symbol_to_ensg(
  gene.symbol = regulons.dorothea$tf,
  genome = "hg38"
)
regulons.dorothea$target <- MethReg:::map_symbol_to_ensg(
  gene.symbol = regulons.dorothea$target,
  genome = "hg38"
)
split_tibble <- function(tibble, col = 'col') tibble %>% split(., .[, col])
regulons.dorothea.list <- regulons.dorothea %>% na.omit() %>% 
  split_tibble('tf') %>% 
  lapply(function(x){x[[3]]})

## ---- message = FALSE, results = "hide", eval = FALSE---------------------------------------------
#  library(GSVA)
#  rnaseq.tf.es.gsva <- gsva(
#    expr = gene.exp.chr21.se %>% SummarizedExperiment::assay(),
#    gset.idx.list = regulons.dorothea.list,
#    method = "gsva",
#    kcdf = "Gaussian",
#    abs.ranking = TRUE,
#    min.sz = 5,
#    max.sz = Inf,
#    parallel.sz = 1L,
#    mx.diff = TRUE,
#    ssgsea.norm = TRUE,
#    verbose = TRUE
#  )

## ----size = 'tiny'--------------------------------------------------------------------------------
sessionInfo()