[parashardhapola]

Summary

Pathway scores turn each sample's expression into one number per gene set, so that
pathways can be compared across samples, drawn as a heatmap and tested with the same
designs and contrasts as genes. This adds the four methods of the R package GSVA 2.0.0
(Haenzelmann, Castelo and Guinney 2013) with its options and defaults: GSVA itself,
ssGSEA (Barbie et al. 2009), combined z-scores (Lee et al. 2008) and PLAGE (Tomfohr, Lu
and Kepler 2005), from a samples x genes matrix or a cell type of a cube, with the
heatmap GSVA's guide draws. GSVA is GPL-licensed, so the methods are written from the
papers and GSVA's documentation, and what these leave open (the kernel's table, ties,
constant genes, missing values, ssGSEA's weights and normalization, PLAGE's sign) was
found by running GSVA as a black box. GSVA's method agrees with GSVA to the last bit
on continuous values and on counts of moderate size, and the other methods to rounding
(at most 1e-13 of the largest score).

Changes

• cytepool.pathway_scores.gsva(data, gene_sets, *, gene_ids, sample_ids, group, values, min_cells, expression_filter, kernel, min_samples_without_kernel, exponent, score_type, min_size, max_size, n_threads):
  GSVA's method, with the Gaussian, Poisson or no kernel (kernel="auto" chooses as
  GSVA chooses), tau (exponent) and GSVA's three scores (score_type:
  "difference", maxDiff = TRUE; "kuiper", absRanking = TRUE;
  "max_deviation", maxDiff = FALSE).
• cytepool.pathway_scores.ssgsea(..., exponent=0.25, normalize=True, ...),
  combined_z_scores(...) and plage(...): the other three methods, with the same
  data, set and size arguments.
• cytepool.pathway_scores.PathwayScores: scores (samples x sets, read-only, so that
  fit_linear_model(scores.scores, Design(scores.samples, ...)) takes them as they
  are), samples (the sample table; a cube's sample columns included), gene_sets (a
  GeneSets of the sets scored with the genes each used, keeping the collection's
  description and source), method, settings (with the kernel "auto" chose);
  n_genes and table() (one row per set and sample).
• Data: an array with gene_ids (and optional sample_ids), or a Cube with group,
  whose samples with at least min_cells cells and some counts are scored over the
  genes filter_by_expression keeps, as TMM-normalized log-CPM (values="log_cpm",
  edgeR's cpm(log = TRUE)) or as counts (values="counts", the Poisson kernel for an
  integer layer). Gene ids may come from a features column (gene_ids="symbol").
• cytepool.plotting.plot_pathway_scores(scores, *, sets, color_by, standardize, ax):
  the heatmap of sets x samples on a diverging scale symmetric about zero (blue,
  neutral grey, red; the red arm matched to the blue ramp's lightness at each step),
  samples grouped by a sample variable or labels under a strip of the categorical
  hues, white lines between groups, and standardize as R's heatmap(scale = "row"). The drawing is a private _draw_heatmap(ax, values, rows, columns, label, groups), ready for heatmaps of gene expression.
• Documentation: API page; capability matrix row; a section "Pathway scores per sample"
  in "Gene-set tests" with a choice bullet; rows in "Choosing a method", the overview,
  the glossary and "How to cite"; GSVA among the tools validated against in "Related
  software" (2.6.14 is now Artistic-2.0; 2.0.0 is GPL (>= 2)); four references
  (DOIs checked against Crossref) on the Bibliography page; the Validation page
  (summary row, totals, section, timings, commands); the README's status rows, plan
  and a known limit.
• tests/test_pathway_scores.py, fixtures in tests/data/pathway_scores/ (85 KB) and
  validation/pathway_scores.py; the masked-input inventory gains the four functions
  (for an array and for a cube), the figure and the record.

Design notes (found by probing GSVA 2.0.0)

• The kernel. Each gene's estimate in sample j is (1/n) sum_k K(x_ij, x_ik), summed
  over the samples in their order. The Gaussian kernel has bandwidth sd/4 (probing
  other bandwidths ranked thousands of genes differently). GSVA does not evaluate Phi
  exactly: an exact Phi ranked 1,839 of 30,000 gene-sample pairs differently, and a
  scan of genes whose estimates differ by 1e-5 in one argument showed flat steps of
  0.001. GSVA reads Phi from a table at |t| rounded down to a multiple of 0.001, and
  takes one minus it for negative t; with that table and sequential sums, every rank
  matched on six data sets (normal, few or many samples, rounded values, an outlier
  sample, log counts). The Poisson kernel is P(Y <= x_ij), Y ~ Poisson(x_ik + 1/2);
  without a kernel, the share of samples at or below x_ij.
• The statistic. Genes are ranked by -log((1 - F)/F), whose rounding merges
  estimates a few ulps apart. Grouping tied ranks in the far tail of the kernel showed
  which of 13 candidate forms GSVA computes: only this one reproduced its tie groups, and with it every rank of 47,592 low-count gene-sample pairs matched
  (given R's ppois). Equal statistics take their ranks in reverse gene order (R's
  ties.method = "last"); weights are |p/2 - r|^tau; the walk is GSVA's
  (gsvaEnrichment returns it), and with maxDiff = FALSE equal deviations give the
  negative one.
• kcdf = "auto" chooses the Poisson kernel for R's integer storage, not for
  whole numbers stored as doubles, and no kernel from kcdfNoneMinSampleSize samples
  on. CytePool reads an integer dtype as R's integer storage.
• Genes and sets. GSVA's method, z-scores and PLAGE drop genes whose values are all
  equal before matching (so sizes count only varying genes); ssGSEA keeps them. Sets are
  matched, deduplicated and filtered by minSize/maxSize after matching; empty sets
  are dropped; the output keeps the input order; a set of every gene, or of one gene
  at rank p/2 (weight zero), gives NA in GSVA's method.
• ssGSEA. The score is the closed form of the walk's sum, sum w a / sum w -
  sum_{other} a / (p - k), with a the rank from the bottom (ties in gene order) and the
  weight the average rank truncated to an integer, to the power alpha: a one-sample
  probe with ties determined the truncation (weights 3, 5, 8 for tie blocks at ranks
  2-4, 5-6 and 7-10). normalize = TRUE divides by the range of all scores.
• z-scores and PLAGE standardize with R's scale (sd with n - 1); PLAGE is
  svd(Z[set, ])$v[, 1], the sign R's SVD happens to give.
• Missing values. GSVA's code recognises R's NA, not other NaN payloads (NaN
  passed from Python crashed GSVA; the validation converts it). use = "na.rm" leaves
  missing values out of each gene's kernel and each sample's ranking; ssGSEA matches it
  exactly; z-scores propagate NA to the sets holding the gene; PLAGE stops.
• Exactness. With the same table and summation order, ranks, and so scores, are
  GSVA's to the last bit on continuous values and on counts of moderate size. Genes
  whose estimates are equal in exact arithmetic, as genes with the same low counts in
  different samples are, are ordered by the last bit of R's ppois (correctly rounded
  for only 51% of a grid of counts up to 120, so no other implementation reproduces
  it): on
  counts of 0 to 6 scores differ by up to 2.2e-3 of the largest. The validation
  reports this as information.

Deliberate differences

• PLAGE's sign is chosen, not arbitrary: the loadings (first left singular vector)
  sum to a positive number, so that scores rise with the set's mean standardized
  expression; when they sum to zero, the largest loading, the first of equal ones, is
  positive. GSVA's sign agreed for 42% of the validation's sets. A set of one gene has
  its standardized values scaled to length one, where GSVA gives 1 in every sample.
• GSVA's method with missing values: GSVA's na.rm walk falls by 1/(p - k_j) with
  p counting genes missing in the sample, so it ends above zero (and scores a set of
  every gene); CytePool's falls by 1/(p_j - k_j) and ends at zero. Samples without
  missing values agree exactly. GSVA's default use = "everything" returns values
  such as -1.5e45 with this method, and a sample without values stops it; CytePool
  gives that sample NaN scores.
• ssGSEA with a set of every gene: NaN for that set and normalization over the
  others, where GSVA stops.
• kernel="auto" with an integer dtype and 200 samples or more uses no kernel, as
  documented; GSVA 2.0.0 stops with an error. The Poisson kernel refuses values that
  are not counts (GSVA stops on negative ones and scores fractional ones).
• Values are samples x genes, as everywhere in CytePool; GSVA takes genes x samples.

Validation

validation/pathway_scores.py against GSVA 2.0.0 (R 4.4.3), on simulated
log-expression values (2,000 genes x 24 samples, with constant, equal, whole-number and
outlier genes), negative binomial counts (stored as integers and as decimals), low
Poisson counts, 400 genes x 210 samples, rounded values, 2% missing values and a cube
of one cell type, with 66 sets of 1 to 300 genes (a fifth outside the data, a set of
one gene, of constant genes, listed twice, identical sets, every gene). Differences are
relative to the largest absolute score of a case; NaN must stand where GSVA gives NA.

Check family	Largest difference	Tolerance	Outcome
GSVA's method: Gaussian kernel, tau 0.5 with maxDiff FALSE, tau 2 with absRanking, no kernel, sizes 10 to 100, counts as integers (Poisson) and as decimals (Gaussian), 210 samples (no kernel)	0	1e-10	8 pass
ssGSEA: defaults, alpha 1 unnormalized, alpha 0 on tied values	1.8e-15	1e-10	3 pass
Combined z-scores; PLAGE (sets of 2+ genes, up to sign)	8.7e-14	1e-10	2 pass
Missing values: GSVA's method (samples without any), ssGSEA, z-scores	2.0e-15	1e-10	3 pass
A cube's cell type, as log-CPM and as counts	0	1e-10	2 pass
plot_pathway_scores: 1,584 standardized cells vs R's heatmap(scale = "row")	3.2e-16	1e-10	1 pass
PLAGE: share of sets with GSVA's sign; Poisson kernel on low counts	0.42; 2.2e-3	n/a	2 information
GSVA's method with missing values; PLAGE's one-gene sets; ssGSEA with a set of every gene	4.8e-2; 1.6; GSVA stops	n/a	3 deliberate differences

24 checks, 0 failed (19 pass, 2 information, 3 deliberate differences); 10.7 s, 949
MiB peak (R included), 2 threads. Benchmark (--benchmark, 500 sets of 15 to 500
genes, 50 samples x 15,000 genes; GSVA on one thread, CytePool on 2):

Method	GSVA (s)	CytePool (s)
GSVA, Gaussian kernel	4.13	0.83
GSVA, Poisson kernel	17.61	7.32
ssGSEA	2.21	0.22
Combined z-scores	0.67	0.15
PLAGE	1.95	0.96

RSCRIPT=$REF/r-gsea/bin/Rscript PYTHONPATH=. python validation/pathway_scores.py
RSCRIPT=$REF/r-gsea/bin/Rscript PYTHONPATH=. python validation/pathway_scores.py --write-fixtures
RSCRIPT=$REF/r-gsea/bin/Rscript PYTHONPATH=. python validation/pathway_scores.py --benchmark

Tests

tests/test_pathway_scores.py, 49 tests, 1.2 s warm, about 23 s cold (compiling the
kernels and, for one test, the linear models):

• GSVA's scores from the fixtures (300 genes, 21 sets) in 12 cases, to 1e-10: every
  kernel and score type, ssGSEA with and without normalization, z-scores, PLAGE up to
  sign, missing values.
• Values derived by hand: GSVA's three scores of a two-gene set on three samples
  without a kernel (ties, a zero weight, the max-deviation tie rule); ssGSEA with tied
  values (GSVA's -1.9166... included); combined z-scores of two genes; PLAGE of
  correlated, anti-correlated (the sign fallback) and single genes.
• Invariants: the order of the sets and of the genes; relabelled samples; shifts and
  power-of-two scalings of a gene for the kernels; monotone transforms for ssGSEA and
  without a kernel; affine invariance and oddness of z-scores; PLAGE's unit length,
  zero mean and sign against z-scores; thread counts; score ranges.
• Matching, sizes, constant genes, genes without ids, the table, missing values, a
  cube's group against the same log-CPM given as an array, the scores of a cube in
  Design, fit_linear_model and moderated_statistics, the docstring example, and
  error messages.
• The heatmap read back from its artists: cells, group order, strip hues, boundaries,
  symmetric scale, standardization, NaN, errors.

Checks: full suite 4,373 passed, 18 skipped (159 s); the new tests with
OPENBLAS_CORETYPE=Sandybridge, 4 threads and NUMBA_CPU_NAME=generic pass; ruff
check and format, mypy, make -C docs api, make -C docs html (Node.js 24) and
tests/test_documentation.py pass.