scDist — Perturbation Distance Between Conditions (per cell type)

Overview

scDist estimates the distance between the two condition means of a cell type in
high-dimensional gene-expression space. For each cell type you get one number
(Dist.), a 95% credible interval and a p-value — i.e. an answer to "which cell
types are most perturbed by the condition?", not "which genes changed?".

How it differs from differential expression (DE):

Internally: per cell type → PCA (irlba::prcomp_irlba, d PCs) on that cell type's
cells → for each PC fit lmer(y ~ condition + (1|random)) → the coefficient of the
condition per PC, shrunk with ashr, gives a posterior over the vector of
coefficients → Dist. = median of the posterior distance (L2 norm of the
condition effect vector across PCs); p-value is a Monte-Carlo test (10^5 reps) of
the summed Wald statistic.

Upstream (not on CRAN or Bioconductor — install from GitHub):

• Repo: <https://github.com/phillipnicol/scDist> (default branch master, package version 1.1.5)
• Install: Rscript -e 'devtools::install_github("phillipnicol/scDist")'
• Cite: Nicol PB, Paulson D, Qian G, Liu XS, Irizarry RA, Sahu AD (2024)
  Robust identification of perturbed cell types in single-cell RNA-seq data,
  Nature Communications 15:7610. <https://doi.org/10.1038/s41467-024-51649-3>

Verified signature (man/scDist.Rd + R/scDist.R on master, 2026-09):

scDist(normalized_counts, meta.data, fixed.effects, random.effects = c(),
       clusters, d = 20, truncate = FALSE, min.count.per.cell = 20, weights = NULL)

Inputs

scDist() takes two objects only — a matrix and a metadata table. There is no
Seurat method inside the package; the user converts the Seurat object themselves.

1. normalized_counts — a matrix (or object coercible with as.matrix),
   genes on rows, cells on columns (the function transposes internally: the
   source comment is "Normalized counts currently in cells x genes" right before a
   t()). rownames (gene names) are required — the function stop()s
   without them. In practice pass Seurat's scale.data: NormalizeData() +
   ScaleData() output, which is already genes × cells, log-normalized and
   centered. Do not pass raw counts and do not pass a plain library-size-
   normalized matrix if a scale.data slot exists.
2. meta.data — a data.frame with one row per cell, in the same cell order
   as the matrix columns. It must contain the three columns referenced below.

The three column arguments

Ordering matters for fixed.effects: the design is
y ~ fixed.effects[1] + ... + (1|random.effects[1]) + ... and the code reads the
condition coefficient as summary(fit)$coefficients[2, ]. Always put the
condition first; extra fixed effects are covariates (design.null drops only
fixed.effects[1]).

Steps

Reusable script: scripts/scdist_perturbation_distance.R.txt — fully command-line
parameterised, no hard-coded paths. Defaults reproduce the original CD4T run
(d = 20, fixed.effects = "response", random.effects = "sample",
clusters = "subCluster1").

1. Inspect the object before running (cheap, avoids silent skips):

srt <- readRDS("CD4T.seu.globalC.rds")
colnames(srt@meta.data)                       # confirm column names
table(srt$subCluster1)                        # cell-type sizes vs d
table(srt$response, srt$sample)               # design: is every subject in both arms?
dim(GetAssayData(srt, assay = "RNA", layer = "scale.data"))   # Seurat v5

2. Build the two inputs (exactly as in the original user script):

libme <- list(
  Y         = as.matrix(srt@assays[["RNA"]]@scale.data),  # genes x cells (Seurat v4)
  meta.data = as.data.frame(srt@meta.data)
)
dim(libme$Y)          # sanity: ncol(Y) == nrow(meta.data)

3. Run scDist for one comparison, or once per condition subset:

out <- scDist(normalized_counts = libme$Y,
              meta.data        = libme$meta.data,
              d                = 20,                    # PCs used
              fixed.effects    = "response",            # condition of interest
              random.effects   = c("sample"),           # absorb subject effect
              clusters         = "subCluster1")         # cell-type column

4. Take the results out — out$results is a data frame whose rownames are
   the cell-type labels:

out$results
#                  Dist. 95% CI (low) 95% CI (upper)      p.val
# CD4_Tn_CCR7          0            0      0.5700683 0.24810752
# CD4_Tcm_IL7R         0            0      0.0000000 0.08592914
# ...
out$design        # y ~ response + (1 | sample)   -> confirms what was fitted
out$gene.names    # genes actually used (= nrow of the input matrix)
out$vals[["CD4_Tcm_IL7R"]]$beta      # condition coefficient per PC
out$vals[["CD4_Tcm_IL7R"]]$D.hat     # unshrunk estimate (vs out$results$Dist. = posterior median)

5. Save and plot:

saveRDS(out, file = "scDist.cd4t.Response.PRE.rds")

scDist::DistPlot(out)        # distance + 95% CI per cell type (built in)
scDist::FDRDistPlot(out)     # distance vs -log10(FDR) (built in)
scDist::distGenes(out, cluster = "CD4_Tcm_IL7R")   # genes driving that cell type's shift

6. Typical PRE/POST/ALL workflow (the original analysis): run the script three
   times (all cells, time == "pre", time == "post"), then merge the result
   tables — merge() appends suffixes .pre / .post, turning Dist. into
   Dist..pre and p.val into p.val.pre. A scatter of -log10(p.val.pre) vs
   -log10(p.val.post), with ggrepel::geom_text_repel() labels for cell types
   significant in either condition, is the published figure:

out.list <- list(PRE  = readRDS("scDist.cd4t.Response.PRE.rds"),
                 POST = readRDS("scDist.cd4t.Response.POST.rds"),
                 ALL  = readRDS("scDist.cd4t.Response.ALL.rds"))
out <- merge(out.list$PRE$results  |> tibble::rownames_to_column("sub"),
             out.list$POST$results |> tibble::rownames_to_column("sub"),
             by = "sub", suffixes = c(".pre", ".post"))
out <- merge(out, out.list$ALL$results |> tibble::rownames_to_column("sub"), by = "sub")

ggplot2::ggplot(out) +
  ggplot2::geom_point(ggplot2::aes(x = -log10(p.val.pre), y = -log10(p.val.post),
                                   size = -log10(p.val), col = sub), alpha = 0.5) +
  ggrepel::geom_text_repel(
    data = subset(out, p.val.pre < 0.05 | p.val.post < 0.05),
    ggplot2::aes(x = -log10(p.val.pre), y = -log10(p.val.post), label = sub, col = sub),
    size = 4, max.overlaps = 20, show.legend = FALSE, min.segment.length = 1) +
  ggplot2::theme_classic() + ggplot2::labs(x = "PRE", y = "POST") +
  ggplot2::ggtitle("Response")

7. Command-line run of the bundled script (does steps 2–5 in one go):

# all cells; original parameters are the defaults
Rscript scripts/scdist_perturbation_distance.R.txt \
    --input-rds CD4T.seu.globalC.rds --output-dir out --prefix cd4t_response_ALL

# one condition at a time (run twice, then merge the two results CSVs)
Rscript scripts/scdist_perturbation_distance.R.txt --input-rds CD4T.seu.globalC.rds \
    --subset-col time --subset-value pre  --output-dir out --prefix cd4t_PRE
Rscript scripts/scdist_perturbation_distance.R.txt --input-rds CD4T.seu.globalC.rds \
    --subset-col time --subset-value post --output-dir out --prefix cd4t_POST

Parameters

Rscript scripts/scdist_perturbation_distance.R.txt --help prints the full usage.

Outputs

In-memory object returned by scDist() (verified against a real run's .rds):

out$results     data.frame, rownames = cell-type labels, columns:
                  "Dist."          posterior median distance (NOT D.hat)
                  "95% CI (low)"   posterior 2.5% quantile
                  "95% CI (upper)" posterior 97.5% quantile
                  "p.val"          Monte-Carlo p-value (resolution ~1e-5)
                cell types that were skipped are ABSENT from this table
out$vals[[ct]]  per cell type: D.hat, D.se, W, beta (per PC), beta_sd, lambda,
                scores, data, dfs, p.sum, p.max, D.post.med, D.post.lb,
                D.post.ub, loadings, beta.hat
out$design      the fitted formula, e.g. y ~ response + (1 | sample)
out$gene.names  genes used (== nrow of the input matrix)

Files written by the script — --output-dir is the output root; every file
keeps its name and lands in one of three fixed sub-directories, and the script
prints a content summary of the three directories at the end:

• <out>/Rdata/<prefix>.rds — the full scDist object (feeds DistPlot / distGenes / plotBetas later)
• <out>/Stat/<prefix>.results.csv — cell_type, Dist., 95% CI (low), 95% CI (upper), p.val, FDR
  (FDR added by the script with p.adjust(method="fdr"); p.val is the raw value)
• <out>/Stat/<prefix>.vals_summary.csv — per cell type n_cells, D.hat, D.se, W, p.sum, p.max, D.post.med/lb/ub
• <out>/Img/<prefix>.DistPlot.{pdf,png}, <out>/Img/<prefix>.FDRDistPlot.{pdf,png} — --make-plots TRUE
• <out>/Stat/<prefix>.log.txt — parameters, design echo, warnings, saved paths

Layout in short: Img/ = figures, Stat/ = tables + log, Rdata/ = saved objects.

Reading the numbers: larger Dist. = the cell type's transcriptome shifted
further between the two conditions. Distances are in the units of the PCA score
space of that cell type, so compare cell types within one run, not across runs
(different gene sets / cell counts give different scales). Dist. = 0 with a
non-significant p-value = "no evidence of a shift", not "cells are identical".

Pitfalls

• Use scale.data, not counts. scDist runs PCA and then lmer on the values.
  Raw counts make the distance dominated by library size / mean expression and
  break the Gaussian mixed-model assumptions. The original script's
  srt@assays[["RNA"]]@scale.data is exactly right; on Seurat v5 use
  GetAssayData(srt, assay = "RNA", layer = "scale.data") (passing slot= is
  defunct in SeuratObject 5). scale.data only exists after ScaleData(); if it
  was run on variable features only, out$gene.names is just those genes (1500 in
  the original CD4T object). Multi-layer (v5) objects may need JoinLayers()
  first, otherwise the layer is empty — the script detects an empty slot and stops
  with that hint.
• d vs cell-type size. d (default 20) is the number of PCs. scDist raises
  min.count.per.cell to at least d and then nexts over any cell type with
  fewer cells, silently dropping it from out$results. With small
  sub-clusters you must lower d (e.g. 5–10) or the cell type simply will not
  appear. Always diff names(out$vals)/rownames(out$results) against
  unique(meta[[clusters]]) to see which cell types were skipped (the script
  warns up front and reports the design). Too-large d also over-fits well before
  it errors (inferred from the source, not documented upstream).
• random.effects must be able to absorb the individual differences — that is
  the whole point of the method. Pass the subject/patient/sample column, not the
  condition. If the same subjects are not present in both conditions (or, worse,
  every subject is single-condition), condition and subject are confounded and the
  distance measures subject identity. Check table(condition, subject) before
  running; the script warns when a subject appears in only one condition.
• Never put clusters into random.effects and never put the condition into
  random.effects: the condition must stay a fixed effect (fixed.effects[1]),
  otherwise the whole estimated distance is absorbed by the random structure.
• Run one cell type per call — you already do. clusters is a column, not a
  value: scDist loops internally. Do not pre-subset the object to one cell
  type (--subset-* is for conditions/timepoints, not for cell types).
• Set a seed. p-values come from a 10^5-replicate Monte-Carlo loop (and ashr
  posterior sampling), so re-running an unchanged analysis can move the smallest
  p-values slightly; --seed (default 1) makes runs reproducible.
• A cell type is skipped, not analysed, when any fixed/random effect has a single
  level inside it (check_effects()), e.g. a sub-cluster present in only one
  of the two conditions. The script warns about this too.
• Per cell type, one lmer per PC: with many cell types and many cells the run
  is the expensive part (each fit is a mixed model). Pre-filtering genes to HVGs
  (as the original object was: 1500 genes) keeps it fast.
• Direction is not reported. scDist gives a magnitude; use
  distGenes(out, cluster=...) / featureImportance / plotBetas to see
  which genes moved and in which direction. For "up vs down" claims, follow up
  with a pseudobulk DE.
• DistPlot/FDRDistPlot need ggplot2 ≥ 3.x and ggpubr; distGenes needs
  ggrepel. Plotting failures must not lose the results — the script writes the
  CSV/RDS first and treats plot errors as warnings.
• fixed.effects order is load-bearing: only fixed.effects[1] is tested
  (design.null <- makeDesign(fixed.effects[-1], ...)); anything after it is a
  covariate. Passing "age,response" computes the age effect while looking like a
  response analysis.

Environment

• R ≥ 4.0 (package requirement; verified on R 4.5.3 in this container).
• Packages: scDist (GitHub only — not on CRAN/Bioconductor), Seurat /
  SeuratObject, plus scDist's own imports: lmerTest, lme4, irlba, ashr,
  data.table, ggplot2, ggpubr, ggrepel (≥0.9.6), ggraph.
  sscVis is only used in the original script for a colour palette
  (sscVis:::auto.colSet(16, name="Paird")) — cosmetic, not required for the
  analysis.
• Install: Rscript -e 'devtools::install_github("phillipnicol/scDist")'
  (or remotes::install_github). No CRAN/Bioconductor equivalent exists.
• The bundled script uses base-R argument parsing only — it has no
  optparse/argparse dependency, so --help works even before scDist is
  installed.
• Status in this container (checked 2026-09, R 4.5.3): installed — Seurat 5.5.1,
  lme4 2.0.6, irlba 2.3.7, data.table 1.18.4, ggplot2 4.0.3, ggpubr 0.6.3,
  ggrepel 0.9.8, ggraph 2.2.2, dplyr 1.1.4, tibble 3.3.1;
  not installed — scDist, lmerTest, ashr (install scDist and let its
  dependencies resolve lmerTest/ashr).

> Packaging note. Scripts and tables in this skill carry a .txt suffix on their
> original extension (e.g. .R.txt, .tsv.txt) because the SkillHub upload allowlist
> does not accept .R or .tsv. This is cosmetic: Rscript is extension-agnostic, so
> Rscript scripts/<name>.R.txt ... behaves exactly like .R, and tab-separated data
> parses identically from .tsv.txt. Keep the suffix when repackaging.