Get started in Python

DenMark measures how gene expression depends on local cell density in single-cell resolution spatial transcriptomics data. Cell locations form a point pattern and the counts of a gene are its marks. After the cells are binned on a regular grid over the tissue, the DenMark model for grid cell $g$ with area $|g|$ inside the tissue, cell count $Y(g)$ and summed expression $M(g)$ is

$$\begin{aligned} Y(g) &\sim \text{Poisson}\{\lambda_1(g)\}, & \log\lambda_1(g) &= \log|g| + \beta_0 + a_{11}w_1(g),\\ M(g)\mid Y(g) &\sim \text{Poisson}\{\lambda_2(g)\}, & \log\lambda_2(g) &= \beta_1 + \log\{Y(g)+0.05\} + a_{21}w_1(g) + a_{22}w_2(g), \end{aligned}$$

where $w_1$ and $w_2$ are independent unit-variance Gaussian processes with Matern-3/2 covariance, approximated by Hilbert space basis functions (HSGP). With $a_{11} = \sigma_1$, $a_{21} = \sigma_2\rho$ and $a_{22} = \sigma_2\sqrt{1-\rho^2}$, the parameter $\rho$ (Astar[2,1] in Stan) is the density-expression correlation.

This page runs the workflow on the bundled toy data, where the true correlations are known, so every estimate can be checked against the truth. The Python package denmark builds the same Stan data as the R package and fits the same Stan models through cmdstanpy; Get started in R walks through the same steps in R, and R and Python side by side lists the function names of both. The code below is shown without its output, and the package tests run it as written. To try it without installing anything, open the same workflow in Google Colab:

Installation

DenMark fits its models with Stan through cmdstanpy, which needs a working CmdStan installation, and computes the tissue boundary with shapely. Install them once:

pip install "denmark[boundary,plot] @ git+https://github.com/StaGill/DenMark#subdirectory=python"
python -m cmdstanpy.install_cmdstan

denmark needs Python 3.9 or later. The extra boundary adds shapely (2.0 or later) for the tissue boundary and plot adds matplotlib for the figures. install_cmdstan downloads and builds CmdStan once per machine and needs a C++ toolchain.

The toy data

dm.datasets.load_toy() returns a tissue section simulated by denmark_simulate() in R: 7,449 cells in an irregular, lobed section inside a 6 x 4 window, one row per cell, with the cell type and the counts of three genes. The 3,449 Tumor cells are the population analyzed; their true density-expression correlations are 0.6 for geneA, 0 for geneB and -0.6 for geneC. The 4,000 Other cells only delineate the tissue.

import matplotlib.pyplot as plt
import pandas as pd
import denmark as dm

toy = dm.datasets.load_toy()
print(toy.head())
print(toy["celltype"].value_counts())
tumor = toy[toy["celltype"] == "Tumor"].drop(columns="celltype")
truth = {"geneA": 0.6, "geneB": 0.0, "geneC": -0.6}

fig, ax = plt.subplots(figsize=(7, 3.6))
points = ax.scatter(tumor["x"], tumor["y"], c=tumor["geneA"], s=2, cmap="magma")
ax.set_aspect("equal")
fig.colorbar(points, ax=ax, label="geneA")

Tissue boundary

A tissue section is rarely a rectangle. If the grid covered the bounding box of the cells, grid cells outside the tissue would hold no cells because there is no tissue there, and the model would read them as regions of low cell density. dm.alpha_shape() delineates the tissue as an alpha shape: it triangulates the cell locations and keeps every Delaunay triangle whose circumradius is below $\alpha$, by default three times the median Delaunay edge length. The boundary is computed from all cells of the section, not only the Tumor cells, because the region where a Tumor cell could have been observed is the whole tissue. The result holds the polygon (.polygon), the alpha (.alpha) and the tissue area (.area).

tissue = dm.alpha_shape(toy["x"], toy["y"])
print(tissue.alpha, tissue.area)

Gridding

dm.buildgridpp() assigns every cell to one grid cell and returns the cell counts $Y(g)$ and the summed expression $M(g)$. Grid cells are ordered with x varying fastest. With boundary, the grid spans the extent of the tissue, and each grid cell is classified as interior (entirely inside the tissue), boundary (partly inside) or exterior, with its area inside the tissue in the column area.

gridded = dm.buildgridpp(10, tumor, gene="geneA", boundary=tissue)
print(gridded["grid"].head())
print(gridded["grid"]["type"].value_counts())
fig = dm.plot_boundary(toy, tissue, grid=gridded, highlight=toy["celltype"] == "Tumor")

Exterior grid cells are dropped from the model, and boundary grid cells enter with their area inside the tissue, so the area offset $\log|g|$ differs between grid cells. dm.plot_grid() shows the gridded counts inside the tissue.

fig = dm.plot_grid(gridded)

dm.prepare_data() builds the Stan input that dm.DenMark() uses: it drops the exterior grid cells, centers the region of interest at the origin, uses the grid centroids as HSGP coordinates and sets the HSGP domain to boundfactor times the half-ranges. With a boundary, log_grid_area holds one value per grid cell in the model.

prep = dm.prepare_data(tumor, gene="geneA", grid_size=10, basis=10, boundfactor=1.5,
                       boundary=tissue)
stan_data = prep["stan_data"]
print({k: stan_data[k] for k in ("n", "mstar", "L", "mrange")})
print(stan_data["log_grid_area"][:5])

Fitting the model

dm.DenMark() grids the data, compiles the Stan model on first use and samples; the compiled model is cached in a user directory (set by the environment variable DENMARK_CACHE_DIR, by default ~/.cache/denmark), so later sessions reuse it. With boundary, it fits the variable-area model DenMark-vararea.stan, the same model with one area offset per grid cell; without it, the bounding box of the cells is gridded with a common grid area. The settings below keep the page fast; for real data, use the defaults (four chains of 1,000 warm-up and 1,000 sampling iterations and basis=25) and check the diagnostics.

fit = dm.DenMark(tumor, gene="geneA", boundary=tissue, grid_size=10, basis=10,
                 boundfactor=1.5, chains=2, iter_warmup=500, iter_sampling=500,
                 seed=2024, quiet=True)
print(fit)

Posterior summaries

fit.summary() reports the intercepts, the length-scales ell[1] and ell[2], $\rho$ and the field standard deviations, with R-hat and effective sample sizes.

print(fit.summary()[["variable", "mean", "sd", "q2.5", "q97.5", "rhat", "ess_bulk"]])

fit.draws() returns one row per posterior draw, including $\rho$. The red line marks the true value.

draws = fit.draws()
fig, ax = plt.subplots(figsize=(5, 3))
ax.hist(draws["rho"], bins=40, density=True, color="0.8")
ax.axvline(0, linestyle="--", color="0.3")
ax.axvline(truth["geneA"], color="#B2182B")
ax.set_xlim(-1, 1)
ax.set_xlabel("rho")
ax.set_ylabel("Posterior density")

Latent fields

fit.fields() decomposes, for every grid cell in the model, the log cell intensity $\beta_0 + a_{11}w_1$ and the log expression per cell $\beta_1 + a_{21}w_1 + a_{22}w_2$ into the density field $a_{11}w_1$, the density-shared part of expression $a_{21}w_1$ and the expression-specific part $a_{22}w_2$.

fields = fit.fields()
print(fields[["grid_index", "Y", "M", "eta1", "a11w1", "eta2", "a21w1", "a22w2"]].head())
fig = dm.plot_fields(fit, fields=fields)

Spatial correlation ranges

fit.decay() gives the posterior spatial correlation of the density and expression fields as a function of distance, and the practical range, the distance at which the correlation falls to 0.05.

decay = fit.decay()
print(decay["ranges"])
fig = dm.plot_decay(fit, decay=decay)

Scanning genes for density-correlated genes

dm.scan_genes() fits DenMark to each gene and calls it a positively (PCG) or negatively (NCG) density-correlated gene when the 95% credible interval of $\rho$ lies above or below 0, and not a DCG otherwise.

scan = dm.scan_genes(tumor, genes=["geneA", "geneB", "geneC"], boundary=tissue,
                     grid_size=10, basis=10, boundfactor=1.5, chains=2,
                     iter_warmup=500, iter_sampling=500, seed=2024, verbose=False)
print(scan[["gene", "rho_mean", "rho_lower", "rho_upper", "class", "max_rhat"]])
fig = dm.plot_scan(scan)

For a full panel, n_jobs=4 fits four genes at a time, with the same results as the sequential scan; every fit also runs its chains in parallel, so plan for n_jobs times chains cores. On a cluster, run one gene per SLURM array task with dm.scan_genes(..., genes=[gene], save_dir=...), then collect the table with the same call on all genes, which loads the saved results instead of refitting; the R package ships a template of this workflow (system.file("slurm", package = "DenMark")).

Checking the estimates against the truth

With true values of 0.6, 0 and -0.6, the expected calls are PCG for geneA, not DCG for geneB and NCG for geneC. The table compares each posterior mean and 95% credible interval with the true $\rho$; covered shows whether the interval contains it. Reading density-correlated genes explains what a call can and cannot tell you, and how to tell a change within the focal cells from a signal of their neighborhood.

check = pd.DataFrame({"gene": scan["gene"], "truth": scan["gene"].map(truth),
                      "estimate": scan["rho_mean"], "lower": scan["rho_lower"],
                      "upper": scan["rho_upper"], "class": scan["class"]})
check["covered"] = (check["lower"] <= check["truth"]) & (check["truth"] <= check["upper"])
print(check)

Using your own data

DenMark needs one row per cell with the coordinates in columns x and y and raw counts in one column per gene. Compute the tissue boundary from all cells of the section, then analyze the cells of one population: filter uninformative genes, fit one gene, then scan the panel; with save_dir, each gene is saved as it finishes and an interrupted scan resumes where it stopped.

cells = pd.read_csv("cells.csv")   # all cells: x, y, celltype and gene counts
# or from an AnnData object (scanpy, or sdata.tables["table"] from spatialdata-io):
# cells = dm.from_anndata(adata, celltype_key="celltype")
tissue = dm.alpha_shape(cells["x"], cells["y"])
focal = cells[cells["celltype"] == "Tumor"].drop(columns="celltype")
genes = dm.filter_genes(focal.drop(columns=["x", "y"]))
fit = dm.DenMark(focal[["x", "y", genes[0]]], gene=genes[0], boundary=tissue,
                 grid_size=15, basis=25, boundfactor=1.5)
dcg = dm.scan_genes(focal[["x", "y"] + genes], boundary=tissue, grid_size=15, basis=25,
                    save_dir="denmark_scan")

dm.annotate(adata, dcg) writes the estimates, their intervals and the calls back into the AnnData object’s gene table (adata.var). fit.check_hsgp() compares the posterior length-scales with a rule of thumb for the number of basis functions and the boundary factor, and dm.save_result() keeps a fit for later sessions.