[parashardhapola]

Summary

Combine several cohorts of one disease, and find out which of their changes replicate:
two short guides under "Other designs and data", executed notebooks with the figures
users know from metafor and from the paper that introduced SumRank. "Combine studies"
aligns six cohorts' cubes, analyses each on its own and combines them by
meta-analysis, or stacks them into one cube for a joint model. "Replication across
studies" measures how well the cohorts replicate each other and which genes
replicate, with combined p-values, partial conjunction tests, π₁, SumRank and
prediction in cohorts left out. The cohorts are simulated, so every page says which
genes each method finds rightly: a new dataset of docs/datasets.py,
simulated_cohorts, six case-control cohorts with changes shared by all and changes
of one cohort alone.

Changes

• docs/source/guides/combine_studies.ipynb, "Combine studies" (8 code cells, 4
  figures, 6 s): the six cohort cubes as their sources wrote them; CubeCollection
  with repeated donor ids, genes that two cohorts did not measure (one cube lacks them,
  one writes zeros and marks them in measured) and cell types aligned by Cell
  Ontology id, with group_mapping rolling a cohort's naive and memory CD4 T cells up
  to their parent term; limma_voom on every aligned_cube, with
  plot_significant_genes for each cohort; align_study_results telling genes not
  measured from genes filtered; differential_expression_meta_analysis with
  plot_volcano; plot_forest of a gene changed in every cohort and of a gene
  changed in one cohort only, with the fixed-effect estimate beside the random-effects
  one; CubeCollection.stack and select with a joint ~ study + condition model;
  the four approaches against the truth; Choices.
• docs/source/guides/replication_across_studies.ipynb, "Replication across studies" (7
  code cells, 4 figures, 19 s, about half of it the permutations): every cohort's calls
  against the truth; differential_expression_replication with plot_replicated_genes
  (genes found in k cohorts); plot_study_agreement for π₁, with π₀ at λ = 0.5 and why,
  and the sign agreement; combine_p_values (Fisher) and partial_conjunction (u = 2
  and 3) against the replication result's weighted Stouffer p-values and the
  meta-analysis; SumRank of all six cohorts with plot_sumrank; the best 4 of 6 summed,
  with Irwin-Hall p-values and with 30 permutations; held_out_prediction for several
  proportions, with plot_held_out_auc; Choices.
• docs/datasets.py: simulated_cohorts() and its registry entry (a cube file per
  cohort and the true log2 fold changes, seed 20261028, fingerprints recorded).
• docs/source/myst.yml: both pages under "Other designs and data", after Bulk
  RNA-seq.
• docs/coverage_exemptions.toml: plot_forest, CubeCollection.select,
  plot_volcano and plot_significant_genes leave [planned]; the pages use them.
• .github/workflows/docs.yml: both pages added to check-outputs; they need no
  download.
• docs/source/reference/datasets.md: a row for simulated_cohorts.
• docs/source/index.md: the questions on replication and on pooled studies link the
  guides.
• docs/source/statistics/meta_analysis.md and statistics/cubes.md: one sentence
  each linking the guides.
• tests/test_documentation.py: simulated_cohorts added to the registered datasets.

Design notes

• Two pages, not one. Everything on one page would have run to about fifteen code
  cells. "Combine studies" answers how to get one result from several cohorts;
  "Replication across studies" answers which results hold up. Both read the same
  cohorts, so the cohorts are a registered dataset rather than a simulation repeated in
  each page.
• The simulation. Six cohorts of 12 to 60 donors (one sample each), four blood cell
  types, 2,000 genes, pseudobulk counts drawn directly. The disease changes 100
  monocyte genes in every cohort (80 alike, 20 by amounts that vary between cohorts),
  each case by a gamma-distributed activity of mean 1, so that cases differ in severity
  and the AUCs of cohorts left out stay below 1; each cohort also changes 10 monocyte
  genes of its own in every case, as a confounder of that cohort would. Effect sizes
  and donor variation make single cohorts find 17 to 62 of the 100 shared genes, and
  most genes are found by one cohort only, as the literature reports for real studies.
  The cohorts differ as real ones do: donor ids repeat, cohort_3's cube lacks 200
  genes, cohort_5 wrote zeros for 150 genes it did not measure and marks them in
  measured, cell types carry each cohort's own names with Cell Ontology ids,
  cohort_4 annotated its CD4 T cells as naive and memory cells, and cohort_6 has no NK
  cells. The ontology ids, and CL:0000624 as the parent of the two CD4 subsets, were
  checked in the EBI Ontology Lookup Service.
• What the pages show, all from executed output. The random-effects meta-analysis
  finds 80 of the 100 shared monocyte genes and none of the 60 single-cohort genes;
  the fixed effect 91 and 23; the stacked joint model 76 and 18, of the 82 shared genes
  it can test. Fisher's combination calls 36 single-cohort genes, weighted Stouffer 19,
  the partial conjunction tests none. SumRank of all six finds 55 shared genes and no
  single-cohort gene; summing the best 4 of 6, Irwin-Hall p-values call 16 unchanged
  genes (the function warns, and the page prints the warning), 30 permutations none.
  Mean AUC in cohorts left out: SumRank 0.92, meta-analysis 0.87.
• π₁ at λ = 0.5. qvalue's default π₀, a smoother read at λ = 0.95, gives π₀ = 1 on
  the few dozen p-values of the genes one cohort finds (π₁ of 0 for some pairs; qvalue
  2.38.0, run as a black box on the same 41 p-values, gives the same), so the page
  uses null_proportion={"lambdas": 0.5} and prints both ranges.
• SumRank's figure from Irwin-Hall p-values. With every cohort summed they are
  valid, and the Manhattan plot spreads as the authors' does; permutation p-values tie
  at their floor, and their labels pile up (see below). Permutations appear where they
  are needed, for the best 4 of 6, with 30 permutations, about 7 s of the page.
• Layout. The site shows about 80 characters of code and output, and scales
  figures to the width of the text: every line of code and printed output is at most
  80 characters, two-panel figures are stacked, and the reprs with longer first lines
  are left out in favour of short summaries.
• Literature numbers, each checked in the open-access full text (PMC12344038): over
  85% of the DEGs of one Alzheimer's dataset not found in any of 16 others; every study
  summed for Parkinson's disease, Huntington's disease and schizophrenia (6, 4 and 3
  studies), 65% and 55% for Alzheimer's disease (17) and COVID-19 (16), chosen by
  cross-validation; for Alzheimer's disease, mean AUCs in datasets left out of 0.78 for
  SumRank and for merged datasets, 0.74 for the inverse-variance meta-analysis and
  0.67 for single datasets, from Table 1.
• The aligned cell types keep their Cell Ontology ids as names in figures and tables;
  the first page prints each cohort's names of its cell types, and the cohorts' names
  under CL:0000624.

Checks

• make -C docs html (Node.js 24): builds; both pages render with their API links,
  citations and the backlinks of the API reference. Checked by eye in a headless
  browser at a width of 1300 pixels.
• make -C docs check-outputs for every page the workflow lists (quickstart,
  differential_expression, guides/export_cubes, guides/bulk_rna_seq,
  guides/combine_studies, guides/replication_across_studies): outputs current.
• uv run pytest tests/test_documentation.py: 76 passed.
• docs/coverage_report.py: 251 public names, 181 used or explained, 53 exempt, 17
  planned; no problems.
• ruff check and ruff format --check on cytepool tests validation docs, and
  mypy cytepool: clean.
• docs/datasets.py verify simulated_cohorts simulated_bulk simulated_pbmc: as
  registered.
• Full suite, on top of the PCA and release pull requests: passes, as part of the
  4,926 tests of the stack's last pull request; the pages that CI checks are current
  when their datasets are made from an empty cache.

Bugs found

None in the library's results. One stale docstring: the module description of
cytepool.replication still lists held-out prediction under "Not done", although
cytepool.held_out_prediction implements it:
python -c "import cytepool.replication as r; print(r.__doc__.split('Not done:')[1][:300])".

Awkward APIs, for the API and language pull request

• Analysing every study takes the same ten lines. Each guide loops over
  collection.studies, takes aligned_cube, builds a Design from its samples and
  runs limma_voom, collecting a dict. A collection-level call (or limma_voom
  accepting a collection and one formula) would remove it.
• Aligned cell types are named by ontology id. Aligning by cell_ontology_id
  names the groups CL:0001054 and so on in every table and figure, with no way to
  give them readable names (an agreed label column, or a group_names argument).
• Signatures differ across the functions of several studies.
  differential_expression_meta_analysis(results) takes results alone, while
  align_study_results, differential_expression_replication and
  differential_expression_sumrank take (collection, results, group, contrast), and
  held_out_prediction(collection, results, contrast, variable, case, control) asks
  for variable, case and control, which the designs and the contrast already
  imply.
• π₁'s default is unsteady where it is used. differential_expression_replication
  estimates π₀ by qvalue's smoother by default, which gives 0 or 1 on the few dozen
  genes a study finds; λ = 0.5 by default, or a warning below some number of genes,
  would serve users better.
• plot_sumrank labels. With permutation p-values many genes tie at the smallest
  p-value; the labels go to the first of them in alphabetical order and overlap, and
  close neighbours overlap with Irwin-Hall p-values too. Breaking ties by the
  statistic and moving labels apart, as plot_volcano does, would fix it.
• plot_forest needs the studies' results passed again besides the
  meta-analysis, and marks a study that did not measure the gene "not tested", like
  one whose filter left it out.
• Combining p-values by hand needs
  align_study_results(...)["p_value"].reshape(n_studies, n_features); a method that
  returns the studies x genes arrays would help.
• Proportions for SumRank read awkwardly for few studies
  (proportion_of_studies=4 / 6); a number of studies summed would be plainer.
• Renaming widens counts. aligned_cube renames groups through Cube.regroup,
  which returns integer layers as int64 even for a group that is only renamed, so
  aligned and stacked cubes of uint32 counts take twice the memory.
• Long reprs and messages. The first lines of the reprs of
  DifferentialExpressionResult, MetaAnalysisResult, ReplicationResult (one line
  of about 170 characters) and SumRankResult, and SumRank's warning, run past the 80
  characters that the site's output cells show; the result reprs also print feature
  counts without thousands separators, where the cube and collection reprs use them.

Merging with the other guide pull requests

• THIRD_PARTY_NOTICES.md says the documentation reads "two simulated" datasets; with
  simulated_cohorts it is three. The file is left as it is here, because the guides
  of "Compare conditions" change the lines next to it; the merged file should say
  "three simulated".
• docs/coverage_exemptions.toml: the guides of "Compare conditions" also remove
  plot_volcano and plot_significant_genes from [planned]; the merged file has
  neither, nor plot_forest and CubeCollection.select.
• .github/workflows/docs.yml: both pull requests append pages to the same
  check-outputs line; the merged line lists all of them.