Matching case study I: CTCF occupancy
Eric S. Davis
05/19/2023
In this vignette we demonstrate generating covariate-matched, null-hypothesis GRanges using the matchRanges() function to test for the occupancy of CCCTC-binding factor (CTCF) at chromatin loop anchors.
Background
One of the fundamental principles of chromatin-looping suggests that most loops are bound at both ends by the CTCF transcription factor (TF). CTCF-bound loops can be formed by loop-extrusion, where the ring-like cohesin complex extrudes chromatin until stopped by bound CTCF. By this mechanism, we expect most loop anchors will be bound by CTCF.
While we could test this hypothesis by simple overlap or permutation testing, these approaches fail to account for non-uniformly distributed covariate genomic features. For example, loop anchors are commonly bound by CTCF and located in open chromatin regions. We can use matchRanges() to test for CTCF occupancy at loop anchors controlling for open chromatin regions.
Here, we generate a set of null-hypothesis GRanges to more rigorously test CTCF occupancy at loop anchors independently from open chromatin regions. We use the hg19_10kb_bins dataset from the nullrangesData package, which contains ranges for every 10Kb bin along the genome with CTCF, DNase, and loop feature annotations from GM12878 (see ?nullrangesData::hg19_10kb_bins).
Matching with matchRanges()
Before we generate our null ranges, let’s take a look at our example dataset:
## GRanges object with 303641 ranges and 5 metadata columns:
## seqnames ranges strand | n_ctcf_sites ctcfSignal
## <Rle> <IRanges> <Rle> | <numeric> <numeric>
## [1] chr1 1-10000 * | 0 0
## [2] chr1 10001-20000 * | 0 0
## [3] chr1 20001-30000 * | 0 0
## [4] chr1 30001-40000 * | 0 0
## [5] chr1 40001-50000 * | 0 0
## ... ... ... ... . ... ...
## [303637] chrX 155230001-155240000 * | 0 0.00000
## [303638] chrX 155240001-155250000 * | 0 0.00000
## [303639] chrX 155250001-155260000 * | 1 4.09522
## [303640] chrX 155260001-155270000 * | 0 0.00000
## [303641] chrX 155270001-155270560 * | 0 0.00000
## n_dnase_sites dnaseSignal looped
## <factor> <numeric> <logical>
## [1] 0 0.00000 FALSE
## [2] 0 5.03572 FALSE
## [3] 0 0.00000 FALSE
## [4] 0 0.00000 FALSE
## [5] 0 0.00000 FALSE
## ... ... ... ...
## [303637] 0 8.42068 FALSE
## [303638] 0 4.08961 FALSE
## [303639] 0 6.00443 FALSE
## [303640] 0 8.07179 FALSE
## [303641] 0 0.00000 FALSE
## -------
## seqinfo: 23 sequences from hg19 genomematchRanges() works by selecting a set of covariate-matched controls from a pool of options based on an input focal set of interest. Here, we define focal as bins that contain a loop anchor, pool as bins that don’t contain a loop anchor, and covar as DNase signal and number of DNase sites per bin:
library(nullranges)
## Match ranges
set.seed(123)
mgr <- matchRanges(focal = bins[bins$looped],
pool = bins[!bins$looped],
covar = ~dnaseSignal + n_dnase_sites)
mgr## MatchedGRanges object with 13979 ranges and 5 metadata columns:
## seqnames ranges strand | n_ctcf_sites ctcfSignal
## <Rle> <IRanges> <Rle> | <numeric> <numeric>
## [1] chr3 107520001-107530000 * | 0 0.00000
## [2] chr5 143180001-143190000 * | 0 0.00000
## [3] chr19 18190001-18200000 * | 0 0.00000
## [4] chr7 30210001-30220000 * | 0 0.00000
## [5] chr2 238530001-238540000 * | 1 5.35384
## ... ... ... ... . ... ...
## [13975] chr14 68060001-68070000 * | 1 6.38626
## [13976] chr2 191700001-191710000 * | 0 0.00000
## [13977] chr9 13220001-13230000 * | 0 0.00000
## [13978] chr4 40300001-40310000 * | 1 4.82893
## [13979] chr16 24670001-24680000 * | 0 0.00000
## n_dnase_sites dnaseSignal looped
## <factor> <numeric> <logical>
## [1] 2 10.20343 FALSE
## [2] 2 10.97785 FALSE
## [3] 3+ 12.26576 FALSE
## [4] 1 11.58243 FALSE
## [5] 3+ 9.43115 FALSE
## ... ... ... ...
## [13975] 3+ 12.50117 FALSE
## [13976] 1 9.95555 FALSE
## [13977] 0 7.17512 FALSE
## [13978] 3+ 14.49135 FALSE
## [13979] 2 11.79441 FALSE
## -------
## seqinfo: 23 sequences from hg19 genomeWhen the focal and pool arguments are GRanges objects, matchRanges() returns a MatchedGRanges object. The MatchedGRanges class extends GRanges, so all of the same operations can be applied:
library(GenomicRanges)
library(plyranges)
library(ggplot2)
## Summarize ctcfSignal by n_ctcf_sites
mgr %>%
group_by(n_ctcf_sites) %>%
summarize(ctcfSignal = mean(ctcfSignal)) %>%
as.data.frame() %>%
ggplot(aes(x = n_ctcf_sites, y = ctcfSignal)) +
geom_line() +
geom_point(shape = 21, stroke = 1, fill = 'white') +
theme_minimal() +
theme(panel.border = element_rect(color = 'black',
fill = NA))
Here, we utilize the plyranges package which provides a set of “tidy” verbs for manipulating genomic ranges for a seamless and integrated genomic analysis workflow.
Assessing quality of matching
We can get a quick summary of the matching quality with overview():
## MatchedGRanges object:
## set N dnaseSignal.mean dnaseSignal.sd n_dnase_sites.0
## focal 13979 10.0 1.9 2341
## matched 13979 10.0 1.9 2327
## pool 289662 7.9 2.7 222164
## unmatched 275683 7.8 2.7 219837
## n_dnase_sites.1 n_dnase_sites.2 n_dnase_sites.3+ ps.mean ps.sd
## 4829 2353 4456 0.130 0.072
## 5089 2519 4044 0.130 0.071
## 34826 13627 19045 0.042 0.061
## 29737 11108 15001 0.037 0.057
## --------
## focal - matched:
## dnaseSignal.mean dnaseSignal.sd n_dnase_sites.0 n_dnase_sites.1
## 0.038 -0.021 14 -260
## n_dnase_sites.2 n_dnase_sites.3+ ps.mean ps.sd
## -170 410 3e-04 0.00042For continuous covariates (such as dnaseSignal), overview() shows the mean and standard deviation between each matched set. For categorical covariates, such as n_dnase_sites, overview() reports the number of observations per category and matched set. The bottom section shows the mean and s.d (or n, for factors) difference between focal and matched sets.
overview() also summarizes the propensity scores for each set to give a quick idea of overall matching quality.
Visualizing matching results
Let’s visualize overall matching quality by plotting propensity scores for the focal, pool, and matched sets:
From this plot, it is clear that the matched set is much closer to the focal set than the pool set.
We can ensure that covariate distributions have been matched appropriately by using the covariates() function to extract matched covariates along with patchwork and plotCovarite to visualize all distributions: