Add GC reclustering presets and Leiden clone embedding review figures

.gitignore

@@ -13,6 +13,7 @@ data/
 *.h5
 manuscript.ipynb
 .coverage
+internal_validation/
 
 # lineage-tree working files: regenerate with case_studies/lineage_trees.py + run_gctree.sh
 case_studies/results/*/lineage/trees/

CHANGELOG.md

@@ -0,0 +1,56 @@
+# Changelog
+
+## 4.1.0 — 2026-10-09
+
+### New optional analyses
+
+- StateDensityModel provides fixed, clone-balanced reference regions, context-aware
+  leave-own-clone-out backgrounds, count-based shrinkage, conditional posterior
+  intervals, separate simultaneous sampling intervals, regional contrasts and
+  clone-balanced gene/module descriptions.
+- Query counts can be accumulated in batches. Numeric NPZ model serialization
+  does not require pickle.
+- FrozenExpressionModel fits receptor-excluded HVGs, scaling and PCA on reference
+  cells only, then projects held-out cells without refitting.
+
+### Correctness and reproducibility
+
+- Kernel profile reconstruction preserves the random-feature basis after bandwidth
+  estimation. Legacy models are matched to saved features or rejected; modified
+  cells, labels or coordinates invalidate the cached fit. Older AnnData files that
+  omit None-valued context parameters are supported.
+- The minimum capture count now solves the actual weighted reliability function.
+- AIRR heavy/light masks remain aligned after UMI sorting. The NumPy affine-gap
+  global alignment fallback penalizes leading gaps.
+- Memory null matches never borrow another donor. Insufficient matched pairs and
+  unidentified noise-corrected ratios are reported explicitly.
+- Degenerate Gaussian split nulls do not yield false p=0 results. Split tests remain
+  approximate because they do not repeat the full model-selection procedure.
+- Additional checks reject ambiguous IDs, missing strata, invalid weights/counts,
+  and undefined zero-variation or empty-data inference.
+
+### Migration and scientific interpretation
+
+The existing mean/kernel pipeline remains the default. Recompute downstream
+programme stability, assignments and memory when upgrading a 4.0 analysis:
+those routines may have reconstructed a different random-feature basis.
+The initial fitted profile formulas are unchanged. No automatic conversion of
+an unreproducible legacy fit is attempted.
+
+The optional state model is not a validated replacement for every endpoint.
+Simulations retain prior-misspecification failures; whole-mouse reporter tests
+often favour RNA means. Its data_weight is not the legacy profile reliability,
+and posterior direction probabilities are not p-values. See
+[the algorithm review](docs/ALGORITHM_REVIEW_2026-10-09_zh.md) and
+[state-distribution documentation](docs/STATE_DENSITY.md).
+
+### Paper and reproducibility assets
+
+Twelve case studies were rerun and checked for identical cells, clone membership,
+capture counts and display eligibility. Publication display coordinates are
+retained; corrected statistics and programme annotations are regenerated.
+The LARRY future-state comparison now includes RNA mean plus variance.
+Figure 4B is the original cell SPRING layout; Figure 4C is Threadfin, and the
+map comparison alone does not establish method superiority. Figure 2's reporter
+and SHM interpretation is supported by conditional and within-family checks,
+with capture and physical-library limitations stated.

README.md

@@ -7,36 +7,34 @@
 [![Python >= 3.10](https://img.shields.io/badge/python-%3E%3D3.10-blue.svg)](https://www.python.org)
 
 A B-cell clone is a family of cells descended from one ancestor that share a
-B-cell receptor. Paired single-cell sequencing tells you, for every cell, which
-clone it belongs to and what it is doing. Threadfin uses the **clone as the
+B-cell receptor. Paired single-cell sequencing supports receptor-based family
+assignment and measures each captured cell's expression state. Threadfin uses the **clone as the
 unit of analysis**.
 
-This matters most where the biology is a cycle. In a germinal centre, B cells
-divide and mutate their receptor, test it, and are then either sent back for
-another round or leave as plasma or memory cells. Nothing in that cycle is a
-beginning or an end, so ordering single cells along a pseudotime asks a
-question the biology does not answer. A clone is different: all of its cells
-descend from one ancestor, so a clone has a history even when its cells do not
-have an order. Threadfin therefore compares clones with one another, and asks
-which clones look as though they were recently selected and which look as
-though they were sent back to divide again.
+In a germinal centre, B cells alternate between selection and division,
+while clonally related cells can occupy different states. Threadfin compares
+the distributions of states captured within sequence-defined families. Measured
+division reporters, sorting gates and antigen probes help interpret those
+distributions. A clone map shows similarities between families' observed
+expression profiles; neighbouring points do not establish shared ancestry,
+future fate or memory-cell re-entry into a GC.
 
 It answers four questions:
 
-1. **Is B-cell state inherited within clones?** How much of the variation in
+1. **Which states recur within clones?** How much of the variation in
    gene expression is explained by which clone a cell belongs to, beyond when
    and where the cells were sampled.
-2. **Which clones behave alike?** Clones are grouped only when the split
-   between groups is real; otherwise they are reported as a continuum.
+2. **Which clones behave alike?** Reliable clones are grouped when matched-null
+   tests and resampling support separation; otherwise they are reported as a continuum.
 3. **What explains the differences between clones?** Clone-level tests against
    anything you measured: antigen binding, isotype, mutation load, division
    history, sort gate, infection, tissue, time.
 4. **Do clones keep their state** over time, across tissues, or across the
    compartments of a germinal centre?
 
-Every result comes with an honest measure of uncertainty, every test counts
-clones rather than cells, and the output is a hypothesis about clones that
-experiments still have to confirm.
+Uncertainty and profile reliability accompany the analysis. Clones supply
+profile units; donor- or mouse-level designs are needed for population and
+treatment claims. The output supports testable biological hypotheses.
 
 ---
 
@@ -49,6 +47,29 @@ pip install "git+https://github.com/Reynold4k/Threadfin.git"
 Python >= 3.10. Scanpy, Leiden and Harmony are installed automatically; no GPU
 and no R are needed.
 
+Clean Linux wheel installations were checked on Python 3.10, 3.11 and 3.12,
+including the full test suite, tutorial and simulated quick start. The legacy
+v3 reclustering repeatability check remains an expected failure on Python 3.12.
+See the [compatibility results](docs/COMPATIBILITY.md) and
+[English input-error guide](docs/TROUBLESHOOTING.md).
+
+## GC manuscript and figures
+
+The [GC manuscript](paper/MANUSCRIPT_draft_v2.md),
+[Figure 1](paper/figure_plan/Figure_1.png),
+[figure legends](paper/FIGURE_LEGENDS.md) and
+[native-method comparison](paper/METHOD_COMPARISON.md) distinguish receptor
+family membership from captured cell-state similarity. Figures 1–5 and
+Supplementary Figures 1–6 and 8 are available. The NP native comparison is complete;
+the RBD BiGCN run and Figure 6/S7 remain pending. The manuscript makes this
+status explicit and reports the RNA-centroid baseline's stronger NP readout.
+
+The [manuscript export guide](paper/EXPORT.md) describes template-style DOCX
+and numbered-reference/RIS exports. The
+[literature review and reviewer questions](paper/LITERATURE_REVIEW_REVIEWER_QUESTIONS_2026-10-06.md)
+record scientific improvements for future work. No additional biological
+analyses are required to build those documents.
+
 ## Quick start
 
 ```python
@@ -122,7 +143,7 @@ print(result.summary())
    split between two states is not mistaken for one sitting in between.
 3. **Coherence test.** The share of transcriptional variance explained by
    clone identity is compared with clones shuffled *within samples*, so
-   differences between samples are never mistaken for clonality.
+   the baseline accounts for differences between sampling contexts.
 4. **Programmes.** Reliable clones are grouped by similarity of their
    profiles. Every split between groups must pass a significance test against
    a single-group model, so clones that only vary along a continuum are
@@ -158,7 +179,7 @@ tf.tl.profile_association(adata, "isotype")               # 5. label vs clone st
 tf.tl.association_test(adata, "isotype")                  #    label vs programmes
 tf.tl.programme_markers(adata)                            #    genes, clones as replicates
 tf.tl.clonal_memory(adata, "timepoint")                   #    memory over time
-tf.tl.gene_heritability(adata)                            #    which genes are clonally inherited
+tf.tl.gene_heritability(adata)                            #    within-clone gene-expression resemblance
 ```
 
 Figures: `tf.pl.clone_map`, `tf.pl.coherence`, `tf.pl.stability`,
@@ -167,43 +188,28 @@ Figures: `tf.pl.clone_map`, `tf.pl.coherence`, `tf.pl.stability`,
 
 ## What Threadfin found in public data
 
-Eight published datasets, human and mouse, were re-analysed with the same
-script (`case_studies/`). A few of the findings:
-
-**Germinal centres have no beginning and no end**, so ordering single cells
-along a pseudotime asks a question the biology does not answer. Comparing
-*clones* works instead, because a clone has a history even when its cells do
-not have an order. In two mouse experiments with a model antigen, where cells
-were sorted by how often they had divided and by germinal-centre zone:
-
-* clone identity explained 9-15% of B-cell state, far more than shuffled
-  clones (p = 0.002), but never most of it: a clone biases what its cells do;
-* clones did **not** fall into distinct programmes - they differ along a
-  continuum;
-* how many times a clone's cells had divided was the strongest explanation of
-  how clones differ (18% and 8%), with the dark/light-zone sort comparable
-  (9%), while antigen binding and mutation load explained much less;
-* about half of what distinguishes a clone was still recognisable when its
-  cells were caught in a different part of the cycle.
-
-**In a vaccinated human lymph node**, spike-binding clones were concentrated in
-the germinal-centre and one antibody-secreting group of clones (odds ratios
-around 3 across donors). Clones kept their state over months (memory index
-0.26), but the same clone's cells in blood and lymph node did not resemble each
-other at all - where a cell is matters more than which clone it came from.
-
-**Seven days after influenza vaccination**, the antibody-secreting burst came
-from class-switched, mutated clones, as expected if it is recall of existing
-memory rather than a new response.
-
-These are associations in observational data: Threadfin produces hypotheses
-about clones that need experiments to confirm. The step-by-step walkthroughs
-are in [`report/PUBLIC_DATASETS_REPORT.md`](report/PUBLIC_DATASETS_REPORT.md)
-and, for the germinal centre,
-[`report/GERMINAL_CENTRE_CASE_STUDIES.md`](report/GERMINAL_CENTRE_CASE_STUDIES.md);
-the analysis code is in [`case_studies/`](case_studies/); how Threadfin relates
-to other tools is in [`docs/GAP_ANALYSIS.md`](docs/GAP_ANALYSIS.md); the
-manuscript figure plan is in [`paper/figure_plan/`](paper/figure_plan/).
+The [current manuscript](paper/MANUSCRIPT_draft_v2.md) centres on GC biology:
+
+* **Measured GC division history:** NP-OVA and RBD reporter cohorts show
+  clone-associated state variation, interpreted using division/FACS labels
+  and independently measured BCR mutation load.
+* **Plasmodium GC/output relationships:** some same-mouse families contain
+  GC, plasmablast or memory-like cells, including exact paired heavy/light
+  matches. GC–PB sharing shows no consistent extra enrichment after controlling
+  isotype composition. It supports candidate selection rather than directed
+  fate or memory re-entry claims.
+* **Longitudinal human vaccination:** repeated same-donor GC-containing families
+  test persistence across dates; an expression map alone cannot identify a
+  parent cell.
+* **One non-GC validation:** verified marrow/blood donor identities and pure
+  PC/memory sorts reveal exact heavy/light receptor sharing across gates.
+
+[Main figures and focused supplements](paper/figure_plan/FIGURE_PLAN.md)
+carry these questions. Other [datasets already tested](paper/figure_plan/tested_datasets/)
+are retained as exploratory coverage. The [biological review](paper/BIOLOGICAL_INFERENCE_REVIEW_zh.md)
+explains the evidence, controls and limitations. Earlier reports in `report/`
+record the original analysis and are superseded by the current manuscript
+where sampling identities or interpretations have changed.
 
 ## Upgrading from v3
 
@@ -221,3 +227,38 @@ A manuscript is in preparation.
 ## License
 
 [MIT](LICENSE) © 2026 Chen Satoshi (Reynold4k)
+
+## Centroid reclustering presets
+
+The legacy `tf.clonotype_recluster` path now accepts `preset="cohesive"`
+(the unchanged default settings), `"continuous"` and `"discrete"`, with explicit
+UMAP/graph overrides and saved effective parameters. See
+[reclustering controls](docs/RECLUSTERING.md). These starting configurations
+remain separate from v4 profile reliability and programme validation.
+
+[Figure 2B](paper/figure_plan/Figure_2.png) now uses audited real GSE246382
+clone embedding of 49 same-mouse families, coloured by unsupervised Leiden
+`clone_cluster` with a captured-size legend. The
+[standalone preview](paper/figure_plan/review/GSE246382_clone_embedding.png)
+uses the notebook's distance-row UMAP and Scanpy-graph recipe, available through
+`embedding_mode="distance_profiles", cluster_on="embedding"`. Measured
+GC/PC gates and Myc RNA annotate the resulting three clusters.
+[Supplementary 8](paper/figure_plan/Supplementary_8.png) retains the historical
+Top2a image separately and adds the new cell/marker evidence. The
+[preset comparison](paper/figure_plan/review/GSE246382_reclustering_presets.png)
+shows identical families under three fixed settings. These views describe
+model-antigen GC state bias without measuring direction or fate; the GC
+selection/output interpretation does not extend to non-GC data. Source and
+execution limits are in the
+[provenance manifest](paper/figure_plan/assets/legacy_gc_provenance.json).
+
+## Optional sampling-aware clone composition
+
+Threadfin 4.1 adds frozen-reference state distributions and train-only RNA
+projection. It also fixes kernel feature reconstruction used by programme
+bootstraps and snapshot analysis. Existing mean/kernel profile defaults remain
+available. See [the API and its limitations](docs/STATE_DENSITY.md) and
+[the algorithm review and validation](docs/ALGORITHM_REVIEW_2026-10-09_zh.md).
+
+New comparisons retain negative results: simple RNA means often predict the
+held-out division reporter better, and shrinkage can bias pure-state clones.

case_studies/benchmark_representation_runtime.py

@@ -0,0 +1,58 @@
+"""Isolated-process representation timing on identical synthetic inputs.
+
+Includes reference fitting and transform for density; includes default
+smoothing/RFF/PCA for kernel. These outputs solve different representation
+tasks, so speed is not evidence of biological accuracy or Clonotrace parity.
+"""
+import json
+from pathlib import Path
+import resource
+import subprocess
+import sys
+import time
+
+import numpy as np
+import pandas as pd
+
+OUT = Path(__file__).resolve().parent / "results" / "algorithm_revision"
+
+
+def worker(method, n_cells):
+    import anndata as ad
+    import threadfin as tf
+    rng = np.random.default_rng(17)
+    n_clones = n_cells//10
+    clone_codes = np.repeat(np.arange(n_clones), 10)
+    x = rng.normal(size=(n_cells, 30)) + rng.normal(size=(n_clones, 30))[clone_codes]
+    ids = np.array([f"c{i}" for i in clone_codes])
+    a = ad.AnnData(np.zeros((n_cells, 1)))
+    a.obs["clone_id"] = ids
+    a.obsm["X_pca"] = x
+    start = time.perf_counter()
+    if method == "density":
+        model = tf.StateDensityModel.fit(x, ids, n_states=16)
+        result = model.transform_batches((x[i:i+2048], ids[i:i+2048], None)
+                                         for i in range(0, len(x), 2048))
+        dimension = result.proportions.shape[1]
+    else:
+        tf.tl.clone_profiles(a, representation=method, verbose=False)
+        dimension = a.uns["threadfin"]["profiles"]["features"].shape[1]
+    return {"method": method, "n_cells": n_cells, "n_clones": n_clones,
+            "seconds": time.perf_counter()-start,
+            "peak_process_rss_mib": resource.getrusage(resource.RUSAGE_SELF).ru_maxrss / 1024,
+            "dimensions": dimension, "seed": 17}
+
+
+if __name__ == "__main__":
+    if len(sys.argv) > 1:
+        print(json.dumps(worker(sys.argv[1], int(sys.argv[2]))))
+    else:
+        OUT.mkdir(parents=True, exist_ok=True)
+        rows = []
+        for n in (3000, 15000, 50000):
+            for method in ("mean", "kernel", "density"):
+                output = subprocess.check_output([sys.executable, __file__, method, str(n)], text=True)
+                row = json.loads(output.strip().splitlines()[-1])
+                rows.append(row)
+                print(row, flush=True)
+                pd.DataFrame(rows).to_csv(OUT / "representation_runtime.csv", index=False)