---
title: "methylTFRAnnotationMm10"
author: "Irem B. Gunduz"
date: "`r Sys.Date()`"
output: BiocStyle::html_document
vignette: >
  %\VignetteIndexEntry{methylTFRAnnotationMm10}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

# Overview

`methylTFRAnnotationMm10` provides the genome annotations that
[methylTFR](https://github.com/EpigenomeInformatics/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.

# Installation

```r
if (!requireNamespace("BiocManager", quietly = TRUE)) {
    install.packages("BiocManager")
}
BiocManager::install("methylTFRAnnotationMm10")

# methylTFR, which uses these annotations
BiocManager::install("methylTFR")
```

# Example data used in this vignette

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.

```{r mock-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)
```

# Available motif sets

`availableMotifSets()` lists the motif sets that can be passed to
`getTFbindsites()` and `getGCfreq()`:

```{r motif-sets}
availableMotifSets()
```

# Binding sites: `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.

```{r tfbs}
tf_bindsites <- getTFbindsites("altius")
tf_bindsites
```

The number of binding sites per motif, and the width of each motif's
window:

```{r tfbs-summary}
lengths(tf_bindsites)
vapply(width(tf_bindsites), unique, integer(1))
```

# GC frequency tables: `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.

```{r gcfreq}
gcfreqs <- getGCfreq("altius")
names(gcfreqs)
dim(gcfreqs$MA0001.1)
round(gcfreqs$MA0001.1, 2)
colSums(gcfreqs$MA0001.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.

# Genome-wide GC distribution: `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.

```{r genome-gc}
gc_dist <- getGenomeGC()
gc_dist
table(gc_dist$GC_bin)
```

The quintile boundaries are stored with the object, so that the genome
table and the motif tables always use the same scale:

```{r genome-gc-breaks}
metadata(gc_dist)$gc_breaks
```

# Using the annotations with methylTFR

The three objects are passed to `methylTFR::run_methyltfr()`. With real
data this looks like:

```r
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.

# Reading the real data from AnnotationHub

Once the local-directory option is removed, the accessors go to
AnnotationHub. The records can also be browsed directly:

```r
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.

# Cleaning up

Restore the option and remove the synthetic files:

```{r cleanup}
options(old_opt)
unlink(mock_dir, recursive = TRUE)
getOption("methylTFRAnnotationMm10.datadir")
```

# Session info

```{r session-info}
sessionInfo()
```