[parashardhapola]

Summary

The guides showed where the public API makes users repeat work: refitting a cell type
to draw its trend figures, filtering tables by hand, rebuilding a design to test
pathway scores, converting a network to gene sets, writing the same loop over the
studies of a collection, naming aligned cell types by ontology ids. This pull request
adds what removes that work, fixes eight problems the guides found, and makes π₁ of
the replication analysis steady on the few genes a study finds. It is the first of
two API pull requests: it adds names and options and changes no existing call, so
that renames and moves, which the second makes, are checked separately.

Every existing call keeps its result. An equivalence check ran master and this branch
in separate processes on the same inputs, call by call: every public entry point with
a valid call, the inputs of every validation script, edge cases, and master's whole
test suite with every public call it makes recorded. It compared every output element
by element: 54,215 calls in 14,804 cases, and master's 4,931 tests in both trees.
Every difference it found is one of the intended differences written down, each with
the rule that checks it, before its results were read, except 43 in calls that draw
permutations without a seed, which differ between any two runs and agree bit for bit
once seeded.

Changes

Results that keep what they analysed:

• DifferentialExpressionResult keeps the cube and design it analysed (new
  properties; the objects given, not copies), and fit(group=None) -> GroupFit
  computes one group's fit again from them with the recorded settings: bit for bit
  what the analysis computed, on the same machine. The last group fitted is kept. The
  constructor takes cube=None, design=None as keywords, both or neither, so that
  results built from their parts as before keep working; a pickled result carries
  them.
• GroupFit(group, samples, features, log2_cpm, weights, voom, moderated), in
  cytepool.differential_expression: the samples fitted, in order, with their design
  variables, n_cells, library_size, norm_factor and any prior or quality
  weights; the features tested; voom's or limma-trend's log2 CPM; voom's precision
  weights or limma-trend's sample weights; the VoomResult (None for limma-trend);
  and the ModeratedFit or dream's moderated tests.
• plot_mean_variance(result, *, group=None, ax=None) and
  plot_residual_variances(fit, *, group=None, ax=None) also take a result and draw
  its group's fit; limma-trend has no voom trend and dream no ModeratedFit, and
  each says so.
• provenance of PathwayScores, Activities and PrincipalComponents computed
  from a cube: the method, the software versions, the cube's fingerprint, digests and
  record, the layer read and the settings; None for values given as an array. The
  three records now pickle, the provenance with them (they did not pickle before).
• table() keys may be left out when a result has one: DifferentialExpressionResult
  and MetaAnalysisResult.table(group=None, contrast=None),
  CompositionResult.table(contrast=None),
  VariancePartitioningResult.table(group=None); with several, the message lists
  them.
• DESeq2Results.table(*, feature_ids=None), QuasiLikelihoodTest.table(*, feature_ids=None) and ModeratedFit.table(coefficient=None, *, feature_ids=None):
  the tests as tables with the columns of limma_voom's (feature_id,
  log2_fold_change, standard_error, ..., adjusted_p_value), ids feature_0,
  feature_1, ... by default.

Tables and cubes:

• Table.filter(mask=None, /, **values): rows where the mask holds and each named
  column equals its value (result.samples.filter(group_id="B cells", used=True)); a
  missing value matches missing ones; a value of the wrong kind, or a list, is refused
  with the mask that does it. Calls that pass mask= by keyword keep working.
  Table.sort(by, *, descending=False): stable, missing values last.
• Cube.with_columns(*, samples=None, groups=None, features=None): columns added to,
  or replaced in, a cube's tables (not the ids), the layers kept as they are (a lazy
  cube reads nothing), with a provenance record of the columns added and replaced.
  Cube.collapse_samples(by=None) also pools by several columns (ids joined by _,
  refused when two combinations join alike) and, with None, every sample into one,
  the values plot_dotplot draws.
• to_anndata(cube, *, x_layer=None): x_layer="counts" copies that layer into X
  too.
• plot_dotplot(..., vmin=None, vmax=None, ax=None): colour limits, so that two dot
  plots share a scale; the docstring says that the colour bar and the size key are
  the figure's other Axes.

Several studies:

• CubeCollection(cubes, *, align_groups_by=None, group_mapping=None, group_names=None, within_study=()): group_names={"CL:0001054": "CD14 monocytes"}
  names aligned groups; the names become their ids and the aligned column keeps the
  ontology ids. The record holds the names only when given, so that existing records
  are unchanged.
• cytepool.study_results (new module): analyse_each_study(collection, formula, *, analysis=limma_voom, random=None, levels=None, contrasts=None, **options) -> StudyResults builds each study's design from its aligned cube's samples and runs
  limma_voom, limma_trend or dream_voom; an error names its study.
  StudyResults(collection, results), a read-only mapping of study to result, has
  collection, studies, designs, summary (every study's summary with a study
  column), table(group=None, contrast=None) (what align_study_results gives, which
  stays) and matrix(column, group=None, contrast=None) (studies x features, float64,
  NaN where a study did not test a gene).

Scores, networks and filters:

• cytepool.score_tests.limma(data, design, *, contrast=None, feature_ids=None, sample_ids=None, weights=None, trend=False, robust=False) -> Table (new module):
  limma's lmFit, contrasts.fit, eBayes and topTable for one contrast of a
  design, on PathwayScores, Activities of samples or an array; the design's rows
  are matched to the scored samples by sample_id, so the design built for the cube
  serves, and a design that those samples do not determine is fitted on its column
  space, as limma_voom fits a group. Columns: the set, source or feature id,
  estimate, standard_error, average, t, p_value, adjusted_p_value,
  log_odds.
• Network.to_gene_sets() -> GeneSets: each source's targets in edge order, without
  weights, with the network's name, version, licence, file and digest.
• expression_filter={"sample_groups": "condition"}: edgeR's filterByExpr(group = ...) with a column of the samples, in limma_voom, limma_trend, dream_voom,
  variance_partitioning, the pathway scores and the activities of a cube, so that
  two models of the same samples can test the same genes; other settings combine with
  it. The design remains the default.

Fixes:

• Design.check warns about numeric ids (numeric_ids): a numeric variable of whole
  numbers, each shared by two samples or more, at least three distinct and not
  consecutive, whose name ends in an id word (donor, patient, subject, individual,
  participant, sample, batch, plate, lane, run, id, or their plurals, or a number or
  code after one of them), such as Kang's donors 101, 107 and 1015, which enter a
  model as one slope.
• gsea tells the sets that drew a few random scores of their sign, for which more
  permutations may help, from those that drew none, for which they are unlikely to.
• plot_significant_genes widens its limits so that the count of the longest bar
  stays inside the frame. It measures last, once its legend beside the plot has
  narrowed it, and again after the new limits change the tick labels, so that the
  counts stay inside panels side by side in a constrained layout.
• plot_sumrank breaks ties of p-values by the statistic, keeps its labels inside the
  frame (a label that would cross a side moves in) and apart (one that would overlap
  another moves up), with a line to the gene of a label that moved, and draws its
  legend beside the plot, clear of the labels along its top.
• A legend of a few entries (at most four) sits inside a scatter plot where it covers
  no point, line, bar or label, at the size the figure is laid out at, so that panels
  side by side keep their width; where every place inside would cover something, or
  the legend is longer, it goes beside the plot, as before. This applies to the MD,
  volcano, mean-variance, residual-variance, dispersion and quasi-dispersion plots and
  to MDS and PCA plots; bar charts keep their legends beside.
• plot_mds and plot_pca draw their samples with one helper.
• Integer layers keep their dtype where nothing is summed: groups only renamed by
  Cube.regroup, samples only renamed by collapse_samples, and so the aligned and
  stacked cubes of a collection (uint32 counts stayed uint32 counts); merged groups and
  pooled samples are summed in int64, as before. Fingerprints hash integers as int64,
  so they do not change.
• Reprs, Cube.history(), describe_provenance, DesignReport and SumRank's and
  gsea's warnings in lines of at most 80 characters, a long item continuing under its
  value; counts with thousands separators. The one-line reprs of records
  (VoomResult, ModeratedFit, ...) are unchanged.

Default:

• differential_expression_replication estimates π₀ for π₁ at λ = 0.5 by default
  (null_proportion=None is {"lambdas": 0.5}); null_proportion={} gives qvalue's
  defaults, the earlier default. proportion_of_true_nulls keeps qvalue's defaults.

Documentation and inventories:

• API pages score_tests.md and study_results.md, with description front matter,
  in the index after activity and meta_analysis; a sentence on fit in the
  differential expression page; the λ = 0.5 default in "Meta-analysis".
• "Replication across studies" asks for qvalue's smoother explicitly
  (null_proportion={}) where it compares it, and its advice on π₁ names the default.
• README: the status table names the additions.
• docs/coverage_exemptions.toml: the 13 new public names under [planned], each
  with the guide that will use it in the second pull request.
• tests/test_masked_inputs.py classifies every new callable and argument;
  tests/test_documentation.py knows the two new modules.

Choices, and why

• Reprs within 80 columns, for results and what the guides print. The reprs of
  results, of Cube, CubeCollection, Design, ModelMatrix, Contrasts,
  DesignReport, GeneSets and Network, Cube.history(), describe_provenance
  and the SumRank and gsea warnings wrap at 80 characters, a long item continuing
  under its value (or under the text after "1. "), so that continued lines are not
  read as new items; a name and its dtype stay together. The one-line reprs of records
  (VoomResult, ModeratedFit, TrueNullProportion, ...) stay one line: they are
  read as a line, and wrapping them changed nothing the guides print.
• Integer dtypes, as one rule. The rule sits in the derivation every operation
  shares, not in aligned_cube alone: an integer layer keeps its dtype whenever no
  two entries are summed, so a regroup that renames, a collapse_samples that pools
  nothing and the aligned and stacked cubes all keep it, and anything summed is int64,
  as before.
• A result's cube and design are keyword arguments of the constructor, both or
  neither, so that results built from their parts, as the meta-analysis tests build
  them, keep working; fit says why it cannot run without them.
• GroupFit.samples holds the design's sample columns as well as the computed
  ones, so that a figure of a group's samples can be coloured by any design variable.
• Trend figures of a result refuse what a method does not have: limma-trend fits
  no voom trend, and dream's moderated tests are not a ModeratedFit; each message
  names the figure that does apply.
• limma on scores names its last column log_odds, the name the next pull
  request gives limma's B in every table (b_statistic today), takes only contrasts
  of the design, and refuses activities of contrasts, which are tested statistics
  already.
• sample_groups works wherever expression_filter does: the three analyses of a
  cube, variance partitioning, and the scores and activities of a cube.
• with_columns records the names of the columns added and replaced, not their
  values: JSON cannot hold every value a column can (NaN, datetimes), and the new
  cube's fingerprint covers the values. collapse_samples of several columns joins
  their values with _ and refuses two combinations that join alike.
• Table.filter takes one value per column; a list is refused with the mask that
  does it (numpy.isin), since "equal to a list" is ambiguous for text. Calls that
  pass mask= by keyword keep working.
• Provenance of scores is None for values given as an array or for statistics
  (whose result keeps its own record), and is kept outside the records' fields, so
  their equality, dataclasses.replace and reprs are unchanged. Pickling carries it
  beside the fields, with the settings as a dict, which pickle can write where it
  cannot write their read-only view; a record without one stays without.
• Tables of tests take the columns of limma_voom's tables, with log_odds for
  limma's B as in limma, and ids feature_0, feature_1, ... when none are given,
  as from_bulk names them.
• Numeric ids are recognised by the last word of a column's name (digits aside),
  or a number or code after an id word: donor_id, subjectID, Batches and
  run_number say ids, while donor_age is an age and n_donors a count, which the
  id words matched anywhere in a name would flag.
• gsea's warnings are two, one per case, so that each gives one piece of advice;
  ten times the permutations gave a set of the first kind its p-value and left one of
  the second without, in the test that checks them.
• Legends inside a plot are checked, not guessed. matplotlib's "best" place
  minimises what a legend covers but counts a label as one point, so it may cover a
  label to spare three points. A legend goes inside only where it covers nothing at
  all, measured once the figure is laid out, and only when it has at most four
  entries; MDS and PCA legends of several styled variables, and SumRank's, whose
  labels run along the top, go beside. A figure that the caller lays out again later,
  by adding panels, can still move points under a legend placed inside.
• Counts and labels are placed last, once everything the function adds is there,
  and again when the new limits change the layout, so that what is measured is what
  is drawn.
• group_names enter the collection's record only when given, so that the
  records of existing collections, and the provenance of their aligned and stacked
  cubes, are unchanged.

Validation

Equivalence check. Master (cd5dbe9, before #67, which changes only the
documentation build) and this branch ran in separate processes, each
in its own environment, on the same inputs: every public entry point with its valid
call, the inputs of every validation script built from the same seeds with the calls
the scripts make (the replication analysis also with every option at its default),
edge cases, and master's whole test suite with every public call recorded. Every
output was compared element by element, NaN equal to NaN and any other difference
measured in units in the last place. The intended differences, each with the rule that
checks it, were written down before any result was read; a second run kept every
call's arguments, to show where the inputs of a call differ.

• 54,215 calls in 14,804 cases, 1,342,730 output parts with 192,180,438 elements;
  13,317 differences, each explained:

  Intended difference	Differences	How it was checked
  π₁ by default	12	in both trees each replication without null_proportion was also made with {} and with {"lambdas": 0.5}: master's default equals master's {} call and the new default the new {"lambdas": 0.5} call, every value; those calls agree across the trees but for their reprs
  integer dtypes kept	37, and 18 texts that also name the dtype	the same shape and values, master's int64 against the input's integer dtype; reprs and descriptions of lazy layers equal once dtype names are made alike
  80 columns and thousands separators	4,091	the same words once lines are joined and separators removed, and every new line at most 80 characters or a single word
  gsea's two warnings	2	the sets warned about sum to master's count; other warnings identical
  new properties (cube, design, provenance)	3,355	parts that only the new tree has, under those names; in two tests that unpickle a result, the calls that unpickling its cube and design makes, the other calls pairing up identically
  master's tests of the changed texts, dtypes and inventories	27 tests, and 5 calls in them	each failing test listed beforehand, failing on the change it meets
  refusal messages	6	the same exception, its message naming what is now accepted (DifferentialExpressionResult, sample_groups)
  the arguments of calls that take such objects	5,748	with the arguments kept, every argument that differs differs as a text, a dtype or a new property does above; Table.filter's mask, now positional-only, is passed by position and compared as master's keyword

  No call of master's inputs met the other intended differences (the numeric-ids
  warning, legend warnings, calls that master refused and the new tree accepts). The
  other 43 differences are calls that draw permutations without rng (gsea of a
  result, weighted_aggregate of a cube or a result, and their edge cases): their
  scores and p-values come from fresh entropy and differ between two runs of one tree,
  as the check's repeated calls found in both trees. Made again with rng=1 in both
  trees, the same 10 calls give 46 output arrays, all identical.

• Master's test suite: 4,931 passed in master; 4,904 passed and the 27 listed failed
  in this branch.

• Timing: of 1,703 calls of 1 ms or more, the new / old ratio has quantiles 0.51, 1.00
  and 1.30 (1%, 50%, 99%); 483 that looked slower were timed again back to back, and
  none is still slower by 20% and 0.3 ms.

• Not compared, in either tree: the inputs of validation/sumrank.py, whose script
  imports scikit-learn, which the environments lack (SumRank's tests, entry points and
  edge cases are compared); and two write_gmt calls to a path that exists, in both
  trees alike.

The changes to the figures and to pickling came after the check; they change no
computed value, and the tests below check them.

Validation rerun. validation/combined_p_values_and_replication.py gained a fourth
replication setting, null_proportion left out, compared with pi0est(lambda = 0.5)
(qvalue 2.38.0): 9 values of π₁, largest difference 0 (qvalue stops for the other 3
pairs, where its estimate is 0). The setting with null_proportion={} gives qvalue's
defaults, which stop for all 12 pairs here, as before, and the smoother is compared
with pi0est's defaults on 4 sets of p-values (1.3e-13). 55 checks, 0 failed;
docs/source/validation.md updated to the four settings.

Pages. Fourteen pages print or draw something these changes reach and were
executed again: the quick start, "AnnData and MuData", "Bulk RNA-seq", "Combine
studies", "Composition", "Export cubes", "Gene sets", "Changed genes per cell type",
"Normalization and quality", "Pathway activity", "Pool cells", "Repeated measures",
"Replication across studies" (whose code changed) and "Save and trace cubes". Every
code cell's printed text was compared with master's: every number the same; 16 cells
differ only in line breaks and thousands separators (the reprs of cubes, results,
networks and gene sets, design checks and histories); one prints gsea's new warning
("Gene sets"); the replication guide prints the same π₁ lines as before, now asking
for qvalue's smoother with null_proportion={}. Twelve figures on six pages differ
from master's, and each was looked at: the counts of plot_significant_genes inside
every frame, also in the two-panel figures; SumRank's labels inside the frame, its
legend beside; legends inside where they cover nothing (two volcano plots, the
residual variances, the PCA of one variable) and beside where they would cover a
point, a line or a label (the guides' MD and mean-variance plots, a volcano plot next
to an MD plot, the MDS and PCA plots of several variables). The other page, "Why
pseudobulk", passes make -C docs check-outputs.

Tests

• New tests, 64 functions in 23 files, each checked against something independent of
  the code under test: DifferentialExpressionResult.fit against the analysis's own
  tables for limma-voom, quality weights, limma-trend and dream, and against voom,
  fit_linear_model, contrasts_fit and moderated_statistics run by hand on a
  rank-deficient group; limma against the same building blocks with a design built
  from the scored samples, for scores, activities and arrays, with weights, trends and a
  donor without samples; analyse_each_study against the loop it replaces, for the
  three analyses; StudyResults.matrix against each study's own table placed at its
  features; Network.to_gene_sets against the comprehension the guides wrote and a
  CollecTRI file; sample_groups against filter_by_expression(group=...), and dream,
  limma-voom and limma-trend testing the same genes; with_columns and
  collapse_samples of several columns or of every sample against aggregation of the
  simulated cells; the DESeq2 and edgeR tables against the R columns the fixtures hold;
  Table.filter and sort against masks and stable argsorts; π₁ by default against
  Storey's estimate computed by hand; the numeric-ids warning on names and values that
  do and do not say ids; gsea's warnings, with ten times the permutations helping one
  set and not the other; pickling of scores, activities and principal components, with
  and without provenance.
• The figures are checked after the figure's final draw (saved to a buffer): every
  count inside its frame in two panels side by side with the eight cell types and
  counts of the Kang guide (constrained layout, shared groups, a legend each), which
  fails on the previous placement; every SumRank label inside its frame in the
  replication guide's two-panel layout, with genes at both ends of the alphabetical
  order; and legends inside plots covering no point, line vertex or label, going
  beside where points, a label or a line fill every place, or where they have more
  than four entries.
• tests/test_masked_inputs.py classifies every new callable and argument;
  tests/test_documentation.py knows the two new modules; master's tests of the
  changed texts, dtypes and inventories are updated (the 27 the check lists).
• uv run pytest -q --durations=15: 5074 passed, 18 skipped (zstd needs Python 3.14).
• Rounding changed (OPENBLAS_CORETYPE=Sandybridge, 4 threads,
  NUMBA_CPU_NAME=generic, a fresh cache): 5074 passed, 18 skipped.
• ruff check and ruff format --check of cytepool tests validation docs: clean;
  mypy cytepool: no issues in 83 files.
• make -C docs api (Sphinx with -W): clean. make -C docs html with Node 24.22.1:
  built; the only warning is Node's deprecation of url.parse in the site tool.
• make -C docs check-outputs on every page not executed again: current.
• tests/test_documentation.py with the docs group: 85 passed.