Skip to contents

Fit DenMark to one gene

Source: R/fit.R, R/postprocess.R
DenMark.Rd

Grids the data, builds the HSGP inputs, compiles the Stan model (once per session), sets valid initial values and samples with cmdstanr.

Usage

DenMark(
  dataset,
  gene = NULL,
  grid_size = 10,
  basis = 25,
  boundfactor = 1.5,
  bbox = NULL,
  chains = 4,
  parallel_chains = chains,
  iter_warmup = 1000,
  iter_sampling = 1000,
  adapt_delta = 0.95,
  max_treedepth = 10,
  seed = 2024,
  refresh = 400,
  init = NULL,
  rho0 = 0.3,
  stan_file = NULL,
  output_dir = NULL,
  keep_fit = TRUE,
  quiet = FALSE,
  L1 = NULL,
  L2 = NULL,
  boundary = NULL,
  min_area_frac = 0,
  ...
)

# S3 method for class 'denmark'
print(x, digits = 3, ...)

Arguments

dataset

Data frame with columns x, y and the expression column(s), one row per cell; expression values are raw counts.

gene

Name of the expression column (default: the first column other than x and y).

grid_size

Grid cells per axis, or c(nx, ny).

basis

HSGP basis functions per axis, or c(m1, m2).

boundfactor

HSGP boundary factor (>= 1).

bbox

Optional region of interest c(xmin, xmax, ymin, ymax).

chains, parallel_chains, iter_warmup, iter_sampling, adapt_delta, max_treedepth, seed, refresh

Passed to the cmdstanr $sample() method.

init

NULL (valid defaults from denmark_init()) or anything accepted by cmdstanr.

rho0

Starting value of rho.

stan_file

Optional path to a custom Stan file with the parameters of DenMark.stan; with a boundary it must declare vector[n] log_grid_area.

output_dir

Optional folder for the CmdStan CSV files (created if needed).

keep_fit

Keep the CmdStanMCMC object (set FALSE in large scans).

quiet

Suppress the sampler progress output.

L1, L2

Optional half-ranges of the region of interest; only checked against the data (see denmark_prepare()).

boundary

Optional tissue polygon from denmark_alpha_shape(); grid cells outside the tissue are dropped and boundary grid cells use their clipped area (variable-area model).

min_area_frac

See denmark_prepare().

...

For DenMark(), further arguments passed to $sample(); ignored by print().

x

A denmark object.

digits

Number of significant digits in the printed summary table.

Value

An object of class denmark: a list with the cmdstanr fit (fit), gene, stan_data, grid, info, the sampler settings, the posterior summary (see denmark_summary()), the sampler diagnostics, the run time in minutes, the HSGP check hsgp (see denmark_check_hsgp()) and, for backward compatibility with DenMarkmain(), rhat, Ymat and Mmat. With keep_fit = FALSE the fit is dropped and the draws of the key parameters are kept in draws_cache.

Details

The DenMark model (DenMark.stan) has two dependent latent fields, log lambda1(g) = log|g| + beta0 + a11 w1(g) and log lambda2(g) = beta1 + log(Y(g) + 0.05) + a21 w1(g) + a22 w2(g); rho = Astar[2,1] is the density-expression correlation.

A tissue section is rarely a rectangle. With boundary (an alpha shape of all cells of the section, from denmark_alpha_shape()), the grid cells outside the tissue are dropped and every remaining grid cell enters with its area inside the tissue, through the variable-area model DenMark-vararea.stan (the same model with one log_grid_area per grid cell). Without boundary, the whole bounding box is gridded with a common grid area.

The print method reports rho, the sampler diagnostics and the posterior summary.

See also

denmark_scan() to fit many genes.

Examples

toy <- denmark_toy_data()
tissue <- denmark_alpha_shape(toy$x, toy$y)
tumor <- subset(toy, celltype == "Tumor", select = -celltype)
fit <- DenMark(tumor, gene = "geneA", boundary = tissue,
               grid_size = 10, basis = 10, chains = 2,
               iter_warmup = 300, iter_sampling = 300,
               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.343, 0.865]
#>   max R-hat 1.011 | min ESS bulk 241 | divergences 0 | time 0.51 min
#> 
#>    variable     mean    sd   q2.5    q50 q97.5 rhat ess_bulk
#>       beta0  5.60133 0.310  4.933 5.6110 6.177 1.00      583
#>       beta1 -0.00687 0.558 -1.349 0.0732 0.839 1.00      508
#>      ell[1]  0.60501 0.150  0.378 0.5820 0.952 1.01      321
#>      ell[2]  0.77199 0.264  0.408 0.7264 1.469 1.01      352
#>  Astar[2,1]  0.63323 0.136  0.343 0.6526 0.865 1.01      292
#>    sigma[1]  0.99615 0.206  0.693 0.9505 1.528 1.00      292
#>    sigma[2]  1.33412 0.315  0.857 1.2882 2.090 1.01      241