# Case study I: CTCF occupancy

#### Eric S. Davis

#### 11/07/2021

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:

    library(nullrangesData)

    ## Load example data
    bins <- hg19_10kb_bins()

    bins

    ## 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 genome

`matchRanges()` 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]     chr8   16360001-16370000      * |            0    0.00000
    ##       [2]     chr1 212260001-212270000      * |            0    0.00000
    ##       [3]     chrX     9570001-9580000      * |            0    0.00000
    ##       [4]     chr9   96840001-96850000      * |            0    0.00000
    ##       [5]    chr16   31480001-31490000      * |            1    4.82341
    ##       ...      ...                 ...    ... .          ...        ...
    ##   [13975]    chr12 105060001-105070000      * |            1    4.29979
    ##   [13976]     chr4 103530001-103540000      * |            1    5.67681
    ##   [13977]    chr14   65860001-65870000      * |            0    0.00000
    ##   [13978]     chr2   96940001-96950000      * |            0    0.00000
    ##   [13979]    chr19   38910001-38920000      * |            0    0.00000
    ##           n_dnase_sites dnaseSignal    looped
    ##                <factor>   <numeric> <logical>
    ##       [1]            0      7.40552     FALSE
    ##       [2]            1      9.72456     FALSE
    ##       [3]            3+    10.87185     FALSE
    ##       [4]            1      8.86143     FALSE
    ##       [5]            3+    10.65095     FALSE
    ##       ...           ...         ...       ...
    ##   [13975]            3+    13.90788     FALSE
    ##   [13976]            1      9.89450     FALSE
    ##   [13977]            1      9.83681     FALSE
    ##   [13978]            0      9.66206     FALSE
    ##   [13979]            2     12.13214     FALSE
    ##   -------
    ##   seqinfo: 23 sequences from hg19 genome

When 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()`:

    overview(mgr)

    ## 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            2297
    ##       pool 289662              7.9            2.7          222164
    ##  unmatched 275683              7.8            2.7          219867
    ##  n_dnase_sites.1 n_dnase_sites.2 n_dnase_sites.3+ ps.mean ps.sd
    ##             4829            2353             4456   0.130 0.072
    ##             5187            2487             4008   0.130 0.071
    ##            34826           13627            19045   0.042 0.061
    ##            29639           11140            15037   0.037 0.057
    ## --------
    ## focal - matched:
    ##  dnaseSignal.mean dnaseSignal.sd n_dnase_sites.0 n_dnase_sites.1
    ##             0.037        -0.0061              44            -360
    ##  n_dnase_sites.2 n_dnase_sites.3+ ps.mean   ps.sd
    ##             -130              450 0.00025 0.00081

For 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:

    plotPropensity(mgr, sets = c('f', 'p', 'm'), type = 'ridges')

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:

    library(patchwork)
    plots <- lapply(covariates(mgr), plotCovariate, x=mgr, sets = c('f', 'm', 'p'))
    Reduce('/', plots)

## Compare CTCF sites

Using our matched ranges, we can compare CTCF occupancy in bins that 1) contain a loop anchor (i.e. looped), 2) don’t contain a loop anchor (i.e. unlooped), or 3) don’t contain a loop anchor, but are also matched for the strength and number of DNase sites (i.e. matched). In this case, we calculate CTCF occupancy as the percent of bins that contain CTCF among our 3 sets by using the `focal()` and `pool()` accessor functions:

    ## Percent of bins with CTCF
    g1 <- (sum(focal(mgr)$n_ctcf_sites >= 1) / length(focal(mgr))) * 100
    g2 <- (sum(pool(mgr)$n_ctcf_sites >= 1) / length(pool(mgr))) * 100
    g3 <- (sum(mgr$n_ctcf_sites >= 1) / length(mgr)) * 100

    ## Visualize
    barplot(height = c(g1, g2, g3),
            names = c('looped\n(focal)', 'unlooped\n(pool)', 'unlooped\n(matched)'),
            ylab = "CTCF occupied bins (%)",
            col = c('#1F78B4', '#33A02C', '#A6CEE3'),
            main = 'CTCF occupancy',
            border = NA,
            las = 1)

# Session information

    sessionInfo()

    ## R version 4.1.1 (2021-08-10)
    ## Platform: x86_64-w64-mingw32/x64 (64-bit)
    ## Running under: Windows Server x64 (build 17763)
    ##
    ## Matrix products: default
    ##
    ## locale:
    ## [1] LC_COLLATE=C
    ## [2] LC_CTYPE=English_United States.1252
    ## [3] LC_MONETARY=English_United States.1252
    ## [4] LC_NUMERIC=C
    ## [5] LC_TIME=English_United States.1252
    ##
    ## attached base packages:
    ## [1] grid      stats4    stats     graphics  grDevices utils     datasets
    ## [8] methods   base
    ##
    ## other attached packages:
    ##  [1] patchwork_1.1.1             plyranges_1.14.0
    ##  [3] nullrangesData_1.0.0        ExperimentHub_2.2.0
    ##  [5] AnnotationHub_3.2.0         BiocFileCache_2.2.0
    ##  [7] dbplyr_2.1.1                ggplot2_3.3.5
    ##  [9] plotgardener_1.0.1          nullranges_1.0.1
    ## [11] InteractionSet_1.22.0       SummarizedExperiment_1.24.0
    ## [13] Biobase_2.54.0              MatrixGenerics_1.6.0
    ## [15] matrixStats_0.61.0          GenomicRanges_1.46.0
    ## [17] GenomeInfoDb_1.30.0         IRanges_2.28.0
    ## [19] S4Vectors_0.32.2            BiocGenerics_0.40.0
    ##
    ## loaded via a namespace (and not attached):
    ##   [1] colorspace_2.0-2              rjson_0.2.20
    ##   [3] ellipsis_0.3.2                ggridges_0.5.3
    ##   [5] mclust_5.4.8                  XVector_0.34.0
    ##   [7] farver_2.1.0                  bit64_4.0.5
    ##   [9] interactiveDisplayBase_1.32.0 AnnotationDbi_1.56.1
    ##  [11] fansi_0.5.0                   mvtnorm_1.1-3
    ##  [13] cachem_1.0.6                  knitr_1.36
    ##  [15] jsonlite_1.7.2                speedglm_0.3-3
    ##  [17] Rsamtools_2.10.0              png_0.1-7
    ##  [19] shiny_1.7.1                   BiocManager_1.30.16
    ##  [21] compiler_4.1.1                httr_1.4.2
    ##  [23] assertthat_0.2.1              Matrix_1.3-4
    ##  [25] fastmap_1.1.0                 later_1.3.0
    ##  [27] prettyunits_1.1.1             htmltools_0.5.2
    ##  [29] tools_4.1.1                   gtable_0.3.0
    ##  [31] glue_1.5.0                    GenomeInfoDbData_1.2.7
    ##  [33] dplyr_1.0.7                   rappdirs_0.3.3
    ##  [35] Rcpp_1.0.7                    jquerylib_0.1.4
    ##  [37] vctrs_0.3.8                   Biostrings_2.62.0
    ##  [39] strawr_0.0.9                  rtracklayer_1.54.0
    ##  [41] xfun_0.28                     stringr_1.4.0
    ##  [43] mime_0.12                     lifecycle_1.0.1
    ##  [45] restfulr_0.0.13               XML_3.99-0.8
    ##  [47] zlibbioc_1.40.0               MASS_7.3-54
    ##  [49] scales_1.1.1                  hms_1.1.1
    ##  [51] promises_1.2.0.1              parallel_4.1.1
    ##  [53] RColorBrewer_1.1-2            yaml_2.2.1
    ##  [55] curl_4.3.2                    memoise_2.0.0
    ##  [57] yulab.utils_0.0.4             sass_0.4.0
    ##  [59] stringi_1.7.5                 RSQLite_2.2.8
    ##  [61] BiocVersion_3.14.0            highr_0.9
    ##  [63] BiocIO_1.4.0                  filelock_1.0.2
    ##  [65] BiocParallel_1.28.0           rlang_0.4.12
    ##  [67] pkgconfig_2.0.3               bitops_1.0-7
    ##  [69] pracma_2.3.3                  evaluate_0.14
    ##  [71] lattice_0.20-45               purrr_0.3.4
    ##  [73] GenomicAlignments_1.30.0      ks_1.13.2
    ##  [75] labeling_0.4.2                bit_4.0.4
    ##  [77] tidyselect_1.1.1              plyr_1.8.6
    ##  [79] magrittr_2.0.1                R6_2.5.1
    ##  [81] generics_0.1.1                DelayedArray_0.20.0
    ##  [83] DBI_1.1.1                     pillar_1.6.4
    ##  [85] withr_2.4.2                   KEGGREST_1.34.0
    ##  [87] RCurl_1.98-1.5                tibble_3.1.6
    ##  [89] crayon_1.4.2                  KernSmooth_2.23-20
    ##  [91] utf8_1.2.2                    rmarkdown_2.11
    ##  [93] progress_1.2.2                data.table_1.14.2
    ##  [95] blob_1.2.2                    digest_0.6.28
    ##  [97] xtable_1.8-4                  httpuv_1.6.3
    ##  [99] gridGraphics_0.5-1            munsell_0.5.0
    ## [101] ggplotify_0.1.0               bslib_0.3.1