Overview

The following vignette describes the nullranges implementation of the block bootstrap with respect to a genomic segmentation. See the main nullranges vignette for an overview of the idea of bootstrapping, or below diagram, and there is additionally a vignette on block boostrapping without respect to segmentation-Unsegmented block bootstrap.

As proposed by Bickel et al. (2010), nullranges contains an implementation of a block bootstrap, such that features are sampled from the genome in blocks. The original block bootstrapping algorithm is implemented in a python software called Genome Structure Correlation, GSC.

Our description of the bootRanges methods is described in Mu et al. (2022).

In a segmented block bootstrap, the blocks are sampled and placed within regions of a genome segmentation. That is, for a genome segmented into states 1,2,…,S, blocks from state s will be used to tile the ranges of state s in each bootstrap sample. The process can be visualized in (A), a block with length \(L_b\) is \(\color{brown}{\text{randomly}}\) selected from state “red” and move to a \(\color{brown}{\text{tile}}\) block across chromosome within same states.

An example workflow of bootRanges used in combination with plyranges is diagrammed in (B), and can be summarized as:

  1. Compute overlap between GRanges of feature \(x\) and GRanges of feature \(y\) to assess overlap or to derive other statistic of interest
  2. bootRanges() with optional segmentation and exclude to create a bootRanges object \(y'\)
  3. Compute overlap between GRanges of feature \(x\) and \(y'\) to derive bootstrap distribution of statistics
  4. A \(z\) test can be performed for testing the null hypothesis that there is no true biological enrichment of the original data (that bootstrap data often has as high an enrichment as the observed data)

The segmented block bootstrap has two options, either:

  • Perform a de-novo segmentation of the genome using feature density, e.g. gene density
  • Use exiting segmentation (e.g. ChromHMM, etc.) downloaded from AnnotationHub or external to Bioconductor (BED files imported with rtracklayer)

In this vignette, we give an example of segmenting the hg38 genome by Ensembl gene density, create bootstrapped peaks and evaluate overlaps for observed peaks and bootstrap peaks, then we profile the time to generate a single block bootstrap sample. Finally, we use a toy dataset to visualize what a segmented block bootstrap sample looks like with respect to a genome segmentation.

A finally consideration is whether the blocks should scale proportionally to the segment state length, with the default setting of proportionLength=TRUE. When blocks scale proportionally, blockLength provides the maximal length of a block, while the actual block length used for a segmentation state is proportional to the fraction of genomic basepairs covered by that state. This option is visualized on toy data at the end of this vignette.

Import excluded regions

\(\color{brown}{\text{To avoid placing bootstrap features into regions of the genome that don’t typically have features}}\). We import excluded regions including ENCODE-produced excludable regions(Amemiya, Kundaje, and Boyle 2019), telomeres from UCSC, centromeres (Commo 2022). For easy use, pre-combined excludable regions is stored in ExperimentHub. These steps using excluderanges package (Dozmorov et al. 2022) are included in nullrangesData in the inst/scripts/make-segmentation-hg38.R script.

Pre-built segmentations

nullranges has generated pre-built segmentations for easy use by following below section Segmentation by gene density. Either pre-built segmentations using CBS or HMM methods with \(L_s=2e6\) considering excludable regions can be selected from ExperimentHub.

Segmentation by gene density

First we obtain the Ensembl genes (Howe et al. 2020) for segmenting by gene density. We obtain these using the ensembldb package (Rainer, Gatto, and Weichenberger 2019).

We perform some processing to align the sequences (chromosomes) of g with our excluded regions and our features of interest (DNase hypersensitive sites, or DHS, defined below).

## 
##  chr1  chr2  chr3  chr4  chr5  chr6  chr7  chr8  chr9 chr10 chr11 chr12 chr13 
##  5194  3971  3010  2505  2868  2863  2867  2353  2242  2204  3235  2940  1304 
## chr14 chr15 chr16 chr17 chr18 chr19 chr20 chr21 chr22  chrX  chrY 
##  2224  2152  2511  2995  1170  2926  1386   835  1318  2359   523

CBS segmentation

We first demonstrate the use a CBS segmentation as implemented in DNAcopy (Olshen et al. 2004).

We load the nullranges and plyranges packages, and patchwork in order to produce grids of plots.

We subset the excluded ranges to those which are 500 bp or larger. The motivation for this step is to avoid segmenting the genome into many small pieces due to an abundance of short excluded regions. Note that we re-save the excluded ranges to exclude.

Here, and below, we need to specify plyranges::filter as it conflicts with filter exported by ensembldb.

## Analyzing: Sample.1

Note here, the default ranges plot gives whole genome and every fractured bind regions represents state transformations happens. However, some transformations within small ranges cannot be visualized, e.g 1kb. If user want to look into specific ranges of segmentation state, the region argument is flexible to support.

Alternatively: HMM segmentation

Here we use an alternative segmentation implemented in the RcppHMM CRAN package, using the initGHMM, learnEM, and viterbi functions.

## Finished at Iteration: 111 with Error: 9.22278e-06

Segmented block bootstrap

We use a set of DNase hypersensitivity sites (DHS) from the ENCODE project (ENCODE 2012) in A549 cell line (ENCSR614GWM). Here, for speed, we work with a pre-processed data object from ExperimentHub, created using the following steps:

  • Download ENCODE DNase hypersensitive peaks in A549 from AnnotationHub
  • Subset to standard chromosomes and remove mitochondrial DNA
  • Use a chain file from UCSC to lift ranges from hg19 to hg38
  • Sort the DHS features to be bootstrapped

These steps are included in nullrangesData in the inst/scripts/make-dhs-data.R script.

For speed of the vignette, we restrict to a smaller number of DHS, filtering by the signal value. We also remove metadata columns that we don’t need for the bootstrap analysis. Consider, when creating bootstrapped data, that you will be creating an object many times larger than your original features, so \(\color{brown}{\text{filtering and trimming}}\) extra metadata can help make the analysis more efficient.

## see ?nullrangesData and browseVignettes('nullrangesData') for documentation
## loading from cache
## [1] 6214
## 
##  chr1  chr2  chr3  chr4  chr5  chr6  chr7  chr8  chr9 chr10 chr11 chr12 chr13 
##  1436   252   108    30   148    51   184   146   155   443   436   526    20 
## chr14 chr15 chr16 chr17 chr18 chr19 chr20 chr21 chr22  chrX  chrY 
##   197   265   214   715    20   649   142    31    19    17    10

Now we apply a segmented block bootstrap with blocks of size 500kb, to the peaks. Here we show generation of 50 iterations of a block bootstrap followed by a typical overlap analysis using plyranges (Lee, Cook, and Lawrence 2019). (We might normally do 100 iterations or more, depending on the granularity of the bootstrap distribution that is needed.)

## BootRanges object with 310726 ranges and 4 metadata columns:
##            seqnames            ranges strand |        id     block  iter
##               <Rle>         <IRanges>  <Rle> | <integer> <integer> <Rle>
##        [1]     chr1     242791-242940      * |       347         5     1
##        [2]     chr1     256031-256180      * |       348         5     1
##        [3]     chr1     391535-391684      * |      5301         8     1
##        [4]     chr1     421046-421195      * |      5302         8     1
##        [5]     chr1     438186-438335      * |      5303         8     1
##        ...      ...               ...    ... .       ...       ...   ...
##   [310722]     chrY 27090908-27091057      * |      2133     12441    50
##   [310723]     chrY 27194968-27195117      * |      2134     12441    50
##   [310724]     chrY 27224188-27224337      * |      2135     12441    50
##   [310725]     chrY 27234153-27234302      * |      2136     12441    50
##   [310726]     chrY 27789879-27790028      * |      2116     12442    50
##            blockLength
##                  <Rle>
##        [1]      500000
##        [2]      500000
##        [3]      500000
##        [4]      500000
##        [5]      500000
##        ...         ...
##   [310722]      500000
##   [310723]      500000
##   [310724]      500000
##   [310725]      500000
##   [310726]      500000
##   -------
##   seqinfo: 24 sequences from hg38 genome

What is returned here? The bootRanges function returns a bootRanges object, which is a simple sub-class of GRanges. The iteration (iter) and the block length (blockLength) are recorded as metadata columns, accessible via mcols. We return the bootstrapped ranges as GRanges rather than GRangesList, as the former is more compatible with downstream overlap joins with plyranges, where the iteration column can be used with group_by to provide per bootstrap summary statistics, as shown below.

Note that we use the exclude object from the previous step, which does not contain small ranges. If one wanted to also avoid generation of bootstrapped features that overlap small excluded ranges, then omit this filtering step (use the original, complete exclude feature set).

Assessing quality of bootstrap samples

We can examine properties of permuted y over iterations, and compare to the original y. To do so, we first add the original features as iter=0. Then compute summaries:

## # A tibble: 6 × 2
##   iter      n
##   <fct> <int>
## 1 0      6214
## 2 1      6150
## 3 2      6276
## 4 3      6249
## 5 4      6293
## 6 5      6365

We can also look at distributions of various aspects, e.g. here the inter-feature distance of features, across a few of the bootstraps and the original feature set y.

## Picking joint bandwidth of 0.198

Derive statistics of interest - Use with plyranges

Suppose we have a set of features x and we are interested in evaluating the \(\color{brown}{\text{enrichment of this set with the DHS}}\). We can calculate for example the sum observed number of overlaps for features in x with DHS in whole genome (or something more complicated, e.g. the maximum log fold change or signal value for DHS peaks within a maxgap window of x).

## [1] 1.28

We can repeat this with the bootstrapped features using a group_by command, a summarize, followed by a second group_by and summarize. It may help to step through these commands one by one to understand what the intermediate output is.

Note that we need to use tidyr::complete in order to fill in combinations of x and iter where the overlap was 0.

The above code, first grouping by x_id and iter, then subsequently by iter is general and allows for more complex analysis than just mean overlap (e.g. how many times an x range has 1 or more overlap, what is the mean or max signal value for peaks overlapping ranges in x, etc.).

If one is interested in assessing \(\color{brown}{\text{feature-wise}}\) statistics instead of \(\color{brown}{\text{genome-wise}}\) statistics, eg.,the mean observed number of overlaps per feature or mean base pair overlap in x, one can also group by both (block,iter). 10,000 total blocks may therefore be sufficient to derive a bootstrap distribution, avoiding the need to generate many bootstrap genomes of data.

Finally we can plot a histogram. In this case, as the x features were arbitrary, our observed value falls within the distribution of mean overlap with bootstrapped data.

For more examples of combining bootRanges from nullranges with plyranges piped operations, see the relevant chapter in the tidy-ranges-tutorial book.

Timing on DHS peaks

Here, we test the speed of the various options for bootstrapping (see below for visualization of the difference).

## Unit: milliseconds
##     expr      min       lq     mean   median       uq      max neval cld
##     prop 96.14043 207.9214 227.2704 227.7603 292.6067 333.9900    10   a
##  no_prop 93.80061 187.6896 263.2685 261.0181 294.0343 508.7864    10   a

Visualizing toy bootstrap samples

Below we present a toy example for visualizing the segmented block bootstrap. First, we define a helper function for plotting GRanges using plotgardener (Kramer et al. 2021). A key aspect here is that we color the original and bootstrapped ranges by the genomic state (the state of the segmentation that the original ranges fall in).

Create a toy genome segmentation:

We can visualize with our helper function:

Now create small ranges distributed uniformly across the toy genome:

Not scaling by segmentation

We can visualize block bootstrapped ranges when the blocks do not scale to segment state length:

Scaling by segmentation

This time the blocks scale to the segment state length. Note that in this case blockLength is the maximal block length possible, but the actual block lengths per segment will be smaller (proportional to the fraction of basepairs covered by that state in the genome segmentation).

Note that some ranges from adjacent states are allowed to be placed within different states in the bootstrap sample. This is because, during the random sampling of blocks of original data, a block is allowed to extend beyond the segmentation region of the state being sampled, and features from adjacent states are not excluded from the sampled block.

Session information

## 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] microbenchmark_1.4.9        purrr_0.3.5                
##  [3] ggridges_0.5.4              tidyr_1.2.1                
##  [5] EnsDb.Hsapiens.v86_2.99.0   ensembldb_2.22.0           
##  [7] AnnotationFilter_1.22.0     GenomicFeatures_1.50.0     
##  [9] AnnotationDbi_1.60.0        patchwork_1.1.2            
## [11] plyranges_1.18.0            nullrangesData_1.3.0       
## [13] ExperimentHub_2.6.0         AnnotationHub_3.6.0        
## [15] BiocFileCache_2.6.0         dbplyr_2.2.1               
## [17] ggplot2_3.3.6               plotgardener_1.4.0         
## [19] nullranges_1.4.0            InteractionSet_1.26.0      
## [21] SummarizedExperiment_1.28.0 Biobase_2.58.0             
## [23] MatrixGenerics_1.10.0       matrixStats_0.62.0         
## [25] GenomicRanges_1.50.0        GenomeInfoDb_1.34.0        
## [27] IRanges_2.32.0              S4Vectors_0.36.0           
## [29] BiocGenerics_0.44.0        
## 
## loaded via a namespace (and not attached):
##   [1] RcppHMM_1.2.2                 lazyeval_0.2.2               
##   [3] splines_4.2.1                 BiocParallel_1.32.0          
##   [5] TH.data_1.1-1                 digest_0.6.30                
##   [7] yulab.utils_0.0.5             htmltools_0.5.3              
##   [9] fansi_1.0.3                   magrittr_2.0.3               
##  [11] memoise_2.0.1                 ks_1.13.5                    
##  [13] Biostrings_2.66.0             sandwich_3.0-2               
##  [15] prettyunits_1.1.1             jpeg_0.1-9                   
##  [17] colorspace_2.0-3              blob_1.2.3                   
##  [19] rappdirs_0.3.3                xfun_0.34                    
##  [21] dplyr_1.0.10                  crayon_1.5.2                 
##  [23] RCurl_1.98-1.9                jsonlite_1.8.3               
##  [25] survival_3.4-0                zoo_1.8-11                   
##  [27] glue_1.6.2                    gtable_0.3.1                 
##  [29] zlibbioc_1.44.0               XVector_0.38.0               
##  [31] strawr_0.0.9                  DelayedArray_0.24.0          
##  [33] scales_1.2.1                  mvtnorm_1.1-3                
##  [35] DBI_1.1.3                     Rcpp_1.0.9                   
##  [37] xtable_1.8-4                  progress_1.2.2               
##  [39] gridGraphics_0.5-1            bit_4.0.4                    
##  [41] mclust_6.0.0                  httr_1.4.4                   
##  [43] RColorBrewer_1.1-3            speedglm_0.3-4               
##  [45] ellipsis_0.3.2                pkgconfig_2.0.3              
##  [47] XML_3.99-0.12                 farver_2.1.1                 
##  [49] sass_0.4.2                    utf8_1.2.2                   
##  [51] DNAcopy_1.72.0                ggplotify_0.1.0              
##  [53] tidyselect_1.2.0              labeling_0.4.2               
##  [55] rlang_1.0.6                   later_1.3.0                  
##  [57] munsell_0.5.0                 BiocVersion_3.16.0           
##  [59] tools_4.2.1                   cachem_1.0.6                 
##  [61] cli_3.4.1                     generics_0.1.3               
##  [63] RSQLite_2.2.18                evaluate_0.17                
##  [65] stringr_1.4.1                 fastmap_1.1.0                
##  [67] yaml_2.3.6                    knitr_1.40                   
##  [69] bit64_4.0.5                   KEGGREST_1.38.0              
##  [71] mime_0.12                     pracma_2.4.2                 
##  [73] xml2_1.3.3                    biomaRt_2.54.0               
##  [75] compiler_4.2.1                filelock_1.0.2               
##  [77] curl_4.3.3                    png_0.1-7                    
##  [79] interactiveDisplayBase_1.36.0 tibble_3.1.8                 
##  [81] bslib_0.4.0                   stringi_1.7.8                
##  [83] highr_0.9                     lattice_0.20-45              
##  [85] ProtGenerics_1.30.0           Matrix_1.5-1                 
##  [87] vctrs_0.5.0                   pillar_1.8.1                 
##  [89] lifecycle_1.0.3               BiocManager_1.30.19          
##  [91] jquerylib_0.1.4               data.table_1.14.4            
##  [93] bitops_1.0-7                  httpuv_1.6.6                 
##  [95] rtracklayer_1.58.0            R6_2.5.1                     
##  [97] BiocIO_1.8.0                  promises_1.2.0.1             
##  [99] KernSmooth_2.23-20            codetools_0.2-18             
## [101] MASS_7.3-58.1                 assertthat_0.2.1             
## [103] rjson_0.2.21                  withr_2.5.0                  
## [105] GenomicAlignments_1.34.0      Rsamtools_2.14.0             
## [107] multcomp_1.4-20               GenomeInfoDbData_1.2.9       
## [109] parallel_4.2.1                hms_1.1.2                    
## [111] rmarkdown_2.17                shiny_1.7.3                  
## [113] restfulr_0.0.15

References

Amemiya, Haley M, Anshul Kundaje, and Alan P Boyle. 2019. “The ENCODE Blacklist: Identification of Problematic Regions of the Genome.” Scientific Reports 9 (1): 9354. https://doi.org/10.1038/s41598-019-45839-z.

Bickel, Peter J., Nathan Boley, James B. Brown, Haiyan Huang, and Nancy R. Zhang. 2010. “Subsampling Methods for Genomic Inference.” The Annals of Applied Statistics 4 (4): 1660–97. https://doi.org/10.1214/10-{AOAS363}.

Commo, Frederic. 2022. RCGH: Comprehensive Pipeline for Analyzing and Visualizing Array-Based Cgh Data. https://github.com/fredcommo/rCGH.

Dozmorov, Mikhail G., Eric Davis, Wancen Mu, Stuart Lee, Tim Triche, Douglas Phanstiel, and Michael Love. 2022. Excluderanges. https://github.com/mdozmorov/excluderanges.

ENCODE. 2012. “An integrated encyclopedia of DNA elements in the human genome.” Nature 489 (7414): 57–74. https://doi.org/10.1038/nature11247.

Howe, Kevin L, Premanand Achuthan, James Allen, Jamie Allen, Jorge Alvarez-Jarreta, M Ridwan Amode, Irina M Armean, et al. 2020. “Ensembl 2021.” Nucleic Acids Research 49 (D1): D884–D891. https://doi.org/10.1093/nar/gkaa942.

Kramer, Nicole E, Eric S Davis, Craig D Wenger, Erika M Deoudes, Sarah M Parker, Michael I Love, and Douglas H Phanstiel. 2021. “Plotgardener: Cultivating precise multi-panel figures in R.” bioRxiv. https://doi.org/10.1101/2021.09.08.459338.

Lee, Stuart, Dianne Cook, and Michael Lawrence. 2019. “Plyranges: A Grammar of Genomic Data Transformation.” Genome Biology 20 (1): 4. https://doi.org/10.1186/s13059-018-1597-8.

Mu, Wancen, Eric S. Davis, Stuart Lee, Mikhail G. Dozmorov, Douglas H. Phanstiel, and Michael I. Love. 2022. “BootRanges: Flexible Generation of Null Sets of Genomic Ranges for Hypothesis Testing.” bioRxiv. https://doi.org/10.1101/2022.09.02.506382.

Olshen, A. B., E. S. Venkatraman, R. Lucito, and M. Wigler. 2004. “Circular binary segmentation for the analysis of array-based DNA copy number data.” Biostatistics 5 (4): 557–72.

Rainer, Johannes, Laurent Gatto, and Christian X Weichenberger. 2019. “ensembldb: an R package to create and use Ensembl-based annotation resources.” Bioinformatics 35 (17): 3151–3. https://doi.org/10.1093/bioinformatics/btz031.