feat: add_pairdist and pair_granges functions

DESCRIPTION

@@ -52,6 +52,7 @@ RoxygenNote: 7.3.1
 Collate: 
     'AllGenerics.R'
     'AllClasses.R'
+    'add_pairdist.R'
     'anchor.R'
     'annotate.R'
     'arrange.R'
@@ -79,6 +80,7 @@ Collate:
     'methods-PinnedGInteractions.R'
     'methods-show.R'
     'mutate.R'
+    'pair-granges.R'
     'pin.R'
     'plyinteractions.R'
     'reexports-dplyr.R'

NAMESPACE

@@ -72,6 +72,7 @@ S3method(tally,GroupedGInteractions)
 S3method(tbl_vars,GInteractions)
 S3method(unanchor,AnchoredPinnedGInteractions)
 S3method(ungroup,GroupedGInteractions)
+export(add_pairdist)
 export(anchor)
 export(anchor_3p)
 export(anchor_5p)
@@ -110,6 +111,7 @@ export(join_overlap_left)
 export(join_overlap_left_directed)
 export(mutate)
 export(n_groups)
+export(pair_granges)
 export(pin)
 export(pin_anchors1)
 export(pin_anchors2)

R/add_pairdist.R

@@ -0,0 +1,47 @@
+#' Appends distance between interaction anchors 
+#' 
+#' Appends distance between interaction anchors, using 
+#' `InteractionSet::pairdist`
+#'
+#' @param x The query GInteractions
+#' @param type A character string specifying the type of distance to compute. Can take values of "mid", "gap", "span", "diag" or "intra".
+#' @param colname name of column to hold pair distance values
+#' 
+#' @return The GInteractions with an additional column containing the 
+#'   distance between each pair of anchors.
+#' 
+#' @rdname add-pairdist
+#' 
+#' @export
+#'
+#' @examples
+#' gi <- read.table(text = "
+#' chr1 100 200 chr1 5000 5100 bedpe_example1 30 + -
+#' chr1 1000 5000 chr2 3000 3800 bedpe_example2 100 + -",
+#' col.names = c(
+#'   "seqnames1", "start1", "end1", 
+#'   "seqnames2", "start2", "end2", "name", "score", "strand1", "strand2")
+#' ) |> as_ginteractions()
+#' 
+#' add_pairdist(gi)
+#' @export
+
+add_pairdist <- function(x, type = 'mid', colname = 'pairdist') {
+
+    if (colname %in% names(GenomicRanges::mcols(x))){
+        stop(paste0(colname, " already exists in destination metadata"))
+    }
+
+    if (is.null(GenomicRanges::mcols(x))){
+        # handle IRanges NULL adding X column of NA's
+        meta <- S4Vectors::DataFrame("distance" = NA_integer_)
+        names(meta) <- colname
+        GenomicRanges::mcols(x) <- meta
+    } else {
+        GenomicRanges::mcols(x)[[colname]] <-  NA_integer_
+    }
+
+    GenomicRanges::mcols(x)[[colname]] <- InteractionSet::pairdist(x, type)
+
+    x
+}

R/pair-granges.R

@@ -0,0 +1,32 @@
+#' Pairwise combination of a GRanges object
+#' 
+#' Create a GInteractions object from a GRanges object, 
+#' containing all possible entry pairs
+#'
+#' @param x A GRanges object
+#' 
+#' @return A GInteractions object
+#' 
+#' @rdname pair-granges
+#' 
+#' @export
+#'
+#' @examples
+#' gr <- read.table(text = "
+#' chr1 100 200
+#' chr1 5000 5100
+#' chr1 1000 5000
+#' chr2 3000 3800",
+#' col.names = c(
+#'   "seqnames", "start", "end" 
+#' )) |> plyranges::as_granges()
+#' 
+#' pair_granges(gr)
+#' @export
+
+pair_granges <- function(x) {
+
+    combs <- combn(length(x), 2)
+    InteractionSet::GInteractions(combs[1,], combs[2,], gr)
+
+}

_pkgdown.yml

@@ -28,6 +28,7 @@ reference:
   contents: 
   - as_ginteractions
   - ginteractions-getters
+  - pair-granges
 - title: "`dplyr` core verbs"
   contents: 
   - starts_with("dplyr-")
@@ -43,9 +44,10 @@ reference:
   - ginteractions-count-overlaps
   - ginteractions-filter-overlaps
   - ginteractions-join-overlap-left
-- title: "Annotating GInteractions"
+- title: "Enriching GInteractions"
   contents: 
   - ginteractions-annotate
+  - add-pairdist
 - title: "Pinning GInteractions"
   contents: 
   - pin

tests/testthat/test-enrich.R

@@ -0,0 +1,84 @@
+test_that("enrich functions work", {
+
+    gi <- read.table(text = "
+        chr1 100 200 chr1 5000 5100 bedpe_example1 30 + -
+        chr1 1000 5000 chr2 3000 3800 bedpe_example2 100 + -",
+        col.names = c(
+            "seqnames1", "start1", "end1", 
+            "seqnames2", "start2", "end2", "name", "score", "strand1", "strand2"
+        )
+    ) |> as_ginteractions()
+
+    add_pairdist(gi) |> expect_identical(
+        new("GInteractions", anchor1 = 1:2, anchor2 = 3:4, regions = new("GRanges", 
+            seqnames = new("Rle", values = structure(1:2, levels = c("chr1", 
+            "chr2"), class = "factor"), lengths = c(3L, 1L), elementMetadata = NULL, 
+                metadata = list()), ranges = new("IRanges", start = c(100L, 
+            1000L, 5000L, 3000L), width = c(101L, 4001L, 101L, 801L), 
+                NAMES = NULL, elementType = "ANY", elementMetadata = NULL, 
+                metadata = list()), strand = new("Rle", values = structure(1:2, levels = c("+", 
+            "-", "*"), class = "factor"), lengths = c(2L, 2L), elementMetadata = NULL, 
+                metadata = list()), seqinfo = new("Seqinfo", seqnames = c("chr1", 
+            "chr2"), seqlengths = c(NA_integer_, NA_integer_), is_circular = c(NA, 
+            NA), genome = c(NA_character_, NA_character_)), elementMetadata = new("DFrame", 
+                rownames = NULL, nrows = 4L, elementType = "ANY", elementMetadata = NULL, 
+                metadata = list(), listData = structure(list(), names = character(0))), 
+            elementType = "ANY", metadata = list()), NAMES = NULL, elementMetadata = new("DFrame", 
+            rownames = NULL, nrows = 2L, elementType = "ANY", elementMetadata = NULL, 
+            metadata = list(), listData = list(name = c("bedpe_example1", 
+            "bedpe_example2"), score = c(30L, 100L), pairdist = c(4900L, 
+            NA))), metadata = list())
+        )
+
+    add_pairdist(gi, colname = 's') |> expect_identical(
+        new("GInteractions", anchor1 = 1:2, anchor2 = 3:4, regions = new("GRanges", 
+            seqnames = new("Rle", values = structure(1:2, levels = c("chr1", 
+            "chr2"), class = "factor"), lengths = c(3L, 1L), elementMetadata = NULL, 
+                metadata = list()), ranges = new("IRanges", start = c(100L, 
+            1000L, 5000L, 3000L), width = c(101L, 4001L, 101L, 801L), 
+                NAMES = NULL, elementType = "ANY", elementMetadata = NULL, 
+                metadata = list()), strand = new("Rle", values = structure(1:2, levels = c("+", 
+            "-", "*"), class = "factor"), lengths = c(2L, 2L), elementMetadata = NULL, 
+                metadata = list()), seqinfo = new("Seqinfo", seqnames = c("chr1", 
+            "chr2"), seqlengths = c(NA_integer_, NA_integer_), is_circular = c(NA, 
+            NA), genome = c(NA_character_, NA_character_)), elementMetadata = new("DFrame", 
+                rownames = NULL, nrows = 4L, elementType = "ANY", elementMetadata = NULL, 
+                metadata = list(), listData = structure(list(), names = character(0))), 
+            elementType = "ANY", metadata = list()), NAMES = NULL, elementMetadata = new("DFrame", 
+            rownames = NULL, nrows = 2L, elementType = "ANY", elementMetadata = NULL, 
+            metadata = list(), listData = list(name = c("bedpe_example1", 
+            "bedpe_example2"), score = c(30L, 100L), s = c(4900L, 
+            NA))), metadata = list())
+        )
+
+    gr <- read.table(text = "
+        chr1 100 200
+        chr1 5000 5100
+        chr1 1000 5000
+        chr2 3000 3800",
+        col.names = c(
+        "seqnames", "start", "end" 
+    )) |> plyranges::as_granges()
+
+    pair_granges(gr) |> expect_identical(
+        new("GInteractions", anchor1 = c(1L, 1L, 1L, 3L, 3L, 2L), anchor2 = c(3L, 
+            2L, 4L, 2L, 4L, 4L), regions = new("GRanges", seqnames = new("Rle", 
+            values = structure(1:2, levels = c("chr1", "chr2"), class = "factor"), 
+            lengths = c(3L, 1L), elementMetadata = NULL, metadata = list()), 
+            ranges = new("IRanges", start = c(100L, 1000L, 5000L, 3000L
+            ), width = c(101L, 4001L, 101L, 801L), NAMES = NULL, elementType = "ANY", 
+                elementMetadata = NULL, metadata = list()), strand = new("Rle", 
+                values = structure(3L, levels = c("+", "-", "*"), class = "factor"), 
+                lengths = 4L, elementMetadata = NULL, metadata = list()), 
+            seqinfo = new("Seqinfo", seqnames = c("chr1", "chr2"), seqlengths = c(NA_integer_, 
+            NA_integer_), is_circular = c(NA, NA), genome = c(NA_character_, 
+            NA_character_)), elementMetadata = new("DFrame", rownames = NULL, 
+                nrows = 4L, elementType = "ANY", elementMetadata = NULL, 
+                metadata = list(), listData = structure(list(), names = character(0))), 
+            elementType = "ANY", metadata = list()), NAMES = NULL, elementMetadata = new("DFrame", 
+            rownames = NULL, nrows = 6L, elementType = "ANY", elementMetadata = NULL, 
+            metadata = list(), listData = structure(list(), names = character(0))), 
+            metadata = list()
+        )
+    )
+})