Two-stage niche differential expression

Estimates each core's (sample's) niche slope from a one-step response anchored on a SpaNorm fit, pools those slopes per patient, then contrasts the patient slopes between conditions with a variance-moderated limma fit. Because condition is assigned per patient, patients are the experimental units; the one-stage fitSpiDE() path treats cells as replicates of that contrast, which is anti-conservative when between-patient variation in the niche-expression relationship is real.

Usage

twoStageSpiDE(
  spe,
  condition,
  sigma,
  index = NULL,
  niche = NULL,
  patient = NULL,
  patient.covariates = character(),
  assay = "counts",
  cell_type = "cell_type",
  sample_id = "sample_id",
  name = "Niche",
  min.cells = 30L,
  fdr = 0.05,
  stage1 = c("spanorm", "ols", "nb"),
  epsilon = NULL,
  winsor = 4,
  lambda.a = 0,
  maxit.psi = 2,
  pool.psi = TRUE,
  backend = "cpu",
  verbose = TRUE
)

Arguments

spe

a SpatialExperiment with niche reducedDims (see buildNiches()).

condition

a character, the colData column of the tested condition, constant within patient.

sigma

a numeric, the bandwidth (one value; slopes are per bandwidth).

index, niche

character vectors restricting index / niche cell types. Use them: the full space is usually too large to detect anything. Stage-1 slopes are estimated jointly across niche columns (see stage1), so restricting niche changes the estimated slope for every remaining niche, not merely which rows of the results table are reported.

patient

a character, the colData column identifying the patient a sample (core) belongs to. Defaults to NULL, in which case each sample is its own patient (no pooling). Multiple samples per patient are pooled by precision (stage 2) before the condition contrast is taken; condition must still be constant within patient.

patient.covariates

a character vector of patient-level colData columns to adjust for in stage 2 (e.g. sex, stage, histology). These are exactly the covariates fitSpiDE() rejects under random != "none".

assay, cell_type, sample_id, name

column/assay names.

min.cells

minimum cells for a sample to contribute a slope (default 30). Samples below it are dropped for that index type and the count is reported in diagnostics$inclusion, rather than contributing an unstable slope.

fdr

the target false discovery rate for the reported table.

stage1

one of "spanorm" (the default), "ols" or "nb", the model stage-1 slopes are estimated from.

"spanorm" fits an NB GLM on the raw counts per (sample, index) subset, with design [1, niche columns] and the stored SpaNorm::SpaNorm() fit's library-size and batch linear predictor as a fixed offset (read from metadata(spe)$SpaNorm; a clear error names SpaNorm::SpaNorm() when none is found). The rationale: SpaNorm models biology only to anchor its normalisation – the fitted biology term is a smooth catch-all, not modelled biology – so spiDE models all biology from scratch and assumes only that the library-size (and within-sample batch, e.g. field-of-view) effects are right. Fixing their predictor at coefficient 1 is that assumption made literal: unlike a fitted library-size covariate, an offset cannot absorb depth-correlated biology. (Before spiDE 0.99.16 this option instead used a one-step "addback" working response linearised at the fitted mean; results from that construction are not comparable and should be re-computed.)

"ols" regresses log-CPM directly on the niche columns with unit weights. It needs no stored SpaNorm fit and no dispersion estimate, so it is the fallback when one is not available (e.g. data(toySpiDE)).

"nb" fits a fresh SpaNorm::fitNB per (sample, index) subset carrying [1, log-library-size, niche columns]. Note the library-size term is a fitted covariate here, so it can absorb a depth-correlated effect; prefer "spanorm" when a SpaNorm fit is available. Both NB paths price per-GENE (dispersion estimation and IRLS setup), not per cell.

epsilon

Deprecated and ignored. The former "addback"/"residual" choice applied to the pre-0.99.16 working-response construction, which no longer exists; supplying the argument raises a warning.

winsor, lambda.a, maxit.psi, backend

forwarded to fitNB by the "nb" path.

pool.psi

logical (default TRUE): for stage1 = "spanorm", estimate one per-gene dispersion per SAMPLE (design: cell-type means plus the ls/batch offset) and supply it to every (sample, index) subset fit, instead of re-estimating per subset. Cuts the dominant stage-1 cost (dispersion estimation, measured at ~half of each subset fit on the real cohort) and estimates psi from all of a sample's cells rather than one type's few. The pooling design omits the niche columns, which errs slightly conservative. Requires SpaNorm with fitNB(psi=); on an older SpaNorm the option is silently inert and per-subset estimation is used.

verbose

a logical.

Value

a SpiDEResults with the tidy results table populated (unchanged schema: gene, ct_index, ct_niche, coef, t, p.niche, fdr.niche, DirectionNiche, bandwidth.max). @fits is empty (no per-bandwidth GLM fit exists for this estimator). Diagnostics are attached at r@diagnostics, a list of three tables: r2 (the niche columns' R2 against the SpaNorm biology AND ls bases – a basis column distinguishes them – per sample x index, "spanorm" stage1 only; high biology overlap is expected, high LS overlap warns that the LS field may absorb the niche signal), inclusion (the per-index patient inclusion table, with a warning when min.cells dropout is associated with condition), and tau2 (the DerSimonian-Laird between-patient variance per index x niche used to weight the stage-2 contrast).

Details

On a patient-label permutation of a 55-patient CosMx cohort, where every triplet is null by construction, fitSpiDE() with random = "intercept" returned ~576 false calls per replicate (realized FDR 1.00 against a nominal 0.05). These next two numbers are historical, measured on this function's predecessor (the nbresid/Welch two-stage estimator, since replaced by the SpaNorm-anchored joint estimator documented here), not on the current implementation; re-measurement on the current estimator is pending. That predecessor returned raw type-I 0.036 and zero false calls, ~100x faster.

It does not solve multiplicity. Also measured on that same predecessor estimator, pending re-measurement here: on a spiked plasmode the true triplets were enriched 5x at $\alpha = 0.05$ and 65x at $\alpha = 0.001$, yet nothing survived BH across a full-panel space of ~1.8 million triplets, and the hierarchical cascade did not help (combining a gene's ~137 mostly-null triplets dilutes a single true effect). Restrict index, niche and the gene set to a pre-specified hypothesis: roughly 4,000 tests is the order at which a $p \approx 10^{-5}$ effect survives.

Examples

data(toySpiDE)
spe <- buildNiches(toySpiDE, sigma = 30)
res <- twoStageSpiDE(spe, condition = "condition", sigma = 30,
                     min.cells = 10, stage1 = "ols", verbose = FALSE)
head(results(res))
#>   gene ct_index ct_niche       coef        t      p.niche  fdr.niche
#> 1  G12        A        B -0.6822765 -3.87396 0.0003693287 0.04431945
#>   DirectionNiche bandwidth.max
#> 1           Down            30

Usage

Arguments

Value

Details

See also

Examples