When performing statistical analysis on any set of genomic ranges it is often important to compare focal sets to null sets that are carefully matched for possible covariates that may influence the analysis. To address this need, the nullranges
package implements matchRanges()
, an efficient and convenient tool for selecting a covariate-matched set of null hypothesis ranges from a pool of background ranges within the Bioconductor framework.
In this vignette, we provide an overview of matchRanges()
and its associated functions. We start with a simulated example generated with the utility function makeExampleMatchedDataSet()
. We also provide an overview of the class struture and a guide for choosing among the supported matching methods. To see matchedRanges()
used in real biological examples, visit the Case study I: CTCF occupancy, and Case study II: CTCF orientation vignettes.
For a description of the method, see Davis et al. (2022).
matchRanges
references four sets of data: focal
, pool
, matched
and unmatched
. The focal
set contains the outcome of interest (Y=1
) while the pool
set contains all other observations (Y=0
). matchRanges
generates the matched
set, which is a subset of the pool
that is matched for provided covariates (i.e. covar
) but does not contain the outcome of interest (i.e Y=0
). Finally, the unmatched
set contains the remaining unselected elements from the pool
. The diagram below depicts the relationships between the four sets.
matchRanges
uses propensity scores to perform subset selection on the pool
set such that the resulting matched
set contains similar distributions of covariates to that of the focal
set. A propensity score is the conditional probability of assigning an element (in our case, a genomic range) to a particular outcome (Y
) given a set of covariates. Propensity scores are estimated using a logistic regression model where the outcome Y=1
for focal
and Y=0
for pool
, over the provided covariates covar
. The resulting propensity scores are used to select matches using one of three available matching options: “nearest”, “rejection”, or “stratified” with or without replacement. For more information see the section on [Choosing the method parameter] below.
matchRanges()
We will use a simulated data set to demonstrate matching across covarying features:
## GRanges object with 10500 ranges and 3 metadata columns:
## seqnames ranges strand | feature1 feature2 feature3
## <Rle> <IRanges> <Rle> | <logical> <numeric> <character>
## [1] chr1 1-100 * | TRUE 2.87905 c
## [2] chr1 2-101 * | TRUE 3.53965 c
## [3] chr1 3-102 * | TRUE 7.11742 c
## [4] chr1 4-103 * | TRUE 4.14102 a
## [5] chr1 5-104 * | TRUE 4.25858 c
## ... ... ... ... . ... ... ...
## [10496] chr1 10496-10595 * | FALSE 1.23578 b
## [10497] chr1 10497-10596 * | FALSE 1.69671 a
## [10498] chr1 10498-10597 * | FALSE 6.11140 a
## [10499] chr1 10499-10598 * | FALSE 2.21657 d
## [10500] chr1 10500-10599 * | FALSE 5.33003 b
## -------
## seqinfo: 1 sequence from an unspecified genome; no seqlengths
Our simulated dataset has 3 features: logical feature1
, numeric feature2
, and character/factor feature3
. We can use matchRanges()
to compare ranges where feature1
is TRUE
to ranges where feature1
is FALSE
, matched by feature2
and/or feature3
:
set.seed(123)
mgr <- matchRanges(focal = x[x$feature1],
pool = x[!x$feature1],
covar = ~feature2 + feature3)
mgr
## MatchedGRanges object with 500 ranges and 3 metadata columns:
## seqnames ranges strand | feature1 feature2 feature3
## <Rle> <IRanges> <Rle> | <logical> <numeric> <character>
## [1] chr1 4373-4472 * | FALSE 8.959578 d
## [2] chr1 9740-9839 * | FALSE 0.959336 e
## [3] chr1 7755-7854 * | FALSE 2.107003 c
## [4] chr1 8266-8365 * | FALSE 6.231860 d
## [5] chr1 4298-4397 * | FALSE 6.955316 c
## ... ... ... ... . ... ... ...
## [496] chr1 2443-2542 * | FALSE 1.12276 b
## [497] chr1 2455-2554 * | FALSE 3.38518 c
## [498] chr1 1285-1384 * | FALSE 1.58546 c
## [499] chr1 10137-10236 * | FALSE 9.39272 c
## [500] chr1 6119-6218 * | FALSE 10.22412 c
## -------
## seqinfo: 1 sequence from an unspecified genome; no seqlengths
The resulting MatchedGRanges
object is a set of null hypothesis ranges selected from our pool
of options that is the same length as our input focal
ranges and matched for covar
features 2 and 3. These matched ranges print and behave just as normal GRanges
would:
## MatchedGRanges object with 500 ranges and 3 metadata columns:
## seqnames ranges strand | feature1 feature2 feature3
## <Rle> <IRanges> <Rle> | <logical> <numeric> <character>
## [1] chr1 511-610 * | FALSE 5.545186 c
## [2] chr1 513-612 * | FALSE 2.221684 b
## [3] chr1 534-633 * | FALSE 1.563458 b
## [4] chr1 565-664 * | FALSE 0.932659 c
## [5] chr1 577-676 * | FALSE 3.256908 c
## ... ... ... ... . ... ... ...
## [496] chr1 10377-10476 * | FALSE 0.795032 c
## [497] chr1 10380-10479 * | FALSE 0.977984 b
## [498] chr1 10409-10508 * | FALSE 3.662119 c
## [499] chr1 10455-10554 * | FALSE 6.815473 c
## [500] chr1 10483-10582 * | FALSE 3.724147 c
## -------
## seqinfo: 1 sequence from an unspecified genome; no seqlengths
We can change the type
argument of makeExampleMatchedDataSet
to input data.frames, data.tables, DataFrames, GRanges and GInteractions objects - all of which work as inputs for matchRanges
. These produce either MatchedDataFrame
, MatchedGRanges
, or MatchedGInteractions
objects. For more information about the Matched
class structure and available methods, see the Class structure section below or the help documentation for each class, ?MatchedDataFrame
, ?MatchedGRanges
, or ?MatchedGInteractions
.
matchRanges()
uses propensity scores to select matches using one of three available matching options: “nearest”, “rejection”, or “stratified” with or without replacement. For more information see the section on [Choosing the method parameter] below.
We can assess the quality of Matched
classes with overview()
, plotCovariate()
, and plotPropensity()
. overview()
provides a quick assessment of overall matching quality by reporting the mean and standard deviation for covariates and propensity scores of the focal, pool, matched, and unmatched sets. For factor, character, or logical covariates (e.g. categorical covariates) the N per set (frequency) is returned. It also reports the mean difference in focal-matched sets:
## MatchedGRanges object:
## set N feature2.mean feature2.sd feature3.a feature3.b feature3.c
## focal 500 4.1 1.9 66 157 206
## matched 500 4.5 2.7 34 160 234
## pool 10000 6.0 3.4 4248 3121 1117
## unmatched 9500 6.1 3.5 4214 2961 883
## feature3.d feature3.e ps.mean ps.sd
## 49 22 0.100 0.076
## 53 19 0.110 0.078
## 992 522 0.045 0.051
## 939 503 0.041 0.047
## --------
## focal - matched:
## feature2.mean feature2.sd feature3.a feature3.b feature3.c feature3.d
## -0.42 -0.84 32 -3 -28 -4
## feature3.e ps.mean ps.sd
## 3 -0.0057 -0.0019
Visualizing propensity scores can show how well sets were matched overall:
The distributions of features can be visualized in each set with plotCovariate()
:
Since these functions return ggplots, patchwork
can be used to visualize all covariates like this:
library(patchwork)
plots <- lapply(covariates(mgr), plotCovariate, x=mgr, sets = c('f', 'm', 'p'))
Reduce('/', plots)
By default, continuous features are plotted as density line plots while categorical features are plotted as stacked bar plots. All sets are also shown by default. Defaults can be overridden by setting the type
and sets
arguments.
Custom plots can be made by extracting data from the Matched
object:
## id feature2 feature3 ps set
## 1: 1 2.879049 c 0.21095908 focal
## 2: 1 3.539645 c 0.19210984 focal
## 3: 1 7.117417 c 0.11193396 focal
## 4: 1 4.141017 a 0.01771986 focal
## 5: 1 4.258575 c 0.17308581 focal
## ---
## 20496: 0 1.235781 b 0.08945367 unmatched
## 20497: 0 1.696712 a 0.02707977 unmatched
## 20498: 0 6.111404 a 0.01255772 unmatched
## 20499: 0 2.216575 d 0.07578989 unmatched
## 20500: 0 5.330029 b 0.04535856 unmatched
Attributes of the Matched
object can be extracted with the following accessor functions:
## [1] "feature2" "feature3"
## [1] "rejection"
## [1] FALSE
Each set can also be extracted with the following accessor functions:
## [1] "GRanges object with 500 ranges and 3 metadata columns"
## [1] "GRanges object with 10000 ranges and 3 metadata columns"
## [1] "GRanges object with 500 ranges and 3 metadata columns"
## [1] "GRanges object with 9500 ranges and 3 metadata columns"
The indices()
function can be used to find the original indices for each set. For example, indices(x, set="matched")
will supply the indices from the pool
set that corresponds to the matched
set. In fact, matched(x)
is a convenient wrapper around pool(x)[indices(x, set='matched')
:
## [1] TRUE
method
parameterThere are currently 3 available methods for selecting a matched set:
Nearest-neighbor matching with replacement
Rejection sampling with/without replacement
Stratified sampling with/without replacement
Currently, nearest-neighbor matching without replacement is not implemented, but stratified sampling without replacement is a suitable substitute.
Attempts to find the nearest neighbor for each range by using a rolling-join (as implemented in the data.table
package) between focal
and pool
propensity scores.
set.seed(123)
mgr <- matchRanges(focal = x[x$feature1],
pool = x[!x$feature1],
covar = ~feature2 + feature3,
method = 'nearest',
replace = TRUE)
nn <- overview(mgr)
plotPropensity(mgr)
This method is best if you have a very large dataset because it is usually the fastest matching method. However, because sampling is done with replacement the user should be careful to assess the number of duplicate ranges pulled. This can be done using the indices()
function:
## Total number of duplicated indices
length(which(duplicated(indices(mgr))))
sum(table(indices(mgr)) > 1) # used more than once
sum(table(indices(mgr)) > 2) # used more than twice
sum(table(indices(mgr)) > 3) # used more than thrice
## [1] 59
## [1] 51
## [1] 8
## [1] 0
Duplicate ranges can be pulled since this method selects the closest matching propensity-score in the focal set to each range in the pool set. It is important to inspect the duplicates when using this method particularly when there are very few well-matching options to select from in your pool set to ensure your matched set has a diverse set of ranges.
Nearest neighbor matching without replacement is not currently supported due to its computational complexity. However, stratified sampling without replacement is an acceptable alternative.
Uses a probability-based approach to select options in the pool
that distributionally match the focal
set based on propensity scores. The rejection sampling method first generates kernal-density estimates for both the focal and pool sets. Then a scale factor is determined by finding the point at which the difference in focal and pool densitites is maximized. This scale factor is then applied such that the pool distribution covers the focal distribution at all points. Random sampling is then conducted, with probability of accepting a pool range into the matched set given by the ratio between the height of the density and the scaled (covering) density. If method
or replace
is not supplied, the default values are rejection sampling without replacement.
set.seed(123)
mgr <- matchRanges(focal = x[x$feature1],
pool = x[!x$feature1],
covar = ~feature2 + feature3,
method = 'rejection',
replace = FALSE)
rs <- overview(mgr)
plotPropensity(mgr)
Rejection sampling is the fastest available matching method for sampling without replacement. Therefore, it is ideal to use on large datasets when sampling without replacement is important. However, this method can be unstable, particularly when the pool set is not much larger than the focal set. In those cases, the best method to use is stratified sampling.
Performs iterative sampling on increasingly large bins of data. focal
and pool
propensity scores are binned by their value with high granularity, options are randomly selected (with or without replacement) within each bin and subsequently removed from the pool of available options. This procedure is repeated, decreasing the number of bins (and increasing bin size) until the number of selected matches is equal to the focal set. While matches are being found in each bin the bins stay small. However, as the number of bins with no matches increases the algorithm expands bin size faster, which maintains matching quality while decreasing run-time.
set.seed(123)
mgr <- matchRanges(focal = x[x$feature1],
pool = x[!x$feature1],
covar = ~feature2 + feature3,
method = 'stratified',
replace = FALSE)
ss <- overview(mgr)
plotPropensity(mgr)
For very large data sets, users might notice a slight increase in run time compared to the other methods. Stratified sampling tends to work very well for discrete data, and often produces the best matches even on continuous data:
## Extract difference in propensity scores
## between focal and matched sets
fmps <- sapply(c(nn, rs, ss), `[[`, "quality")
c('nearest', 'rejection', 'stratified')[which.min(fmps)]
## [1] "stratified"
Since matchRanges()
automatically constructs the relevant classes, this section is not essential for using any of the nullranges
package functionality. Instead, this section serves as a guide for developers who wish to extend these classes or those more interested in S4 implementation details.
matchRanges()
acts as a constructor, combining a Matched
superclass - which contains the matching results - with either a DataFrame
(data.frame
/data.table
), GRanges
, or GInteractions
superclass. This results in the MatchedDataFrame
, MatchedGRanges
, or MatchedGInteractions
subclasses.
Internally, each Matched
subclass uses a “delegate” object of the same type to assign its slots. The delegate object used is the matched
set. Therefore, the resulting Matched*
object behaves as a combination of both its superclasses - with access to methods from both.
For example, using matchRanges()
on GRanges
objects assigns a GRanges
delegate object which is used to populate GRanges-specific slots. This results in a MatchedGRanges
object, with access to both Matched
functions (e.g. plotCovariate
) as well as normal GRanges
methods (e.g.s seqnames
, resize
, etc…).
## R version 4.2.1 (2022-06-23)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.5 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.16-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.16-bioc/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] grid stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] patchwork_1.1.2 plyranges_1.18.0
## [3] nullrangesData_1.3.0 ExperimentHub_2.6.0
## [5] AnnotationHub_3.6.0 BiocFileCache_2.6.0
## [7] dbplyr_2.2.1 ggplot2_3.3.6
## [9] plotgardener_1.4.0 nullranges_1.4.0
## [11] InteractionSet_1.26.0 SummarizedExperiment_1.28.0
## [13] Biobase_2.58.0 MatrixGenerics_1.10.0
## [15] matrixStats_0.62.0 GenomicRanges_1.50.0
## [17] GenomeInfoDb_1.34.0 IRanges_2.32.0
## [19] S4Vectors_0.36.0 BiocGenerics_0.44.0
##
## loaded via a namespace (and not attached):
## [1] colorspace_2.0-3 rjson_0.2.21
## [3] ellipsis_0.3.2 ggridges_0.5.4
## [5] mclust_6.0.0 XVector_0.38.0
## [7] farver_2.1.1 bit64_4.0.5
## [9] interactiveDisplayBase_1.36.0 AnnotationDbi_1.60.0
## [11] fansi_1.0.3 mvtnorm_1.1-3
## [13] codetools_0.2-18 cachem_1.0.6
## [15] knitr_1.40 jsonlite_1.8.3
## [17] speedglm_0.3-4 Rsamtools_2.14.0
## [19] png_0.1-7 shiny_1.7.3
## [21] BiocManager_1.30.19 compiler_4.2.1
## [23] httr_1.4.4 assertthat_0.2.1
## [25] Matrix_1.5-1 fastmap_1.1.0
## [27] cli_3.4.1 later_1.3.0
## [29] prettyunits_1.1.1 htmltools_0.5.3
## [31] tools_4.2.1 gtable_0.3.1
## [33] glue_1.6.2 GenomeInfoDbData_1.2.9
## [35] dplyr_1.0.10 rappdirs_0.3.3
## [37] Rcpp_1.0.9 jquerylib_0.1.4
## [39] vctrs_0.5.0 Biostrings_2.66.0
## [41] strawr_0.0.9 rtracklayer_1.58.0
## [43] xfun_0.34 stringr_1.4.1
## [45] mime_0.12 lifecycle_1.0.3
## [47] restfulr_0.0.15 XML_3.99-0.12
## [49] zlibbioc_1.44.0 MASS_7.3-58.1
## [51] scales_1.2.1 hms_1.1.2
## [53] promises_1.2.0.1 parallel_4.2.1
## [55] RColorBrewer_1.1-3 yaml_2.3.6
## [57] curl_4.3.3 memoise_2.0.1
## [59] yulab.utils_0.0.5 sass_0.4.2
## [61] stringi_1.7.8 RSQLite_2.2.18
## [63] BiocVersion_3.16.0 highr_0.9
## [65] BiocIO_1.8.0 filelock_1.0.2
## [67] BiocParallel_1.32.0 rlang_1.0.6
## [69] pkgconfig_2.0.3 bitops_1.0-7
## [71] pracma_2.4.2 evaluate_0.17
## [73] lattice_0.20-45 purrr_0.3.5
## [75] GenomicAlignments_1.34.0 ks_1.13.5
## [77] labeling_0.4.2 bit_4.0.4
## [79] tidyselect_1.2.0 magrittr_2.0.3
## [81] R6_2.5.1 generics_0.1.3
## [83] DelayedArray_0.24.0 DBI_1.1.3
## [85] pillar_1.8.1 withr_2.5.0
## [87] KEGGREST_1.38.0 RCurl_1.98-1.9
## [89] tibble_3.1.8 crayon_1.5.2
## [91] KernSmooth_2.23-20 utf8_1.2.2
## [93] rmarkdown_2.17 progress_1.2.2
## [95] data.table_1.14.4 blob_1.2.3
## [97] digest_0.6.30 xtable_1.8-4
## [99] httpuv_1.6.6 gridGraphics_0.5-1
## [101] munsell_0.5.0 ggplotify_0.1.0
## [103] bslib_0.4.0
Davis, Eric S., Wancen Mu, Stuart Lee, Mikhail G. Dozmorov, Michael I. Love, and Douglas H. Phanstiel. 2022. “MatchRanges: Generating Null Hypothesis Genomic Ranges via Covariate-Matched Sampling.” bioRxiv. https://doi.org/10.1101/2022.08.05.502985.