getTFbindsites()getGCfreq()getGenomeGC()methylTFRAnnotationMm10 provides the genome annotations that
methylTFR needs to
compute GC-bias-corrected transcription factor (TF) deviation scores on
mm10. It contains three kinds of resources, each with its own accessor:
| Accessor | Returns | Help page |
|---|---|---|
getTFbindsites() |
TF binding sites (a GRangesList) |
?tf_bindsites |
getGCfreq() |
GC frequency tables, one per motif | ?motif_gcfreq |
getGenomeGC() |
GC content in 30 nt windows | ?genomewide_GC |
The GC tables are what make the deviation score bias-corrected. A motif whose binding sites sit in GC-rich sequence overlaps CpG islands more often than one that does not, and CpG islands are hypomethylated. Without a correction, that difference in sequence composition would look like a difference in TF activity.
The resources are several gigabytes in total, so they are hosted on AnnotationHub instead of being shipped with the package. The first call to an accessor downloads the resource; later calls read it from the local AnnotationHub cache.
if (!requireNamespace("BiocManager", quietly = TRUE)) {
install.packages("BiocManager")
}
BiocManager::install("methylTFRAnnotationMm10")
# methylTFR, which uses these annotations
BiocManager::install("methylTFR")Downloading the real resources would take too long for a vignette, so
this vignette works on small synthetic objects that have exactly the
structure of the real ones. Every accessor can read .rds files from a
local directory instead of AnnotationHub when the option
methylTFRAnnotationMm10.datadir is set. We write the synthetic objects to a
temporary directory and point the option at it. With the option unset
(the default), the same calls below would download the real data.
suppressPackageStartupMessages(library(GenomicRanges))
library(methylTFRAnnotationMm10)
mock_dir <- file.path(tempdir(), "methylTFRAnnotationMm10_vignette")
dir.create(mock_dir, showWarnings = FALSE)
# Binding sites: a GRangesList with one GRanges per motif. Each range is
# a motif match extended by 200 bp on either side.
mock_tfbs <- GRangesList(
MA0001.1 = GRanges(
c("chr1:10001-10411", "chr1:52001-52411", "chr2:3001-3411"),
strand = c("+", "-", "+")
),
MA0002.1 = GRanges(c("chr1:20001-20415", "chr3:8001-8415"),
strand = c("-", "+")
)
)
saveRDS(mock_tfbs, file.path(mock_dir, "altius_tf_bindsites.rds"))
# GC frequency tables: one 5-row matrix per motif. Rows are GC
# quintiles, columns are positions of a 30 nt window along the binding
# site, and every column sums to 1. We use 6 positions to keep it short.
set.seed(1)
mock_freq <- function(n_pos) {
m <- matrix(runif(5 * n_pos), nrow = 5)
sweep(m, 2, colSums(m), "/")
}
mock_gcfreq <- list(MA0001.1 = mock_freq(6), MA0002.1 = mock_freq(6))
saveRDS(mock_gcfreq, file.path(mock_dir, "altius_motif_gcfreq.rds"))
# Genome-wide GC: non-overlapping 30 nt windows with their GC fraction
# and GC quintile (1 = lowest, 5 = highest).
mock_gc <- GRanges("chr1", IRanges(start = seq(1, by = 30, length.out = 10),
width = 30))
mock_gc$GC_bias <- round(runif(10, 0.2, 0.7), 2)
breaks <- quantile(mock_gc$GC_bias, probs = seq(0, 1, 0.2))
mock_gc$GC_bin <- findInterval(mock_gc$GC_bias, breaks,
rightmost.closed = TRUE)
metadata(mock_gc) <- list(gc_breaks = breaks, bin_scope = "genome",
step = 30L)
saveRDS(mock_gc, file.path(mock_dir, "genomewide_GC_mm10.rds"))
old_opt <- options(methylTFRAnnotationMm10.datadir = mock_dir)
list.files(mock_dir)## [1] "altius_motif_gcfreq.rds" "altius_tf_bindsites.rds"
## [3] "genomewide_GC_mm10.rds"availableMotifSets() lists the motif sets that can be passed to
getTFbindsites() and getGCfreq():
availableMotifSets()## [1] "altius" "cisbpv2" "jaspar2020"getTFbindsites()getTFbindsites() returns a GRangesList named by motif ID, with one
GRanges of binding sites per motif. Each range covers the motif match
plus 200 bp on either side: the footprint window over which methylTFR
reads methylation. All ranges of a motif therefore have the same width.
tf_bindsites <- getTFbindsites("altius")
tf_bindsites## GRangesList object of length 2:
## $MA0001.1
## GRanges object with 3 ranges and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## [1] chr1 10001-10411 +
## [2] chr1 52001-52411 -
## [3] chr2 3001-3411 +
## -------
## seqinfo: 3 sequences from an unspecified genome; no seqlengths
##
## $MA0002.1
## GRanges object with 2 ranges and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## [1] chr1 20001-20415 -
## [2] chr3 8001-8415 +
## -------
## seqinfo: 3 sequences from an unspecified genome; no seqlengthsThe number of binding sites per motif, and the width of each motif’s window:
lengths(tf_bindsites)## MA0001.1 MA0002.1
## 3 2vapply(width(tf_bindsites), unique, integer(1))## MA0001.1 MA0002.1
## 411 415getGCfreq()getGCfreq() returns a named list with one matrix per motif, using the
same motif IDs as getTFbindsites(). Each matrix has five rows, one per
genome-wide GC quintile (row 1 lowest GC, row 5 highest). Each column is
one position of a 30 nt window sliding along the binding site. An entry
is the fraction of the motif’s binding sites whose window at that
position falls into that GC quintile, so every column sums to one. In the
real data a matrix has several hundred columns.
gcfreqs <- getGCfreq("altius")
names(gcfreqs)## [1] "MA0001.1" "MA0002.1"dim(gcfreqs$MA0001.1)## [1] 5 6round(gcfreqs$MA0001.1, 2)## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.11 0.28 0.09 0.15 0.43 0.19
## [2,] 0.16 0.30 0.08 0.21 0.10 0.01
## [3,] 0.25 0.21 0.31 0.29 0.30 0.19
## [4,] 0.39 0.20 0.17 0.11 0.06 0.44
## [5,] 0.09 0.02 0.35 0.23 0.12 0.17colSums(gcfreqs$MA0001.1)## [1] 1 1 1 1 1 1methylTFR combines these fractions with the average methylation it observes in each GC quintile to get the methylation a motif would show from GC content alone, and reports the deviation from that expectation.
getGenomeGC()getGenomeGC() returns a GRanges of non-overlapping 30 nt windows.
GC_bias is the GC fraction of each window and GC_bin its genome-wide
GC quintile. methylTFR overlaps each methylation call with these windows
to assign it to a GC bin. The package serves a single assembly, so no
argument is needed; getGenomeGC("mm10") is accepted too.
gc_dist <- getGenomeGC()
gc_dist## GRanges object with 10 ranges and 2 metadata columns:
## seqnames ranges strand | GC_bias GC_bin
## <Rle> <IRanges> <Rle> | <numeric> <integer>
## [1] chr1 1-30 * | 0.66 5
## [2] chr1 31-60 * | 0.35 2
## [3] chr1 61-90 * | 0.43 3
## [4] chr1 91-120 * | 0.37 2
## [5] chr1 121-150 * | 0.53 4
## [6] chr1 151-180 * | 0.33 1
## [7] chr1 181-210 * | 0.44 3
## [8] chr1 211-240 * | 0.58 4
## [9] chr1 241-270 * | 0.24 1
## [10] chr1 271-300 * | 0.64 5
## -------
## seqinfo: 1 sequence from an unspecified genome; no seqlengthstable(gc_dist$GC_bin)##
## 1 2 3 4 5
## 2 2 2 2 2The quintile boundaries are stored with the object, so that the genome table and the motif tables always use the same scale:
metadata(gc_dist)$gc_breaks## 0% 20% 40% 60% 80% 100%
## 0.240 0.346 0.406 0.476 0.592 0.660The three objects are passed to methylTFR::run_methyltfr(). With real
data this looks like:
library(methylTFR)
deviations <- run_methyltfr(
"samples.tsv", # sample annotation
"samples_dir", # directory with methylation calls
tf_bindsites = getTFbindsites("jaspar2020"),
gcfreqs = getGCfreq("jaspar2020"),
gc_dist = getGenomeGC(),
filetype = "EPP"
)See the methylTFR vignettes for a complete analysis.
Once the local-directory option is removed, the accessors go to AnnotationHub. The records can also be browsed directly:
library(AnnotationHub)
ah <- AnnotationHub()
query(ah, "methylTFRAnnotationMm10")If you have downloaded the .rds files by hand, set
options(methylTFRAnnotationMm10.datadir = "/path/to/files") (or the
environment variable METHYL_TFRANNOTATION_Mm10_DIR) to read them from disk.
Restore the option and remove the synthetic files:
options(old_opt)
unlink(mock_dir, recursive = TRUE)
getOption("methylTFRAnnotationMm10.datadir")## NULLsessionInfo()## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.5 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.24-bioc/R/lib/libRblas.so
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: America/New_York
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] methylTFRAnnotationMm10_0.99.13 GenomicRanges_1.65.4
## [3] Seqinfo_1.3.2 IRanges_2.47.5
## [5] S4Vectors_0.51.10 BiocGenerics_0.59.12
## [7] generics_0.1.4 BiocStyle_2.41.0
##
## loaded via a namespace (and not attached):
## [1] cli_3.6.6 knitr_1.52 rlang_1.3.0
## [4] xfun_0.62 otel_0.3.0 jsonlite_2.0.0
## [7] htmltools_0.5.9 sass_0.4.10 rmarkdown_2.32
## [10] evaluate_1.0.5 jquerylib_0.1.4 fastmap_1.2.0
## [13] yaml_2.3.12 lifecycle_1.0.5 bookdown_0.48
## [16] BiocManager_1.30.27 compiler_4.6.1 digest_0.6.39
## [19] R6_2.6.1 bslib_0.12.0 tools_4.6.1
## [22] cachem_1.1.0