Skip to contents

Fit DenMark gene by gene and call density-correlated genes (DCG)

Source: R/scan.R
denmark_scan.Rd

A gene is a positively density-correlated gene (PCG) when the cred_level credible interval of rho lies above 0, a negatively density-correlated gene (NCG) when it lies below 0, and "not DCG" otherwise.

Usage

denmark_scan(
  data,
  genes = NULL,
  save_dir = NULL,
  cred_level = 0.95,
  quiet = TRUE,
  verbose = TRUE,
  workers = 1,
  ...
)

Arguments

data

Data frame with columns x, y and one column per gene (cells in rows).

genes

Genes to scan (default: all non-coordinate columns).

save_dir

Optional folder; one .rds file per gene is written, and genes already present are loaded instead of refitted (resumable scans).

cred_level

Credible level used to call DCGs.

quiet

Suppress the sampler output, passed to DenMark().

verbose

Report one message per gene with rho and its credible interval.

workers

Number of genes fitted at the same time, each in its own R session (needs the future.apply package and an installed DenMark). Every gene keeps its seed, so the results equal those of the sequential scan. Each fit also runs its chains in parallel (parallel_chains, by default chains), so a scan uses up to workers * chains cores. For a cluster, system.file("slurm", package = "DenMark") holds a template that fits one gene per SLURM array task.

...

Further arguments passed to DenMark() (grid_size, basis, ...).

Value

Data frame with one row per gene, sorted by decreasing posterior mean of rho: gene, rho_mean, rho_median, rho_lower, rho_upper, prob_positive, total_counts, max_rhat, min_ess_bulk, divergences, time_min, error (NA unless the fit failed), class (PCG, NCG or not DCG) and converged (maximum R-hat below 1.01 and no divergent transitions). The credible level is stored in attr(, "cred_level").

Examples

toy <- denmark_toy_data()
tissue <- denmark_alpha_shape(toy$x, toy$y)
tumor <- subset(toy, celltype == "Tumor", select = -celltype)
scan <- denmark_scan(tumor, genes = c("geneA", "geneC"), boundary = tissue,
                     grid_size = 10,
                     basis = 10, chains = 2, iter_warmup = 300,
                     iter_sampling = 300, seed = 2024, verbose = FALSE)
scan[, c("gene", "rho_mean", "rho_lower", "rho_upper", "class")]
#>    gene   rho_mean  rho_lower  rho_upper class
#> 1 geneA  0.6332270  0.3432107  0.8651198   PCG
#> 2 geneC -0.5917726 -0.8330993 -0.3213833   NCG