Getting started with bgnormR
Introduction
bgnorm (Kharbanda and Tubelleza et al. 2026) is a generative statistical framework for background correction, normalisation, and quality control of multiplex spatial proteomics data, including data from the Akoya PhenoCycler-Fusion (formerly CODEX), Cell DIVE, IMC, and CosMx platforms.
The framework models log2-transformed fluorescence intensities as a three-component Gaussian Mixture Model (GMM):
| Component | Biological interpretation |
|---|---|
| 1 | Background (empty space, instrument noise) |
| 2 | Non-specific binding / autofluorescence |
| 3 | True biological signal |
Three-component signal model: each pixel intensity X is the sum of independent latent sources U1 (background), U2 (non-specific), and U3 (signal).
A closed-form deconvolution step isolates the signal component from background sources. An optional quantile normalisation step (bgnormQ) unifies the dynamic range across markers, samples, and 3-D tissue slices, enabling a single fixed positivity threshold across the dataset.
The Jensen-Shannon Divergence (JSD) between components 2 and 3 provides an automated quality control metric: higher JSD indicates better staining quality.
Example dataset
The package includes a cropped head-and-neck cancer (HNC) tissue
section imaged on an Akoya PhenoCycler-Fusion platform. The file ships
at inst/extdata/PA_HNC_sample.ome.tiff and contains five
markers across a 550 × 800 pixel region.
path <- system.file("extdata", "PA_HNC_sample.ome.tiff", package = "bgnormR")
img <- read_qptiff(path)
img
#> QPTIFFImage
#> Dimensions: 550 x 800 (H x W)
#> Channels : 5
#> Names : CD20, CD3e, CD8, PanCK, Vimentin
#> Format : ome_tiff
#> Levels : 1read_qptiff() parses the channel names from the embedded
metadata, so names(img) returns the protein names
directly:
names(img) # five-plex panel
#> [1] "CD20" "CD3e" "CD8" "PanCK" "Vimentin"
dim(img) # height × width × channels
#> [1] 550 800 5The same function also reads standard OME-TIFF and
OME-Zarr (OME-NGFF) images, auto-detecting the
container from the path (OME-Zarr requires the Rarr
package). Rich per-channel and slide metadata is organised into an
OME-style QPTIFFMetadata object, available via
metadata() and summarised as a tidy table by
channel_table():
meta <- metadata(img)
qpi_format(meta) # e.g. "fusion_paged", "ome_tiff", "ome_zarr"
#> [1] "ome_tiff"
channel_table(meta)[, c("name", "fluorophore", "exposure_time_us")]
#> name fluorophore exposure_time_us
#> 1 CD20 CY5 150000
#> 2 CD3e ATTO550 150000
#> 3 CD8 CY5 150000
#> 4 PanCK AF750 150000
#> 5 Vimentin CY5 150000The native, Java-free TIFF/QPTIFF reader and writer implemented here was translated from the bioio-tifffile fork by Rafael Tubelleza.
Pixel-level normalisation
Raw intensity maps
Before normalisation, visualise the raw 16-bit intensities to confirm
the image loaded correctly. plot_qptiff() renders one panel
per requested channel.
Running bgnorm_pixels
bgnorm_pixels() fits a three-component GMM to each
channel independently. sample_prop controls the fraction of
non-zero pixels used for fitting — the default 0.1 (10 %) is sufficient
for full-resolution Akoya images while keeping runtime short.
res <- bgnorm_pixels(img, sample_prop = 0.1)
res
#> QPTIFFImage
#> Dimensions: 550 x 800 (H x W)
#> Channels : 5
#> Names : CD20, CD3e, CD8, PanCK, Vimentin
#> Format : ome_tiff
#> Levels : 1
#> bgnorm : yes (5 channel(s))res is a QPTIFFImage whose pixel values are
the background-adjusted log2-intensities. Per-channel model
parameters are stored as an attribute and retrieved with
bgnorm_results():
br_panck <- bgnorm_results(res)[["PanCK"]]
print(br_panck)
#> BgnormResult (pixel-level)
#> n = 440000
#> Component means: 0.03 0.163 1.865
#> JSD (QC metric): 0.8809
#> BIC (G=2): -1.68 BIC (G=3): -1.4
#> No signal detected: FALSE
#> Quantile normalised: FALSE
#> Tissue positivity: 67.5%summary() shows the full component table:
Intensity distributions with GMM overlay
plot_distributions() draws the raw
log2-intensity histogram alongside the fitted GMM component
densities.
The three overlaid curves correspond to the background (blue), non-specific (orange), and signal (red) components. Markers with well-separated signal components (PanCK, Vimentin) show a clearly distinct red peak at high intensity.
Background class assignment
plot_pixel_classes() colours each pixel by its MAP
component assignment (background / non-specific / signal).
Adjusted intensity maps
After normalisation, visualise the background-adjusted images. The
adjusted values are stored inside res (the returned
QPTIFFImage), so plot_qptiff(res) displays
them directly — applying the 2^x back-transform so the
display is in linear intensity units.
p_bg <- plot_qptiff(res, markers = c("PanCK", "Vimentin", "CD20"), scale = "marker") +
labs(title = "Adjusted (bgnorm)")
p_bg
Quantile normalisation (bgnormQ)
Setting quantile_norm = TRUE rescales adjusted
intensities so that the 75th percentile of the fitted signal component
equals 1 across all markers. This enables a single positivity threshold
to be applied consistently across markers, samples, and tissue
sections.
res_q <- bgnorm_pixels(img, sample_prop = 0.1, quantile_norm = TRUE)
p_bgq <- plot_qptiff(res_q, markers = c("PanCK", "Vimentin", "CD20"), scale = "marker") +
labs(title = "Adjusted (bgnormQ)")
p_raw
Quality control
JSD summary table
qc_summary() computes the JSD between the non-specific
and signal components for every marker in a single call.
qc_df <- qc_summary(res)
qc_df
#> name jsd prop_signal
#> 1 CD20 0.4650031 0.13427661
#> 2 CD3e 0.4336173 0.11827368
#> 3 CD8 0.6985864 0.06917809
#> 4 PanCK 0.8808861 0.51594629
#> 5 Vimentin 0.7574030 0.44192710JSD values below 0.1 indicate likely failed staining; values between 0.1 and 0.2 warrant manual inspection (Kharbanda and Tubelleza et al. 2026).
Exporting normalised images
write_qptiff() saves a QPTIFFImage as a
multi-page 16-bit TIFF, one page per channel. For bgnorm-adjusted images
the pixel values are back-transformed to a linear scale
(2^x) before writing. Channel names and bgnorm model
parameters are embedded in each page’s ImageDescription tag
so the file remains self-documenting.
out <- file.path(tempdir(), "PA_HNC_bgnorm.tif")
write_qptiff(res, out)
# Read back – channel names are preserved
img_out <- read_qptiff(out)
names(img_out)
#> [1] "CD20" "CD3e" "CD8" "PanCK" "Vimentin"
dim(img_out)
#> [1] 550 800 5Cell-level normalisation
When only per-cell intensity summaries (mean or median intensity per
cell) are available, use bgnorm_cells(), which applies a
two-component GMM (non-specific / signal) rather than the full
three-component pixel model.
Single-marker: bgnorm_cells
set.seed(3)
cell_intensities <- inv_log_transform(
c(rnorm(300, mean = 3, sd = 0.5),
rnorm(100, mean = 7, sd = 0.9)),
cofactor = 150
)
res_cell <- bgnorm_cells(cell_intensities)
print(res_cell)
#> BgnormResult (cell-level)
#> n = 400
#> Component means: 3.022 6.965
#> JSD (QC metric): 0.986
#> No signal detected: FALSE
#> Quantile normalised: FALSE
#> Tissue positivity: 25.1%Integration with SummarizedExperiment
bgnorm_sce() processes a
SummarizedExperiment (or any subclass such as
SingleCellExperiment or SpatialExperiment): it
reads the requested assay, runs bgnorm_cells()
independently for each marker, and stores the adjusted intensities in a
new assay alongside per-marker BgnormResult objects in
metadata().
library(SummarizedExperiment)
set.seed(5)
n_cells <- 200L
n_markers <- 5L
marker_names <- c("CD20", "CD3e", "CD8", "PanCK", "Vimentin")
make_cell_raw <- function(j) {
inv_log_transform(
c(rnorm(round(n_cells * 0.7), 3, 0.5 + 0.1 * j),
rnorm(round(n_cells * 0.3), 7, 0.8 + 0.1 * j)),
cofactor = 150
)
}
raw_mat <- do.call(rbind, lapply(seq_len(n_markers), make_cell_raw))
rownames(raw_mat) <- marker_names
colnames(raw_mat) <- paste0("cell", seq_len(n_cells))
se <- SummarizedExperiment(assays = list(counts = raw_mat))
se <- bgnorm_sce(se, assay.type = "counts", name = "bgnorm")
assayNames(se)
#> [1] "counts" "bgnorm"Per-marker model parameters are accessible via
metadata():
Distributions and QC for cell-level data
plot_distributions() and plot_jsd_heatmap()
both accept a SummarizedExperiment with bgnorm results
stored in metadata().
qc_summary(bgnorm_meta)
#> name jsd prop_signal
#> 1 CD20 0.9892432 0.2964197
#> 2 CD3e 0.9684877 0.2811681
#> 3 CD8 0.9548634 0.2952482
#> 4 PanCK 0.9244953 0.2446019
#> 5 Vimentin 0.8745865 0.2712931Mathematical background
Log2 transform
\[I_{\log} = \log_2\!\left(\frac{I}{c} + 1\right), \quad c = 150\]
Three-component GMM
\[f(X_i) = \sum_{j=1}^{3} \pi_j \, \mathcal{N}(X_i \mid \mu_j, \sigma_j^2)\]
Components are ordered so that \(\mu_1 < \mu_2 < \mu_3\) (background, non-specific, signal).
Background deconvolution
\[X_{\text{adj}}(x) = P(C\!=\!3 \mid x) \left[ (\mu_3 - \mu_2) + \frac{\sigma_3^2 + \sigma_2^2 - 2\min(\sigma_2^2, \sigma_3^2)}{\sigma_3^2} (x - \mu_3) \right]\]
Quantile normalisation factor
\[q_{\text{norm}} = P(C\!=\!3 \mid q_{0.75}) \left[ (\mu_3 - \mu_2) + \frac{\sigma_3^2 + \sigma_2^2 - 2\min(\sigma_2^2,\sigma_3^2)}{\sigma_3^2} (q_{0.75} - \mu_3) \right]\]
Jensen-Shannon Divergence
\[\text{JSD}(P \| Q) = \frac{1}{2} D_{\text{KL}}(P \| M) + \frac{1}{2} D_{\text{KL}}(Q \| M), \quad M = \frac{P + Q}{2}\]
Session information
sessionInfo()
#> R version 4.6.2 (2026-10-29)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so; LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
#> [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
#>
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats4 stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] SummarizedExperiment_1.43.0 Biobase_2.73.2
#> [3] GenomicRanges_1.65.4 Seqinfo_1.3.2
#> [5] IRanges_2.47.5 S4Vectors_0.51.10
#> [7] BiocGenerics_0.59.12 generics_0.1.4
#> [9] MatrixGenerics_1.25.0 matrixStats_1.5.0
#> [11] ggplot2_4.0.3 bgnormR_0.99.5
#> [13] BiocStyle_2.41.0
#>
#> loaded via a namespace (and not attached):
#> [1] gtable_0.3.6 rjson_0.2.23
#> [3] xfun_0.62 bslib_0.12.0
#> [5] lattice_0.23-1 vctrs_0.7.3
#> [7] tools_4.6.2 parallel_4.6.2
#> [9] tibble_3.3.1 pkgconfig_2.0.3
#> [11] Matrix_1.7-6 RColorBrewer_1.1-3
#> [13] S7_0.2.2 lifecycle_1.0.5
#> [15] compiler_4.6.2 farver_2.1.2
#> [17] tiff_0.1-12 codetools_0.2-20
#> [19] htmltools_0.5.9 sys_3.4.3
#> [21] buildtools_1.0.0 sass_0.4.10
#> [23] yaml_2.3.12 pillar_1.11.1
#> [25] jquerylib_0.1.4 BiocParallel_1.47.0
#> [27] SingleCellExperiment_1.35.2 DelayedArray_0.39.8
#> [29] cachem_1.1.0 magick_2.9.1
#> [31] abind_1.4-8 mclust_6.1.3
#> [33] tidyselect_1.2.1 digest_0.6.39
#> [35] dplyr_1.2.1 labeling_0.4.3
#> [37] maketools_1.3.2 fastmap_1.2.0
#> [39] grid_4.6.2 cli_3.6.6
#> [41] SparseArray_1.13.4 magrittr_2.0.5
#> [43] S4Arrays_1.13.2 withr_3.0.3
#> [45] scales_1.4.0 rmarkdown_2.32
#> [47] XVector_0.53.0 otel_0.3.0
#> [49] SpatialExperiment_1.23.0 evaluate_1.0.5
#> [51] knitr_1.52 viridisLite_0.4.3
#> [53] rlang_1.3.0 Rcpp_1.1.2
#> [55] glue_1.8.1 BiocManager_1.30.27
#> [57] xml2_1.6.0 jsonlite_2.0.0
#> [59] R6_2.6.1