run_scNPM: Run Single-Cell Network Propensity Matching (scNPM)

In scNPM: Single-Cell Network Propensity Matching

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 samples under the selected final network prior. At the adaptation limit, the default freezes a recent consensus map and records that stability was not established.

Usage

run_scNPM(
  goal_gene_mt,
  cell_t,
  K,
  ds_seed = 20250102,
  num_iter = 5000,
  num_save = 1000,
  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,
  matching_score = c("pairwise", "global"),
  freeze_policy = c("consensus", "error"),
  freeze_window = 100L,
  progress_every = 60,
  progress_file = NULL,
  network_prior = c("auto", "paired", "independent"),
  paired_discordance = 0.75,
  shared_edge_correlation = 0.995,
  ppi_reference = c("control", "both", "treatment")
)

Arguments

goal_gene_mt	Finite nonnegative numeric matrix or data frame, with genes in rows and cells in columns. At least two genes are required. Row names must be unique gene identifiers; missing names are generated.
cell_t	Finite numeric matrix or data frame, with cells in the same order as columns of 'goal_gene_mt'. The first column is 0 (control) or 1 (treatment); both groups must be present. Remaining columns are numeric covariates, standardized within each pair's training fold (pooled over both groups in global mode). Encode categorical covariates numerically before calling this function.
K	Positive integer number of clusters in each condition, no larger than the number of distinct expression profiles in either condition.
ds_seed	Integer seed between 0 and 2147483646. Reproducibility requires the same inputs, software environment and OpenMP thread count.
num_iter	Positive integer total number of sampler iterations.
num_save	Positive integer number of final iterations retained, no greater than 'num_iter'.
protein_input	Optional PPI data frame, local TSV path, or TSV text. Required columns are 'node1', 'node2' and 'combined_score'. A leading '#' on the header is accepted. Positive scores indicate an interaction; duplicate pairs are combined, self-interactions and genes outside the expression matrix are ignored. 'NULL' assigns the lower prior to all edges.
ssp_v0	Positive spike standard deviation, smaller than 'ssp_v1'.
ssp_v1	Positive slab standard deviation.
ssp_l	Positive diagonal shrinkage parameter; the exponential-prior rate for precision diagonals is 'ssp_l / 2'.
edge_prob_threshold	Confidence threshold in (0.5, 1) for deciding whether an edge receives matching feedback; this is separate from the differential-network cutoff.
edge_prob_adj	Feedback multiplier in (0, 1]. A value of 1 leaves baseline edge probabilities unchanged by matching feedback.
epsilon_theta	Positive fixed HMC step size for latent expression.
num_step_theta	Positive integer number of HMC steps for latent expression.
epsilon_lam	Positive fixed HMC step size for dropout parameters.
num_step_lam	Positive integer number of HMC steps for dropout parameters.
ppi_prior_low	Baseline edge-state factor for edges without a PPI interaction (default 0.05).
ppi_prior_high	Baseline inclusion probability for PPI-supported edges (default 0.5). In the paired model this defines an odds multiplier, not a marginal edge probability. Requires '0 < ppi_prior_low <= ppi_prior_high <= 1 - 1e-6'.
network_covariates	Network features in the propensity model. '"all_reference"' (default) evaluates every cell against all current treatment networks followed by all control networks, using (\theta_i-\mu_r)^T\Omega_r(\theta_i-\mu_r). '"legacy"' uses the single uncentered quadratic form from version 1.0.2. Other cell covariates are retained; matching_method selects the assignment rule.
propensity_ridge	Positive diagonal ridge stabilization for the pooled logistic IRLS fit (default 1e-6, also applied to the intercept). Used only with matching_score = "global". Pairwise regression uses a fixed slope penalty of 0.001 times the total observation weight.
inference_mode	'staged' (default) or 'coupled' (continuous feedback; use global scores and greedy matching for the earlier workflow).
adapt_min	Earliest freezing iteration. NULL uses min(100, floor(adapt_max/2)), at least 2.
adapt_max	Latest adaptation iteration. NULL uses min(3000, num_iter-num_save-final_burnin).
freeze_patience	Consecutive identical complete maps needed to freeze, at least 2. NULL uses min(5, floor(adapt_max/2)), at least 2.
final_burnin	Minimum iterations after disabling feedback and before retention. NULL uses min(1000, floor((num_iter-num_save)/3)), at least 1.
diagnostic_every	Positive interval for final-stage J_aug diagnostics; adaptation diagnostics and matching histories are recorded every iteration.
differential_method	'auto' chooses 'bayes_fdr' for staged fits and descriptive 'frequency' selection for coupled fits. Explicit 'frequency' is available for staged comparisons; 'bayes_fdr' requires staged sampling.
fdr_level	Posterior directional error target in (0,1), default 0.05.
fdr_scope	'frozen' (default) uses the frozen correspondence. 'all_matchings' uses a conservative bound over all control counterparts.
fdr_min_draws	Minimum retained draws for FDR discoveries, at least 20, default 200. Fewer draws return probabilities and no selected discoveries.
matching_method	'lexicographic' (default) minimizes maximum pair SMD, then total SMD, with deterministic ties and no SMD cap. 'greedy' reproduces the previous assignment rule. A frozen map is not reoptimized.
fdr_rule	'delta_threshold' (default) scans all distinct absolute contrast cutoffs and retains the largest complete threshold set whose mean posterior directional error is at most fdr_level. 'posterior_error' retains the version 1.0.4 error-sorted prefix rule.
group_fdr_level	Group-network BFDR target in (0,1), default 0.05. Each condition/subtype is selected separately over its full gene panel.
verbose	Logical; FALSE (default) suppresses informational output. TRUE emits initialization, progress and timing through message conditions, including native sampler progress; suppressMessages() can suppress these. Warnings and errors are unaffected.
matching_score	Pairwise cross-fitted propensity scores (default), or original global scores.
freeze_policy	At the adaptation limit use a consensus map (default), or stop with error.
freeze_window	Number of most recent valid maps used for consensus, default 100.
progress_every	Positive elapsed-time interval in seconds for progress reports (default 60). Reports occur at completed iterations and also at start, freezing, first retained draw and completion. A single long sweep can delay a report. ETA is an approximate extrapolation of elapsed sampling time.
progress_file	Optional CSV path. Append timestamped progress records even when verbose is FALSE; parent directory must exist. Use a separate file for each concurrent fit. NULL (default) writes no progress file.
network_prior	'auto' (default) resolves to 'paired' for staged fits and 'independent' for coupled fits. 'paired' jointly models matched T/C edge states after freezing. 'independent' retains the 1.0.8 network prior.
paired_discordance	Positive factor at most 1 (default 0.75) controlling discordant-state weights. With b=ppi_prior_low and d=paired_discordance*b*(1-b), base weights for 00,10,01,11 are (1-b-d,d,d,b-d). It is not the fraction of differential edges and is separate from adaptation's edge_prob_adj.
shared_edge_correlation	Correlation in the bivariate slab kernel for state 11, in [0,1), default 0.995. Larger values favor closer signed precision coefficients for common edges; 10/01 remain uncorrelated.
ppi_reference	Condition supported by the PPI source: 'control' (default), ‘treatment', or 'both'. Let R be the odds of an edge’s PPI weight divided by the odds of ppi_prior_low. Multiply the base state factor by R for each supported present condition, then normalize. Choose both for condition-neutral PPI evidence; control reproduces the control-derived PPI design. Ignored by the independent network prior.

Details

Network features are centered distances to all 2*K current references. Pairwise matching standardizes predictors within each training fold; global matching uses pooled scaling. Default differential selection bounds the estimated posterior directional error for the frozen counterpart. This is model-based Bayesian error selection, not a finite-sample frequentist FDR guarantee. The coupled mode uses descriptive frequency contrasts; combine global scores and greedy assignment for the earlier matching behavior. 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()'):

runtime	MCMC runtime in minutes (difftime object)
K	Number of cell clusters used in analysis
G	Number of genes (rows in goal_gene_mt)
N	Number of cells (columns in goal_gene_mt)
lam0_est	Retained-draw means of lambda0 parameter (length = G)
lam1_est	Retained-draw means of lambda1 parameter (length = G)
mu_treat_est	Retained-draw summaries for treatment group (G x K matrix)
mu_ctrl_est	Retained-draw summaries for control group (G x K matrix)
invcov_treat_est	Retained-draw precision summaries for treatment group (G x G x K array)
invcov_ctrl_est	Retained-draw precision summaries for control group (G x G x K array)
edge_treat_est	Retained-draw edge frequencies for treatment group (G x G x K array)
edge_ctrl_est	Retained-draw edge frequencies for control group (G x G x K array)
group_treat_est	Cluster labels for treatment cells (length = number of treatment cells)
group_ctrl_est	Cluster labels for control cells (length = number of control cells)
pi_treat_est	Retained-draw mixing-weight means for treatment group (length = K)
pi_ctrl_est	Retained-draw mixing-weight means for control group (length = K)
theta_treat_est	Retained-draw latent-expression summaries for treatment group
theta_ctrl_est	Retained-draw latent-expression summaries for control group
theta_est	Retained-draw latent-expression summaries (cells x genes)
matching_result_df	Cluster matching table (treatment_cluster vs matched_control_cluster)
matching_method, smd_cost_matrix, frozen_smd_cost_matrix, matching_objective	Assignment rule, all K by K final-state pair costs, and selected maximum/total SMD. Frozen maps are retained even if final costs change.
smd_results	Standardized Mean Difference (SMD) values for each treatment cluster (computed on the propensity score)
smd_covariate_results	Per-covariate SMD matrix for matched cluster pairs (rows = treatment clusters, columns = the QC covariates in cell_t[, -1] plus the selected network features; NA for unmatched pairs). A value below 0.1 is a conventional balance diagnostic
match_map	Named vector mapping treatment clusters to matched control clusters
qf_global	Legacy uncentered per-cell quadratic form retained for compatibility. It enters the propensity model only in legacy mode.
network_covariates	Raw network-feature matrix from the final matching call. Columns are qf_ref_T1 through qf_ref_TK then qf_ref_C1 through qf_ref_CK in all-reference mode, or qf_network in legacy mode.
network_covariate_mode	Selected network feature mode.
propensity_scaling	Pooled centers and scales for global matching; NULL for pairwise matching, which scales each training fold separately.
global_reference_scaling	Pooled feature scaling retained as a diagnostic in both modes.
propensity_scores	Per-cell scores for the selected pair in pairwise mode; pooled scores in global mode. Pairwise scores from different pairs are not a single global propensity score.
final_group_label	Final MCMC cluster label per cell (0-based: treatment cells 0..K-1, control cells K..2K-1; used by compute_qf_smd())
edge_diff_networks	Per-cluster networks selected by differential_method; solid/dashed indicates treatment/control enrichment. The matrix uses Rao-Blackwell contrasts for bayes_fdr and empirical frequency contrasts for frequency.
inference_mode, differential_method	Resolved inference and selection modes.
sampling_schedule	Freeze time, actual post-switch burn-in and retained interval.
stage_history, matching_history	Per-iteration phase, feedback status and aligned mapping history.
criterion_trace	Recorded augmented likelihood, network prior kernel, SMD loss and J_aug.
joint_edge_probabilities	Final all-pair 00/10/01/11 Rao-Blackwell probabilities, empirical 10/01 frequencies and batch means; slice (k-1)*K+c.
group_fdr	Per-condition/subtype BFDR selection from marginal Rao-Blackwell probabilities; NULL in coupled mode. Contains networks, masks, probabilities, summary and edge tables. Pass a named network to plot_group_network().
differential_fdr	Full directional error table, selected rows, estimated bound, scope and status; NULL for frequency mode.
Result	Retained sampler output, including parameter draws and selected trace summaries
protein_adjacency	Symmetric binary PPI adjacency matrix; diagonal entries are zero.
edge_prior	Per-edge PPI weights; the paired prior converts their odds relative to ppi_prior_low into state multipliers. Diagonal entries are unused.
network_prior, paired_prior	Resolved prior and paired-model settings; paired_prior is NULL for independent fits.
provenance	Package version, input fingerprint, native-library fingerprint, execution environment, times and analysis settings.

Paired network prior

The four states are (treatment, control): 00, 10, 01, 11. Conditional Gaussian kernels have diagonal variances (v0^2,v0^2), (v1^2,v0^2), (v0^2,v1^2), and v1^2 in both conditions for state 11; only state 11 has the specified correlation. All factors belong to one globally normalized prior restricted to positive-definite precision matrices. There are no gene modules or learned module parameters. After freezing, paired column updates integrate one coefficient pair, draw its state, and immediately restore both coefficients. Matched components are relabeled together. The observation-model HMC updates and NPM criterion remain unchanged. Rao-Blackwell probabilities use the joint four-state conditional on each retained draw. Unmatched-counterpart probabilities use within-draw marginal products under this same frozen-pair model, not separate refits with new priors.

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)

scNPM documentation built on Oct. 19, 2026, 5:10 p.m.