TmCalculator (version 1.1.0)

integrate_granges: Integrate a Tm GRanges with multi-omic feature ranges

Description

Combines the output of tm_calculate (a GRanges object with Tm and GC columns) with a second GRanges carrying arbitrary multi-omic metadata (ChIP-seq peaks, ATAC-seq signal, methylation sites, gene annotations, etc.) using one of four positional strategies:

"overlap"

Each tm range is annotated with the aggregated metadata of all feature ranges it directly overlaps.

"nearest"

Each tm range is annotated with the metadata of its single closest feature range, plus an added distance column.

"window"

Each tm range is expanded symmetrically by window_size bp and annotated with aggregated metadata from all features that fall within the expanded window.

"bin"

The genomic space covered by the data is tiled into equal-width bins. Each bin is annotated with the mean tm / GC of overlapping tm ranges and the aggregated feature values - suitable for joint heatmaps and genome-wide correlation analyses.

For strategies "overlap" and "window", when a single Tm range matches multiple features the default behaviour is to summarise: numeric columns are aggregated via agg_fun (default mean), and categorical columns are collapsed to a comma-separated string of unique values.

Usage

integrate_granges(
  gr_tm,
  gr_features,
  strategy = c("overlap", "nearest", "window", "bin"),
  feature_cols = NULL,
  prefix = "",
  window_size = 1000L,
  bin_size = 1e+06,
  agg_fun = mean,
  weight = c("none", "overlap"),
  report_coverage = FALSE,
  min_overlap = 1L,
  ignore_strand = TRUE,
  keep_unmatched = TRUE,
  distance_col = "distance_to_feature"
)

Value

• "overlap", "nearest", "window": A GRanges object with the same ranges as gr_tm (minus unmatched ranges if keep_unmatched = FALSE), with additional metadata columns from gr_features.

• "bin": A new GRanges of genomic bins. Each bin carries Tm_mean, GC_mean (if available), n_tm_ranges, n_features, and one aggregated column per requested feature column.

Arguments

gr_tm

A GRanges object produced by tm_calculate() (or tm_calculate()$gr). Must contain at least a Tm metadata column. A gc column is used automatically when present.

gr_features

A GRanges object with multi-omic feature ranges. All (or a subset of) its metadata columns are transferred / aggregated.

strategy

Character. Integration strategy. One of "overlap" (default), "nearest", "window", or "bin".

feature_cols

Character vector. Names of metadata columns in gr_features to transfer. NULL (default) transfers all metadata columns.

prefix

Character. Prefix prepended to transferred column names to avoid clashes with existing columns in gr_tm. Default: "". Use e.g. "feat_" if there are naming conflicts.

window_size

Integer. Half-width (bp) of the symmetric window added around each Tm range in "window" mode. Default: 1000.

bin_size

Integer. Width (bp) of genomic bins in "bin" mode. Default: 1e6 (1 Mb). Smaller values give finer resolution but sparser coverage.

agg_fun

Function. Applied to numeric feature values when multiple features map to the same Tm range / bin. It is called as agg_fun(values, na.rm = TRUE), so it must accept an na.rm argument; mean, median, sum and max all qualify, whereas function(x) x[1] does not. Ignored for character columns, which are always joined as a comma-separated list of their unique values. Default: mean.

weight

Character. How features are combined within a range. "none" (default) gives every feature that passes the overlap test the same weight, whatever the length of its overlap. "overlap" computes a mean weighted by the number of overlapping base pairs, \(\sum_j w_{ij} v_j / \sum_j w_{ij}\) with \(w_{ij}\) the width of the intersection of range \(i\) and feature \(j\).

The distinction matters wherever a signal changes sharply. A 200 bp window whose first 190 bp are covered at depth 5 and whose last 10 bp are covered at depth 200 has a true mean depth of 14.75; unweighted aggregation returns 102.5, because the 10 bp feature counts as much as the 190 bp one. Raising min_overlap does not fix this, it only reverses the sign of the bias by discarding the short feature entirely.

"overlap" requires agg_fun = mean: an overlap-weighted maximum has no accepted definition, and quietly ignoring agg_fun would return an unweighted value that looks weighted. It does not apply to strategy = "nearest", which performs no aggregation.

report_coverage

Logical. Add a covered_frac column giving the fraction of each range covered by at least one feature, computed after reducing the features so that overlapping ones are not counted twice. A weighted mean normalises by the covered bases, so a value derived from a quarter of a window is indistinguishable from one derived from all of it; this column is what makes the difference visible. Produced by the "overlap" and "window" strategies. Default: FALSE.

min_overlap

Integer. Minimum overlap in base pairs required between a Tm range and a feature range. It is a filter and not a weight: a feature either qualifies or does not, and a qualifying feature counts in full. Applies to the "overlap" strategy only; "window" and "bin" require a single overlapping base pair. Default: 1.

ignore_strand

Logical. If TRUE (default), strand is ignored when finding overlaps / nearest neighbours.

keep_unmatched

Logical. In "overlap" mode only: if TRUE (default) Tm ranges with no overlapping feature are retained with NA in the transferred columns. If FALSE, unmatched Tm ranges are dropped.

distance_col

Character. Name of the distance column added in "nearest" mode. Default: "distance_to_feature".

Aggregating continuous signals

When several features map to the same range, numeric columns are summarised by agg_fun and character columns are joined as their unique values. Which features take part is decided by min_overlap, and how much each one counts is decided by weight. The two are easy to confuse, and the distinction is what determines whether a coverage-like signal is summarised correctly at a boundary.

Let \(R_i\) be the \(i\)-th range of gr_tm, \(F_j\) the \(j\)-th feature, \(x_j\) its value, and \(w_{ij}\) the number of base pairs \(R_i\) and \(F_j\) share. Write \(m\) for min_overlap and \(f\) for agg_fun. The features entering the summary of \(R_i\) are those with \(w_{ij} \ge m\), and

$$v_i = f(\{x_j : w_{ij} \ge m\})$$

with weight = "none", or

$$v_i = \frac{\sum_j w_{ij} x_j}{\sum_j w_{ij}}$$

taken over the same features, with weight = "overlap".

The unweighted form gives a feature that overlaps by one base pair the same influence as one that spans the whole range. This is harmless where a signal is flat and wrong where it steps, which is to say at exon boundaries, peak edges and promoters. Raising min_overlap does not repair it: the threshold is a filter, so a short feature is either counted in full or discarded in full, and the bias changes sign rather than disappearing. The example below shows both failures against a case with a known answer.

The weighted mean normalises by the covered base pairs, not by the width of the range, so a value derived from a quarter of a window is indistinguishable from one derived from all of it. report_coverage = TRUE adds a covered_frac column giving the fraction of each range covered by at least one feature, computed after reduce-ing the features so that overlapping ones are not double counted. Multiply by it to convert a mean over covered bases into a mean over the range.

Weighting is not the default. Enabling it changes numeric output, and existing analyses should stay reproducible unless their author decides otherwise.

Author

Junhui Li

Examples

## Aggregation: a coverage track with a known answer -----------------------
## Two 200 bp windows over a signal that steps sharply inside the first.
##
##   window 1  [  1 .. 200]        window 2  [401 .. 600]
##   depth     [  1 .. 190] = 5
##             [191 .. 400] = 200
##                                           [401 .. 450] = 60
library(GenomicRanges)

win <- GRanges("chr1", IRanges(start = c(1, 401), width = 200),
               Tm = c(70, 72))
cov <- GRanges("chr1", IRanges(start = c(1, 191, 401),
                               end   = c(190, 400, 450)),
               cov = c(5, 200, 60))

## Window 1 truly averages (5 * 190 + 200 * 10) / 200 = 14.75.
integrate_granges(win, cov, strategy = "overlap")$cov
## 102.5  60   the 10 bp feature counts as much as the 190 bp one

integrate_granges(win, cov, strategy = "overlap",
                  weight = "overlap")$cov
## 14.75  60   overlap-weighted mean recovers the true depth

integrate_granges(win, cov, strategy = "overlap", min_overlap = 20L)$cov
## 5      60   the threshold discards the short feature; bias reverses

## Window 2 is only a quarter covered. Weighting cannot show that, because
## it normalises by the covered bases; the coverage column can.
res <- integrate_granges(win, cov, strategy = "overlap",
                         weight = "overlap", report_coverage = TRUE)
res$cov                      # 14.75  60
res$covered_frac             # 1.00   0.25
res$cov * res$covered_frac   # 14.75  15    mean over the whole window

if (FALSE) {
library(GenomicRanges)

# -- Sample data ----------------------------------------------------------
set.seed(42)
gr_tm <- GRanges(
  seqnames = c(rep("chr1", 60), rep("chr2", 30)),
  ranges   = IRanges(
    start = c(sort(sample(1:249e6, 60)),
              sort(sample(1:243e6, 30))),
    width = sample(50:200, 90, replace = TRUE)
  ),
  Tm = runif(90, 55, 85),
  GC = runif(90, 30, 70)
)

gr_features <- GRanges(
  seqnames = c(rep("chr1", 40), rep("chr2", 20)),
  ranges   = IRanges(
    start = c(sort(sample(1:249e6, 40)),
              sort(sample(1:243e6, 20))),
    width = sample(500:5000, 60, replace = TRUE)
  ),
  score      = runif(60, 0, 100),
  peak_type  = sample(c("narrow", "broad"), 60, replace = TRUE),
  signal     = rnorm(60, 5, 2)
)

# Strategy 1: overlap - annotate Tm ranges with overlapping peak features
res_overlap <- integrate_granges(gr_tm, gr_features,
                                  strategy = "overlap")

# Strategy 2: nearest - every Tm range gets its closest peak + distance
res_nearest <- integrate_granges(gr_tm, gr_features,
                                  strategy = "nearest")
head(res_nearest$distance_to_feature)

# Strategy 3: window - 5 kb window around each probe
res_window <- integrate_granges(gr_tm, gr_features,
                                 strategy = "window", window_size = 5000)

# Strategy 4: bin - 500 kb genome bins with mean Tm and aggregated signal
res_bin <- integrate_granges(gr_tm, gr_features,
                              strategy = "bin", bin_size = 5e5)
as.data.frame(res_bin) |> head()

# Use a subset of feature columns and add a prefix
integrate_granges(gr_tm, gr_features,
                  strategy    = "overlap",
                  feature_cols = c("score", "peak_type"),
                  prefix       = "chip_")
}