knitr::opts_chunk$set(fig.width=5, fig.height=5, message=FALSE, warning=FALSE)
In this vignette, we demonstrate the block bootstrap functionality implemented in nullranges. See the main nullranges vignette for an overview of the idea of bootstrapping, or the diagram below.
nullranges contains an implementation of a block bootstrap for genomic data, as proposed by @bickel_2010, such that features (ranges) are sampled from the genome in blocks. The original block bootstrapping algorithm for genomic data is implemented in a python software called Genome Structure Correlation, GSC.
Our description of the bootRanges methods is described in @bootRanges.
Minimal code for running bootRanges()
is shown below. Genome
segmentation seg
and excluded regions exclude
are optional.
eh <- ExperimentHub() ah <- AnnotationHub() # some default resources: seg <- eh[["EH7307"]] # pre-built genome segmentation for hg38 exclude <- ah[["AH107305"]] # Kundaje excluded regions for hg38, see below set.seed(5) # set seed for reproducibility blockLength <- 5e5 # size of blocks to bootstrap R <- 10 # number of iterations of the bootstrap # input `ranges` require seqlengths, if missing see `GenomeInfoDb::Seqinfo` seqlengths(ranges) # next, we remove non-standard chromosomes ... ranges <- keepStandardChromosomes(ranges, pruning.mode="coarse") # ... and mitochondrial genome, as these are too short seqlevels(ranges, pruning.mode="coarse") <- setdiff(seqlevels(ranges), "MT") # generate bootstraps boots <- bootRanges(ranges, blockLength=blockLength, R=R, seg=seg, exclude=exclude) # `boots` can then be used with plyranges commands, e.g. join_overlap_*
The boots
object will contain a column, iter
which marks the
different bootstrap samples that were generated. This allows for tidy
analysis with plyranges, e.g. counting the number of overlapping
ranges, per bootstrap iteration. For more examples of combining
bootRanges
with plyranges operations, see the
tidy ranges tutorial.
Several algorithms are implemented in bootRanges()
, including a
segmented and unsegmented version, where in the former, blocks are
sampled with respect to a particular genome segmentation.
Overall, we recommend the segmented block bootstrap given the
heterogeneity of structure across the entire genome. If the purpose is
block bootstrapping ranges within a smaller set of sequences, such as
motifs within transcript sequence, then the unsegmented algorithm
would be sufficient.
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, \dots, S$, only blocks from state s will be used to sample the ranges of state s in each bootstrap sample. The process can be visualized in diagram panel (A) below, where a block with length $L_b$ is randomly sampled with replacement from state "red" and the features (ranges) that overlap this block are then copied to the first tile (which is in the "red" state). The sampling is allowed across chromosome (as shown here), as long as the two blocks are in the same state.
Note that nullranges provides both functions for generating a genome segmentation from e.g. gene density, described below, as well as a default segmentation for hg38 that can be used directly.
An example workflow of bootRanges()
used in combination with
plyranges [@Lee2019] is diagrammed in panel (B) below, and can be
summarized as:
bootRanges()
with optional
arguments seg
(segmentation) and exclude
(excluded regions as
compiled by @excluderanges) to create a BootRanges object ($y'$)knitr::include_graphics("images/bootRanges.jpeg")
In this vignette, we give an example of segmenting the hg38 genome by
Ensembl gene density, performing bootstrap sampling of peaks ranges,
and evaluating overlaps for observed peaks and bootstrap peaks. We
also provide other examples of statistics that can be computed with
the bootRanges
framework, including a single cell multi-omics
example and a special case of bootstrapping features in one region of
the genome.
Proportional blocks:
A finally consideration is whether the blocks used to generate the
bootstrap samples 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. It is
theoretically motivated to have the blocks scale with the overall
extent of the segment state. However, in practice, if the genome
segmentation states are very heterogeneous in size (e.g. orders of
magnitude differences), then the blocks constructed via the
proportional length method for the smaller segmentation states can be
too short to effectively capture inter-range distances. We therefore
recommend proportional length blocks unless some segmentation states
have a much smaller extent than others, in which case fixed length
blocks can be used. This option is visualized on toy data at the end
of this vignette.
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[@encode_exclude], telomeres from UCSC, centromeres. These, and other excludable sets, are assembled in the excluderanges package [@excluderanges].
suppressPackageStartupMessages(library(AnnotationHub)) ah <- AnnotationHub() # hg38.Kundaje.GRCh38_unified_Excludable exclude_1 <- ah[["AH107305"]] # hg38.UCSC.centromere exclude_2 <- ah[["AH107354"]] # hg38.UCSC.telomere exclude_3 <- ah[["AH107355"]] # hg38.UCSC.short_arm exclude_4 <- ah[["AH107356"]] # combine them suppressWarnings({ exclude <- trim(c(exclude_1, exclude_2, exclude_3, exclude_4)) }) exclude <- sort(GenomicRanges::reduce(exclude))
For most genomic datasets we examine, the density of ranges of interest (e.g. ChIP- or ATAC-seq peaks) is often correlated to other large-scale patterns of other genomic features, such as density of genes. @bickel_2010 therefore proposed the idea of bootstrapping with respect to a segmented genome given known, large-scale genomic structures such as isochores ("larger than 300kb").
A genomic segmentation can be considered if it defines large (e.g. on the order of ∼1 Mb), relatively homogeneous segments with respect to feature density, and the variance of the distribution of the test statistics become stable as block length increases (see @bootRanges Fig 2A).
There are two options for choosing a segmentation, either:
Pre-built segmentations
Given that these genome segmentation evaluations take time and involve consideration of multiple criteria, we have provided our recommended segmentation for hg38. nullranges has generated pre-built segmentations for easy use, which were generated using code outlined below in the Segmentation by gene density section.
Pre-built segmentations using either CBS or HMM methods with $L_s=2e6$ considering excludable regions can be downloaded directly from ExperimentHub. We find that the segmentation and block length (500kb) shown in the case study below could be used for most analyses of hg38.
suppressPackageStartupMessages(library(ExperimentHub)) eh <- ExperimentHub() seg_cbs <- eh[["EH7307"]] # prefer CBS for hg38 seg_hmm <- eh[["EH7308"]] seg <- seg_cbs
Segmentation by gene density
This section describes how we generated the pre-built segmentations, such that users with a different genome can generate a segmentation for their own purposes. First we obtain the Ensembl genes [@ensembl2021] for segmenting by gene density. We obtain these using the ensembldb package [@ensembldb].
suppressPackageStartupMessages(library(ensembldb)) suppressPackageStartupMessages(library(EnsDb.Hsapiens.v86)) edb <- EnsDb.Hsapiens.v86 filt <- AnnotationFilterList(GeneIdFilter("ENSG", "startsWith")) g <- genes(edb, filter = filt)
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).
library(GenomeInfoDb) g <- keepStandardChromosomes(g, pruning.mode = "coarse") # MT is too small for bootstrapping, so must be removed seqlevels(g, pruning.mode="coarse") <- setdiff(seqlevels(g), "MT") # normally we would assign a new style, but for recent host issues # that produced vignette build problems, we use `paste0` ## seqlevelsStyle(g) <- "UCSC" seqlevels(g) <- paste0("chr", seqlevels(g)) genome(g) <- "hg38" g <- sortSeqlevels(g) g <- sort(g) table(seqnames(g))
CBS segmentation
We first demonstrate the use of a CBS segmentation as implemented in DNAcopy [@dnacopy].
We load the nullranges and plyranges packages, and patchwork in order to produce grids of plots.
library(nullranges) suppressPackageStartupMessages(library(plyranges)) library(patchwork)
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.
set.seed(5) exclude <- exclude %>% plyranges::filter(width(exclude) >= 500) L_s <- 1e6 seg_cbs <- segmentDensity(g, n = 3, L_s = L_s, exclude = exclude, type = "cbs") plots <- lapply(c("ranges","barplot","boxplot"), function(t) { plotSegment(seg_cbs, exclude, type = t) }) plots[[1]] plots[[2]] + plots[[3]]
Note here, the default ranges plot shows the whole genome. Some of
the state transitions within small regions cannot be visualized. One
can look into specific regions to observe segmentation states, by
specifying the region
argument.
region <- GRanges("chr16", IRanges(3e7,4e7)) plotSegment(seg_cbs, exclude, type="ranges", region=region)
Alternatively: HMM segmentation
Here we use an alternative segmentation implemented in the RcppHMM
CRAN package, using the initGHMM
, learnEM
, and viterbi
functions.
seg_hmm <- segmentDensity(g, n = 3, L_s = L_s, exclude = exclude, type = "hmm") plots <- lapply(c("ranges","barplot","boxplot"), function(t) { plotSegment(seg_hmm, exclude, type = t) }) plots[[1]] plots[[2]] + plots[[3]]
We use a set of DNase hypersensitivity sites (DHS) from the ENCODE project [@encode] in A549 cell line (ENCSR614GWM). Here, for speed, we work with a pre-processed data object from ExperimentHub, created using the following steps:
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 unrelated metadata
columns that we don't need for the bootstrap analysis.
Because we are interested in signal value for DHS peaks later, we only
keep this column. Consider, when creating bootstrapped data, that you will be
creating an object many times larger than your original data
(i.e. multipled by R
the number of bootstrap iterations), so
filtering down to key ranges and selecting only the relevant
metadata can help make the analysis much more efficient.
suppressPackageStartupMessages(library(nullrangesData)) dhs <- DHSA549Hg38() dhs <- dhs %>% plyranges::filter(signalValue > 100) %>% mutate(id = seq_along(.)) %>% plyranges::select(id, signalValue) length(dhs) table(seqnames(dhs))
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 [@Lee2019].
Note that we have already removed non-standard chromosomes and
mitochondrial chromosome, as these are typically shorter than our
desired blockLength
(see e.g. code in Quick Start above).
set.seed(5) # for reproducibility R <- 50 blockLength <- 5e5 boots <- bootRanges(dhs, blockLength, R = R, seg = seg, exclude=exclude) boots
What is returned here? The bootRanges
function returns a
BootRanges object, which is a simple sub-class of GRanges. The
iteration (iter
) and (optionally) 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).
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:
suppressPackageStartupMessages(library(tidyr)) combined <- dhs %>% mutate(iter=0) %>% bind_ranges(boots) %>% plyranges::select(iter) stats <- combined %>% group_by(iter) %>% summarize(n = n()) %>% as_tibble() head(stats)
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.
suppressPackageStartupMessages(library(ggridges)) suppressPackageStartupMessages(library(purrr)) suppressPackageStartupMessages(library(ggplot2)) interdist <- function(dat) { x <- dat[-1,] y <- dat[-nrow(dat),] ifelse(x$seqnames == y$seqnames, x$start + floor((x$width - 1)/2) - y$start - floor((y$width - 1)/2), NA) } # just looking at first 3 iterations... combined %>% plyranges::filter(iter %in% 0:3) %>% mutate(iter = droplevels(iter)) %>% plyranges::select(iter) %>% as_tibble() %>% nest(data = !iter) %>% mutate(interdist = map(data, interdist)) %>% dplyr::select(iter, interdist) %>% unnest(interdist) %>% mutate(type = ifelse(iter == 0, "original", "boot"), interdist = pmax(interdist, 0)) %>% filter(!is.na(interdist)) %>% ggplot(aes(log10(interdist + 1), iter, fill=type)) + geom_density_ridges(alpha = 0.75) + geom_text(data = head(stats, 4), aes(x=1.5, y=iter, label=paste0("n=",n), fill=NULL), vjust=1.5)
We will now show how to combine bootstrapping with plyranges to
perform statistical enrichment analysis.
The general idea will be to combine the long vector of bootstrapped
ranges, indexed by iter
, with another set of ranges to compute
enrichment. We will explore this idea across a number of case studies
below. In pseudocode, the general outline will be:
# pseudocode for the general paradigm: boots <- bootRanges(y) # make bootstrapped y x %>% join_overlap_inner(boots) %>% # overlaps of x with bootstrapped y group_by(x_id, iter) %>% # collate by x ID and the bootstrap iteration summarize(some_statistic = ...) %>% # compute some summary on metadata as_tibble() %>% # pass to tibble complete( x_id, iter, # for any missing combinations of x ID and iter... fill=list(some_statistic = 0) # ...fill in missing values )
Suppose we have a set of features x
and we are interested in
evaluating the 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
).
x <- GRanges("chr2", IRanges(1 + 50:99 * 1e6, width=1e6), x_id=1:50)
x <- x %>% mutate(n_overlaps = count_overlaps(., dhs)) sum( x$n_overlaps )
We can repeat this with the bootstrapped features using a group_by
command, a summarize
, followed by a second group_by
and summarize
.
If it is your first time working with plyranges, it may help you to
step through these commands one by one, breaking the pipe at
intermediate points, to understand what the intermediate
output is, and how it combines to provide the final statistics.
Note that we need to use tidyr::complete
in order to fill in
combinations of x
and iter
where the overlap was 0.
boot_stats <- x %>% join_overlap_inner(boots) %>% group_by(x_id, iter) %>% summarize(n_overlaps = n()) %>% as_tibble() %>% complete(x_id, iter, fill=list(n_overlaps = 0)) %>% group_by(iter) %>% summarize(sumOverlaps = sum(n_overlaps))
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 feature-wise
statistics instead of 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 sum
number of overlap bootstrapped peaks with $x$.
suppressPackageStartupMessages(library(ggplot2)) ggplot(boot_stats, aes(sumOverlaps)) + geom_histogram(binwidth=5)+ geom_vline(xintercept = sum(x$n_overlaps), linetype = "dashed")
x_obs <- x %>% join_overlap_inner(dhs,maxgap=1e3) sum(x_obs$signalValue ) boot_stats <- x %>% join_overlap_inner(boots,maxgap=1e3) %>% group_by(x_id, iter) %>% summarize(Signal = sum(signalValue)) %>% as_tibble() %>% complete(x_id, iter, fill=list(Signal = 0)) %>% group_by(iter) %>% summarize(sumSignal = sum(Signal))
Still in this case, our observed value falls within the distribution of bootstrapped statistics.
ggplot(boot_stats, aes(sumSignal)) + geom_histogram()+ geom_vline(xintercept = sum(x_obs$signalValue), linetype = "dashed")
This case study provides an example of a more complex statistic of
interest.
Instead of computing simple overlap, or some function of a single
metadata column, e.g. signal
above, a statistic can be computed on a
count matrix from a SummerizedExperiment or SingleCellExperiment.
Here we provide unevaluated code for such an example.
We use a trick to perform computation on a count matrix within plyranges:
we extract count matrix data from a SingleCellExperiment
(e.g. pseudobulk summaries of cell type specific expression or
accessibility) and save it in a GRanges's metadata column as a
NumericList()
format. This allows us to perform the block bootstrap
and quickly compute statistics using plyranges. The case study we
consider here is to assess the correlation of gene expression and
promoter peak accessibility from a dataset of Chromium Single Cell
Multiome ATAC + Gene Expression.
We first load the 10X
Multiome
data, as compiled by the MOFA2 tutorial.
The first two pre-processed steps are included in nullrangesData in the
inst/scripts/make-multiome-data.R
script. Here, we aggregate 14
cell types across cluster-sample groups.
## split sparse count matrix into NumericList # sc_rna <- rna_Granges[-which(rna.sd==0)] %>% # mutate(counts1 = NumericList(asplit(rna.scaled, 1))) %>% sort() # sc_promoter <- promoter_Granges[-which(promoter.sd==0)] %>% # mutate(counts2 = NumericList(asplit(promoter.scaled, 1))) %>% sort() suppressPackageStartupMessages(library(nullrangesData)) data("sc_rna") sc_rna sc_rna$counts1 nct <- length(sc_rna$counts1[[1]]) print(paste("There are", nct, "cell types")) data("sc_promoter") ## bootstrap promoters library(BSgenome.Hsapiens.UCSC.hg38) genome <- BSgenome.Hsapiens.UCSC.hg38 seqlengths(sc_promoter) <- seqlengths(genome)[1:22] # pull chrom lens from USCS bootranges <- bootRanges(sc_promoter, blockLength = 5e5, R=50)
There are two options for downstream pipelines for computing statistics after generating a BootRanges object. One is to use plyranges in downstream analysis as in the previous sections, while another is to create a tidySummarizedExperiment object used to compute the statistics of interest.
Plyranges
We can compute the mean correlation of the all pairs of genes and promoter peaks over all bootstrap iterations and compare to the observed statistic. The histogram will show that our observed value falls far from the distribution of bootstrapped statistics, as we expected for this dataset.
cor_gr <- sc_rna %>% join_overlap_inner(sc_promoter, maxgap=1000) %>% mutate(rho = 1/(nct-1) * sum(counts1 * counts2)) %>% summarise(meanCor = mean(rho)) ## mean correlation distribution cor_boot <- sc_rna %>% join_overlap_inner(bootranges, maxgap=1000) %>% # vectorized code are 10 times faster than cor(counts1, counts2) for plyranges pipeline mutate(rho = 1/(nct-1) * sum(counts1 * counts2)) %>% select(rho, iter) %>% group_by(iter) %>% summarise(meanCor = mean(rho)) %>% as.data.frame() ggplot(cor_boot, aes(meanCor)) + geom_histogram(binwidth=0.01)+ geom_vline(xintercept = cor_gr$meanCor,linetype = "dashed")+ theme_classic()
If one is interested in assessing the significance of gene-promoter pairs
per gene, one can also use plyranges::select
(gene, peak, rho)
and plyranges::group_by
gene.
tidySummarizedExperiment
In order to use tidySummarizedExperiment to compute the same operation, we first must create an RangedSE for each modality, and then we can find the overlaps..
library(tidySummarizedExperiment) library(purrr) # make an SE where each row is an overlap makeOverlapSE <- function(se_rna, se_promoter) { idx <- rowRanges(se_rna) %>% join_overlap_inner(rowRanges(se_promoter),maxgap = 1000) assay_x <- assay(se_rna, "rna")[ idx$gene, ] assay_y <- assay(se_promoter, "promoter")[ idx$peak, ] # this is needed to build the SE rownames(assay_x) <- rownames(assay_y) <- seq_along( idx$gene ) names(idx) <- seq_along( idx$gene ) SummarizedExperiment( assays=list(x=assay_x, y=assay_y), rowRanges=idx ) } # create SE for observed data se_rna <- SummarizedExperiment( assays=list(rna=do.call(rbind,sc_rna$counts1)), rowRanges=sc_rna) se_promoter <- SummarizedExperiment( assays=list(promoter=do.call(rbind,sc_promoter$counts2)), rowRanges=sc_promoter) se <- makeOverlapSE(se_rna, se_promoter) se <- se %>% as_tibble() %>% nest(data = -.feature) %>% mutate(rho = map(data, function(data) data %>% summarize(rho = cor(x, y)) )) %>% unnest(rho) %>% select(-data) print(paste("mean correlation is", mean(se$rho)))
Repeating for the bootstrapping case
# create SE for bootranges data se_promoter_boots <- SummarizedExperiment( assays=list(promoter=do.call(rbind,bootranges$counts2)), rowRanges=bootranges) se_boots <- makeOverlapSE(se_rna, se_promoter_boots) se_boots <- se_boots %>% as_tibble() %>% nest(data = -c(.feature,iter)) %>% # vectorized code similar to cor(counts1, counts2) for tidySE pipeline mutate(rho = map(data, function(data) data %>% summarize(rho = cor(x, y)) )) %>% unnest(rho) %>% select(-data) cor_boot <- se_boots %>% group_by(iter) %>% summarise(meanCor = mean(rho)) ggplot(cor_boot, aes(meanCor)) + geom_histogram(binwidth=0.01)+ geom_vline(xintercept = mean(se$rho),linetype = "dashed")+ theme_classic()
For more examples of combining bootRanges
from nullranges with
plyranges piped operations, see the relevant chapter in the
tidy-ranges-tutorial
book.
Generally, it makes sense to block bootstrap the entire genome at once. This is motivated by the "tidy analysis" paradigm where loops are avoided by stacking data into a longer format. This makes computation more efficient in our case (as a single overlap call can be made with all regions of interest at once, across multiple bootstrap iterations), and it also can simplify code and avoid repetition.
However, in some cases, there is a single region of interest, and it is desired to generate bootstrap data within this one region. For this, we have a convenience function that enables bootstrap computation.
Suppose we have data in the following region of chromosome 1:
suppressPackageStartupMessages(library(nullrangesData)) dhs <- DHSA549Hg38() region <- GRanges("chr1", IRanges(10e6 + 1, width=1e6)) x <- GRanges("chr1", IRanges(10e6 + 0:9 * 1e5 + 1, width=1e4)) y <- dhs %>% filter_by_overlaps(region) %>% select(NULL) x %>% mutate(num_overlaps = count_overlaps(., y))
We can easily bootstrap data just in this region using the following code:
seg <- oneRegionSegment(region, seqlength=248956422) y <- keepSeqlevels(y, "chr1") set.seed(1) boot <- bootRanges(y, blockLength=1e5, R=1, seg=seg, proportionLength=FALSE) boot x %>% mutate(num_overlaps = count_overlaps(., boot))
Here it is important to use proportionLength=FALSE
so that the
blocks will be of the size specified and not smaller (they would
otherwise be scaled down proportional to the fraction of region
compared to the chromosome).
Below we present a toy example for visualizing the segmented block bootstrap. First, we define a helper function for plotting GRanges using plotgardener [@Kramer2022]. 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).
suppressPackageStartupMessages(library(plotgardener)) my_palette <- function(n) { head(c("red","green3","red3","dodgerblue", "blue2","green4","darkred"), n) } plotGRanges <- function(gr) { pageCreate(width = 5, height = 5, xgrid = 0, ygrid = 0, showGuides = TRUE) for (i in seq_along(seqlevels(gr))) { chrom <- seqlevels(gr)[i] chromend <- seqlengths(gr)[[chrom]] suppressMessages({ p <- pgParams(chromstart = 0, chromend = chromend, x = 0.5, width = 4*chromend/500, height = 2, at = seq(0, chromend, 50), fill = colorby("state_col", palette=my_palette)) prngs <- plotRanges(data = gr, params = p, chrom = chrom, y = 2 * i, just = c("left", "bottom")) annoGenomeLabel(plot = prngs, params = p, y = 0.1 + 2 * i) }) } }
Create a toy genome segmentation:
library(GenomicRanges) seq_nms <- rep(c("chr1","chr2"), c(4,3)) seg <- GRanges( seqnames = seq_nms, IRanges(start = c(1, 101, 201, 401, 1, 201, 301), width = c(100, 100, 200, 100, 200, 100, 100)), seqlengths=c(chr1=500,chr2=400), state = c(1,2,1,3,3,2,1), state_col = factor(1:7) )
We can visualize with our helper function:
plotGRanges(seg)
Now create small ranges distributed uniformly across the toy genome:
set.seed(1) n <- 200 gr <- GRanges( seqnames=sort(sample(c("chr1","chr2"), n, TRUE)), IRanges(start=round(runif(n, 1, 500-10+1)), width=10) ) suppressWarnings({ seqlengths(gr) <- seqlengths(seg) }) gr <- gr[!(seqnames(gr) == "chr2" & end(gr) > 400)] gr <- sort(gr) idx <- findOverlaps(gr, seg, type="within", select="first") gr <- gr[!is.na(idx)] idx <- idx[!is.na(idx)] gr$state <- seg$state[idx] gr$state_col <- factor(seg$state_col[idx]) plotGRanges(gr)
Scaling vs. not scaling by segment length
We can visualize block bootstrapped ranges when the blocks do not scale to segment state length:
set.seed(1) gr_prime <- bootRanges(gr, blockLength = 25, seg = seg, proportionLength = FALSE) plotGRanges(gr_prime)
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).
set.seed(1) gr_prime <- bootRanges(gr, blockLength = 50, seg = seg, proportionLength = TRUE) plotGRanges(gr_prime)
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.
sessionInfo()
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.