The 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.

library(plyinteractions)
library(tidyverse)
#> ── Attaching core tidyverse packages ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 2.0.0 ──
#> ✔ forcats   1.0.1     ✔ readr     2.1.6
#> ✔ ggplot2   4.0.1     ✔ stringr   1.6.0
#> ✔ lubridate 1.9.4     ✔ tibble    3.3.1
#> ✔ purrr     1.2.1     ✔ tidyr     1.3.2
#> ── Conflicts ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
#> ✖ lubridate::%within%()   masks IRanges::%within%()
#> ✖ ggplot2::Position()     masks BiocGenerics::Position(), base::Position()
#> ✖ ggplot2::annotate()     masks plyinteractions::annotate()
#> ✖ plyranges::between()    masks dplyr::between()
#> ✖ dplyr::collapse()       masks IRanges::collapse()
#> ✖ dplyr::combine()        masks Biobase::combine(), BiocGenerics::combine()
#> ✖ dplyr::count()          masks matrixStats::count()
#> ✖ dplyr::desc()           masks IRanges::desc()
#> ✖ tidyr::expand()         masks S4Vectors::expand()
#> ✖ dplyr::filter()         masks stats::filter()
#> ✖ dplyr::first()          masks S4Vectors::first()
#> ✖ dplyr::lag()            masks stats::lag()
#> ✖ plyranges::n()          masks dplyr::n()
#> ✖ plyranges::n_distinct() masks dplyr::n_distinct()
#> ✖ purrr::reduce()         masks GenomicRanges::reduce(), IRanges::reduce()
#> ✖ dplyr::rename()         masks S4Vectors::rename()
#> ✖ lubridate::second()     masks S4Vectors::second()
#> ✖ lubridate::second<-()   masks S4Vectors::second<-()
#> ✖ dplyr::slice()          masks IRanges::slice()
#> ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

## Importing it in R
pairs_file <- HiContactsData::HiContactsData('mESCs', 'pairs.gz')
#> see ?HiContactsData and browseVignettes('HiContactsData') for documentation
#> loading from cache
pairs_df <- read.delim(
    pairs_file, sep = "\t", header = FALSE, comment.char = "#", nrows = 1e6
) |> 
    set_names(c(
        "ID", "seqnames1", "start1", 
        "seqnames2", "start2", "strand1", "strand2"
    )) 
pairs <- as_ginteractions(
    pairs_df, end1 = start1, end2 = start2, keep.extra.columns = TRUE
)
pairs
#> GInteractions object with 1000000 interactions and 1 metadata column:
#>             seqnames1   ranges1 strand1     seqnames2   ranges2 strand2 |                  ID
#>                 <Rle> <IRanges>   <Rle>         <Rle> <IRanges>   <Rle> |         <character>
#>         [1]     chr13  17057558       + ---     chr13  17176616       - |       SRR5339749.58
#>         [2]     chr13  68759440       - ---     chr13 113578864       - |      SRR5339749.105
#>         [3]     chr13  47940999       + ---     chr13  48134537       + |      SRR5339749.169
#>         [4]     chr13  80638451       + ---     chr13  80638826       - |      SRR5339749.170
#>         [5]     chr13   4362498       - ---     chr13  96982617       + |      SRR5339749.249
#>         ...       ...       ...     ... ...       ...       ...     ... .                 ...
#>    [999996]     chr13  17722638       - ---     chr13  20561010       - | SRR5339749.45907723
#>    [999997]     chr13  91988792       + ---     chr13  92333140       - | SRR5339749.45907730
#>    [999998]     chr13  98284294       - ---     chr13  98531510       + | SRR5339749.45907742
#>    [999999]     chr13  36053113       + ---     chr13  36076689       - | SRR5339749.45907855
#>   [1000000]     chr13  25585835       + ---     chr13  26197381       - | SRR5339749.45907876
#>   -------
#>   regions: 1936743 ranges and 0 metadata columns
#>   seqinfo: 1 sequence from an unspecified genome; no seqlengths

1 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.

## Prepare DpnII/HinfI-digested genomic fragments
library(Biostrings)
#> Loading required package: XVector
#> 
#> Attaching package: 'XVector'
#> The following object is masked from 'package:purrr':
#> 
#>     compact
#> 
#> Attaching package: 'Biostrings'
#> The following object is masked from 'package:base':
#> 
#>     strsplit
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 = 2), 
    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 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.

## 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
#> GInteractions object with 1000000 interactions and 5 metadata columns:
#>             seqnames1   ranges1 strand1     seqnames2   ranges2 strand2 |                  ID   binID.1   binID.2 n_fragments       group
#>                 <Rle> <IRanges>   <Rle>         <Rle> <IRanges>   <Rle> |         <character> <integer> <integer>   <integer> <character>
#>         [1]     chr13  17057558       + ---     chr13  17176616       - |       SRR5339749.58   9591352   9592012         660          +-
#>         [2]     chr13  68759440       - ---     chr13 113578864       - |      SRR5339749.105   9880169  10124404      244235          --
#>         [3]     chr13  47940999       + ---     chr13  48134537       + |      SRR5339749.169   9762274   9763393        1119          ++
#>         [4]     chr13  80638451       + ---     chr13  80638826       - |      SRR5339749.170   9946878   9946878           0          +-
#>         [5]     chr13   4362498       - ---     chr13  96982617       + |      SRR5339749.249   9521271  10034142      512871          -+
#>         ...       ...       ...     ... ...       ...       ...     ... .                 ...       ...       ...         ...         ...
#>    [999996]     chr13  17722638       - ---     chr13  20561010       - | SRR5339749.45907723   9594967   9610226       15259          --
#>    [999997]     chr13  91988792       + ---     chr13  92333140       - | SRR5339749.45907730  10007449  10009172        1723          +-
#>    [999998]     chr13  98284294       - ---     chr13  98531510       + | SRR5339749.45907742  10041410  10042771        1361          -+
#>    [999999]     chr13  36053113       + ---     chr13  36076689       - | SRR5339749.45907855   9696787   9696912         125          +-
#>   [1000000]     chr13  25585835       + ---     chr13  26197381       - | SRR5339749.45907876   9638511   9641687        3176          +-
#>   -------
#>   regions: 1936743 ranges and 0 metadata columns
#>   seqinfo: 1 sequence from an unspecified genome; no seqlengths

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.

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"
    )
#> Warning: Removed 267493 rows containing missing values or values outside the scale range (`geom_line()`).
#> Warning: Removed 267493 rows containing missing values or values outside the scale range (`geom_point()`).

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.

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) |> 
    rename(threshold = n_fragments)
filters
#> # A tibble: 4 × 2
#>   group threshold
#>   <chr>     <int>
#> 1 ++            1
#> 2 --            1
#> 3 -+            8
#> 4 +-           10

2 Filtering pairs using appropriate thresholds

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"
    ))
#> Joining with `by = join_by(group)`
mcols(annotated_pairs) |>
    as_tibble() |> 
    count(type) |> 
    mutate(n = scales::percent(n/sum(n)))
#> # A tibble: 4 × 2
#>   type      n     
#>   <chr>     <chr> 
#> 1 excluded  1.08% 
#> 2 kept      78.70%
#> 3 religated 0.40% 
#> 4 uncut     19.82%

filtered_pairs <- filter(annotated_pairs, type == 'kept')

3 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.

3.1 Plotting distance law: first try

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.

3.2 Second try: switching to logarithmic scale

Switching to a log scale in ggplot2 is very easy.

p + scale_x_log10() + scale_y_log10() + annotation_logticks()

3.3 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:

# 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()

3.4 With some polishing

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
#> Warning: Removed 41 rows containing missing values or values outside the scale range (`geom_line()`).

4 Reproducibility

R session information:

#> ─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────
#>  setting  value
#>  version  R Under development (unstable) (2025-12-22 r89219)
#>  os       Ubuntu 24.04.3 LTS
#>  system   x86_64, linux-gnu
#>  ui       X11
#>  language (EN)
#>  collate  C
#>  ctype    en_US.UTF-8
#>  tz       America/New_York
#>  date     2026-01-11
#>  pandoc   2.7.3 @ /usr/bin/ (via rmarkdown)
#>  quarto   1.8.25 @ /usr/local/bin/quarto
#> 
#> ─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────
#>  package                      * version   date (UTC) lib source
#>  abind                          1.4-8     2024-09-12 [2] CRAN (R 4.6.0)
#>  AnnotationDbi                  1.73.0    2026-01-11 [2] Bioconductor 3.23 (R 4.6.0)
#>  AnnotationHub                * 4.1.0     2026-01-11 [2] Bioconductor 3.23 (R 4.6.0)
#>  backports                      1.5.0     2024-05-23 [2] CRAN (R 4.6.0)
#>  bibtex                         0.5.1     2023-01-26 [2] CRAN (R 4.6.0)
#>  Biobase                      * 2.71.0    2026-01-11 [2] Bioconductor 3.23 (R 4.6.0)
#>  BiocFileCache                * 3.1.0     2026-01-11 [2] Bioconductor 3.23 (R 4.6.0)
#>  BiocGenerics                 * 0.57.0    2026-01-11 [2] Bioconductor 3.23 (R 4.6.0)
#>  BiocIO                         1.21.0    2026-01-11 [2] Bioconductor 3.23 (R 4.6.0)
#>  BiocManager                    1.30.27   2025-11-14 [2] CRAN (R 4.6.0)
#>  BiocParallel                   1.45.0    2026-01-05 [2] Bioconductor 3.23 (R 4.6.0)
#>  BiocStyle                    * 2.39.0    2026-01-11 [2] Bioconductor 3.23 (R 4.6.0)
#>  BiocVersion                    3.23.1    2026-01-11 [2] Bioconductor 3.23 (R 4.6.0)
#>  Biostrings                   * 2.79.4    2026-01-11 [2] Bioconductor 3.23 (R 4.6.0)
#>  bit                            4.6.0     2025-03-06 [2] CRAN (R 4.6.0)
#>  bit64                          4.6.0-1   2025-01-16 [2] CRAN (R 4.6.0)
#>  bitops                         1.0-9     2024-10-03 [2] CRAN (R 4.6.0)
#>  blob                           1.2.4     2023-03-17 [2] CRAN (R 4.6.0)
#>  bookdown                       0.46      2025-12-05 [2] CRAN (R 4.6.0)
#>  BSgenome                       1.79.1    2026-01-11 [2] Bioconductor 3.23 (R 4.6.0)
#>  BSgenome.Mmusculus.UCSC.mm10   1.4.3     2025-12-23 [2] Bioconductor
#>  bslib                          0.9.0     2025-01-30 [2] CRAN (R 4.6.0)
#>  cachem                         1.1.0     2024-05-16 [2] CRAN (R 4.6.0)
#>  cigarillo                      1.1.0     2026-01-11 [2] Bioconductor 3.23 (R 4.6.0)
#>  cli                            3.6.5     2025-04-23 [2] CRAN (R 4.6.0)
#>  codetools                      0.2-20    2024-03-31 [3] CRAN (R 4.6.0)
#>  crayon                         1.5.3     2024-06-20 [2] CRAN (R 4.6.0)
#>  curl                           7.0.0     2025-08-19 [2] CRAN (R 4.6.0)
#>  DBI                            1.2.3     2024-06-02 [2] CRAN (R 4.6.0)
#>  dbplyr                       * 2.5.1     2025-09-10 [2] CRAN (R 4.6.0)
#>  DelayedArray                   0.37.0    2026-01-11 [2] Bioconductor 3.23 (R 4.6.0)
#>  dichromat                      2.0-0.1   2022-05-02 [2] CRAN (R 4.6.0)
#>  digest                         0.6.39    2025-11-19 [2] CRAN (R 4.6.0)
#>  dplyr                        * 1.1.4     2023-11-17 [2] CRAN (R 4.6.0)
#>  evaluate                       1.0.5     2025-08-27 [2] CRAN (R 4.6.0)
#>  ExperimentHub                * 3.1.0     2026-01-11 [2] Bioconductor 3.23 (R 4.6.0)
#>  farver                         2.1.2     2024-05-13 [2] CRAN (R 4.6.0)
#>  fastmap                        1.2.0     2024-05-15 [2] CRAN (R 4.6.0)
#>  filelock                       1.0.3     2023-12-11 [2] CRAN (R 4.6.0)
#>  forcats                      * 1.0.1     2025-09-25 [2] CRAN (R 4.6.0)
#>  generics                     * 0.1.4     2025-05-09 [2] CRAN (R 4.6.0)
#>  GenomicAlignments              1.47.0    2026-01-11 [2] Bioconductor 3.23 (R 4.6.0)
#>  GenomicRanges                * 1.63.1    2026-01-11 [2] Bioconductor 3.23 (R 4.6.0)
#>  ggplot2                      * 4.0.1     2025-11-14 [2] CRAN (R 4.6.0)
#>  glue                           1.8.0     2024-09-30 [2] CRAN (R 4.6.0)
#>  gtable                         0.3.6     2024-10-25 [2] CRAN (R 4.6.0)
#>  HiContactsData               * 1.13.0    2026-01-08 [2] Bioconductor 3.23 (R 4.6.0)
#>  hms                            1.1.4     2025-10-17 [2] CRAN (R 4.6.0)
#>  htmltools                      0.5.9     2025-12-04 [2] CRAN (R 4.6.0)
#>  httr                           1.4.7     2023-08-15 [2] CRAN (R 4.6.0)
#>  httr2                          1.2.2     2025-12-08 [2] CRAN (R 4.6.0)
#>  InteractionSet               * 1.39.0    2026-01-11 [2] Bioconductor 3.23 (R 4.6.0)
#>  IRanges                      * 2.45.0    2026-01-11 [2] Bioconductor 3.23 (R 4.6.0)
#>  jquerylib                      0.1.4     2021-04-26 [2] CRAN (R 4.6.0)
#>  jsonlite                       2.0.0     2025-03-27 [2] CRAN (R 4.6.0)
#>  KEGGREST                       1.51.1    2026-01-11 [2] Bioconductor 3.23 (R 4.6.0)
#>  knitr                          1.51      2025-12-20 [2] CRAN (R 4.6.0)
#>  labeling                       0.4.3     2023-08-29 [2] CRAN (R 4.6.0)
#>  lattice                        0.22-7    2025-04-02 [3] CRAN (R 4.6.0)
#>  lifecycle                      1.0.5     2026-01-08 [2] CRAN (R 4.6.0)
#>  lubridate                    * 1.9.4     2024-12-08 [2] CRAN (R 4.6.0)
#>  magrittr                       2.0.4     2025-09-12 [2] CRAN (R 4.6.0)
#>  Matrix                         1.7-4     2025-08-28 [3] CRAN (R 4.6.0)
#>  MatrixGenerics               * 1.23.0    2026-01-11 [2] Bioconductor 3.23 (R 4.6.0)
#>  matrixStats                  * 1.5.0     2025-01-07 [2] CRAN (R 4.6.0)
#>  memoise                        2.0.1     2021-11-26 [2] CRAN (R 4.6.0)
#>  otel                           0.2.0     2025-08-29 [2] CRAN (R 4.6.0)
#>  pillar                         1.11.1    2025-09-17 [2] CRAN (R 4.6.0)
#>  pkgconfig                      2.0.3     2019-09-22 [2] CRAN (R 4.6.0)
#>  plyinteractions              * 1.9.2     2026-01-11 [1] Bioconductor 3.23 (R 4.6.0)
#>  plyr                           1.8.9     2023-10-02 [2] CRAN (R 4.6.0)
#>  plyranges                    * 1.31.1    2026-01-11 [2] Bioconductor 3.23 (R 4.6.0)
#>  png                            0.1-8     2022-11-29 [2] CRAN (R 4.6.0)
#>  purrr                        * 1.2.1     2026-01-09 [2] CRAN (R 4.6.0)
#>  R6                             2.6.1     2025-02-15 [2] CRAN (R 4.6.0)
#>  rappdirs                       0.3.3     2021-01-31 [2] CRAN (R 4.6.0)
#>  RColorBrewer                   1.1-3     2022-04-03 [2] CRAN (R 4.6.0)
#>  Rcpp                           1.1.1     2026-01-10 [2] CRAN (R 4.6.0)
#>  RCurl                          1.98-1.17 2025-03-22 [2] CRAN (R 4.6.0)
#>  readr                        * 2.1.6     2025-11-14 [2] CRAN (R 4.6.0)
#>  RefManageR                   * 1.4.0     2022-09-30 [2] CRAN (R 4.6.0)
#>  restfulr                       0.0.16    2025-06-27 [2] CRAN (R 4.6.0)
#>  rjson                          0.2.23    2024-09-16 [2] CRAN (R 4.6.0)
#>  rlang                          1.1.7     2026-01-09 [2] CRAN (R 4.6.0)
#>  rmarkdown                      2.30      2025-09-28 [2] CRAN (R 4.6.0)
#>  Rsamtools                      2.27.0    2026-01-11 [2] Bioconductor 3.23 (R 4.6.0)
#>  RSQLite                        2.4.5     2025-11-30 [2] CRAN (R 4.6.0)
#>  rtracklayer                    1.71.3    2026-01-11 [2] Bioconductor 3.23 (R 4.6.0)
#>  S4Arrays                       1.11.1    2026-01-11 [2] Bioconductor 3.23 (R 4.6.0)
#>  S4Vectors                    * 0.49.0    2026-01-11 [2] Bioconductor 3.23 (R 4.6.0)
#>  S7                             0.2.1     2025-11-14 [2] CRAN (R 4.6.0)
#>  sass                           0.4.10    2025-04-11 [2] CRAN (R 4.6.0)
#>  scales                         1.4.0     2025-04-24 [2] CRAN (R 4.6.0)
#>  Seqinfo                      * 1.1.0     2026-01-11 [2] Bioconductor 3.23 (R 4.6.0)
#>  sessioninfo                  * 1.2.3     2025-02-05 [2] CRAN (R 4.6.0)
#>  SparseArray                    1.11.10   2026-01-11 [2] Bioconductor 3.23 (R 4.6.0)
#>  stringi                        1.8.7     2025-03-27 [2] CRAN (R 4.6.0)
#>  stringr                      * 1.6.0     2025-11-04 [2] CRAN (R 4.6.0)
#>  SummarizedExperiment         * 1.41.0    2026-01-11 [2] Bioconductor 3.23 (R 4.6.0)
#>  tibble                       * 3.3.1     2026-01-11 [2] CRAN (R 4.6.0)
#>  tidyr                        * 1.3.2     2025-12-19 [2] CRAN (R 4.6.0)
#>  tidyselect                     1.2.1     2024-03-11 [2] CRAN (R 4.6.0)
#>  tidyverse                    * 2.0.0     2023-02-22 [2] CRAN (R 4.6.0)
#>  timechange                     0.3.0     2024-01-18 [2] CRAN (R 4.6.0)
#>  tzdb                           0.5.0     2025-03-15 [2] CRAN (R 4.6.0)
#>  utf8                           1.2.6     2025-06-08 [2] CRAN (R 4.6.0)
#>  vctrs                          0.6.5     2023-12-01 [2] CRAN (R 4.6.0)
#>  withr                          3.0.2     2024-10-28 [2] CRAN (R 4.6.0)
#>  xfun                           0.55      2025-12-16 [2] CRAN (R 4.6.0)
#>  XML                            3.99-0.20 2025-11-08 [2] CRAN (R 4.6.0)
#>  xml2                           1.5.1     2025-12-01 [2] CRAN (R 4.6.0)
#>  XVector                      * 0.51.0    2026-01-11 [2] Bioconductor 3.23 (R 4.6.0)
#>  yaml                           2.3.12    2025-12-10 [2] CRAN (R 4.6.0)
#> 
#>  [1] /tmp/RtmpQBUECW/Rinst13ae4754b5501b
#>  [2] /home/biocbuild/bbs-3.23-bioc/R/site-library
#>  [3] /home/biocbuild/bbs-3.23-bioc/R/library
#>  * ── Packages attached to the search path.
#> 
#> ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────