Methodology — how these figures were made (in detail)
1. Question. When a differential-expression study has few replicates, how much would a repeat study of the same size find again (replicability), and how much of the full-data truth does one small study recover (precision/recall)? Design follows Degen & Medo (2025): repeatedly draw small random cohorts from a large RNA-seq data set, run DE, and compare cohort-vs-cohort and cohort-vs-full-data.
2. Data provenance. TCGA breast cancer, Basal vs Luminal A, 126 + 126 tumours, 18,188 genes after min_total_count: 10 filtering — the authors' BSLA set, taken from pdegen/BootstrapSeq (MIT), not re-downloaded. Files: data/BSLA.csv (genes × samples counts), data/BSLA.meta.csv (Condition, TumorPurity, days_to_birth), data/upstream_BSLA.deseq2.lfc1.csv (authors' full-data DESeq2: logCPM, logFC, lfcSE, pvalue, FDR), data/published_degen_medo_2025_results.csv (their 18-dataset table, row BSLA), data/gene_annotation.csv (16,879/18,188 symbols). All parameters in config/config.yaml (seed: 20261012).
3. Preprocessing (both engines). Size-factor normalisation → log2 CPM with prior 0.5 (as the DE engine) for embeddings; raw counts + covariates for modelling. Genes kept per-cohort if summed count ≥ 10.
4. Cohort design (workflow/cohorts.py). N per group in {3, 4, 5, 6, 7, 8, 9, 10, 12, 15} × 100 trials = 1,000 cohorts, drawn without replacement within a cohort, seeded from (seed, N, trial) so reruns are identical. Same-N cohorts are independent draws and can share samples, as in the paper. Pairwise Jaccard uses all 4,950 pairs per N at the primary setting (results/pair_jaccard.csv.gz).
5. DE models. (a) Fast numpy engine (workflow/fast_de.py, limma-trend style): weighted least squares on log2 CPM with purity + age covariates, empirical-Bayes moderated t, formal threshold test H0: |log2FC| ≤ L against L ∈ {0, 0.5, 1.0}, BH-FDR ∈ {1%, 5%}. Runs 1,000 cohorts in ~1 min but is conservative below ~8/group (64%/47%/21% no-calls at N = 3/4/5). (b) PyDESeq2 anchor (workflow/pydeseq2_crosscheck.py): Python port of DESeq2, Wald + formal |log2FC|>1 test with covariates, ~1 min/cohort → 10 cohorts/size at N = 3, 5, 7, 10, 15 only (45 pairs/size; results/pydeseq2_anchor_sets.json). R/DESeq2/edgeR were unavailable at build time.
6. References (“truth”). Two full-252 analyses: authors' DESeq2 (2,039 DEGs at |log2FC|>1, FDR 5%) and own fast-engine full-data call (1,852 DEGs). Primary precision/recall in the gallery use the authors' reference.
- Metrics (
workflow/run_sweep.py). Per cohort:precision = |DEG_cohort ∩ DEG_ref| / |DEG_cohort|(NaN/undefined if no calls, excluded from medians);recall = |∩| / |DEG_ref|;replicability = |A∩B| / |A∪B|(Jaccard; 0 if either list empty), median over all pairs;spearman_truth= Spearman(logFC_cohort, logFC_full). Summaries: median + IQR (bands/boxes), whiskers 1.5×IQR, violins via Gaussian KDE. - Threshold grid. 6 settings (LFC × FDR) in
results/metrics_by_cohort.csv(6,000 rows); primary L=1.0, P=0.05. Volcano/MA sections show 4 of the 6. - Embeddings (
workflow/embeddings.py, seed 7). 2,000 most variable genes → z-score → PCA-50 → t-SNE (perplexity 30, PCA init) + UMAP (n_neighbors 15, min_dist 0.3) on the 50 PCs. Stats: PC1 23.0%, PC1–10 49.8%, silhouettes PCA2 0.587 / t-SNE 0.794 / UMAP 0.803; PC1–purity r = 0.26, PC1–age r = −0.17. Top-40 box genes = full-data DEGs with largest |log2FC|/SE. - 100 figures by section. Volcano (4) from authors'
logFC/pvalue/FDR; MA (4) fromlogCPM/logFC; Heatmap (4) fromtop_gene_expression.csv(log2 CPM, raw + row z-score); PCA/t-SNE/UMAP scatter (10 = 1 overview + 9 projection×colouring) fromembeddings.csv; Box+dots (44 = 4 distributions + 40 genes); Violin (4 distributions); Dot/recurrence (10, one per N) fromgene_recurrence.csv.gz; Bar (7 = 1 published + 4 thresholds + 2 rankings); Line (1 metrics trend); QC/diagnostics (6: 2 p-hist + scree + 2 PC-covariate + library size); Curves (2 PR at N = 5/10); Histograms (4 Jaccard at N = 3/5/10/15). - Units. PCA: standardised log2 CPM units. t-SNE/UMAP: arbitrary units (neighbourhoods only — gaps/sizes not quantitative). Expression: log2 CPM (zeros at panel bottom). Volcano y = −log10 p (capped at 300); MA x = logCPM (log axis).
pip install -r requirements.txt python workflow/run_sweep.py config/config.yaml # ~1 min: 1,000 cohorts python workflow/embeddings.py config/config.yaml # PCA / t-SNE / UMAP + top genes python workflow/validate_engine.py config/config.yaml # report/validation.json python workflow/make_gallery_figures.py config/config.yaml # 70 PNGs python workflow/make_gallery_bioplots.py # +8 volcano/heatmap → 82 python workflow/make_gallery_more.py # +18 MA/QC/curves → 100 python workflow/build_figures_page.py # this page → docs/figures.html
Validation (see dashboard). Full-data fast vs authors' DESeq2: Spearman logFC 0.965, 1,546/2,039 shared (Jaccard 0.659). Authors' N = 5 edgeR cohort: Spearman 0.966; engine 48 vs edgeR 340 DEGs (48 shared — conservative at small N). Fast-vs-PyDESeq2 medians agree within 0.04 at N = 7/10/15. Limitations. One of 18 data sets; no enrichment analysis; GEO sets outside network allowlist; reference = same-sample full-data (precision flattering at large N); anchor = 10 cohorts/size (noisy medians); cohort-size list + BSLA label taken from authors' repo, not paper Table 1.