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.
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():
## [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.
## 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 seqlengths
The number of binding sites per motif, and the width of each motif’s window:
## MA0001.1 MA0002.1
## 3 2
## MA0001.1 MA0002.1
## 411 415
getGCfreq()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.
## [1] "MA0001.1" "MA0002.1"
## [1] 5 6
## [,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.17
## [1] 1 1 1 1 1 1
methylTFR 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.
## 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 seqlengths
##
## 1 2 3 4 5
## 2 2 2 2 2
The quintile boundaries are stored with the object, so that the genome table and the motif tables always use the same scale:
## 0% 20% 40% 60% 80% 100%
## 0.240 0.346 0.406 0.476 0.592 0.660
The 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:
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:
## NULL
## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 26.04.1 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so; LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [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: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] methylTFRAnnotationMm10_0.99.12 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.61 otel_0.2.0 jsonlite_2.0.0
## [7] buildtools_1.0.0 htmltools_0.5.9 maketools_1.3.2
## [10] sys_3.4.3 sass_0.4.10 rmarkdown_2.32
## [13] evaluate_1.0.5 jquerylib_0.1.4 fastmap_1.2.0
## [16] yaml_2.3.12 lifecycle_1.0.5 BiocManager_1.30.27
## [19] compiler_4.6.1 digest_0.6.39 R6_2.6.1
## [22] bslib_0.12.0 tools_4.6.1 cachem_1.1.0