Reproduction vs the paper
Median values for this data set at the primary setting (|log2FC| > 1, 5 % FDR). "Published" is the authors' table; "Reproduction" is PyDESeq2 on 10 random cohorts per size, with precision and recall measured against the authors' full-data DESeq2 result.
Table view
Which genes get found?
Each dot is one gene called differentially expressed in the full data (|log2FC| > 1, 5 % FDR; genes). Height is the share of 100 random cohorts of the chosen size that also called it, using the fast engine. Small cohorts mostly find only the genes with the very largest fold changes.
The same measure across the paper's 18 data sets
Median replicability reported in the authors' results table (their repository), sorted. The marked row is the data set reproduced here; dots show this reproduction's estimates.
Table view
How much do results vary from cohort to cohort?
The same 100 random cohorts per size as above, one dot each. Boxes show the median, quartiles and whiskers (1.5 x IQR); violins show the shape of the spread with the median (white dot) and quartiles (black bar). Replicability uses every cohort pair (4,950 per size), with 100 random pairs drawn as dots. Fast engine (blue), primary setting. Orange squares are the 10 PyDESeq2 cohorts per size that were run (sizes 3, 5, 7, 10, 15 only).
Table view
The tumours themselves
All 252 samples projected from their 2,000 most variable genes (PCA to 50 components, then t-SNE and UMAP on those components). Basal and Luminal A tumours separate cleanly. PCA axes are in standardised log2 CPM units (z-scored expression summed over genes). t-SNE and UMAP coordinates are arbitrary units: they come from an optimisation that preserves neighbourhoods, not distances, so they have no physical scale and only relative position within a panel means anything. Gaps and cluster sizes should not be read quantitatively. Optionally overlay one of the random cohorts used above to see how a small study samples the population.
Single genes: what a large effect looks like
The 40 reference genes with the largest standardised effect (full-data log2FC divided by its standard error), expression as log2 counts per million, one dot per tumour (126 per group). Names are HGNC gene symbols (Ensembl ID where no symbol exists); zero-count samples sit at the bottom of a panel.
Table view
Question
When a differential-expression study has few replicates, how much of what it reports would another study of the same size find again, and how much of the full truth does it recover? Degen & Medo (2025) answer this by repeatedly drawing small random cohorts from large RNA-seq data sets. This site reproduces that design on one of their data sets.
Methodology: what was actually done
- Data. TCGA breast cancer, Basal vs Luminal A, 126 + 126 tumours, 18,188 genes (the authors'
BSLAdata set, taken from their repository, not re-downloaded). - Cohorts. For each N in {3, 4, 5, 6, 7, 8, 9, 10, 12, 15} per group, 100 random cohorts drawn without replacement. Cohorts at the same N are independent draws and can share samples, as in the paper. Every draw is seeded from
(seed, N, trial)so reruns are identical. - Reference ("truth"). The full 252-sample analysis. Two versions: the authors' own full-data DESeq2 result (2,039 genes), and the same pipeline as ours on all samples (1,852 genes).
- Calling genes. Covariates tumour purity and age, formal test of |log2FC| > 1, Benjamini-Hochberg FDR below 5 % (primary setting; sensitivity over no filter / 0.5 / 1.0 and FDR 1 % / 5 % in the filters above).
- Metrics. Precision and recall against the reference. Replicability is the Jaccard index between the gene lists of two cohorts, over all pairs at a given N, summarised by the median (0 if either list is empty). Rank agreement is Spearman correlation of log2FC with full data.
- Two engines. PyDESeq2 (Wald test, formal threshold, covariates) — closest match to the paper, ~1 min/cohort, so 10 cohorts per size (50 cohorts, 45 pairs per size). A numpy engine in the style of limma-trend ran all 1,000 cohorts in ~1 min; used for distributions, gene-level view and threshold sensitivity. It is conservative below ~8 per group.
Findings
- 1. PyDESeq2 lands on the paper's numbers. At N = 5: precision, recall, replicability 0.84 / 0.15 / 0.20 vs published 0.84 / 0.14 / 0.19. At N = 10: 0.91 / 0.21 / 0.33 vs 0.90 / 0.25 / 0.30. With 10 cohorts per size these medians are noisy — agreement is good, not exact.
- 2. Replicability stays low. PyDESeq2 median overlap: 0.14 at N = 3, 0.20 at 5, 0.29 at 7, 0.33 at 10, 0.41 at 15. Two 15-vs-15 studies agree on fewer than half the genes either reports.
- 3. Small studies are precise but incomplete. Precision 0.71 at N = 3, ≥0.84 from N = 5. Recall 0.14 at N = 3, only 0.35 at N = 15 — mostly real genes, but most missed.
- 4. Which genes survive is not random. At N = 10, the 128 reference genes (of 1,852) called by ≥half the cohorts have median |log2FC| 4.78 vs 2.02 for the whole list. Small cohorts find very large effects.
- 5. Not a threshold artefact. Across six settings (no filter / 0.5 / 1.0; FDR 1 % / 5 %), fast-engine replicability at N = 10 stays 0.28–0.35.
- 6. Engine caveat. Fast engine calls no genes in 64 %, 47 %, 21 % of cohorts at N = 3, 4, 5, pinning its replicability near 0 there. At N = 7, 10, 15 it is within 0.04 of PyDESeq2. Its small-N curve should not be read as biology.
Detailed results tables
Primary setting: |log2FC| > 1, 5 % FDR. Precision/recall vs authors' full-data DESeq2 unless noted.
Table 1 — Fast engine, 100 cohorts per size (medians)
| N | Median genes | No calls | Precision | Recall | Replicability | Rank agreement |
|---|---|---|---|---|---|---|
| 3 | 0 | 64 % | 0.86 | 0.00 | 0.00 | 0.61 |
| 4 | 1 | 47 % | 1.00 | 0.00 | 0.00 | 0.72 |
| 5 | 22 | 21 % | 0.96 | 0.01 | 0.02 | 0.75 |
| 6 | 54 | 4 % | 0.96 | 0.03 | 0.13 | 0.80 |
| 7 | 96 | 1 % | 0.96 | 0.05 | 0.21 | 0.82 |
| 8 | 128 | 0 % | 0.95 | 0.06 | 0.24 | 0.84 |
| 9 | 175 | 0 % | 0.95 | 0.08 | 0.28 | 0.85 |
| 10 | 200 | 0 % | 0.95 | 0.09 | 0.32 | 0.86 |
| 12 | 322 | 0 % | 0.95 | 0.15 | 0.37 | 0.89 |
| 15 | 444 | 0 % | 0.95 | 0.21 | 0.44 | 0.91 |
Table 2 — PyDESeq2 anchor vs fast engine, same cohorts
| N | Cohorts | PyDESeq2 genes | PyDESeq2 prec. | PyDESeq2 rec. | PyDESeq2 repl. | Fast repl. |
|---|---|---|---|---|---|---|
| 3 | 10 | 423 | 0.71 | 0.14 | 0.14 | 0.00 |
| 5 | 10 | 365 | 0.84 | 0.15 | 0.20 | 0.02 |
| 7 | 10 | 519 | 0.86 | 0.22 | 0.29 | 0.25 |
| 10 | 10 | 465 | 0.91 | 0.21 | 0.33 | 0.33 |
| 15 | 10 | 778 | 0.91 | 0.35 | 0.41 | 0.43 |
Table 3 — Authors' published values (row BSLA)
| N | Precision | Recall | Replicability | Rank agreement (mean) |
|---|---|---|---|---|
| 5 | 0.84 | 0.14 | 0.19 | 0.59 |
| 10 | 0.90 | 0.25 | 0.30 | 0.73 |
Table 4 — Threshold sensitivity at N = 10 (fast engine, own reference)
| FC threshold | FDR | Ref genes | Median called | Precision | Recall | Replicability |
|---|---|---|---|---|---|---|
| none | 1 % | 12146 | 1100 | 1.00 | 0.09 | 0.30 |
| none | 5 % | 13498 | 2648 | 0.99 | 0.19 | 0.35 |
| |log2FC| > 0.5 | 1 % | 4332 | 276 | 1.00 | 0.06 | 0.29 |
| |log2FC| > 0.5 | 5 % | 4928 | 600 | 0.99 | 0.12 | 0.32 |
| |log2FC| > 1.0 | 1 % | 1630 | 102 | 1.00 | 0.06 | 0.28 |
| |log2FC| > 1.0 | 5 % | 1852 | 200 | 1.00 | 0.11 | 0.32 |
How well do the engines agree with the authors?
R is not available in the environment this was built in, so the paper's DESeq2 and edgeR could not be run. Two stand-ins were used, and both were checked against results the authors published in their repository.
- Fast engine (blue): a numpy re-implementation in the style of limma-trend with a formal fold-change test and the same covariates. It runs 1,000 cohorts in about a minute, but it is conservative below about 8 samples per group, so its small-cohort numbers understate what DESeq2 or edgeR would find.
- PyDESeq2 anchor (orange): the Python port of DESeq2 with covariates and a formal |log2FC| > 1 test, the closest available match to the paper's main pipeline. It takes minutes per cohort here, so only 10 cohorts per size (45 cohort pairs) were run and the ranges are wide.
What this does and does not show
- One data set of 18. Gene-set enrichment, the paper's second analysis, is not reproduced. The GEO data sets are outside this environment's network allowlist.
- The reference ("truth") is the full-data analysis, and the cohorts are subsets of those same samples, as in the paper. Precision at large N is therefore somewhat flattering.
- Replicability is the Jaccard index between the DEG lists of two cohorts, taken over all pairs of random cohorts at a given N and summarised by the median. Cohorts at the same N can share samples. The list of cohort sizes (3 to 15, ten values) is an assumption.
- "Published" values come from the authors' results table in their repository (row
BSLA), which stores N = 5 and N = 10 only. Column meanings were not cross-checked against the paper's Table 1. - Anchor and published points exist for the primary setting only (|log2FC| > 1, FDR 5 %).
Practical reading
For a study of this kind (a heterogeneous tumour comparison), 5 to 10 replicates per group yields gene lists that are mostly correct but cover only about a fifth to a quarter of the full signal, and a repeat study would find only about a fifth to a third of the same genes. The paper's bootstrap tool (pdegen/BootstrapSeq) is the authors' suggested way to estimate this for your own small data set.
Reproduce it
pip install numpy pandas scipy pyyaml pydeseq2 python workflow/run_sweep.py config/config.yaml # ~1 min: 1,000 cohorts python workflow/pydeseq2_crosscheck.py config/config.yaml --sizes 3 5 7 10 15 --k 10 # hours; resumable python workflow/validate_engine.py config/config.yaml python workflow/build_dashboard.py config/config.yaml # writes docs/index.html