betterChromVAR 1.1.12
The chromVAR R package was originally developed by Alicia Schep and colleagues from the Greenleaf lab (Schep et al. 2017). Although originally designed for single-cell ATAC-seq data, it has been shown to be highly sensitive for bulk data as well (Gerbaldo et al. 2024). The aim of the method is to infer, based on the accessibility of motif matches, the relative activity of transcription factors (TFs) in each sample or cell, adjusting for technical biases (GC content and enrichment bias). It is recommended that you read the chromVAR documentation before using this package.
betterChromVAR is first and foremost a considerably faster, analytical re-implementation of the original method (it is also considerably faster than the C++ reimplementation in ArchR (Granja et al. 2021). Contrarily to the original chromVAR, it is entirely deterministic, and achieves much higher efficiency by replacing permutations with an analytical solution and computing expectations and variance at the level of bias bins, instead of in the peak-space (for more detail, see (Wang et al. 2026)).
In addition, betterChromVAR includes a few additions, such as simpler weighted expectations, bias shrinkage, and an ATAC-seq normalization method based on the chromVAR logic. Importantly, however, not all functionalities of chromVAR have been reimplemented here: betterChromVAR chiefly focuses on the key task of efficiently computing motif deviations.
betterChromVARbetterChromVAR is a R package available via the
Bioconductor repository. It can be installed
using the following commands in your R session:
if (!requireNamespace("BiocManager", quietly = TRUE)) {
install.packages("BiocManager")
}
BiocManager::install("betterChromVAR")
## Check that you have a valid Bioconductor installation
BiocManager::valid()betterChromVAR takes two primary inputs: 1) the counts in
peaks across cells or samples, and 2) an annotation of which TFs/motifs
match which peak1 The annotation can be binary or probabilistic (i.e. from 0 to 1).
If using motif matches, however, it is recommended to input binary
matches, as the magnitude of motif scores is generally poorly
correlated to actual binding.. The peak counts should be provided as a
RangedSummarizedExperiment object2 See SummarizedExperiment if you’re not familiar
with those, while the annotation can either
be in that format too or provided as a (sparse) matrix.
There are multiple ways of generating your peak counts; the original
chromVAR package includes such a function
(getCounts()), and the epiwraps
package has some with more options (functions peakCountsFromBAM() and
peakCountsFromFrags()). Similarly, the motifmatchr
package can be used to generate the motif matching annotation. An
important consideration is that the peaks or regions used for the
purpose of this analysis should have similar widths. It is thus highly
recommended, before generating the count and annotation matrices, to
resize your regions (e.g. using
peaks <- resize(peaks, width=300, fix="center")). The exact size can
be something not too large (otherwise the presence/absence of a motif
becomes meaningless) and ideally close to the median size of your
original peaks. In addition, is it advisable to restrict your peaks to
those that are on standard chromosomes3 See keepStandardChromosomes() from the
GenomeInfoDb package.
Here we’ll use dummy data as example:
suppressPackageStartupMessages({
library(SummarizedExperiment)
library(betterChromVAR)
})
attach(getDummyData())
counts## class: RangedSummarizedExperiment
## dim: 500 10
## metadata(0):
## assays(1): counts
## rownames: NULL
## rowData names(1): bias
## colnames(10): sample1 sample2 ... sample9 sample10
## colData names(0):head(motifMatches)## 6 x 5 sparse Matrix of class "dgCMatrix"
## motif1 motif2 motif3 motif4 motif5
## [1,] . . . 1 1
## [2,] . . . . .
## [3,] . . . . .
## [4,] . . 1 . .
## [5,] . . . . .
## [6,] . . . . .In this case, we can see that rowData(counts) already has a bias
column indicating the GC content of the regions:
rowData(counts)## DataFrame with 500 rows and 1 column
## bias
## <numeric>
## 1 0.444343
## 2 0.528662
## 3 0.496576
## 4 0.505714
## 5 0.407049
## ... ...
## 496 0.577362
## 497 0.505315
## 498 0.528168
## 499 0.481967
## 500 0.548394Had this not been the case, we would first need to add this using:
# not run
counts <- addGCBias(counts, genome=my_genome)where my_genome is a BSgenome object or similar (e.g.
an FaFile4 See the Rsamtools package.).
Once we have this, we can launch the computation of the deviations :
dev <- betterChromVAR(counts, motifMatches)
dev## class: SummarizedExperiment
## dim: 5 10
## metadata(0):
## assays(2): deviations z
## rownames(5): motif1 motif2 motif3 motif4 motif5
## rowData names(8): N Neff ... var.pval var.adjPval
## colnames(10): sample1 sample2 ... sample9 sample10
## colData names(0):The resulting dev object is a SummarizedExperiment with the same
columns (and colData) as the original counts object, but with the
motifs as rows (instead of the original peaks). The values can be
interpreted as the relative activity, across samples or cells, of the
corresponding TFs5 Note, however, that similar motifs will be given similar activity
estimates, so that it is often hard to know which of a set of
highly-similar motifs is in fact responsible for the observed
signal.. When the function is used with fpc=FALSE (and
otherwise default arguments) the results are virtually the same as the original
chromVAR (averaging over the noise coming from the random seeds). fpc=TRUE
(default), it includes an additional finite-population correction which
affects motifs matching a large fraction of the peaks
(see the fpc argument of computeDeviationsAnalytic()).
The object contains two assays: the deviations assay contains the
bias-adjusted deviations from the expectation (by default, the average),
i.e. the difference to the expectation divided by the expectation, and
the z assay contains z-scores, i.e. the difference to the expectation
divided by the variance of the expectation.
If your samples/cells have similar library sizes, it is recommended that
you use the z assay for downstream analysis (such as differential TF
activity using limma). If the samples have very different library
sizes, the z scores will be influenced by that, and it might be
preferable to use the deviations assay.
The betterChromVAR() function is actually a wrapper around three
steps, which can also be executed individually for more customization,
or to avoid repeating some computation multiple times.
The first step is creating the bias bins and their pairwise sampling probabilities. This is achieved with:
bg <- getBackgroundBins(counts)## Creating 50*50=2500 bias bins and computing their sampling distancesbg## bcvBackground object with 500 peaks,
## split into 50*50 ( 2500 ) bins.The second step is computing the per-bin background expectations and variances for each sample or cell:
bg <- computeBackgrounds(counts, bg)
bg## bcvBackground object with 500 peaks,
## split into 50*50 ( 2500 ) bins.
## Background data filled for 10 samples.The bg object has now been filled with the additional information. It
can now be used for:
dev <- computeDeviationsAnalytic(counts, background=bg, annotations=motifMatches)
dev## class: SummarizedExperiment
## dim: 5 10
## metadata(0):
## assays(2): deviations z
## rownames(5): motif1 motif2 motif3 motif4 motif5
## rowData names(8): N Neff ... var.pval var.adjPval
## colnames(10): sample1 sample2 ... sample9 sample10
## colData names(0):This last call could be made with different annotations, without needing to recompute the previous steps.
In addition to the enrichment and GC bias taken into account by the
original chromVAR, betterChromVAR supports
an optional third dimension: fragment length bias (see
?getBackgroundBins() for more information). This requires the
compilation of an additional bias component (the log10-transformed mean
or median fragment length per region). If this is desired, on option is to use
the counting functions from the
epiwraps package, which can
provide this information.
If you are using betterChromVAR (or ChromVAR, for that matter) on
single-cell data with multiple cell types, the way to best run it
depends on whether your are interested in differences within cell types,
or between cell types. If interested in differences between cell types,
run it on the entire dataset, specifying the cell types as grouping
argument to betterChromVAR(). In this way, rare and abundant cell
types will be given the same weight in computing the expectation.
If you are interested in differences within cell types (e.g. between samples/conditions), you should instead run the method separately for each cell type. In this way, the background bins will be defined based on the enrichment in the given cell type, leading to better capture of the bias. The drawback, however, is that the values won’t be comparable across cell types. To test across conditions, we also recommend using it on pseudobulk data.
The only step of the process that needs to be done across the entire
dataset is the computing of the expectation, i.e. the mean (or weighted
mean) for each region across all cells, and the creation of the
background bins, which depends on it. Everything else can easily be
executed in chunks (of cells), as is done for multi-threading in the
betterChromVAR() function.
ex <- getExpectation(counts)
bg <- getBackgroundBins(ex, bias=rowData(counts)$bias)## Creating 50*50=2500 bias bins and computing their sampling distances# this is equivalent to bg <- getBackgroundBins(counts)Then we can apply the next steps only on subset of the data:
bg2 <- computeBackgrounds(counts[,1:3], bg, expectation = ex)
dev2 <- computeDeviationsAnalytic(counts[,1:3], bg2, motifMatches)
# this should be identical to what we had run on the whole object:
identical(assay(dev)[,1:3], assay(dev2))## [1] TRUEIf data is on disk, rather than in memory (and you want it to stay that
way), make sure to use this procedure. In recent versions,
getExpectation() should use chunk-based processing. However it is
quite likely that better performance for this step would be obtained
using for instance using the aggregateAcrossCells() function of the
scrapper package.
The variability of each motif across the dataset is stored in the
rowData of the output of betterChromVAR() :
rowData(dev)## DataFrame with 5 rows and 8 columns
## N Neff total variability var.lower var.upper var.pval
## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
## motif1 68 48.1480 34689 1.321291 0.966496 2.24667 0.0358041
## motif2 71 48.1501 36188 0.909367 0.603440 1.70446 0.4788213
## motif3 73 52.2140 42142 1.253822 0.851014 2.29761 0.0928770
## motif4 83 55.0115 48543 1.210001 0.909765 2.00163 0.0681323
## motif5 69 48.5535 35453 1.023475 0.756641 1.72189 0.2641557
## var.adjPval
## <numeric>
## motif1 0.154795
## motif2 0.478821
## motif3 0.154795
## motif4 0.154795
## motif5 0.330195This is comparable to the variability produced by the original chromVAR,
however the confidence intervals are analytic estimates. While these are
generally good enough, it is also possible to obtain bootstrap confidence
intervals (as in the original approach) through the
computeMotifVariability() function.
An alternative way to assess how useful or informative motifs are is how
reliably they capture biological variation. This can be estimated by the
proportion of their variance across cells that is not attributable to sampling
noise. The estimateNoiseKappa() function can be used to assess this:
motifRel <- estimateNoiseKappa(counts, annotation=motifMatches, background=bg)## Counts assumed to not be paired insertions (`pairs=FALSE`).head(motifRel)## kappa reliability reliability.se
## motif1 0.09709886 0.9683240 0.022692364
## motif2 0.07961262 0.9462814 0.063164628
## motif3 0.04245781 0.9843512 0.008054397
## motif4 0.10037171 0.9608653 0.030952176
## motif5 0.07394526 0.9596145 0.031323623While in practice reliability and variability of a motif are correlated, they do not track the same thing, even in this toy example:
cor(rowData(dev)$variability, motifRel$reliability)## [1] 0.7837231Two motifs can have the same reliability while having very different variability, or the converse. Both are therefore relevant information in selecting useful motifs.
The deviations of individual cells are very noisy: for most motifs, almost all
of the cell-to-cell variation is sampling noise (as estimated with
estimateNoiseKappa()). Working at the pseudobulk level addresses this issue;
alternatively, it is possible to use the k nearest neighbors (kNN) of a cell
to ‘smooth’ its estimates. smoothDeviationsKNN() does exactly this.
Note that the kNN used for smoothing should not be based on the motif
deviations, but can for example be computed from some embedding of the peak
counts (e.g. LSI or similar). smoothDeviationsKNN() shrinks the motif
deviations of each cell towards that of its neighbors in an adaptive fashion:
cells whose deviations are more trustworthy (e.g. because the cell has more
fragments) and motifs with higher reliability are smoothed less. The smoothing
assumes that neighboring cells have similar motif activities: effects that are
specific to individual cells and small relative to the noise are partly
smoothed out (the tauScale argument can be used as a knob to adjust how
aggressive the smoothing is).
# for the sake of the example we do a very basic embedding, using a proper
# scATAC-seq embedding is recommended
emb <- prcomp(t(log1p(assay(counts))))$x[,1:3]
knn <- BiocNeighbors::findKNN(emb, k=3)
dev <- smoothDeviationsKNN(dev, knn)
assayNames(dev)## [1] "deviations" "z" "smoothed" "smoothed.z"The smoothed deviations and z-scores are added to the object as the assays. Note that the smoothed values of neighboring cells are not independent observations anymore, and as such they should not be used in tests that treat cells as independent (the kind of tests you should arguably never do anyway).
The general chromVAR approach, and in particular the deterministic
version implemented here, is also amenable to be used to normalize GC-
and enrichment bias out of bulk ATAC-seq data. This is implemented in
the CVnorm() function. In addition, if given a grouping of the samples
(e.g. experimental conditions), the function applies a variance-based
smoothing of the bias correction, inspired by smoothed quantile
normalization (see the qsmooth package or
Hicks et al. (2018) ). In a nutshell, if the bias in a certain background
bin is explained by experimental groups, it will be less corrected than
if it varies across samples of the groups6 Note that this is different from the original
qsmooth approach – see ?CVnorm for more detail..
Example usage:
# we assign arbitrary groups to the samples:
counts$group <- rep(LETTERS[1:2], each=5)
# we run the smoothed CVnorm:
counts <- CVnorm(counts, grouping=counts$group)
# (the normal CVnorm could be run by omitting the grouping)
counts## class: RangedSummarizedExperiment
## dim: 500 10
## metadata(0):
## assays(2): counts corrected
## rownames: NULL
## rowData names(1): bias
## colnames(10): sample1 sample2 ... sample9 sample10
## colData names(1): groupThis adds a corrected assay to the object. Note that although the
assay is corrected for enrichment- and GC-bias (and fragment length bias if
that was provided), it is not corrected for library size differences. Rather,
it is on the original count scale (although not integer anymore), so that it is
amenable to use in downstream count-based analysis methods such as
edgeR.
The scPRINTER python package developed by Ruochi Zhang in the Buenrostro lab includes a GPU-based reimplementation of , with an important modification: rather than cutting the bias space into bins, a continuous, multidimensional bias space space is created from which \(k\) nearest neighbors (kNN) are selected as background. This approach scales nicely to the inclusion of more dimensions of bias. In scPRINTER, the deviations and variance across the \(k\) neighbors is then computed through permutations, as in the classical (although sped up by the GPU).
If a more complex multi-dimensional bias needs to be taken into account, betterChromVAR implements a similar algorithm that combines our analytical framework with the kNN-based approach. An additional feature is that rather than taking the transformed bias space directly to compute distances, we weight the different dimensions by their correlation with peak overdispersion.
The kNN version of betterChromVAR can be used in the following way:
bg <- getBackgroundKNN(counts)
dev <- computeDeviationsFromKNN(counts, bg, motifMatches)## Computing motif expectations and variances## Calculating final deviations and Z-scoresWe note, however, that this approach is considerably slower (although still
substantially faster than the original chromVAR or the ArchR re-implementation),
requires more memory, and, at the single-cell level, produces slightly more
variable deviations across replicate cells (Wang et al. 2026).
We therefore recommend using the more scalable bin-based variant, unless
working with bulk data with fragment length bias (or any other extra bias
dimensions, see the bias argument of getBackgroundKNN()).
## R version 4.6.1 Patched (2026-06-24 r90190)
## Platform: x86_64-apple-darwin20
## Running under: macOS Ventura 13.7.8
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.6-x86_64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.6-x86_64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
##
## locale:
## [1] C/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: America/New_York
## tzcode source: internal
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] betterChromVAR_1.1.12 SummarizedExperiment_1.43.0
## [3] Biobase_2.73.2 GenomicRanges_1.65.4
## [5] Seqinfo_1.3.2 IRanges_2.47.5
## [7] S4Vectors_0.51.10 BiocGenerics_0.59.12
## [9] generics_0.1.4 MatrixGenerics_1.25.0
## [11] matrixStats_1.5.0 BiocStyle_2.41.0
##
## loaded via a namespace (and not attached):
## [1] Matrix_1.7-6 jsonlite_2.0.0
## [3] crayon_1.5.3 compiler_4.6.1
## [5] BiocManager_1.30.27 Rcpp_1.1.2
## [7] DelayedMatrixStats_1.35.0 Biostrings_2.81.9
## [9] parallel_4.6.1 jquerylib_0.1.4
## [11] BiocParallel_1.47.0 yaml_2.3.12
## [13] fastmap_1.2.0 lattice_0.23-1
## [15] R6_2.6.1 XVector_0.53.0
## [17] S4Arrays_1.13.2 knitr_1.52
## [19] DelayedArray_0.39.8 bookdown_0.48
## [21] bslib_0.12.0 BiocNeighbors_2.7.3
## [23] rlang_1.3.0 cachem_1.1.0
## [25] xfun_0.62 sass_0.4.10
## [27] otel_0.2.0 SparseArray_1.13.4
## [29] cli_3.6.6 digest_0.6.39
## [31] grid_4.6.1 sparseMatrixStats_1.25.0
## [33] lifecycle_1.0.5 evaluate_1.0.5
## [35] codetools_0.2-20 abind_1.4-8
## [37] rmarkdown_2.32 tools_4.6.1
## [39] htmltools_0.5.9