diff --git a/.gitignore b/.gitignore index d2a02cb..fc21a93 100644 --- a/.gitignore +++ b/.gitignore @@ -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 diff --git a/CHANGELOG.md b/CHANGELOG.md index 4cd739e..481ea40 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -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 diff --git a/docs/fidelity.md b/docs/fidelity.md index b7d3821..a9ddd46 100644 --- a/docs/fidelity.md +++ b/docs/fidelity.md @@ -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. diff --git a/docs/quickstart.md b/docs/quickstart.md index 7cd9871..152e0a3 100644 --- a/docs/quickstart.md +++ b/docs/quickstart.md @@ -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. diff --git a/tests/test_pbmc_handoff.py b/tests/test_pbmc_handoff.py index 7381cf2..3a7cd7c 100644 --- a/tests/test_pbmc_handoff.py +++ b/tests/test_pbmc_handoff.py @@ -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" diff --git a/tests/test_tutorial_csv_precision.py b/tests/test_tutorial_csv_precision.py index b0b6530..4aa6a19 100644 --- a/tests/test_tutorial_csv_precision.py +++ b/tests/test_tutorial_csv_precision.py @@ -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\(") @@ -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 @@ -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})" diff --git a/tutorials/README.md b/tutorials/README.md index 299de5f..5ccfcbb 100644 --- a/tutorials/README.md +++ b/tutorials/README.md @@ -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 | @@ -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 | diff --git a/tutorials/pbmc3k_tutorial.md b/tutorials/pbmc3k_tutorial.md index d10b241..5e6cc62 100644 --- a/tutorials/pbmc3k_tutorial.md +++ b/tutorials/pbmc3k_tutorial.md @@ -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 @@ -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 diff --git a/tutorials/pbmc3k_tutorial.py b/tutorials/pbmc3k_tutorial.py index c451e61..abbaa47 100644 --- a/tutorials/pbmc3k_tutorial.py +++ b/tutorials/pbmc3k_tutorial.py @@ -45,7 +45,7 @@ ) from truecell.reduction import run_pca from truecell.neighbors import find_neighbors -from truecell.clustering import find_clusters +from truecell.clustering import _res_label, find_clusters from truecell.umap import run_umap from truecell.markers import find_markers, find_all_markers @@ -77,6 +77,34 @@ "n_clusters_expected": 9, } +#: The resolutions this tutorial scans, in the order they are handed to +#: `find_clusters`. Choosing a resolution means running a few and comparing +#: them, so a fidelity claim pinned to one point says less than it appears to. +#: +#: **0.5 must stay last.** Seurat leaves the object on the last resolution in +#: the sequence, and every step after the clustering call — UMAP, markers, +#: annotation, the handoff CSV — is written against 0.5. Reordering this list +#: silently re-points all of them. +RESOLUTION_SWEEP = [0.4, 0.8, 1.2, 0.5] + +#: Clustering agreement, per resolution, as declared ranges rather than prose. +#: +#: These exist because the number they guard **drifted unnoticed**. This file +#: documented ARI 0.938 at resolution 0.5; the measured value is 0.899, and the +#: graph fixes in #67-#71 are the likely cause — they moved cells between +#: clusters, which is exactly what the DE tutorial's `deseq2 top50` band caught +#: at the time (25 -> 22). The guided tutorial had no band, so its headline +#: number went stale in six documents instead. +#: +#: Bounds are wide enough to absorb a Louvain local-optimum shift and narrow +#: enough that a real regression fails. Re-measure before widening. +CLUSTER_BANDS = { + "ARI at 0.4": (0.85, 0.95), + "ARI at 0.5": (0.85, 0.95), + "ARI at 0.8": (0.78, 0.90), + "ARI at 1.2": (0.75, 0.88), +} + def validate(label: str, value, expected=None, atol: float = 0.05) -> None: """Print a validation line. Green check if matches, red ✗ otherwise.""" @@ -244,10 +272,22 @@ def run_tutorial(data_dir: str | None = None) -> None: print(f" Graphs: {list(pbmc.graphs)}") # ----------------------------------------------------------------------- - section("10. Find Clusters (resolution=0.5)") + section("10. Find Clusters (resolutions 0.4 / 0.8 / 1.2 / 0.5)") # ----------------------------------------------------------------------- + # Seurat's `FindClusters(obj, resolution = c(...))` idiom: run several, look + # at them, then pick. Each lands in its own `RNA_snn_res.` column. + # + # 0.5 is given **last** on purpose. Seurat leaves the object on the last + # resolution in the sequence — last as given, not largest — so every step + # below this one sees exactly the partition it saw before this sweep + # existed, and the tutorial's published numbers (9 clusters, ARI 0.938, + # the marker tables) are unaffected. A resolution's partition does not + # depend on the ones before it, so 0.5 here is the same 0.5 as a lone call. t0 = time.time() - find_clusters(pbmc, resolution=0.5, algorithm=1, random_seed=0) + find_clusters(pbmc, resolution=RESOLUTION_SWEEP, algorithm=1, random_seed=0) + for res in RESOLUTION_SWEEP: + col = f"RNA_snn_res.{_res_label(res)}" + print(f" resolution {res:<4} -> {pbmc.meta_data[col].nunique()} clusters") cluster_counts = pbmc.meta_data["seurat_clusters"].value_counts().sort_index() n_clusters = len(cluster_counts) print(f" {n_clusters} clusters found ({time.time() - t0:.1f}s)") @@ -459,6 +499,14 @@ def write_anchors(pbmc, all_markers) -> None: cells[f"PC_{k + 1}"] = emb[:, k] cells.to_csv(FIGURES / "py_cell_meta.csv", index=False) + # Every resolution from the sweep, so the R side can be scored at each rather + # than only at the one the rest of this file happens to use. + res_cols = [f"RNA_snn_res.{_res_label(r)}" for r in RESOLUTION_SWEEP] + pd.DataFrame( + {"cell": list(pbmc.cell_names()), + **{c: md[c].astype(str).to_numpy() for c in res_cols}} + ).to_csv(FIGURES / "py_resolutions.csv", index=False) + rna = pbmc.assays["RNA"] # truecell stores the VST statistics on the assay's meta_data under the names # HVFInfo() uses, so these columns carry straight through to the R side of @@ -582,6 +630,40 @@ def report() -> None: print(f" mean |r| {matched.mean():.6f} min |r| {matched.min():.6f} " f"in order: {'yes' if list(rows) == list(cols) else 'NO — reordered'}") + # ---- 3b. the same comparison across the resolution sweep --------------- + # A single resolution is one point on a curve users routinely scan. Scoring + # only 0.5 leaves open whether agreement is a property of the port or of + # that setting — so every resolution in RESOLUTION_SWEEP is scored here. + py_res_path, r_res_path = (FIGURES / "py_resolutions.csv", + FIGURES / "r_resolutions.csv") + if py_res_path.exists() and r_res_path.exists(): + pr = pd.read_csv(py_res_path, dtype=str).set_index("cell") + rr = pd.read_csv(r_res_path, dtype=str).set_index("cell") + shared_cells = pr.index.intersection(rr.index) + out_of_band: list[tuple] = [] + print("\n Resolution stability — the whole sweep, not just the one " + "the rest of this file uses") + print(f" {'resolution':>11}{'truecell':>10}{'R':>5}{'ARI':>9}" + f"{'concordance':>13}") + print(f" {'-' * 48}") + for res in sorted(RESOLUTION_SWEEP): + col = f"RNA_snn_res.{_res_label(res)}" + if col not in pr.columns or col not in rr.columns: + continue + mr = match_partitions(pr.loc[shared_cells, col], + rr.loc[shared_cells, col]) + flag = " <- the one used below" if res == RESOLUTION_SWEEP[-1] else "" + band = CLUSTER_BANDS.get(f"ARI at {res}") + if band and not (band[0] <= mr["ari"] <= band[1]): + flag = f" OUT OF BAND [{band[0]}, {band[1]}]" + out_of_band.append((res, mr["ari"], band)) + print(f" {res:>11}{mr['n_a']:>10}{mr['n_b']:>5}" + f"{mr['ari']:>9.4f}{mr['concordance']:>13.4f}{flag}") + if out_of_band: + print("\n A clustering ARI outside its declared band is either a " + "regression or a\n band that was never re-measured — see " + "CLUSTER_BANDS. Do not widen it\n without saying why.") + # ---- 4. clusters: match the partitions, then score --------------------- m = match_partitions(p["cluster"], r["cluster"]) print(f"\n Clusters — truecell {m['n_a']}, R {m['n_b']}") diff --git a/tutorials/pbmc3k_verify.R b/tutorials/pbmc3k_verify.R index 789b93f..c28794b 100644 --- a/tutorials/pbmc3k_verify.R +++ b/tutorials/pbmc3k_verify.R @@ -58,7 +58,21 @@ pbmc <- RunPCA(pbmc, npcs = 50, verbose = FALSE) # that belongs to annoy rather than to either implementation. The same trap cost # `pbmc3k_objects_verify.R` a false negative of 182 SNN edges. pbmc <- FindNeighbors(pbmc, dims = 1:10, k.param = 20, nn.method = "rann", verbose = FALSE) -pbmc <- FindClusters(pbmc, resolution = 0.5, algorithm = 1, verbose = FALSE) +# The vector form, matching the Python side's RESOLUTION_SWEEP. Each resolution +# lands in its own `RNA_snn_res.` column, and Seurat leaves the object on the +# **last** one given — last as given, not largest — so 0.5 stays here at the end +# and every step below sees the partition it always saw. +pbmc <- FindClusters(pbmc, resolution = c(0.4, 0.8, 1.2, 0.5), algorithm = 1, + verbose = FALSE) +res_cols <- c("RNA_snn_res.0.4", "RNA_snn_res.0.8", + "RNA_snn_res.1.2", "RNA_snn_res.0.5") +write.csv( + data.frame(cell = colnames(pbmc), + setNames(lapply(res_cols, function(c) as.character(pbmc[[c]][, 1])), + res_cols), + check.names = FALSE, stringsAsFactors = FALSE), + file.path(FIG, "r_resolutions.csv"), row.names = FALSE) +cat("Wrote r_resolutions.csv (", paste(res_cols, collapse = ", "), ")\n") pbmc <- RunUMAP(pbmc, dims = 1:10, verbose = FALSE) cat(sprintf("PBMC 3k: %d cells -> %d after QC, %d clusters\n", n_cells_raw, ncol(pbmc), length(levels(pbmc))))