---
title: "Hi-C arithmetic with plyinteractions"
output:
BiocStyle::html_document:
self_contained: yes
toc: true
toc_float: true
toc_depth: 4
code_folding: show
date: "`r doc_date()`"
package: "`r pkg_ver('plyinteractions')`"
vignette: >
%\VignetteIndexEntry{HiCarithmetic}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r setup, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
crop = NULL ## Related to https://stat.ethz.ch/pipermail/bioc-devel/2020-April/016656.html
)
options(width = 9999)
```
The `r Biocpkg("plyinteractions")` package facilitates data aggregation, for
up to hundreds of thousands and even millions of
genomic interactions. In this vignette, we explore several use cases
which can arise when exploring Hi-C data stored in `pairs` files.
We will use a real-life `pairs` file provided by the `4DN` Consortium. This
file has been generated from processing Hi-C performed in mouse from brain
cell primary culture during neural development (Bonev et al., Cell 2017). Pairs
have been filtered to only those mapped over `chr13`.
```{r importPairs}
library(tidyverse)
library(plyinteractions)
## Importing it in R
pairs_file <- HiContactsData::HiContactsData('mESCs', 'pairs.gz')
pairs_df <- read.delim(
pairs_file, sep = "\t", header = FALSE, comment.char = "#"
) |>
set_names(c(
"ID", "seqnames1", "start1",
"seqnames2", "start2", "strand1", "strand2"
))
pairs <- as_ginteractions(
pairs_df, end1 = start1, end2 = start2, keep.extra.columns = TRUE
)
pairs
```
# Estimating pairs filtering thresholds
We can first *in silico* digest the mouse genome to obtain the coordinates
of each genomic fragment after digestion by **DpnII and HinfI**.
```{r cutGenome}
## Prepare DpnII/HinfI-digested genomic fragments
library(GenomicRanges)
library(Biostrings)
library(plyranges)
genome <- BSgenome.Mmusculus.UCSC.mm10::BSgenome.Mmusculus.UCSC.mm10
cutter <- DNAStringSet(c("GATC", "GANTC")) ## DpnII/HinfI cutting site
fragments <- BiocParallel::bplapply(BPPARAM = BiocParallel::MulticoreParam(workers = 8),
names(genome), function(.x) {
seq <- genome[[.x]]
mids <- lapply(
cutter,
function(cutsite) {
hits <- matchPattern(cutsite, seq, fixed = "subject")
start(hits) + {end(hits) - start(hits)}
}
) |> unlist() |> sort()
GRanges(seqnames = .x, IRanges(
start = c(1, mids), end = c(mids-1, length(seq))
))
}
) |>
set_names(names(genome)) |>
GRangesList() |>
unlist()
fragments$binID <- seq_along(fragments)
```
We can then use the `annotate()` function from `r Biocpkg("plyinteractions")` to recover,
for each interaction, which restriction enzyme fragment each anchor
overlaps with, and how many restriction enzyme cutting sites are found between
them.
```{r annotatePairs}
## Annotate for each anchor set which genomic fragment it overlaps with
annotated_pairs <- pairs |>
plyinteractions::annotate(fragments, by = "binID") |>
mutate(n_fragments = binID.2 - binID.1, group = paste0(strand1, strand2))
annotated_pairs
```
Next, we can plot the distribution of `strand1` and `strand2` cominations
as a function of the number of restriction enzyme cutting sites between
anchors of each interaction.
```{r getDistrib}
df <- annotated_pairs |>
head(n = 1e6) |>
group_by(strand1, strand2, n_fragments) |>
count() |>
as_tibble() |>
mutate(group = paste0(strand1, strand2)) |>
select(group, n_fragments, n)
ggplot(df, aes(x = n_fragments, y = n, group = group, col = group)) +
geom_line() +
geom_point() +
xlim(c(0, 15)) +
annotation_logticks(sides = 'l') +
theme_bw() +
labs(
x = "Number of restriction sites between anchors",
y = "Number of pairs"
)
```
From this distribution, we can see that `--` and `++` pairs have a decreasing
frequency over increasing numbers of cut sites between
anchors of each interaction. These pairs are unambiguous, as the orientation
of each sequenced end can only come from true cutting and religation event,
(except the set of `--` and `++` pairs which have `0` cut sites between
each anchor, which cannot be explained); all these pairs can be kept.
The over-representation of `+-` pairs at short distance likely represent
uncut fragments subsequently sequenced on each end. The under-representation
of `-+` pairs at short distance likely represent self-religated fragments. We
can estimate a threshold for each of these pairs sets by computing the MAD and
expected , as described in
[Cournac et al., 2012](https://bmcgenomics.biomedcentral.com/articles/10.1186/1471-2164-13-436).
```{r getThresholds}
filters <- df |>
filter(n_fragments <= 50) |>
arrange(n_fragments) |>
group_by(n_fragments) |>
mutate(median = median(n)) |>
ungroup() |>
mutate(MAD = median(abs(n - median))) |>
mutate(withinMAD = abs(n - median) <= MAD / 0.67449) |>
filter(withinMAD) |>
slice_head(by = group, n = 1) |>
select(group, n_fragments) |>
dplyr::rename(threshold = n_fragments)
filters
```
# Filtering pairs using appropriate thresholds
```{r filterPairs}
annotated_pairs <- annotated_pairs |>
mutate(threshold = left_join(as_tibble(mcols(annotated_pairs)), filters)$threshold) |>
mutate(type = case_when(
group %in% c('--', '++') & n_fragments < threshold ~ "excluded",
group == '+-' & n_fragments < threshold ~ "uncut",
group == '-+' & n_fragments < threshold ~ "religated",
.default = "kept"
))
mcols(annotated_pairs) |>
as_tibble() |>
count(type) |>
mutate(n = scales::percent(n/sum(n)))
filtered_pairs <- filter(annotated_pairs, type == 'kept')
```
# Computing distance law from pairs
Another typical step when analyzing Hi-C processed data is the modeling of a so-called
"distance law", (a.k.a "P(s)"), which describes the genomic distance-dependent
contact frequency between pairs of genomic loci from a Hi-C experiment.
We can easily recover the distance between the two anchors of each
interaction (noted *s*) and plot the interaction frequency
(noted *P(s)*) as a function of this genomic distance.
## Plotting distance law: first try
```{r Ps1}
dat <- filtered_pairs |>
mutate(s = abs(end2 - start1)) |>
group_by(s) |>
count(name = "n") |>
as_tibble() |>
mutate(Ps = n/sum(n))
p <- ggplot(dat, aes(x = s, y = Ps)) + geom_line()
p
```
This is not very informative, as the distances span several orders of magnitude
in both dimensions.
## Second try: switching to logarithmic scale
Switching to a `log` scale in `r CRANpkg("ggplot2")` is very easy.
```{r Ps2}
p + scale_x_log10() + scale_y_log10() + annotation_logticks()
```
## Third try: aggregating data before plotting
The previous P(s) plot is precise at the base-pair resolution.
We can aggregate counts by binned distances:
```{r Ps3}
# Calculate distance breaks evenly spaced on a log scale (base 1.1)
x <- 1.1^(1:200-1)
lmc <- coef(lm(c(1,1161443398)~c(x[1], x[200])))
bins_breaks <- unique(round(lmc[2]*x + lmc[1]))
bins_widths <- lead(bins_breaks) - bins_breaks
# Bin distances
dat <- filtered_pairs |>
mutate(s = abs(end2 - start1)) |>
mutate(
binned_s = bins_breaks[as.numeric(cut(s, bins_breaks))],
bin_width = bins_widths[as.numeric(cut(s, bins_breaks))]
) |>
group_by(binned_s, bin_width) |>
count(name = "n") |>
as_tibble() |>
mutate(Ps = n / sum(n) / bin_width)
# Plot results
ggplot(dat, aes(x = binned_s, y = Ps)) + geom_line() +
scale_x_log10() + scale_y_log10() + annotation_logticks()
```
## With some polishing
```{r Ps4}
ggplot(dat, aes(x = binned_s, y = Ps)) +
geom_line() +
scale_x_log10(limits = c(1e3, 1e8)) + ## This changes x axis to log scale
scale_y_log10() + ## This changes y axis to log scale
annotation_logticks() + ## This adds log ticks
labs(
x = "Genomic distance (s)",
y = "P(s)",
title = "Distance-dependent genomic frequency P(s) in mESC (chr. 13)"
) + ## This fixes axes titles
theme_bw() ## This changes default plot theme
```
# Reproducibility
`R` session information:
```{r sessioninfo, echo=FALSE}
## Session info
library("sessioninfo")
options(width = 120)
session_info()
```