Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -121,6 +121,8 @@ tutorials/figures/py_hvg.csv
tutorials/figures/r_hvg.csv
tutorials/figures/py_markers.csv
tutorials/figures/r_markers.csv
tutorials/figures/py_resolutions.csv
tutorials/figures/r_resolutions.csv
tutorials/figures/py_anchors.json
tutorials/figures/r_anchors.json
# PBMC 8k subclustering tutorial — global + T/NK compartment handoff
Expand Down
43 changes: 43 additions & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -18,8 +18,51 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0

## [Unreleased]

### Fixed

- **The PBMC 3k guided tutorial documented a clustering ARI of 0.938 while
measuring 0.899**, across six files (`tutorials/README.md` ×2,
`pbmc3k_tutorial.md` ×2, `docs/fidelity.md`, `docs/quickstart.md`). The
associated concordance figures were stale too — 2,554/2,638 cells and 0.968,
against a measured 2,519 and 0.955.

The likely cause is the graph fixes in #67–#71, which moved cells between
clusters — the same drift the DE tutorial's `deseq2 top50` band caught at the
time (25 → 22). This tutorial had **no band on its headline number**, so it
went stale in six documents instead of failing once. Every swept resolution
now carries a declared band, checked by `--report`.

Verified this predates the change: R's vector-form `FindClusters` produces a
partition identical to the single 0.5 call, and the tutorial's four handoff
outputs are byte-identical to what `main` produces.

### Changed

- **The guided clustering tutorial scans four resolutions instead of one.**
Choosing a resolution means running a few and comparing them, so a fidelity
claim pinned to a single setting says less than it appears to. Both sides now
use the vector idiom — `find_clusters(pbmc, resolution=[0.4, 0.8, 1.2, 0.5])`
and `FindClusters(pbmc, resolution = c(...))` — and every resolution is scored
against R:

| resolution | truecell | Seurat | ARI |
|---:|---:|---:|---:|
| 0.4 | 9 | 9 | 0.8958 |
| 0.5 | 8 | 9 | 0.8987 |
| 0.8 | 11 | 11 | 0.8264 |
| 1.2 | 12 | 12 | 0.7995 |

Two readings. The **cluster count matches exactly at 0.4, 0.8 and 1.2**, so
the 8-vs-9 split this tutorial describes is specific to resolution 0.5 rather
than a standing property of the port. And agreement falls as resolution rises,
which is expected: finer partitions put more boundaries in play.

0.5 is given **last** deliberately, since Seurat leaves the object on the last
resolution in the sequence and every step below the clustering call is written
against it. The four handoff outputs are byte-identical to before, and a test
pins the ordering, because reordering the list would silently re-point UMAP,
the markers, the annotation and the handoff without raising.

- **`add_module_score` is now verified against Seurat as an equality, not a
correlation.** Every prior comparison was correlation-based and structurally
had to be: `AddModuleScore` picks control genes with `sample()`, R's RNG is
Expand Down
3 changes: 2 additions & 1 deletion docs/fidelity.md
Original file line number Diff line number Diff line change
Expand Up @@ -78,7 +78,8 @@ Real, understood, and not going away:

**Louvain cluster counts drift by one.** Both tools run the same algorithm at the
same resolution and land on different local optima. On PBMC 3k, truecell finds 8
clusters to Seurat's 9 at ARI 0.938. On ifnb RPCA, Seurat's deeper modularity
clusters to Seurat's 9 at ARI 0.899 — and the count matches exactly at
resolutions 0.4, 0.8 and 1.2, so the 8-vs-9 split is specific to 0.5. On ifnb RPCA, Seurat's deeper modularity
search buys 0.17 % modularity by splitting CD14 Mono along the batch — and
truecell's coarser partition then scores **ARI 0.92 against the annotations to
Seurat's 0.74**. The coarser answer is the better one there.
Expand Down
2 changes: 1 addition & 1 deletion docs/quickstart.md
Original file line number Diff line number Diff line change
Expand Up @@ -122,7 +122,7 @@ variances agree to three decimals; see [Fidelity](fidelity.md#what-actually-diff
clusters: [692, 515, 458, 344, 301, 159, 155, 14]
```

Eight clusters, against Seurat's nine on the same data, at **ARI 0.938** — the
Eight clusters, against Seurat's nine on the same data, at **ARI 0.899** — the
extra one is a 32-cell dendritic-cell population Seurat's deeper modularity
search separates.

Expand Down
61 changes: 61 additions & 0 deletions tests/test_pbmc_handoff.py
Original file line number Diff line number Diff line change
Expand Up @@ -216,3 +216,64 @@ def test_numeric_anchor_names_agree_across_the_two_sides():
f"{module} and {script} disagree about the anchor set: "
f"only python {sorted(py_keys - r_keys)}, only R {sorted(r_keys - py_keys)}"
)


# ---------------------------------------------------------------------------
# The resolution sweep
# ---------------------------------------------------------------------------

def test_the_sweep_ends_on_the_resolution_everything_downstream_uses():
"""0.5 must be last in RESOLUTION_SWEEP.

`find_clusters` leaves the object on the **last** resolution given, and UMAP,
the markers, the cell-type annotation and `py_cell_meta.csv` are all written
against 0.5. Reordering the list silently re-points every one of them, and
nothing downstream would raise — the tutorial would simply start describing a
different clustering while still calling it resolution 0.5.
"""
from tutorials.pbmc3k_tutorial import RESOLUTION_SWEEP

assert RESOLUTION_SWEEP[-1] == 0.5


def test_both_languages_sweep_the_same_resolutions():
"""The list lives in two files; a drift between them is invisible.

If R swept {0.4, 0.7, 0.5} and Python {0.4, 0.8, 0.5}, the comparison would
simply skip the columns that did not line up and report the ones that did —
fewer rows, no error, no sign that anything was missed.
"""
import re

from tutorials.pbmc3k_tutorial import RESOLUTION_SWEEP

r_text = (TUTORIALS / "pbmc3k_verify.R").read_text()
call = re.search(r"FindClusters\(pbmc,\s*resolution\s*=\s*c\(([^)]*)\)", r_text)
assert call, "pbmc3k_verify.R no longer calls FindClusters with a vector"
r_values = [float(v) for v in call.group(1).split(",")]

assert r_values == RESOLUTION_SWEEP, (
f"R sweeps {r_values}, Python sweeps {RESOLUTION_SWEEP}; the comparison "
"would silently score only the overlap")


def test_every_swept_resolution_has_a_declared_band():
"""A resolution scored with no band is one that can drift unnoticed.

Which is not hypothetical here: this tutorial's ARI at 0.5 was documented as
0.938 across six files while measuring 0.899, because there was no band on it.
"""
from tutorials.pbmc3k_tutorial import CLUSTER_BANDS, RESOLUTION_SWEEP

for res in RESOLUTION_SWEEP:
assert f"ARI at {res}" in CLUSTER_BANDS, f"resolution {res} has no band"


def test_the_bands_are_not_vacuous():
"""A band spanning [0, 1] would pass anything and prove nothing."""
from tutorials.pbmc3k_tutorial import CLUSTER_BANDS

for name, (low, high) in CLUSTER_BANDS.items():
assert high > low, f"{name}: empty band"
assert high - low <= 0.20, f"{name}: band too wide to catch a regression"
assert 0.0 <= low and high <= 1.0, f"{name}: ARI is bounded by 1"
17 changes: 14 additions & 3 deletions tests/test_tutorial_csv_precision.py
Original file line number Diff line number Diff line change
Expand Up @@ -26,12 +26,17 @@

TUTORIALS = Path(__file__).parent.parent / "tutorials"

# Reads that carry no float into a comparison. Each is exempt for a stated
# reason; adding to this set should require one.
# `dtype=str` parses no floats at all, so `float_precision` cannot apply — the
# rule is general rather than a per-file carve-out. It covers the DE tutorial's
# hex reader (floats arrive via `float.fromhex`) and any read of purely
# categorical columns, such as the resolution sweep's cluster labels.
GENERALLY_EXEMPT = ("dtype=str",)

# Reads that carry no float into a comparison for a file-specific reason. Each is
# exempt for a stated reason; adding to this set should require one.
EXEMPT = {
("anchors_tutorial.py", "metadata.csv"): "cell labels, no floats compared",
("visium_tutorial.py", "header=header"): "wrapped in len() — a row count",
("pbmc3k_de_tutorial.py", "dtype=str"): "the hex reader; floats via float.fromhex",
}

READ_CSV = re.compile(r"read_csv\(")
Expand Down Expand Up @@ -67,6 +72,8 @@ def _offending_calls(path: Path) -> list[str]:
call = _call_text(source, m.end() - 1)
if "float_precision" in call:
continue
if any(marker in call for marker in GENERALLY_EXEMPT):
continue
if any(name == path.name and marker in call
for (name, marker) in EXEMPT):
continue
Expand All @@ -91,6 +98,10 @@ def test_the_exemptions_still_exist():
If one of the exempt calls is edited away, the entry stops matching anything
and quietly permits a real offender in the same file.
"""
for marker in GENERALLY_EXEMPT:
assert any(marker in p.read_text() for p in TUTORIALS.glob("*.py")), (
f"general exemption {marker!r} no longer matches anything — remove it")

for (name, marker), reason in EXEMPT.items():
path = TUTORIALS / name
assert path.exists(), f"{name} is gone; drop its exemption ({reason})"
Expand Down
5 changes: 3 additions & 2 deletions tutorials/README.md
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@ each pairing **R Seurat** code side-by-side with the equivalent **Python Truecel

| # | Tutorial | Dataset | Key Concepts | Complexity |
|---|----------|---------|--------------|-----------|
| 1 | [PBMC 3k — Guided Clustering](pbmc3k_tutorial.md) | 3,000 PBMCs · 10x Genomics (2016) | QC · Normalization · HVG/VST · PCA · Louvain · UMAP · Markers. **Compared end to end**, both sides running their own pipeline: the same 2,638 barcodes survive QC, 1,998/2,000 variable features shared, PCA matched \|r\| **0.9988**, clusters at **ARI 0.938** (8 vs 9 — one 32-cell DC population), and on the two clusters whose cells match exactly the marker tables are **identical gene sets** agreeing to 4.6e-14 | Beginner |
| 1 | [PBMC 3k — Guided Clustering](pbmc3k_tutorial.md) | 3,000 PBMCs · 10x Genomics (2016) | QC · Normalization · HVG/VST · PCA · Louvain · UMAP · Markers. **Compared end to end**, both sides running their own pipeline: the same 2,638 barcodes survive QC, 1,998/2,000 variable features shared, PCA matched \|r\| **0.9988**, clusters at **ARI 0.899** (8 vs 9 — one 32-cell DC population), and on the two clusters whose cells match exactly the marker tables are **identical gene sets** agreeing to 4.6e-14 | Beginner |
| 2 | [PBMC 8k — Advanced Subclustering](advanced_pbmc8k_subclustering.md) | 8,400 PBMCs · GRCh38 · 10x Genomics | All of Tutorial 1 + subclustering, hierarchical cell-type gating, T/NK annotation. **Both stages compared by barcode**: same 7,475 cells after QC, global clusters at **ARI 0.977**, and the T/NK compartment handed to stage 2 matches at **Jaccard 0.9991** (4,631 of 4,635 cells) — subclusters then at ARI 0.916 and subset labels **98.2%** concordant | Intermediate |
| 3 | [CBMC CITE-seq — Multimodal](multimodal_citeseq.md) | 8,600 CBMCs · RNA + 13 surface proteins | Multi-assay objects · CLR normalization · Protein feature plots · RNA-protein comparison · WNN joint clustering. **Compared per protein and per cell**: CLR to **4.2e-15**, WNN modality weights at Pearson **0.9847** over 8,617 shared barcodes, cell-type labels **99.29%** concordant, all cluster counts identical. Settled the long-open progenitor question — it was a labelling difference, not a WNN one | Advanced |
| 4 | [PBMC 3k — SCTransform](sctransform_vignette.md) | 3,000 PBMCs · 10x Genomics (2016) | Regularized NB normalization · Pearson residuals · `vars.to.regress` · 30-PC workflow · SCT-vs-LogNormalize. The **fitted model is compared per gene** against Seurat's `SCTModel` feature attributes — `detection_rate`/`gmean` to machine precision, intercept and theta at **Spearman 1.0000**, the 3,848 non-overdispersed genes exactly the same set, residual variance at 0.9986; both arms now agree on cluster count (12 and 11) | Advanced |
Expand Down Expand Up @@ -131,7 +131,8 @@ pipeline from the same 10x bytes; nothing is pinned across them.
| The 2,000 variable features | **1,998 shared**; the two swaps are genes 0.03 apart in a LOESS fit, at ranks 1,982–2,000 |
| PCA, the 10 dims clustering uses | matched \|r\| mean **0.9988**, min 0.9946, no reordering |
| kNN graph | **52,760 = 2,638 × 20 on both** |
| Clusters at resolution 0.5 | truecell **8**, Seurat **9** — ARI **0.938**, 2,554/2,638 cells agree |
| Clusters at resolution 0.5 | truecell **8**, Seurat **9** — ARI **0.899**, 2,519/2,638 cells agree |
| Clusters across the sweep | 0.4 → 9 vs 9 (ARI 0.896) · 0.8 → 11 vs 11 (0.826) · 1.2 → 12 vs 12 (0.800). The cluster **count** matches at every resolution except 0.5 |
| The one cluster Seurat has and truecell does not | its 32 DC cells land, **all 32**, in truecell's CD14+ Mono cluster |
| Markers on the two clusters whose cells match exactly | **identical gene sets** (151/151 and 242/242), `avg_log2FC` to 4.9e-15 and 4.6e-14 respectively |

Expand Down
37 changes: 35 additions & 2 deletions tutorials/pbmc3k_tutorial.md
Original file line number Diff line number Diff line change
Expand Up @@ -544,11 +544,44 @@ find_clusters(

> Both use Louvain community detection via `igraph`. At `resolution = 0.5`
> Seurat returns **9** clusters here and truecell **8** — the two runs agree
> about 2,554 of the 2,638 cells (**ARI 0.938**), and the single cluster
> about 2,519 of the 2,638 cells (**ARI 0.899**), and the single cluster
> Seurat has and truecell does not is a 32-cell dendritic-cell population whose
> cells land, **all 32 of them**, in truecell's CD14+ Mono cluster. DCs are
> monocyte-lineage, so this is one borderline split at this resolution, not a
> scattered disagreement.

### The whole sweep, not one point on it

Picking a resolution means running a few and comparing them, so a fidelity claim
pinned to a single setting says less than it looks like it does. The tutorial
runs `find_clusters(pbmc, resolution=[0.4, 0.8, 1.2, 0.5])` — Seurat's own
vector idiom — and scores every one against R:

| resolution | truecell | Seurat | ARI | concordance |
|---:|---:|---:|---:|---:|
| 0.4 | 9 | 9 | 0.8958 | 0.9602 |
| **0.5** | **8** | **9** | **0.8987** | **0.9549** |
| 0.8 | 11 | 11 | 0.8264 | 0.9174 |
| 1.2 | 12 | 12 | 0.7995 | 0.8647 |

Two things worth reading off this. **The cluster count matches exactly at 0.4,
0.8 and 1.2** — the 8-vs-9 split described above is specific to resolution 0.5,
not a standing property of the port. And **agreement falls as resolution rises**
(0.90 → 0.83 → 0.80), which is what you would expect: finer partitions put more
boundaries in play, and each is another chance for the two Louvain runs to land
in different local optima.

0.5 is given **last** on purpose. Seurat leaves the object on the last
resolution in the sequence, so every step below — UMAP, markers, annotation, the
handoff — sees exactly the partition it saw before this sweep existed. The four
outputs this tutorial writes are byte-identical to what it produced without it.

Each ARI now carries a **declared band** (`CLUSTER_BANDS`), for a reason: this
file documented ARI **0.938** at resolution 0.5 while measuring 0.899, across six
documents. The number had drifted — most likely with the graph fixes in #67–#71,
which moved cells between clusters — and with no band on it, nothing failed.
That is the same drift the DE tutorial's `deseq2 top50` band *did* catch at the
time (25 → 22).
>
> Where it comes from is traceable: the two runs keep the **same 2,638
> barcodes** and agree on the per-gene VST means to 4.8e-14 and observed
Expand Down Expand Up @@ -1037,7 +1070,7 @@ on this dataset:
| Variable features | **1,998 of 2,000 shared**, rank Spearman 0.9999 |
| PCA (10 dims) | matched \|r\| mean **0.9988**, min 0.9946, no reordering |
| kNN graph | **52,760 on both** (2,638 × 20) |
| Clusters | 8 vs 9 — **ARI 0.938**, concordance 0.968 (see Step 11) |
| Clusters | 8 vs 9 — **ARI 0.899**, concordance 0.955 (see Step 11) |
| Markers, where the clusters hold identical cells | **identical gene sets** (151/151, 242/242), `avg_log2FC` to 4.9e-15 and 4.6e-14 respectively |

That last row is the one to read carefully. Two clusters — B cells and
Expand Down
Loading
Loading