ELMER.data: Supporting data for the ELMER package

1 Introduction

This document provides an introduction of the ELMER.data, which contains supporting data for ELMER (Yao, L., Shen, H., Laird, P. W., Farnham, P. J., & Berman, B. P. 2015). ELMER is package using DNA methylation to identify enhancers, and correlates enhancer state with expression of nearby genes to identify one or more transcriptional targets. Transcription factor (TF) binding site analysis of enhancers is coupled with expression analysis of all TFs to infer upstream regulators. ELMER.data provide 3 necessary data for ELMER analysis:

1. Probes information: files with DNA methylation platforms metadata retrieved from http://zwdzwd.github.io/InfiniumAnnotation (Zhou, Wanding and Laird, Peter W and Shen, Hui 2016).
2. Probes.motif: motif occurences within \(\pm 250bp\) of probe sites on HM450K/EPIC array aligned against hg19/hg38.
3. Union.enhancer: A comprehensive list of genomic strong enhancers.

1.1 Installing and loading ELMER.data

To install this package, start R and enter

devtools::install_github(repo = "tiagochst/ELMER.data")
library("ELMER.data")
library("GenomicRanges")

2 Contents

2.1 Probes information

Probes information were retrieved from http://zwdzwd.github.io/InfiniumAnnotation (Zhou, Wanding and Laird, Peter W and Shen, Hui 2016).

for(plat in c("450K","EPIC")) {
  for(genome in c("hg38","hg19")) {
    base <- "http://zwdzwd.io/InfiniumAnnotation/current/"
    if(plat == "EPIC") {
      annotation <- paste0(base,"EPIC/EPIC.manifest.rda")
    } else {
      annotation <- paste0(base,"hm450/hm450.manifest.rda")
    }
    if(genome == "hg38") annotation <- gsub(".rda",".hg38.rda", annotation)
    if(!file.exists(basename(annotation))) {
      if(Sys.info()["sysname"] == "Windows") mode <- "wb" else  mode <- "w"
      downloader::download(annotation, basename(annotation), mode = mode)
    }
  }
}

data("EPIC.manifest")
as.data.frame(EPIC.manifest)

data("EPIC.manifest.hg38")
as.data.frame(EPIC.manifest.hg38)

data("hm450.manifest.hg38")
as.data.frame(hm450.manifest.hg38)

data("hm450.manifest")
as.data.frame(hm450.manifest)

2.2 Probes.motif

Probes.motif provides information for motif occurences within\(\pm 250bp\) of probe sites on HM450K/EPIC array. HOMER (Heinz, Sven and Benner, Christopher and Spann, Nathanael and Bertolino, Eric and Lin, Yin C and Laslo, Peter and Cheng, Jason X and Murre, Cornelis and Singh, Harinder and Glass, Christopher K 2010) was used with a p-value < 0.0001 to scan a \(\pm 250bp\) region around each probe on HM450K/EPIC using HOCOMOCO V10 motif position weight matrices (PWMs) which provides transcription factor (TF) binding models for more than 600 human TFs. This data set is used in get.enriched.motif function in ELMER to calculate Odds Ratio of motif enrichments for a given set of probes. This data is storaged in a sparse matrix with wuth 640 columns, there is one matrix for HM450K aligned to hg19, one for HM450K aligned to hg38, one for EPIC aligned to hg19, one for EPIC aligned to hg38. Each row is each probe regions (annotation of the regions used can be found in this repository) and each column is motif from http://hocomoco.autosome.ru/ (Kulakovskiy, Ivan V and Medvedeva, Yulia A and Schaefer, Ulf and Kasianov, Artem S and Vorontsov, Ilya E and Bajic, Vladimir B and Makeev, Vsevolod 2013). The value 1 indicates the occurrence of a motif in a particular probe and 0 means no occurrence.

data("Probes.motif.hg19.450K")
dim(Probes.motif.hg19.450K)

[1] 435391    640

str(Probes.motif.hg19.450K)

Formal class 'ngCMatrix' [package "Matrix"] with 5 slots
  ..@ i       : int [1:38198775] 13 33 40 44 71 101 137 149 161 163 ...
  ..@ p       : int [1:641] 0 19923 54244 87292 122353 161278 212385 323860 464485 601197 ...
  ..@ Dim     : int [1:2] 435391 640
  ..@ Dimnames:List of 2
  .. ..$ : chr [1:435391] "cg07092448" "cg20389653" "cg13452258" "cg07344096" ...
  .. ..$ : chr [1:640] "AHR_HUMAN.H10MO.B" "AIRE_HUMAN.H10MO.C" "ALX1_HUMAN.H10MO.B" "ALX3_HUMAN.H10MO.D" ...
  ..@ factors : list()

data("Probes.motif.hg38.450K")
dim(Probes.motif.hg38.450K)

[1] 435391    640

str(Probes.motif.hg38.450K)

Formal class 'ngCMatrix' [package "Matrix"] with 5 slots
  ..@ i       : int [1:38199819] 3 28 31 36 37 65 66 77 84 85 ...
  ..@ p       : int [1:641] 0 19916 54239 87297 122364 161361 212466 323956 464581 601305 ...
  ..@ Dim     : int [1:2] 435391 640
  ..@ Dimnames:List of 2
  .. ..$ : chr [1:435391] "cg15451548" "cg07212761" "cg24463320" "cg11787186" ...
  .. ..$ : chr [1:640] "AHR_HUMAN.H10MO.B" "AIRE_HUMAN.H10MO.C" "ALX1_HUMAN.H10MO.B" "ALX3_HUMAN.H10MO.D" ...
  ..@ factors : list()

data("Probes.motif.hg19.EPIC")
dim(Probes.motif.hg19.EPIC)

[1] 786296    640

str(Probes.motif.hg19.EPIC)

Formal class 'ngCMatrix' [package "Matrix"] with 5 slots
  ..@ i       : int [1:71168136] 13 40 93 107 110 111 112 134 173 174 ...
  ..@ p       : int [1:641] 0 29135 107569 190350 276086 371235 478173 660941 848126 1056604 ...
  ..@ Dim     : int [1:2] 786296 640
  ..@ Dimnames:List of 2
  .. ..$ : chr [1:786296] "cg18818484" "cg18977839" "cg25234611" "cg11752380" ...
  .. ..$ : chr [1:640] "AHR_HUMAN.H10MO.B" "AIRE_HUMAN.H10MO.C" "ALX1_HUMAN.H10MO.B" "ALX3_HUMAN.H10MO.D" ...
  ..@ factors : list()

data("Probes.motif.hg38.EPIC")
dim(Probes.motif.hg38.EPIC)

[1] 786296    640

str(Probes.motif.hg38.EPIC)

Formal class 'ngCMatrix' [package "Matrix"] with 5 slots
  ..@ i       : int [1:71165854] 44 89 90 109 129 133 137 141 149 169 ...
  ..@ p       : int [1:641] 0 29132 107562 190310 276041 371155 478091 661176 848358 1056846 ...
  ..@ Dim     : int [1:2] 786296 640
  ..@ Dimnames:List of 2
  .. ..$ : chr [1:786296] "cg10817250" "cg06145569" "cg23606775" "cg09283376" ...
  .. ..$ : chr [1:640] "AHR_HUMAN.H10MO.B" "AIRE_HUMAN.H10MO.C" "ALX1_HUMAN.H10MO.B" "ALX3_HUMAN.H10MO.D" ...
  ..@ factors : list()

The following code was used to create the objects:

getInfiniumAnnotation <- function(plat = "450K", genome = "hg38"){
    message("Loading object: ",file)
    newenv <- new.env()
    if(plat == "EPIC" & genome == "hg19") data("EPIC.manifest", package = "ELMER.data",envir=newenv)
    if(plat == "EPIC" & genome == "hg38") data("EPIC.manifest.hg38", package = "ELMER.data",envir=newenv)
    if(plat == "450K" & genome == "hg19") data("hm450.manifest", package = "ELMER.data",envir=newenv)
    if(plat == "450K" & genome == "hg38") data("hm450.manifest.hg38", package = "ELMER.data",envir=newenv)
    annotation <- get(ls(newenv)[1],envir=newenv)   
    return(annotation)  
}
# To find for each probe the know motif we will use HOMER software (http://homer.salk.edu/homer/)
# Step:
# 1 - get DNA methylation probes annotation with the regions
# 2 - Make a bed file from it
# 3 - Execute section: Finding Instance of Specific Motifs from http://homer.salk.edu/homer/ngs/peakMotifs.html to the HOCOMOCO TF motifs
# Also, As HOMER is using more RAM than the available we will split the files in to 100k probes.
# Obs: for each probe we create a winddow of 500 bp (-size 500) around it. This might lead to false positives, but will not have false negatives.
# The false posives will be removed latter with some statistical tests.
TFBS.motif <- "http://hocomoco.autosome.ru/final_bundle/HUMAN/mono/HOCOMOCOv10_HUMAN_mono_homer_format_0.0001.motif"
if(!file.exists(basename(TFBS.motif))) downloader::download(TFBS.motif,basename(TFBS.motif))
for(plat in c("450K","EPIC")){
    for(gen in c("hg19","hg38")){
      
      file <- paste0(plat,gen,".txt")
      print(file)
      if(!file.exists(file)){
        # STEP 1
        gr <- getInfiniumAnnotation(plat = plat,genome =  gen)
        
        # This will remove masked probes. They have poor quality and might be arbitrarily positioned (Wanding Zhou)
        print(table(gr$MASK.general))
        gr <- gr[!gr$MASK.general]
        print(table(gr$MASK.general))
        
        df <- data.frame(seqnames=seqnames(gr),
                         starts=as.integer(start(gr)),
                         ends=end(gr),
                         names=names(gr),
                         scores=c(rep(".", length(gr))),
                         strands=strand(gr))
        step <- 50000 # nb of lines in each file 50K was selected to not explode RAM
        n <- nrow(df)
        for(j in 0:floor(n/step)){
          # STEP 2
          file.aux <- paste0(plat,gen,"_",j,".bed")
          if(!file.exists(gsub(".bed",".txt",file.aux))){
            end <- ifelse(((j + 1) * step) > n, n,((j + 1) * step))
            write.table(df[((j * step) + 1):end,], file = file.aux, col.names = F, quote = F,row.names = F,sep = "\t")
            
            # STEP 3
            cmd <- paste("source ~/.bash_rc; annotatePeaks.pl" ,file.aux, gen, "-m", basename(TFBS.motif), "-size 500 -cpu 12 >", gsub(".bed",".txt",file.aux))
            system(cmd)
          }
        }
        # We will merge the results from each file into one
        peaks <- NULL
        for(j in 0:floor(n/step)){
          aux <-  readr::read_tsv(paste0(plat,gen,"_",j,".txt"))
          colnames(aux)[1] <- "PeakID"
          if(is.null(peaks)) {
            peaks <- aux
          } else {
            peaks <- rbind(peaks, aux)
          }
        }
        readr::write_tsv(peaks,path=file,col_names = TRUE)
        print("DONE!")
        gc()
      }
  }
}

# This code will read the table with the motifs, save it as a sparce matrix
# and save all as a .rda that will be placed in ELMER.data
for(plat in c("450K","EPIC")){
    for(gen in c("hg19","hg38")){
        file <- paste0(plat,gen,".txt")
        motifs <- readr::read_tsv(file)
        # From 1 to 21 we have annotations then we have 640 motifs
        matrix <- Matrix::Matrix(0, nrow = nrow(motifs), ncol = 640,sparse = TRUE)
        colnames(matrix) <- gsub(" Distance From Peak\\(sequence,strand,conservation\\)","",colnames(motifs)[-c(1:21)])
        rownames(matrix) <- motifs$PeakID
        matrix[!is.na(motifs[,-c(1:21)])] <- 1
        matrix <- as(matrix, "nsparseMatrix")
        assign(paste0("Probes.motif.",gen,".",plat),matrix) # For each probe that there is a bind to the motif, we will add as 1 in the matrix. (Are hg38, hg19,450k,epic matrices equal?)
        rm(matrix)
        rm(motifs)
        gc()
    }
}
save(Probes.motif.hg19.450K, file = "Probes.motif.hg19.450K.rda", compress = "xz")
save(Probes.motif.hg38.450K, file = "Probes.motif.hg38.450K.rda", compress = "xz")
save(Probes.motif.hg19.EPIC, file = "Probes.motif.hg19.EPIC.rda", compress = "xz")
save(Probes.motif.hg38.EPIC, file = "Probes.motif.hg38.EPIC.rda", compress = "xz")

data("Probes.motif.hg19.450K")
as.data.frame(as.matrix(Probes.motif.hg19.450K[1:20,1:20]))

3 Session Information

---

sessionInfo()

R version 3.4.0 (2017-04-21)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Sierra 10.12.4

Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats4    parallel  stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] Biostrings_2.44.1   XVector_0.16.0      IRanges_2.10.2      S4Vectors_0.14.3    BiocGenerics_0.22.0 ELMER.data_2.0.0   

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.11               compiler_3.4.0             GenomeInfoDb_1.12.1        base64enc_0.1-3            bitops_1.0-6               tools_3.4.0                zlibbioc_1.22.0            digest_0.6.12              jsonlite_1.5              
[10] evaluate_0.10              memoise_1.1.0              lattice_0.20-35            Matrix_1.2-10              DelayedArray_0.2.6         commonmark_1.2             yaml_2.1.14                GenomeInfoDbData_0.99.0    withr_1.0.2               
[19] stringr_1.2.0              roxygen2_6.0.1             xml2_1.1.1                 knitr_1.16                 devtools_1.13.2            rprojroot_1.2              grid_3.4.0                 Biobase_2.36.2             R6_2.2.1                  
[28] rmarkdown_1.5              magrittr_1.5               backports_1.1.0            htmltools_0.3.6            matrixStats_0.52.2         GenomicRanges_1.28.3       SummarizedExperiment_1.6.3 BiocStyle_2.4.0            stringi_1.1.5             
[37] RCurl_1.95-4.8            

4 References

Heinz, Sven and Benner, Christopher and Spann, Nathanael and Bertolino, Eric and Lin, Yin C and Laslo, Peter and Cheng, Jason X and Murre, Cornelis and Singh, Harinder and Glass, Christopher K. 2010. “Simple Combinations of Lineage-Determining Transcription Factors Prime Cis-Regulatory Elements Required for Macrophage and B Cell Identities.”

Kulakovskiy, Ivan V and Medvedeva, Yulia A and Schaefer, Ulf and Kasianov, Artem S and Vorontsov, Ilya E and Bajic, Vladimir B and Makeev, Vsevolod. 2013. “HOCOMOCO a Comprehensive Collection of Human Transcription Factor Binding Sites Models.”

Yao, L., Shen, H., Laird, P. W., Farnham, P. J., & Berman, B. P. 2015. “Inferring Regulatory Element Landscapes and Transcription Factor Networks from Cancer Methylomes.”

Zhou, Wanding and Laird, Peter W and Shen, Hui. 2016. “Comprehensive Characterization, Annotation and Innovative Use of Infinium DNA Methylation BeadChip Probes.”