[r] Fix iterative peak merging for isolated peaks (#338)

* [r] Fix iterative peak merging for isolated peaks

Author: koefoeden

r/R/atac_utils.R

@@ -296,25 +296,30 @@ merge_peaks_iterative <- function(peaks) {
   peaks <- tibble::as_tibble(peaks) %>%
     dplyr::distinct(chr, start, end, .keep_all = TRUE)
 
-  overlaps <- range_overlaps(peaks, peaks)
-  # Initialize keeper set with any non-overlapping peaks
-  keeper_sets <- list(setdiff(seq_len(nrow(peaks)), overlaps$from))
-
-  # Maintain invariant: overlaps contains only overlap pairs where we have not
-  #    permanently kept or discarded either of the elements in the pair
-  overlaps <- dplyr::filter(overlaps, from < to)
-  while (nrow(overlaps) > 0) {
-    # Add peaks with no higher-ranked overlap to the keeper set
-    keeper_set <- setdiff(overlaps$from, overlaps$to)
+  # Only keep one direction for each overlap. With peaks ordered by priority,
+  # lower row numbers have higher priority.
+  overlaps <- range_overlaps(peaks, peaks) %>%
+    dplyr::filter(from < to)
+
+  # Track peak indices which are neither conclusively excluded or included from a keeper set
+  remaining <- seq_len(nrow(peaks))
+  keeper_sets <- list()
+  while (length(remaining) > 0) {
+
+    # Add peaks with no higher-ranked overlap to the keeper set. This includes
+    # peaks that are isolated from the start or become isolated after discards.
+    keeper_set <- setdiff(remaining, overlaps$to)
     keeper_sets <- c(keeper_sets, list(keeper_set))
+
     # Mark peaks that overlap directly with the keeper set
-    discard_set <- overlaps$to[overlaps$from %in% keeper_set] %>% unique()
+    discard_set <- overlaps$to[overlaps$from %in% keeper_set] %>%
+      unique()
+    remove_set <- c(keeper_set, discard_set)
+    remaining <- setdiff(remaining, remove_set)
+
     # Discard all overlap information on the keeper and discard sets
     overlaps <- overlaps %>%
-      dplyr::filter(
-        !(to %in% discard_set), !(from %in% discard_set),
-        !(from %in% keeper_set)
-      )
+      dplyr::filter(!(from %in% remove_set), !(to %in% remove_set))
   }
   return(peaks[sort(unlist(keeper_sets)), ])
 }

r/tests/testthat/test-atac_utils.R

@@ -84,6 +84,13 @@ test_that("tile_ranges works", {
 })
 
 test_that("merge_peaks works", {
+  peaks0 <- tibble::tibble(
+    chr = "chr1",
+    start = as.integer(c(1, 5)),
+    end = start + 2L
+  )
+  expect_identical(merge_peaks_iterative(peaks0), peaks0)
+
   peaks1 <- tibble::tibble(
     chr = "chr1",
     start = as.integer(1:10),
@@ -96,6 +103,18 @@ test_that("merge_peaks works", {
   )
   expect_identical(merge_peaks_iterative(peaks1), res1)
 
+  peaks1_chain <- tibble::tibble(
+    chr = "chr1",
+    start = as.integer(1:3),
+    end = start + 2L
+  )
+  res1_chain <- tibble::tibble(
+    chr = "chr1",
+    start = as.integer(c(1, 3)),
+    end = start + 2L
+  )
+  expect_identical(merge_peaks_iterative(peaks1_chain), res1_chain)
+
   peaks2 <- tibble::tibble(
     chr = "chr1",
     start = as.integer(c(1, 22, 11, 22, 2, 3, 13, 1, 4, 21, 12)),
@@ -740,4 +759,4 @@ test_that("qc_scATAC handles genes near the start of a chromosome", {
 
   qc <- qc_scATAC(fragments, genes, blacklist)
   expect_identical(sort(qc$cellName), sort(cellNames(fragments)))
-})
+})