Package {scNPM}

Single-Cell Network Propensity Matching

Description

Estimate and compare heterogeneous gene conditional-dependence networks across two conditions. A latent Gaussian mixture and probit observation mechanism represent cellular heterogeneity and zero observations. Protein interaction information supplies edge-specific prior weights, and network propensity matching guides stochastic updates of uncertain edges.

Details

Matching uses centered quadratic forms against all current treatment and control reference networks (2 * K features), with pooled predictor scaling. Default staged inference freezes a stable complete matching, disables feedback and discards a post-switch burn-in before retaining baseline-model draws. Group networks use per-subtype BFDR selection from marginal slab probabilities. Directional edge selection controls an estimated posterior expected FDP, with a conservative all-counterparts option; it is not a uniform frequentist FDR guarantee. Coupled mode retains version 1.0.3 behavior.

The augmented criterion trace records likelihood, network-prior and SMD contributions. It is diagnostic and need not increase. HMC uses fixed step sizes, leapfrog counts and Metropolis acceptance correction.

Author(s)

Maintainer: Tao Liu 25210690190@m.fudan.edu.cn

Authors:

• Tao Liu 25210690190@m.fudan.edu.cn

See Also

run_scNPM, compute_qf_smd

Cell-Level Metadata (AMD Example Data)

Description

A data frame with per-cell metadata (quality-control metrics computed with scater, sample provenance and the PR_or_MR region indicator) for the 400 cells of the AMD example subset (200 PR, 200 MR). Rows are aligned with the columns of goal_gene_mt. The first example covariate column in run_scNPM is typically built from PR_or_MR (1 = MR, treatment; 0 = PR, control).

Usage

data(top50_goal_gene_mt)

Format

A data frame with 400 rows and 52 columns: cell barcodes (Barcode), cell type labels (Labels), the region indicator PR_or_MR, donor identifier (individual), tissue, spatial coordinates (x, y), and scater quality-control metrics (columns prefixed scater_qc.: library size, detected features, percentage of counts in top features, mitochondrial percentage, etc.).

Source

Derived from a public single-cell RNA-seq study of human retinas (see the scNPM reference publication).

Column Indices of MR Cells (AMD Example Data)

Description

Integer vector of column indices (in goal_gene_mt/mt) of the cells belonging to the MR group, used as the treatment arm in the example analysis.

Usage

data(top50_goal_gene_mt)

Format

An integer vector of length 200.

Source

Derived from the PR_or_MR metadata of the example dataset.

Column Indices of PR Cells (AMD Example Data)

Description

Integer vector of column indices (in goal_gene_mt/mt) of the cells belonging to the PR group, used as the control arm in the example analysis.

Usage

data(top50_goal_gene_mt)

Format

An integer vector of length 200.

Source

Derived from the PR_or_MR metadata of the example dataset.

Compute the SMD of the Network Quadratic Form for Matched Cluster Pairs

Description

In all-reference mode (version 1.0.3), this helper returns one row per matched pair and reference network, using the exact final network features. The additional column 'network_feature' identifies each reference. The legacy vector interface below remains supported for older results.

Usage

compute_qf_smd(result_list)

Arguments

Details

The legacy network quadratic form qf_i = \theta_i' \Omega_{k(i)} \theta_i enters the propensity score model as an additional covariate. This helper reports the matched-pair SMD of the quadratic form using the same formula as 'smd_covariate_results':

SMD = |mean_t - mean_c| / sqrt((var_t + var_c) / 2)

Two modes are supported (returned in the attribute 'mode'):

• **exact**: 'result_list' contains 'qf_global' and 'final_group_label' (scNPM >= 1.0.1). The helper uses the quadratic forms and labels from the internal final matching call.

• **posterior-approx**: for results saved by older versions (with the heavy 'Result' element stripped), the quadratic form is recomputed from stored sampling summaries ('theta_est', 'invcov_treat_est', 'invcov_ctrl_est', 'group_est'). The legacy mode name is retained for compatibility; these summaries need not approximate the final matching state closely.

Value

A data frame with columns 'treat_cluster', 'matched_control_cluster' and 'qf_smd', rounded to four decimal places. All-reference results additionally include 'network_feature' and contain K * (2 * K) rows; their mode attribute is 'exact-all-reference'. Unmatched or empty pairs have NA. Singleton samples have variance zero. At the numerical zero-variance boundary, equal means give zero and unequal means give infinity, using the matching routine's tolerances.

Examples

# Small synthetic example (exact mode), runs in well under 5 seconds
set.seed(1)
G <- 5; K <- 2; nT <- 30; nC <- 30
invT <- array(diag(G), c(G, G, K)); invC <- array(diag(G), c(G, G, K))
thT <- matrix(rnorm(nT * G), nT, G)
thC <- matrix(rnorm(nC * G, mean = 0.5), nC, G)
qfT <- rowSums((thT %*% invT[, , 1]) * thT)
qfC <- rowSums((thC %*% invC[, , 2]) * thC)
rl <- list(
  K = K,
  matching_result_df = data.frame(treatment_cluster = 1:K,
                                  matched_control_cluster = c(2, 1)),
  qf_global = c(qfT, qfC),
  final_group_label = c(rep(0, nT), rep(K + 1, nC))
)
compute_qf_smd(rl)

Target Gene Names (AMD Example Data)

Description

Character vector with the symbols of the 50 target genes (rows of goal_gene_mt).

Usage

data(top50_goal_gene_mt)

Format

A character vector of length 50 (e.g., "FTH1", "SAG", "TF", "CLU", "PDE6G", "GLUL", ...).

Source

Selected in the scNPM reference publication.

Top-50 Gene Expression Matrix (AMD Example Data)

Description

Sparse UMI count matrix (genes x cells) for the 50 target genes used in the AMD (age-related macular degeneration) example analysis. This is a stratified random subset (200 PR control and 200 MR treatment cells) of the full 20091-cell dataset, drawn with a fixed seed. Loaded via data(top50_goal_gene_mt) together with the companion objects cell, goal_gene, mt, cols_to_PR and cols_to_MR.

Usage

data(top50_goal_gene_mt)

Format

A dgTMatrix with 50 rows (genes) and 400 columns (cells); entries are raw UMI counts.

Source

Derived from a public single-cell RNA-seq study of human retinas (see the scNPM reference publication).

Full Gene Expression Matrix (AMD Example Data)

Description

Sparse UMI count matrix (all detected genes x cells) from which the top-50 target gene matrix goal_gene_mt was extracted. Provided so users can derive alternative gene subsets. This object covers the same 400 cells as goal_gene_mt; genes detected in fewer than 3 of these cells were removed.

Usage

data(top50_goal_gene_mt)

Format

A dgTMatrix with 12530 rows (genes) and 400 columns (cells); entries are raw UMI counts.

Source

Derived from a public single-cell RNA-seq study of human retinas (see the scNPM reference publication).

Plot Differential Gene Regulatory Network (Combined/Separate Mode)

Description

Plots differential networks (treatment - control) with flexible output: - Combined mode: solid/dashed edges in a single plot (default) - Separate mode: solid/dashed edges as two independent plots

Usage

plot_diff_network(
  diff_mat,
  gene_names,
  solid_threshold = 0.5,
  dashed_threshold = -0.5,
  solid_color = "#ADD8E6",
  dashed_color = "#ADD8E6",
  output_mode = c("combined", "separate"),
  ...
)

Arguments

Value

- If output_mode = "combined": ggplot object (single plot with solid/dashed edges) or NULL (no edges) - If output_mode = "separate": list of ggplot objects (solid/dashed) or NULL for empty edges

Examples

# Self-contained example with a small random differential matrix (< 5 sec)
set.seed(1)
gn <- paste0("Gene_", 1:10)
dm <- matrix(rnorm(100), 10, 10, dimnames = list(gn, gn))
dm <- (dm + t(dm)) / 2
diag(dm) <- 0

# Mode 1 - Combined (solid+dashed in one plot, default colors: light blue)
p_combined <- plot_diff_network(
  diff_mat = dm,
  gene_names = gn,
  output_mode = "combined"
)

# Mode 2 - Separate (solid/dashed as two plots; custom colors optional)
p_separate <- plot_diff_network(
  diff_mat = dm,
  gene_names = gn,
  output_mode = "separate",
  solid_color = "#2ECC71",
  dashed_color = "#F39C12",
  label_size = 3
)

Plot a Selected Group Network

Description

Uses the saved BFDR mask from [select_group_edges()] or a staged [run_scNPM()] fit. No probability rethresholding occurs in BFDR mode.

Usage

plot_group_network(
  edge_mat,
  gene_names = NULL,
  group_type = c("treatment", "control"),
  edge_color = "#ADD8E6",
  ...,
  selection = c("bfdr", "frequency"),
  frequency_threshold = 0.5
)

Arguments

Value

A ggplot object or NULL for an empty selected network.

Examples

# Synthetic retained posterior summaries for one matched pair.
# These values illustrate selection; no MCMC is run in this example.
gn <- paste0("G", 1:3)
p <- matrix(c(0,.99,.97,.99,0,.60,.97,.60,0), 3,
            dimnames = list(gn, gn))
z <- array(0, c(3,3,1))
p10 <- array(p, c(3,3,1)); p00 <- array(1-p, c(3,3,1))
diag(p00[, , 1]) <- 0
example_fit <- list(
  K = 1L, G = 3L, gene_names = gn, inference_mode = "staged",
  matching_result_df = data.frame(treatment_cluster=1L, matched_control_cluster=1L),
  sampling_schedule = list(feedback_disabled_for_retained=TRUE, retained_draws=500L),
  joint_edge_probabilities = list(p00=p00, p10=p10, p01=z, p11=z)
)
group <- select_group_edges(example_fit)
plot_group_network(group$networks$treatment_1)

Core Network Plotting Function

Description

Generates publication-quality gene regulatory network plots from adjacency matrices.

Usage

plot_network_core(
  adj_matrix,
  mode = "undirected",
  layout_fun = igraph::layout_nicely,
  node_id_col = "gene",
  size_var = "degree",
  color_var = "degree",
  size_range = c(1, 5),
  color_palette = c("steelblue", "white", "orangered"),
  label_size = 2,
  edge_color = "#ADD8E6",
  edge_alpha = 0.6,
  edge_size = 0.3,
  curvature = 0,
  edge_linetype = "solid",
  legend_position = "right"
)

Arguments

Value

ggplot object (network plot) or NULL if no edges exist.

Run Single-Cell Network Propensity Matching (scNPM)

Description

Estimate condition-specific gene networks using staged stochastic inference. The default first adapts networks and matching, aligns component labels to fixed initial mean anchors, and freezes a complete mapping after sustained stability. It then disables counterpart feedback, completes post-switch burn-in and retains final baseline-model samples. Failure to stabilize by the adaptation limit stops the run without producing final discoveries.

Usage

run_scNPM(
  goal_gene_mt,
  cell_t,
  K,
  ds_seed = 20250102,
  num_iter = 1000,
  num_save = 500,
  protein_input = NULL,
  ssp_v0 = 0.02,
  ssp_v1 = 1,
  ssp_l = 1,
  edge_prob_threshold = 0.6,
  edge_prob_adj = 0.75,
  epsilon_theta = 0.2,
  num_step_theta = 20,
  epsilon_lam = 0.01,
  num_step_lam = 10,
  ppi_prior_low = 0.05,
  ppi_prior_high = 0.5,
  network_covariates = c("all_reference", "legacy"),
  propensity_ridge = 1e-06,
  inference_mode = c("staged", "coupled"),
  adapt_min = NULL,
  adapt_max = NULL,
  freeze_patience = NULL,
  final_burnin = NULL,
  diagnostic_every = 10L,
  differential_method = c("auto", "bayes_fdr", "frequency"),
  fdr_level = 0.05,
  fdr_scope = c("frozen", "all_matchings"),
  fdr_min_draws = 200L,
  matching_method = c("lexicographic", "greedy"),
  fdr_rule = c("delta_threshold", "posterior_error"),
  group_fdr_level = 0.05,
  verbose = FALSE
)

Arguments

Details

Network features are centered distances to all 2*K current references, with pooled scaling; the legacy uncentered feature remains available. Default differential selection bounds the estimated posterior directional error over every possible control counterpart. This is model-based Bayesian error selection, not a finite-sample frequentist FDR guarantee. The coupled mode reproduces version 1.0.3 and uses descriptive frequency contrasts. The J_aug trace uses an augmented-data likelihood, a state-wise network prior kernel and total SMD. It is a diagnostic, not an optimized objective or the observed-data marginal likelihood; stochastic decreases are allowed.

Value

A named list containing all key analysis results (can be saved with 'save()'):

Examples

# Small demonstration of staged inference; use longer runs for an analysis.
set.seed(7)
x <- matrix(rpois(6 * 40, lambda = 4), nrow = 6,
            dimnames = list(paste0("G", 1:6), NULL))
cells <- cbind(condition = rep(0:1, each = 20))
ppi <- data.frame(node1 = "G1", node2 = "G2", combined_score = 0.9)
fit <- run_scNPM(x, cells, K = 1, num_iter = 60, num_save = 20,
                 adapt_min = 5, adapt_max = 15, freeze_patience = 3,
                 final_burnin = 20, fdr_min_draws = 20,
                 num_step_theta = 2, num_step_lam = 2,
                 epsilon_theta = 0.01, epsilon_lam = 0.01,
                 protein_input = ppi)
fit$matching_result_df
fit$group_fdr$summary
compute_qf_smd(fit)

Select Differential Edges by Posterior Directional Error

Description

Uses the final baseline-model Rao-Blackwell probabilities from a staged run. The selection unit is one treatment cluster and unordered gene pair; all K clusters form a single family. The most supported direction under the frozen map is reported. For 'all_matchings', its error is bounded by the largest posterior error over all K possible control counterparts. This is a conservative bound for every fixed bijection, not a posterior distribution over matchings and not a frequentist FDR guarantee.

Usage

select_differential_edges(
  result,
  fdr_level = 0.05,
  scope = c("frozen", "all_matchings"),
  min_draws = 200L,
  rule = c("delta_threshold", "posterior_error")
)

Arguments

Details

The selected set has mean estimated directional error at most 'fdr_level'. Monte Carlo approximation, label identification and model assumptions still apply. Edge states refer to continuous spike/slab allocations, not exact zero precision entries. Batch MCSE is a diagnostic, not a confidence bound or convergence certificate.

Value

A list with threshold (strict |Delta P| > threshold), threshold_curve, rule, and the full edge table, selected rows, estimated error bound, target, scope, sample count and status. Each table row has an 'error_bound' for that edge and a 'frozen_error' for its reported pairing. 'estimated_error_bound' is the mean bound across selected rows.

Select Group Networks by Posterior Expected False Discovery Proportion

Description

Reconstructs marginal Rao-Blackwell slab probabilities from final-stage joint probabilities: p10 + p11 for treatment and p01 + p11 for control. Each condition and subtype is a separate family of G*(G-1)/2 unordered pairs. Probabilities are sorted decreasingly, complete tie blocks are evaluated, and the largest prefix with mean error 1-p at most fdr_level is retained. Selection uses p >= cutoff; an empty set has cutoff Inf and BFDR zero.

Usage

select_group_edges(result, fdr_level = 0.05, min_draws = 200L)

Arguments

Details

Selection is performed on the full gene panel before any subnetwork is extracted. Returned masks should be shared by plotting and sensitivity calculations. BFDR applies to each original network, not to a subsequent union over subtypes. Edge states are continuous spike/slab allocations. This estimates posterior expected FDP under the baseline model; it is not a uniform frequentist FDR guarantee. Existing staged 1.0.5 fits with saved joint probabilities can be processed without rerunning the sampler.

Value

A list with 'networks', 'summary', 'table', 'selected', 'probabilities', 'masks', 'target', 'status' and draw counts. Each named network (for example 'treatment_1') contains 'probability_matrix', 'mask', 'cutoff', 'estimated_bfdr', 'threshold_curve', 'table' and 'selected'. 'probabilities' and 'masks' each contain treatment/control G x G x K arrays. 'probability_matrix' carries a saved selection attribute for plotting.

Examples

# Synthetic retained posterior summaries for one matched pair.
# These values illustrate selection; no MCMC is run in this example.
gn <- paste0("G", 1:3)
p <- matrix(c(0,.99,.97,.99,0,.60,.97,.60,0), 3,
            dimnames = list(gn, gn))
z <- array(0, c(3,3,1))
p10 <- array(p, c(3,3,1)); p00 <- array(1-p, c(3,3,1))
diag(p00[, , 1]) <- 0
example_fit <- list(
  K = 1L, G = 3L, gene_names = gn, inference_mode = "staged",
  matching_result_df = data.frame(treatment_cluster=1L, matched_control_cluster=1L),
  sampling_schedule = list(feedback_disabled_for_retained=TRUE, retained_draws=500L),
  joint_edge_probabilities = list(p00=p00, p10=p10, p01=z, p11=z)
)
group <- select_group_edges(example_fit)
group$summary