Numeric handoffs for PBMC 3k and PBMC 8k, and fix find_all_markers - #61
Merged
Conversation
…nd_all_markers Tutorials 1 and 2 were the last two in the suite compared entirely by eye: every R panel in the PBMC 3k walkthrough links a canonical satijalab.org image, so nothing failed if the numbers behind those pictures drifted. Both verify scripts now write a numeric handoff and both tutorials gained a `--report` that compares it. Neither side is pinned to the other — each runs its own pipeline from the same 10x bytes, which is what makes the comparison worth anything. Per-cell rows are keyed by barcode and per-gene rows by symbol so cluster *numbering*, which is arbitrary in both tools, can be matched one-to-one on overlap before anything is scored. PBMC 3k vs Seurat 5.5.1: the same 2,638 barcodes survive QC (nCount/nFeature exact, percent.mt 5.3e-15), VST means to 4.8e-14, 1,998/2,000 variable features shared, PCA over the 10 clustering dims at |r| 0.9988, kNN 52,760 on both, clusters at ARI 0.938 (8 vs 9). PBMC 8k: the same 7,475 barcodes, global clusters at ARI 0.977, and the T/NK compartment handed to stage two matches at Jaccard 0.9991 — 4,631 of 4,635 cells are the same barcodes. That is the load-bearing number: everything in stage two is conditioned on it, so a compartment of the right size drawn from the wrong clusters would pass a count check and make every later number incomparable. Found one real defect. `find_all_markers` was missing Seurat's `return.thresh = 1e-2` and returned rows Seurat does not. The handoff pinned it decisively: two PBMC 3k clusters came out with identical cell membership on both sides, and on those two shanuz returned 190 and 383 genes against Seurat's 151 and 242 — every extra row a gene the two agreed about to 4.9e-15 and that Seurat simply does not return. Applying the filter reproduces Seurat's gene sets exactly, 151/151 and 242/242. Rows are now also ordered by p_val then descending avg_log2FC, matching `order(gde$p_val, -gde[, 2])`; without the tie-break "the top N markers" was decided by incoming row order, because the strongest markers tie at p = 0 (40-302 genes per cluster here). Cluster labels iterate numerically, so eleven or more clusters no longer come back as 0, 1, 10, 11, 2. Two stale claims corrected: tutorials/README recorded "Clusters at resolution 0.5: 9" and "DC recovered" for PBMC 3k, and the vignette said "9 clusters in both R and Python". shanuz finds 8; the missing cluster is a 32-cell dendritic population whose cells land, all 32, in its CD14+ Mono cluster. PBMC 8k runs the opposite way — shanuz splits Seurat's merged Platelet/DC cluster — so neither run is uniformly finer, and the docs now say so. Both R references pin nn.method = "rann", the same trap that cost pbmc3k_objects_verify.R a false negative of 182 SNN edges. On PBMC 8k it is worse than cosmetic: the graph decides the global clusters, which decide which cells enter the subclustering stage. 18 new tests (806 -> 824), each mutation-tested. One caught a decorative guard: the tie-break test survived removing the tie-break because its fixture produced no ties, and had to be rebuilt around six perfectly-separating markers. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Closes the last gap the fresh-install audit left open: tutorials 1 and 2
were the only two in the suite still compared entirely by eye. Every R panel in
the PBMC 3k walkthrough links a canonical satijalab.org image, so nothing
failed if the numbers behind those pictures drifted.
Both verify scripts now write a numeric handoff, and both tutorials gained a
--report. Neither side is pinned to the other — each runs its ownpipeline from the same 10x bytes, which is what makes the comparison worth
anything. (
pbmc3k_dimreduc_verify.Ris the opposite experiment: it does pinPython's cells and features, to isolate the post-PCA machinery.)
Per-cell rows are keyed by barcode and per-gene rows by symbol, and cluster
comparisons match the two partitions one-to-one on overlap before scoring —
cluster numbering is arbitrary in both tools, so without that a pure
relabelling reads as total disagreement.
PBMC 3k, end to end
variance.standardized2.6e-2 relativeavg_log2FCto 4.9e-15PBMC 8k, both stages
The compartment Jaccard is the load-bearing number, and the reason the handoff
dumps barcodes rather than counts. Everything in stage two is conditioned on
which cells stage one selected, so a compartment of the right size drawn from
the wrong clusters would still produce subclusters, still produce a marker
table, and pass a count check while making every later number incomparable.
The defect this found
find_all_markerswas missing Seurat'sreturn.thresh = 1e-2and returnedrows Seurat does not. The handoff pinned it the decisive way: two PBMC 3k
clusters came out with identical cell membership on both sides, and on those
two shanuz returned 190 and 383 genes against Seurat's 151 and 242 — every
extra row a gene the two agreed about numerically (max
avg_log2FCdifference4.9e-15) and that Seurat simply does not return. Applying the same filter
reproduces Seurat's gene sets exactly: 151/151 and 242/242.
Two more fixes in the same function:
p_valascending thenavg_log2FCdescending,matching
order(gde$p_val, -gde[, 2]). Not cosmetic — Wilcoxon p-values tieat exactly 0 for the strongest markers (40–302 genes per cluster here), so
"the top N markers", which every tutorial prints, was decided by incoming row
order.
back as 0, 1, 10, 11, 2, … Invisible below eleven, which is why every tutorial
in the suite missed it.
return_thresh=Nonerestores the old unfiltered table.Stale docs corrected
tutorials/README.mdrecorded "Clusters at resolution 0.5: 9 ✅" and "All 6canonical cell types recovered … DC ✅"; the vignette said "9 clusters in
both R and Python". shanuz finds 8, and the cluster Seurat has that it does
not is a 32-cell dendritic population whose cells land — all 32 — in
shanuz's CD14+ Mono cluster. The provenance is now traced in the vignette from
the LOESS fit through the 2 swapped variable features and the 286 moved SNN
edges to that one boundary.
PBMC 8k runs the opposite way: Seurat's 100-cell cluster 11 holds both
Platelet and DC and shanuz splits it in two. Same borderline population, other
side of the threshold — so neither run is uniformly finer, and the docs now say
so rather than implying a verdict.
Also
nn.method = "rann"— the annoy trap that costpbmc3k_objects_verify.Ra false negative of 182 SNN edges. On PBMC 8k it isworse than cosmetic: the graph decides the global clusters, which decide which
cells enter the subclustering stage.
their cells, instead of letting the Hungarian matching drop them off the table.
_assign_cell_typesfixed to match the R port.Verification
the tie-break test survived removing the tie-break, because its fixture
produced no ties. Rebuilt around six perfectly-separating markers (which share
a U statistic and so a p-value) with an anti-vacuity assertion; mutation now
caught.
ruff check shanuz72 and whole-repo 222, both unchanged frommainmeasured with the branch stashed.
🤖 Generated with Claude Code