Get started in R

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 vignette runs the workflow on the bundled toy data, where the true correlations are known, so every estimate can be checked against the truth.

Installation

DenMark is installed from the StaGILL R-universe. It fits its models with Stan through cmdstanr, which needs a working CmdStan installation, and computes the tissue boundary with sf. Install them once:

install.packages(c("DenMark", "sf", "patchwork"),
                 repos = c("https://stagill.r-universe.dev",
                           "https://stan-dev.r-universe.dev", getOption("repos")))
cmdstanr::install_cmdstan()

The toy data

denmark_toy_data() returns a tissue section simulated by denmark_simulate(): 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.

library(DenMark)
library(ggplot2)

toy <- denmark_toy_data()
head(toy)
#>          x        y celltype geneA geneB geneC
#> 1 4.958866 2.517796    Tumor     3     1     1
#> 2 3.398667 1.109772    Tumor     2     1     4
#> 3 2.060565 2.566543    Tumor     1     3     3
#> 4 1.463342 2.365366    Tumor     3     0     4
#> 5 5.085220 2.566559    Tumor    10     1     1
#> 6 2.149059 1.976934    Tumor     0     8     1
table(toy$celltype)
#> 
#> Other Tumor 
#>  4000  3449
tumor <- subset(toy, celltype == "Tumor", select = -celltype)
truth <- c(geneA = 0.6, geneB = 0, geneC = -0.6)

ggplot(tumor, aes(x, y, color = geneA)) +
  geom_point(size = 0.6) +
  scale_color_viridis_c(option = "magma") +
  coord_equal() +
  theme_classic()

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. denmark_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.

tissue <- denmark_alpha_shape(toy$x, toy$y)
c(alpha = attr(tissue, "alpha"), area = as.numeric(sf::st_area(tissue)))
#>     alpha      area 
#> 0.1068794 8.5385318

Gridding

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 <- buildgridpp(grid_size = 10, dataset = tumor, gene = "geneA",
                       boundary = tissue)
head(gridded$grid)
#>   grid_index i j         x         y  Y  M       area     type
#> 1          1 1 1 0.8713599 0.7530424  0  0 0.00000000 exterior
#> 2          2 2 1 1.3814337 0.7530424  0  0 0.00000000 exterior
#> 3          3 3 1 1.8915075 0.7530424 12 18 0.02704043 boundary
#> 4          4 4 1 2.4015813 0.7530424 29 63 0.09246043 boundary
#> 5          5 5 1 2.9116551 0.7530424 31 78 0.11468348 boundary
#> 6          6 6 1 3.4217289 0.7530424 14 13 0.09659054 boundary
table(gridded$grid$type)
#> 
#> boundary exterior interior 
#>       36       18       46
plot_denmark_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. plot_denmark_grid() shows the gridded counts inside the tissue.

plot_denmark_grid(gridded)

denmark_prepare() builds the Stan input that 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 <- denmark_prepare(tumor, gene = "geneA", grid_size = 10, basis = 10,
                        boundfactor = 1.5, boundary = tissue)
str(prep$stan_data[c("n", "mstar", "L", "log_grid_area", "mrange")])
#> List of 5
#>  $ n            : int 82
#>  $ mstar        : int 100
#>  $ L            : num [1:2] 3.83 1.91
#>  $ log_grid_area: num [1:82] -3.61 -2.38 -2.17 -2.34 -4.68 ...
#>  $ mrange       : num 0.933

Fitting the model

DenMark() grids the data, compiles the Stan model once per session and samples. 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 vignette 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 <- 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)
fit
#> DenMark fit  |  gene geneA
#>   grid 10 x 10 (82 grids, 3449 cells)  |  basis 10  |  boundfactor 1.50
#>   rho (density-expression correlation): 0.633  [95% CrI 0.355, 0.855]
#>   max R-hat 1.007 | min ESS bulk 503 | divergences 0 | time 0.92 min
#> 
#>    variable     mean    sd   q2.5   q50 q97.5 rhat ess_bulk
#>       beta0  5.59004 0.358  4.837 5.622 6.216 1.00     1267
#>       beta1 -0.00594 0.516 -1.198 0.062 0.845 1.00      605
#>      ell[1]  0.60924 0.144  0.388 0.588 0.943 1.00      638
#>      ell[2]  0.76688 0.253  0.400 0.727 1.348 1.00      597
#>  Astar[2,1]  0.63314 0.133  0.355 0.651 0.855 1.01      514
#>    sigma[1]  1.00114 0.196  0.713 0.971 1.487 1.00      510
#>    sigma[2]  1.34816 0.300  0.899 1.309 2.053 1.00      503

Posterior summaries

denmark_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.

denmark_summary(fit)[, c("variable", "mean", "sd", "q2.5", "q97.5",
                         "rhat", "ess_bulk")]
#>     variable         mean        sd       q2.5     q97.5      rhat  ess_bulk
#> 1      beta0  5.590040466 0.3579061  4.8366445 6.2162162 1.0002666 1267.4047
#> 2      beta1 -0.005944119 0.5159840 -1.1983909 0.8449484 1.0009315  605.2964
#> 3     ell[1]  0.609235540 0.1438940  0.3881686 0.9426524 1.0026561  637.9042
#> 4     ell[2]  0.766881690 0.2531255  0.3997421 1.3484433 0.9996954  596.6830
#> 5 Astar[2,1]  0.633138810 0.1332939  0.3548368 0.8550668 1.0072087  514.2703
#> 6   sigma[1]  1.001143618 0.1963631  0.7132929 1.4869231 1.0011645  510.1265
#> 7   sigma[2]  1.348159905 0.2996388  0.8986182 2.0530672 1.0004490  502.7491

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

draws <- denmark_draws(fit)
ggplot(draws, aes(rho)) +
  geom_density(fill = "grey80") +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_vline(xintercept = truth[["geneA"]], color = "#B2182B") +
  xlim(-1, 1) +
  labs(x = expression(rho), y = "Posterior density") +
  theme_classic()

Latent fields

denmark_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 <- denmark_fields(fit)
head(fields[, c("grid_index", "Y", "M", "eta1", "a11w1", "eta2", "a21w1", "a22w2")])
#>   grid_index  Y  M     eta1      a11w1       eta2      a21w1      a22w2
#> 1          3 12 18 6.078127  0.4880869 0.50803381  0.4180511 0.09592686
#> 2          4 29 63 5.718201  0.1281600 0.81180341  0.1143760 0.70337154
#> 3          5 31 78 5.417882 -0.1721584 0.91561175 -0.1400866 1.06164252
#> 4          6 14 13 5.038056 -0.5519845 0.04724763 -0.4590707 0.51226249
#> 5          7  1  1 5.449936 -0.1401043 0.04575002 -0.1123506 0.16404478
#> 6         12  8 15 6.300830  0.7107892 0.74736645  0.6047756 0.14853492
plot_denmark_fields(fit, fields = fields)

Spatial correlation ranges

denmark_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 <- denmark_decay(fit)
decay$ranges
#>        field     mean   median    lower    upper
#> 1    density 1.668618 1.609570 1.063144 2.581804
#> 2 expression 2.066007 1.952302 1.308894 3.449081
plot_denmark_decay(fit, decay = decay)

Scanning genes for density-correlated genes

denmark_scan() 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 <- denmark_scan(tumor, genes = c("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)
scan[, c("gene", "rho_mean", "rho_lower", "rho_upper", "class", "max_rhat")]
#>    gene   rho_mean  rho_lower  rho_upper   class max_rhat
#> 1 geneA  0.6331388  0.3548368  0.8550668     PCG 1.007209
#> 2 geneB  0.1377871 -0.2435719  0.4730265 not DCG 1.007989
#> 3 geneC -0.5949906 -0.8136031 -0.3324967     NCG 1.005691

plot_denmark_scan(scan)
#> `height` was translated to `width`.

For a full panel, workers = 4 fits four genes at a time, each in its own R session, with the same results as the sequential scan; every fit also runs its chains in parallel, so plan for workers times chains cores. On a cluster, the template in system.file("slurm", package = "DenMark") fits one gene per SLURM array task and then collects the results with denmark_scan().

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 <- data.frame(gene = scan$gene, truth = truth[scan$gene],
                    estimate = scan$rho_mean, lower = scan$rho_lower,
                    upper = scan$rho_upper, class = scan$class,
                    row.names = NULL)
check$covered <- check$lower <= check$truth & check$truth <= check$upper
check
#>    gene truth   estimate      lower      upper   class covered
#> 1 geneA   0.6  0.6331388  0.3548368  0.8550668     PCG    TRUE
#> 2 geneB   0.0  0.1377871 -0.2435719  0.4730265 not DCG    TRUE
#> 3 geneC  -0.6 -0.5949906 -0.8136031 -0.3324967     NCG    TRUE

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 <- read.csv("cells.csv")   # all cells: x, y, celltype and gene counts
# or from a Seurat object (Xenium, MERSCOPE or CosMx) with cell types in its metadata:
# cells <- denmark_from_seurat(obj, celltype = "celltype")
tissue <- denmark_alpha_shape(cells$x, cells$y)
focal <- subset(cells, celltype == "Tumor", select = -celltype)
genes <- denmark_filter_genes(focal[, setdiff(names(focal), c("x", "y"))])
fit <- DenMark(focal[, c("x", "y", genes[1])], gene = genes[1],
               boundary = tissue, grid_size = 15, basis = 25, boundfactor = 1.5)
dcg <- denmark_scan(focal[, c("x", "y", genes)], boundary = tissue,
                    grid_size = 15, basis = 25, save_dir = "denmark_scan")

denmark_annotate(obj, dcg) writes the estimates, their intervals and the calls back into the Seurat object’s gene metadata. denmark_check_hsgp() compares the posterior length-scales with a rule of thumb for the number of basis functions and the boundary factor, and denmark_save() keeps a fit for later sessions.