Methodology
Everything below describes the state of this repository: inputs, environment, every parameter, every processing step with library and function names, every output field, every dashboard figure, and every known deviation from Mereu, Lafzi et al., Nat Biotechnol 38:747–755 (2020) (paper, open preprint 630087v1, authors' R code).
- Study design (the paper)
- Input data: accessions, files, formats, quirks
- Environment: language, libraries, versions
- Configuration reference (every key)
- Pipeline scripts, in order
- Analysis internals: formulas and algorithms
- Outputs: results.json schema and friends
- Dashboard: every tab and figure (R1–R9)
- Metric definitions
- Paper vs replication map (Figs. 1–6)
- Limitations and known artefacts
- Generalizability: tissues over species
- Reproduce it / extend it
- Results snapshot (live data)
1. Study design (the paper)
A single reference sample was prepared in one batch (250,000 viable cells per 0.5 ml aliquot), cryopreserved, and distributed to the partners. Design composition (% of viable cells):
Human PBMC 60% · mouse colon 30% · HEK293T 6% (human, RFP-labelled) · NIH3T3 3% (mouse, GFP-labelled) · MDCK 1% (dog, TurboFP650-labelled). Mixed species allow doublet exclusion; spiked-in lines expose batch/capture bias. Expected species mix per protocol (paper Fig. 1): human 63–95%, mouse 4–34%, dog 1–9%.
13 sc/snRNA-seq protocols (~3,000 cells each; 5 plate-based, 7 microfluidic — 4 droplet, 1 nanowell, 2 integrated fluidic circuit C1HT-small/medium — plus 1 single-nucleus experiment), plus a large Chromium V2 reference (30,807 human + 19,749 mouse cells) used only to define cell types and markers (Table S2), never to score protocols. Downsampling standard: 20,000 reads/cell (lowest first quartile) via the zUMIs downsampling function, plus stepwise 25/50/75% fractions and equal cell numbers. Stratified cell types for human: HEK293T, monocytes, B cells; mouse: secretory and transit-amplifying (TA) colon cells.
2. Input data: accessions, files, formats, quirks
2.1 Where the data live
GSE133549 is a SuperSeries with no files of its own; count matrices live in 15 per-protocol SubSeries (plus GSE141469 for 10X V3), retrieved 2026-10-13 into data/raw/ and processed the same day. Base URL pattern:
https://ftp.ncbi.nlm.nih.gov/geo/series/{PREFIX}nnn/{ACC}/suppl/
# e.g. PREFIX=GSE133, ACC=GSE133540 ; GSE141469 uses PREFIX=GSE141
Each SubSeries ships 4 gzipped TSVs: {ACC}_{PROTO}_human_exp_mat.tsv.gz, …_human_metada.tsv.gz (note the authors' truncation, single "t"), and the mouse pair. This repo processes the human matrices only (largest table per SubSeries, picked automatically — §5.2).
2.2 Processed matrices (data/processed/)
| File | SubSeries | Genes (rows) | Cells (cols) | Note |
|---|---|---|---|---|
| 10X_2x5K_250K.tsv.gz | GSE133535 | 38,418 | 10,000 | Chromium V2, deep |
| 10X_8x10K_25K.tsv.gz | GSE133536 | 39,689 | 80,000 | Chromium V2, shallow; largest file (~147 MB gz) |
| C1HT_medium.tsv.gz | GSE133537 | 41,790 | 5,600 | Read counts (non-UMI) |
| C1HT_small.tsv.gz | GSE133538 | 39,588 | 6,400 | Read counts (non-UMI) |
| CELseq2.tsv.gz | GSE133539 | 41,289 | 2,880 | UMI |
| Dropseq.tsv.gz | GSE133540 | 33,921 | 6,000 | UMI |
| ICELL8.tsv.gz | GSE133541 | 34,271 | 3,928 | Read counts (non-UMI) |
| MARSseq.tsv.gz | GSE133542 | 34,408 | 2,304 | UMI |
| QUARTZseq2.tsv.gz | GSE133543 | 36,595 | 3,070 | UMI |
| gmSCRBseq.tsv.gz | GSE133544 | 34,733 | 3,455 | UMI |
| SMARTseq2.tsv.gz | GSE133545 | 44,120 | 2,992 | Read counts (non-UMI) |
| SingleNuclei.tsv.gz | GSE133546 | 33,455 | 6,000 | snRNA-seq, UMI |
| ddSEQ.tsv.gz | GSE133547 | 33,895 | 4,000 | UMI |
| inDrop.tsv.gz | GSE133548 | 28,880 | 2,988 | UMI |
| 10X_V3.tsv.gz | GSE141469 | 39,202 | 9,035 | Chromium V3; postdates the paper, unranked |
2.2b Mouse track
The same 15 SubSeries ship mouse colon matrices (…_mouse_exp_mat.tsv.gz), prepared with PROCESSED_DIR=data/processed_mouse SPECIES=mouse into an identical per-protocol layout and analysed with config.mouse.yaml (same parameters; output site/data/results_mouse.json + .js fallback). The dashboard's Human/Mouse toggle switches datasets without reload. Gene IDs are ENSMUSG… (same mito-prefix artefact); paper Fig. 6 ranks shown on mouse views are human-derived.
2.3 Format quirks (load-bearing)
- No corner cell: the header line holds N barcodes while data lines hold a gene ID + N counts.
pandas.read_csv(..., index_col=0)resolves this correctly (index name = first barcode, gene IDs as index); naiveusecols=[0]reads mistake the first count column for row names. - Gene identifiers are ENSEMBL IDs with versions (e.g.
ENSG00000000003.14), not symbols — so theMT-/MT.mitochondrial prefix match finds nothing (§11). - Integer UMI/read counts (verified: all values integral) — raw, unnormalised, as the pipeline assumes.
- Metadata has no cell-type labels. Columns are QC/mapping metrics (
nTReads, nExonReads, nIntronReads, nIntergenicReads, nUnmappedReads, nAmbiguityReads, nMultimapReads, nUMIs, nGenes, Species, Library, Donor, nVars, is_cell_control, total_features_by_counts, …). Hencelabel_coverage = 0and ARI/NMI/silhouette are skipped until an externaldata/metadata.csv(index = cell barcode, columncelltype) is supplied; the authors' labeller is matchSCore2.
3. Environment: language, libraries, versions
Python 3.9.6 (/usr/bin/python3), isolated virtualenv .venv built by bash run.sh setup (python3 -m venv .venv, then pip install -r requirements.txt). Deliberately no compiled single-cell stack (no scanpy/seurat) — every step below uses only:
| Library (version) | Used for (exact API) |
|---|---|
| numpy 2.0.2 | RNG (default_rng(seed)), subsampling (rng.choice), medians/percentiles, log/standardisation arithmetic |
| pandas 2.3.3 | Chunked TSV reading (read_csv(..., chunksize=1000/2000)), metadata reindexing, composition value_counts |
| scipy 1.13.1 | Sparse cells×genes matrices (scipy.sparse.csr_matrix, vstack), 10x MTX reading (scipy.io.mmread) |
| scikit-learn 1.6.1 | PCA via TruncatedSVD(n_components, random_state), kNN graph (kneighbors_graph(k, include_self=False)), t-SNE (TSNE(n_components=2, perplexity, init="pca", random_state)), metrics (adjusted_rand_score, normalized_mutual_info_score, silhouette_score); UMAP only if umap-learn present (see next row) |
| umap-learn 0.5.12 (optional) | UMAP(n_components=2, random_state) on the 30-PC embedding. Import is guarded: without it the pipeline still runs and stores umap: null. |
| networkx 3.2.1 | Graph (from_scipy_sparse_array on the symmetrised kNN graph) + Louvain (community.louvain_communities(G, resolution, seed)) |
| PyYAML 6.0.3 | yaml.safe_load of config.yaml / config.synthetic.yaml |
Stdlib only otherwise (json, glob, os, sys, datetime, tarfile, shutil). Every analysis run embeds its software environment in results.json: environment (pipeline/run_analysis.py:_env: Python version plus numpy/pandas/scipy/scikit-learn/networkx/PyYAML/umap-learn versions and the UMAP-availability flag), so the table above can be audited per run instead of trusted from documentation. Dev-only: pypdf (paper text extraction) and Pillow (docs) are not needed to run the pipeline. Known numerical note: TruncatedSVD's randomized solver prints RuntimeWarning overflow/invalid lines on these matrices; output is unaffected (all protocols complete, metrics finite).
4. Configuration reference (every key)
config.yaml (real data) and config.synthetic.yaml (smoke test; smaller min_genes/max_cells/tsne_cells/repeats) share these keys:
| Key (default) | Meaning |
|---|---|
metadata (data/metadata.csv) | Label file; index = cell barcode, needs column celltype_col. Absent → label metrics skipped (current state). |
celltype_col (celltype) | Label column name inside metadata. |
matrix_dir (data/processed) | One file or 10x folder per protocol; basename minus extension = protocol name. |
orientation (auto) | genes_by_cells | cells_by_genes | auto (auto: rows>cols ⇒ genes×cells, then transposed). |
protocols ([]) | Empty = every entry in matrix_dir; else an allow-list. |
output (site/data/results.json) | JSON written at the end; a same-named .js fallback is auto-written beside it. |
mito_prefixes (["MT-","MT."]) | Case-insensitive gene-symbol prefixes for the mito filter — inert on ENSEMBL IDs. |
| min_genes (200) / max_mito_pct (20) / min_cells (3) | QC gates: per-cell genes ≥ 200, mito % ≤ 20, per-gene expressing cells ≥ 3. |
| equalize_cells (true) / max_cells (5000) | Every protocol subsampled without replacement to min(QC-pass counts, 5000). |
| n_hvg (2000) / n_pcs (30) / knn (15) / resolution (1.0) | Top-variance genes → SVD components (capped at genes−1) → kNN degree → Louvain granularity. |
| stability_repeats (5) / tsne_cells (1500) / seed (0) | 80%-subsample repeats; t-SNE cell cap; master seed for RNG, SVD, Louvain, t-SNE. |
5. Pipeline scripts, in order
run.sh — one entry point
setup | download | prepare | inspect | analyse | test | serve. serve runs python -m http.server -d site 8000 (serve site/, never the repo root; never double-click the HTML — file:// blocks data loading).
pipeline/download_geo.sh list|get — fetch from NCBI
list prints supplementary files per SubSeries; get ACC… downloads them to data/raw/{ACC}/. Fixed in this repo to skip non-data links (an hhs.gov href used to abort multi-accession downloads under set -euo pipefail).
pipeline/prepare_data.py ACC… — raw → processed
Extracts *.tar*, picks the largest table file per SubSeries (always the human matrix here), copies to data/processed/{protocol}{ext} via a hardcoded accession→protocol MAP (GSE133535→10X_2x5K_250K … GSE141469→10X_V3). Prints chosen vs alternative files — verify the choice.
pipeline/inspect_data.py — diagnose before analysis
Per file: rows/cols, numeric-column fraction, integer-count check (False ⇒ already normalised ⇒ assumptions break), orientation guess, metadata barcode overlap. Its row-name display misreads corner-less TSVs (shows nan); the analyser is unaffected.
pipeline/run_analysis.py [config] — the analysis
Sections §6.1–§6.8: chunked load → QC → equalise → normalise/HVG/PCA → kNN+Louvain → score → results.json (+ auto .js fallback). Preserves previously curated discrepancies across re-runs. Guards: protocols with <50 cells or <100 genes post-QC are skipped into errors, not crashed on.
pipeline/qc_distributions.py [results.json] [--dtype float32] — pre-QC distributions
Second chunked pass (chunksize=2000, float64 accumulators) over the same matrices: per-cell genes/UMIs percentiles, genes seen in ≥1/≥3 cells, top-50/200 count shares, top-5 gene IDs. Merged into each protocol's qc_dist (powers Quality Figs. R6–R8), including a ≤1000-cell per-cell sample (cell_sample{genes,umis}) backing the R13 violins. --dtype float32 halves peak chunk memory (≈2× faster on large matrices) with identical rounded outputs; default float64 preserves existing results bit-for-bit.
pipeline/similarity.py [results.json] — cross-protocol similarity
Two chunked passes per matrix: (1) per-cell gene counts → QC mask (≥ min_genes) → fixed ≤100-cell sample (seed 0); (2) per-gene sums over QC cells + sampled presence submatrix. Then: gene intersection across protocols (common genes) → mean CP10K per gene → log1p → Pearson matrix; cumulative union genes over 1–100 cells, mean of 3 shuffles. Writes d["similarity"] (+ .js).
pipeline/make_synthetic.py + tests/test_pipeline.sh — smoke test
Generates 3 synthetic protocols (1200 genes; depths 4000/1500/700; 400 cells each) into data/synthetic/, runs the analyser with config.synthetic.yaml → site/data/demo.json, asserts 3 protocols and ARI>0.5. Output is never real data; the dashboard's ?demo mode loads it.
6. Analysis internals: formulas and algorithms
6.1 Load
Text matrices stream in 1000-row gene blocks into a CSR cells×genes matrix (float32). 10x folders: mmread(matrix.mtx).T with features.tsv/genes.tsv + barcodes.tsv. Orientation: auto treats input as genes×cells when rows>cols and transposes. Memory-safe: peak ≈ one block (e.g. 1000×80,000 float32 ≈ 320 MB for the largest 10X set).
6.2 QC
Per cell: n_genes = #(X>0), umis = ΣX, mito% = 100·ΣX_mito/max(umis,1) (0 everywhere here — §2.3). Keep (n_genes ≥ 200) & (mito% ≤ 20); then drop genes seen in <3 cells. Medians of the survivors are the reported median_genes/median_umis/median_mito_pct.
6.3 Normalise → HVG → PCA
6.4 kNN graph → Louvain
kneighbors_graph(PCs, k=min(15, n−1), include_self=False), symmetrised by element-wise maximum with its transpose, converted with from_scipy_sparse_array, partitioned by louvain_communities(resolution=1.0, seed=0); community index = cluster label; count = n_clusters.
6.5 Equalisation, labels, stability
Each protocol is subsampled (seeded, no replacement) to n_use = min(5000, min QC-pass counts) — currently 1,451 (inDrop is smallest). Reference labels (when present) reindex by barcode; ARI/NMI/silhouette (silhouette on ≤2000 cells) need >50% labelled cells. Stability (always computed): 5× draw 80% of cells, re-embed and re-cluster from counts, ARI vs full-data labels restricted to the subsample; report mean ± SD. t-SNE: ≤1500 cells, perplexity=min(30,(n−1)/3), init="pca", seed 0, coordinates rounded to 2 dp.
6.5b Stored embeddings
For the t-SNE cell subset (t, ≤1500 cells): PCA scatter (P[t,:2] + the full 30-component explained_variance_ratio_), UMAP coordinates (or null), and cellstat (cluster + per-cell genes/umis from the post-QC arrays). These power the Explorer views and the R14/R15 cluster boxes. t-SNE perplexity = min(30, (n−1)/3), init "pca", seed 0.
6.6 Cluster markers
On the equalised cells, log1p(CP10K) values over the top-n_hvg variable genes: per Louvain community, score = (meanin − meanout) × fraction expressing; top 5 kept. The union set (≤5 genes × communities) is stored with per-cluster means/fractions for the dashboard dot plot (Fig. R10). Deterministic — uses no RNG, so equalisation, stability draws and t-SNE are bit-identical with or without this block. Names stay as-matrix (ENSEMBL IDs); pipeline/map_symbols.py adds symbols when resolvable.
6.7 Synthetic data (test only)
4 cell types × 1200 genes; type-signal genes amplified 8×; depths 4000/1500/700 via Poisson sampling; 400 cells/protocol; MT-0…12 synthetic mito genes.
7. Outputs: results.json schema and friends
| Path | Content |
|---|---|
site/data/results.json | {generated, config, cells_per_protocol, protocols[], discrepancies[], errors[]}. Per protocol: protocol, cells_pass_qc, median_genes, median_umis, median_mito_pct, n_clusters, label_coverage, optional ari/nmi/silhouette/composition, stability_mean/sd, qc_dist{input_cells, genes_total_1cell/3cells, genes_q[7], umis_q[7], top50/200_share, top5_genes} (percentiles = min,P10,P25,median,P75,P90,max), tsne{x,y,cluster,label}. |
markers (per protocol) | {genes[union], top{cluster:[idx]}, clusters{cluster:{means[], fracs[]}}} — top-5 cluster-associated HVGs + the matrix behind the R10 dot plot; optional symbols{ensg: symbol}. |
pca / umap / cellstat | t-subset coordinates + PC variances + per-cell cluster/genes/umis (Explorer views, R14/R15 boxes). |
similarity | {genes_common, qc_note, order[], correlation[15×15], cumulative{protocol:[100]}} — powers Similarity Figs. R11–R12. |
site/data/results.js | Auto-written alongside: window.RESULTS = {...} — offline fallback so the page works over file://. |
site/data/demo.json / demo.js | Same schema, synthetic; loaded only with ?demo (plus embedded fallback). |
site/index.html | The dashboard (9 tabs, figures R1–R9, §8). Loads data/results.json (tries data/ → ./data/ → site/data/), then the embedded fallback, else a diagnostics box. |
8. Dashboard: every tab and figure (R1–R59 + P1–P13)
| Tab | Content |
|---|---|
| Overview | Hero (live complexity order + stability range), KPI cards, sortable metrics table, CSV download, auto-computed key findings, data dictionary. |
| Quality | R1a/R1b log-scale median bars (genes, UMIs) · R2 log–log bubble plot (area ∝ QC-pass cells; labels de-overlapped with leader lines, edge-flipped) · R3 input→QC-pass retention bars · genes-per-1K-UMIs efficiency · R6/R7 pre-QC per-cell distribution plots (whiskers min–max, box IQR, thick median, thin P10/P90) · R8 top-50/200 concentration bars · R13a/b KDE violin distributions (rotated labels) (≤1000 QC-pass cells, Silverman bandwidth, IQR bar + median dot). |
| Clustering | R4a stability mean ± SD (axis fixed 0.5–1.0) · R4b Louvain community counts · honest empty-state for ARI/NMI/silhouette · per-protocol cluster-size distributions from the t-SNE subset · R14/R15 Tukey box plots of genes/UMIs per cluster (same subset, protocol selector) · R16a–c embedding gallery (t-SNE/PCA/UMAP side-by-side, same cells) · R16d PCA scree with 80% mark · R17 communities-vs-genes scatter with Pearson r · R53a/b custom-Louvain vs scanpy-Leiden paired bars (partition count k; label-agreement ARI) · R55 NMI parity · R56 silhouette parity · R57 marker Jaccard ranking · R58 union marker counts. |
| Benchmark | R5 Fig. 6-style dot plot (5 min–max metrics, vertical headers, orange arithmetic-mean composite + paper ranks) with value table · R9 paper-vs-replication rank bump chart. |
| Generalizability | R18 human-vs-mouse median-genes scatter (log–log + identity, Pearson r) · R19 human-vs-mouse rank bump · R20 stability human-vs-mouse paired dots · R27 median-counts concordance · R28 stability gap (mouse − human) · R29 human-vs-mouse ARI · R30 per-species composite-rank bump. Joint 15 protocols; needs both JSONs. |
| Signatures | R21 depth-adjusted complexity (residual genes vs depth trend) · R22 cell heterogeneity (CV of genes/UMIs from cellstat) · R23 cluster evenness (k vs effective k = exp entropy) · R24 accuracy vs precision (ARI vs stability, size ∝ k) · R25 species capture (stacked Human/Mouse/Doublet/Dog) · R26 marker definition (union markers vs strength, size ∝ k) · R31 count efficiency · R32 capture skew (mean/median) · R33 rarefaction increment (50→100 cells) · R34 k vs stability · R35 silhouette vs ARI · R36 NMI vs ARI · R37 stability SD · R38 cluster dominance · R39 total PCA variance · R40 markers per cluster · R41 marker strength ranking · R42 marker-set Jaccard heatmap · R43 expression hubness · R44 rarefaction endpoint · R45 early saturation · R46 doublet rate · R47 spike-in recovery · R48 silhouette ranking · R49 rank instability · R50 PC1 dominance. |
| Paper | Paper-figure analogues from pipeline/paper_figures.py (results.json: paper): P1 species composition (Fig. 1) · P2 read mapping by location (Fig. 1) · P3 depth-equalised genes via binomial thinning (Fig. 2d–e) · P4 dropout vs expression (Fig. 2f) · P5 joint HVG-PCA (Fig. 3a–b) · P6 kBET-style mixability pre-integration (Fig. 5g–h) · P7 Harmony-corrected joint embedding (Fig. 5a–d; harmonypy, protocol = batch) · P8 mixability post-Harmony (Fig. 5g–h) · P9 Harmony wall per protocol, shared axes, species colours · P10 pre-vs-post embedding pair, protocol colours · P11 Harmony wall, protocol colours · P12 rejection pre → post per protocol · P13 Harmony wall, genes per cell. Species-mix labels stand in for cell types throughout. 72 figures in total (R1–R59 + P1–P13). |
| Datasets | GEO accession catalogue plus download centre (JSON + full-metrics/markers/composition/correlation/cumulative/QC CSVs, generated client-side) and full data tables: metrics, label agreement, composition, QC percentiles, embedded software environment. All 15 SubSeries with input/QC-pass/retention/status + paper-rank pills. |
| Protocols | Per-protocol deep-dive cards: full metric table with Δ-vs-best, auto interpretation, paper-rank badge, cluster distribution. |
| Markers | R10 dot plot (union top markers × clusters; size = fraction, shade = gene-normalised mean) with protocol selector. |
| Similarity | R11 correlation heatmap (common genes, composite-rank order) · R12 cumulative-gene curves with end labels. |
| Explorer | Canvas scatter with view switch (t-SNE / PCA / UMAP) and colour switch (cluster / genes / UMIs / label); PCA captions carry PC1/PC2 variance; per-cluster counts and hover tooltips; labels all NA without metadata. R51 t-SNE wall (all protocols, cluster colours) · R52 PCA wall (all protocols + variance shares) · R54 scanpy-Leiden UMAP wall (all protocols, shared axes) · R59 scanpy wall, species colours. |
| Glossary | 22 plain-language term definitions. |
| Discrepancies | Curated paper-vs-replication cards (survive pipeline re-runs by design). |
Every figure card carries a numbered caption and a collapsible methods/provenance note; every SVG has a Download button and <title> tooltips. Scales/axes/ranges are computed from live data (no magic constants), except the R4a axis which is fixed to [0.5, 1.0] by definition of the stability range.
9. Metric definitions and reporting standard
| Metric | Estimator, range, and reporting |
|---|---|
| Median genes / counts | Median over QC-pass cells at full depth (paper: 20K reads cell−1 per cell type). Log10 axes; protocols ordered descending. UMI and read-count chemistries estimate distinct quantities. |
| Retention | QC-pass / input cells (min_genes ≥ 200; mito gate inoperative on ENSEMBL IDs). |
| ARI / NMI / silhouette | ARI ∈ [−1,1] (≈0, chance); NMI ∈ [0,1]; silhouette ∈ [−1,1] (≈ paper ASW). Require >50% labelled cells; silhouette on ≤2,000 cells. Current labels are species-level, so ARI indexes species agreement, not cell-type accuracy. |
| Stability | Mean ± SD (n = 5 seeded 80% re-embeddings/re-clusterings) of ARI(full, subsample). Axis fixed [0.5,1.0]. Precision, not accuracy. |
| Composite | Per-metric min–max over the processed set; arithmetic mean (paper: harmonic mean of 6 at 20K reads — degenerate under few-protocol min–max). Set-relative; rescaling on protocol addition. |
| Concentration / heterogeneity / evenness | Top-50/200 count shares; CV = SD/mean across cells (R22); effective k = exp(Shannon entropy), evenness = H/log k (R23). |
| Residuals (R21) | OLS residual of log10 genes on log10 counts across protocols; + exceeds depth expectation. |
Legend standard: every figure states the display, estimator, n, transform, ordering, and error representation. Claims are scoped to the processed set and depth standardisation.
10. Paper vs replication map (Figs. 1–6)
| Paper step | Paper | This repo | Match |
|---|---|---|---|
| Primary processing | Joint-genome mapping, per-species counts | Starts from GEO count matrices | Not reproduced |
| Cell-type labels | Reference + matchSCore2 (p ≥ 0.5) | Needs external labels in metadata.csv | Not implemented |
| Downsampling | 20K reads/cell and equal cell number | Equal cell number only | Partial |
| Fig. 2 genes/dropout | By cell type, vs depth | Quality tab: log medians, bubble plot, retention, per-cell distributions, concentration | Partial |
| Fig. 3 technical effects | HVG PCA, correlation, markers, cumulative genes | Distributions, concentration + cluster markers (R10 dot plot); no HVG-PCA panel | Partial |
| Fig. 4 clustering | Cluster accuracy + ASW | Stability ± SD, community counts, size distributions, t-SNE explorer | Approximate |
| Fig. 5 integration | Separation vs mixing; projection | — | Not implemented |
| Fig. 6 overall score | Harmonic mean of 6 metric scores | Benchmark tab: dot plot + arithmetic-mean composite + rank bump chart | Approximate |
"Cluster accuracy" is undefined in the accessible text; ARI/NMI are stand-ins. Subsample stability is an added metric, not from the paper.
11. Threats to validity and known artefacts
- Construct. Full-depth medians conflate depth with efficiency (mitigated by R21 residualisation); UMI and read counts mixed; mito gate void on ENSEMBL IDs; Louvain/30-PC/res-1.0 stack differs from Seurat/8-PC/res-0.6; variance HVGs differ from paper; current labels are species-level, so label agreement ≠ cell-type accuracy.
- Internal. Single master seed (0); 5 stability repeats; min–max composites are set-relative and rescale on protocol addition.
- External. One batch, two tissues, no 20K-read standardisation — rankings are conditional. Cross-tissue checks (R18–R20) bound, but do not remove, this limit.
- Artefacts. Species doublets retained in matrix totals; TruncatedSVD RuntimeWarnings are cosmetic (outputs finite and verified); software versions are embedded per run (
results.json: environment);paper.mappingshares are null where GEO QC supplements lack read-location columns (P2 omits those protocols), and all JSON output is strictly parseable (no NaN).
11b. Generalizability: more tissues, not more species
Human + mouse is sufficient for protocol benchmarking — the ranking signal comes from protocol chemistry, not species count. For better generalizability, add more tissues / complexities / depths with the same 2 species (e.g. heterogeneous vs homogeneous tissue, high- vs low-RNA-content cells, full depth vs downsampled to 20K reads/cell), not new species.
New species add confounders without changing the ranking: annotation completeness, mito-gene prefixes (mito_prefixes in config.yaml), ortholog mapping, and dissociation bias. Add a species only to test species-mixing doublet detection or cross-species mapping.
Implemented as dashboard tab Generalizability (loads results.json + results_mouse.json jointly; 15 shared protocols): R18 median-genes human-vs-mouse scatter (log–log, Pearson r ≈ 0.84) · R19 human-vs-mouse rank bump (parallel = generalizable) · R20 stability human (dot) vs mouse (square) per protocol.
12. Reproduce it / extend it
bash run.sh setup && source .venv/bin/activate bash run.sh test # synthetic smoke test (~1 min), asserts ARI>0.5 bash pipeline/download_geo.sh list bash run.sh download GSE133540 GSE133535 GSE133539 # or all 15 accessions bash run.sh prepare GSE133540 GSE133535 GSE133539 bash run.sh inspect # check shapes, orientation, integer counts # optional: data/metadata.csv (index = barcode, column celltype) bash run.sh analyse # → site/data/results.json (+ .js fallback; preserves discrepancies) python pipeline/qc_distributions.py # → adds qc_dist (powers Figs. R6–R8) python pipeline/map_symbols.py # → adds marker symbols (needs network) python pipeline/similarity.py # → adds similarity (powers Figs. R11–R12) bash run.sh serve # mouse track: PROCESSED_DIR=data/processed_mouse SPECIES=mouse python pipeline/prepare_data.py [ACC…] # python pipeline/run_analysis.py config.mouse.yaml && python pipeline/qc_distributions.py site/data/results_mouse.json # python pipeline/map_symbols.py site/data/results_mouse.json && python pipeline/similarity.py site/data/results_mouse.json # http://localhost:8000
Resources: raw downloads ≈ 750 MB for all 15 human matrices (largest: 10X 8×10K, ~147 MB gz, 80,000 cells); full 15-protocol analysis is CPU-bound (single core) and takes on the order of an hour; site/data/results.json is a few MB (t-SNE coordinates dominate).
Adding a protocol: download its SubSeries → prepare (extend MAP if new) → analyse → qc_distributions.py → map_symbols.py → similarity.py; add its paper rank to PAPER_RANK and its CAT→protocol entry to the datasets MAP in site/index.html; all figures/scales/captions adapt automatically.
13. Results snapshot (live data)
Read live from site/data/results.json — always current, no duplicated numbers.
Metric illustrations (study data — joint HVG-PCA sample, n=1,500 cells, 24,514 common genes, this run):
Shared PC1/PC2 axes across panels (PC1 9.6%, PC2 5.9% of HVG variance). Rendered live from site/data/results.json: paper.joint — serve site/ to draw (see §12). Panel values describe the current run.