12. Cell-state TF activity with chromVAR¶

Bolero's score model (next two pages) conditions accessibility prediction on a chromVAR TF-activity score — a per-cell-state number saying how active a transcription-factor family (e.g. AP-1) is. This page shows how that score is computed: a GPU chromVAR run that turns a single-cell peak × cell ATAC matrix into a motif × cell matrix of activity deviations (bias-corrected, background-subtracted z-scores).

Everything here is bolero-native — bolero.tl.chromvar for the chromVAR math and bolero.tl.motif for motif matching — so no external chromVAR/scPrinter install is needed.

This is a data-lake notebook (like the multi-dataset and variant-effect pages): it reads a dataset's cell-level peak matrices from bolerodata.DATASETS and runs on a GPU (cupy). It is rendered from baked outputs; the docs build never re-runs it.

We demo one dataset (ChromiumPBMC). The final section explains the cross-dataset chromVAR that actually feeds the score model, and why it differs.

Setup¶

In [1]:

import pathlib

import anndata
import numpy as np
import pandas as pd

from bolero import Genome
from bolero.tl.chromvar.chromvar import (
    get_peak_bias,  # per-peak GC content + background nucleotide frequency
    scan_peak_motifs,  # peak x motif match matrix (via bolero.tl.motif.Motifs)
    sample_bg_peaks,  # GC/depth-matched background peaks
    compute_deviations,  # the GPU chromVAR deviation z-scores
)
from bolerodata import DATASETS

import pathlib import anndata import numpy as np import pandas as pd from bolero import Genome from bolero.tl.chromvar.chromvar import ( get_peak_bias, # per-peak GC content + background nucleotide frequency scan_peak_motifs, # peak x motif match matrix (via bolero.tl.motif.Motifs) sample_bg_peaks, # GC/depth-matched background peaks compute_deviations, # the GPU chromVAR deviation z-scores ) from bolerodata import DATASETS

The input: a dataset's cell-level peak matrices¶

bolerodata.DATASETS points at each dataset's peak AnnData — one file per sample, each a sparse cell × peak count matrix. chromVAR is run on these raw single cells (not on meta-cells / pseudobulks).

In [2]:

dataset_name = "ChromiumPBMC"
dataset = DATASETS[dataset_name]
genome = Genome(dataset.genome)

adata_paths = sorted(pathlib.Path(dataset.peak_adata_path).glob("*.h5ad"))
print(dataset_name, "genome:", dataset.genome, "| n sample files:", len(adata_paths))

# peek at one: cells x peaks, peaks named "chrom:start-end"
_a = anndata.read_h5ad(adata_paths[0], backed="r")
print("example shape (cells x peaks):", _a.shape)
print("example peaks:", list(_a.var_names[:3]))

dataset_name = "ChromiumPBMC" dataset = DATASETS[dataset_name] genome = Genome(dataset.genome) adata_paths = sorted(pathlib.Path(dataset.peak_adata_path).glob("*.h5ad")) print(dataset_name, "genome:", dataset.genome, "| n sample files:", len(adata_paths)) # peek at one: cells x peaks, peaks named "chrom:start-end" _a = anndata.read_h5ad(adata_paths[0], backed="r") print("example shape (cells x peaks):", _a.shape) print("example peaks:", list(_a.var_names[:3]))

ChromiumPBMC genome: hg38 | n sample files: 4
example shape (cells x peaks): (10599, 227031)
example peaks: ['chr1:816081-816581', 'chr1:817090-817590', 'chr1:821042-821542']

1 · Peak read depth, and peak filtering¶

chromVAR needs a per-peak read depth (the total accessibility of each peak across all cells). We sum every sample's counts, then drop the lowest and highest peaks (very rare or blacklist-like) by quantile. This filtered depth doubles as the peak-level expectation (expectation_var) chromVAR subtracts as background.

In [3]:

peak_total_sum = None
for path in adata_paths:
    adata = anndata.read_h5ad(path)
    peak_sum = np.asarray(adata.X.sum(axis=0)).reshape(-1)
    peak_total_sum = peak_sum if peak_total_sum is None else peak_total_sum + peak_sum
peak_total_sum = pd.Series(peak_total_sum, index=adata.var_names)

# keep the informative middle of the depth distribution
min_cov, max_cov = peak_total_sum.quantile([0.05, 0.997])
peak_total_sum = peak_total_sum[(peak_total_sum > min_cov) & (peak_total_sum < max_cov)]
print(f"{peak_total_sum.size:,} peaks kept after depth filtering")

peak_total_sum = None for path in adata_paths: adata = anndata.read_h5ad(path) peak_sum = np.asarray(adata.X.sum(axis=0)).reshape(-1) peak_total_sum = peak_sum if peak_total_sum is None else peak_total_sum + peak_sum peak_total_sum = pd.Series(peak_total_sum, index=adata.var_names) # keep the informative middle of the depth distribution min_cov, max_cov = peak_total_sum.quantile([0.05, 0.997]) peak_total_sum = peak_total_sum[(peak_total_sum > min_cov) & (peak_total_sum < max_cov)] print(f"{peak_total_sum.size:,} peaks kept after depth filtering")

214,888 peaks kept after depth filtering

2 · Build the shared peak data — bias, background peaks, motif matches¶

chromVAR corrects each observed motif accessibility against a set of background peaks matched for GC content and read depth. So once, for the whole peak set, we compute:

1. get_peak_bias → per-peak GC content (var["gc_content"]) and the overall background base frequency (uns["bg_freq"]), read straight from the genome FASTA.
2. sample_bg_peaks → for every peak, a panel of niterations background peaks with matched GC + depth.
3. scan_peak_motifs → a peak × motif binary match matrix (varm["motif_match"]) using bolero's built-in motif scanner in motifmatchr mode. The JASPAR 2024 vertebrate PFMs are fetched and cached automatically (bolero.tl.motif.jaspar.get_jaspar_motif_file); columns are JASPAR matrix ids.

We store the filtered depth as expectation_var — the background expectation compute_deviations will use.

In [4]:

# load one sample just to carry the peak-level annotations; obs does not matter here
example_adata = anndata.read_h5ad(adata_paths[0])
example_adata = example_adata[
    :, example_adata.var_names.isin(peak_total_sum.index)
].copy()
peak_total_sum = peak_total_sum.reindex(example_adata.var_names)

# 1. GC content + background base frequency (fills var["gc_content"], uns["bg_freq"])
get_peak_bias(example_adata, genome)

# 2. GC/depth-matched background peaks
peak_bg = sample_bg_peaks(
    reads_per_peak=peak_total_sum.values,
    gc_bias=example_adata.var["gc_content"].values,
    niterations=250,
    w=0.1,
    bs=50,
)

# 3. peak x motif matches (fills varm["motif_match"], uns["motif_name"])
scan_peak_motifs(example_adata, genome, pvalue=5e-5, n_jobs=16)

# the per-peak background expectation chromVAR subtracts
example_adata.var["expectation_var"] = peak_total_sum

# bundle everything the per-sample chromVAR runs need
peak_data = {
    "uns": example_adata.uns,
    "varm.motif_match": example_adata.varm["motif_match"],
    "varm.bg_peaks": peak_bg,
    "var": example_adata.var,
}
print("motif_match:", peak_data["varm.motif_match"].shape, "(peaks x motifs)")
print(
    "n motifs:",
    len(peak_data["uns"]["motif_name"]),
    "| example motif ids:",
    list(peak_data["uns"]["motif_name"][:3]),
)

# load one sample just to carry the peak-level annotations; obs does not matter here example_adata = anndata.read_h5ad(adata_paths[0]) example_adata = example_adata[ :, example_adata.var_names.isin(peak_total_sum.index) ].copy() peak_total_sum = peak_total_sum.reindex(example_adata.var_names) # 1. GC content + background base frequency (fills var["gc_content"], uns["bg_freq"]) get_peak_bias(example_adata, genome) # 2. GC/depth-matched background peaks peak_bg = sample_bg_peaks( reads_per_peak=peak_total_sum.values, gc_bias=example_adata.var["gc_content"].values, niterations=250, w=0.1, bs=50, ) # 3. peak x motif matches (fills varm["motif_match"], uns["motif_name"]) scan_peak_motifs(example_adata, genome, pvalue=5e-5, n_jobs=16) # the per-peak background expectation chromVAR subtracts example_adata.var["expectation_var"] = peak_total_sum # bundle everything the per-sample chromVAR runs need peak_data = { "uns": example_adata.uns, "varm.motif_match": example_adata.varm["motif_match"], "varm.bg_peaks": peak_bg, "var": example_adata.var, } print("motif_match:", peak_data["varm.motif_match"].shape, "(peaks x motifs)") print( "n motifs:", len(peak_data["uns"]["motif_name"]), "| example motif ids:", list(peak_data["uns"]["motif_name"][:3]), )

Sampling nearest neighbors

NNDescent (2500, 2)

motif_match: (214888, 879) (peaks x motifs)
n motifs: 879 | example motif ids: ['MA0002.3', 'MA0003.5', 'MA0004.1']

3 · Run chromVAR per sample¶

compute_deviations runs on the GPU. For each observed motif it compares the cells' accessibility at motif-bearing peaks against the matched background panel, and returns a deviation z-score per (cell, motif). We run it per sample and concatenate.

In [5]:

def run_single_adata(path, peak_data):
    """Run chromVAR on one sample's peak AnnData; return its deviation AnnData."""
    adata = anndata.read_h5ad(path)
    # attach the shared peak annotations; keep only peaks in the filtered set
    for k, v in peak_data["var"].items():
        adata.var[k] = v
    adata = adata[:, ~adata.var["expectation_var"].isna().values].copy()
    adata.uns = peak_data["uns"]
    adata.varm["motif_match"] = peak_data["varm.motif_match"]
    adata.varm["bg_peaks"] = peak_data["varm.bg_peaks"]
    return compute_deviations(adata, chunk_size=20000, device="cuda")


dev_per_sample = [run_single_adata(path, peak_data) for path in adata_paths]
chromvar = anndata.concat(dev_per_sample)
print("chromVAR deviations:", chromvar.shape, "(cells x motifs)")

def run_single_adata(path, peak_data): """Run chromVAR on one sample's peak AnnData; return its deviation AnnData.""" adata = anndata.read_h5ad(path) # attach the shared peak annotations; keep only peaks in the filtered set for k, v in peak_data["var"].items(): adata.var[k] = v adata = adata[:, ~adata.var["expectation_var"].isna().values].copy() adata.uns = peak_data["uns"] adata.varm["motif_match"] = peak_data["varm.motif_match"] adata.varm["bg_peaks"] = peak_data["varm.bg_peaks"] return compute_deviations(adata, chunk_size=20000, device="cuda") dev_per_sample = [run_single_adata(path, peak_data) for path in adata_paths] chromvar = anndata.concat(dev_per_sample) print("chromVAR deviations:", chromvar.shape, "(cells x motifs)")

/large_storage/zhoulab/hanliu/pkg/liuhlab/bolero/src/bolero/tl/chromvar/chromvar.py:344: FutureWarning: The method varm_keys is deprecated and will be removed in the future. Use varm instead of varm_keys. (e.g. `k in adata.varm` or `adata.varm.keys() | {'u'}`)
  "bg_peaks" in adata.varm_keys()

Computing expectation reads per cell and peak...
precomputed expectation_var

Computing deviations...

Finish computing deviations...

/large_storage/zhoulab/hanliu/pkg/liuhlab/bolero/src/bolero/tl/chromvar/chromvar.py:344: FutureWarning: The method varm_keys is deprecated and will be removed in the future. Use varm instead of varm_keys. (e.g. `k in adata.varm` or `adata.varm.keys() | {'u'}`)
  "bg_peaks" in adata.varm_keys()

Computing expectation reads per cell and peak...
precomputed expectation_var

Computing deviations...

Finish computing deviations...

/large_storage/zhoulab/hanliu/pkg/liuhlab/bolero/src/bolero/tl/chromvar/chromvar.py:344: FutureWarning: The method varm_keys is deprecated and will be removed in the future. Use varm instead of varm_keys. (e.g. `k in adata.varm` or `adata.varm.keys() | {'u'}`)
  "bg_peaks" in adata.varm_keys()

Computing expectation reads per cell and peak...
precomputed expectation_var

Computing deviations...

Finish computing deviations...

/large_storage/zhoulab/hanliu/pkg/liuhlab/bolero/src/bolero/tl/chromvar/chromvar.py:344: FutureWarning: The method varm_keys is deprecated and will be removed in the future. Use varm instead of varm_keys. (e.g. `k in adata.varm` or `adata.varm.keys() | {'u'}`)
  "bg_peaks" in adata.varm_keys()

Computing expectation reads per cell and peak...
precomputed expectation_var

Computing deviations...

Finish computing deviations...

chromVAR deviations: (45270, 879) (cells x motifs)

In [6]:

out_path = f"{dataset_name}.chromvar.h5ad"
chromvar.write_h5ad(out_path)
print("wrote", out_path)

out_path = f"{dataset_name}.chromvar.h5ad" chromvar.write_h5ad(out_path) print("wrote", out_path)

wrote ChromiumPBMC.chromvar.h5ad

4 · What chromVAR gives you¶

The result is a cell × motif matrix of activity deviations. Positive means the cells' open chromatin is enriched for that motif above the GC/depth-matched background — i.e. the factor is active. Below, the most variable motifs across cells are exactly the TFs whose activity differs most between cell states.

In [7]:

dev = chromvar.to_df()
top = dev.std().sort_values(ascending=False).head(10)
print("most variable TF motifs (by deviation std across cells):")
top

dev = chromvar.to_df() top = dev.std().sort_values(ascending=False).head(10) print("most variable TF motifs (by deviation std across cells):") top

most variable TF motifs (by deviation std across cells):

Out[7]:

MA0081.3    7.018835
MA0462.3    6.966119
MA1634.2    6.966119
MA1928.2    6.966119
MA1988.2    6.966119
MA0835.3    6.966119
MA0080.7    6.941144
MA0478.2    6.926793
MA0102.5    6.862599
MA1128.2    6.825898
dtype: float32

Cross-dataset chromVAR — what the score model actually uses¶

The demo above is a single-dataset chromVAR: peaks, background depth, and the expectation_var all come from ChromiumPBMC alone. That is perfect for looking at TF activity within one dataset — but the scores are not comparable across datasets, because every dataset would use its own peaks and its own background. To train one score model over the whole Bolero-10M corpus, the conditioning score has to live on a common scale across datasets. So the score model is trained on a cross-dataset chromVAR that differs in three ways:

1. A shared, species-level ensemble peak set. Instead of each dataset's own peaks, all datasets on a genome are quantified on one consensus peak set built by pooling and re-calling peaks across every dataset of that species (a fragment pile-up → summit calling → clean-up pipeline). Every dataset is now measured on the same coordinates.

2. A shared background expectation. The per-peak depth (expectation_var) is pooled across all datasets of the genome (and lightly rescaled), and the background-peak sampling + motif scan are done once per genome. Every dataset's chromVAR then subtracts the same background — which is exactly why compute_deviations reads a precomputed expectation_var (the generalization you saw in step 2) instead of re-deriving it from each dataset's own counts.

3. Pseudobulk-level, standardized scores. chromVAR is run on pseudobulks (not raw cells), and for each TF family a representative JASPAR motif is chosen (AP-1 → MA0478.2) and its deviation is z-scored. The result is a single number per pseudobulk, comparable across every dataset — the __shared_data__ value the score model conditions on.

We don't reproduce that whole multi-genome pipeline here; the point is that it is the same bolero.tl.chromvar machinery applied on shared peaks + a shared background so the resulting AP-1 (or other TF) scores are dataset-comparable. The next page picks up from those precomputed per-pseudobulk scores.