LLM cell type annotation: ranking GPT-4 against classical tools with beam

Goal

Hou and Ji (2024) scored GPT-4 and GPT-3.5 against three reference-based annotators (CellMarker2.0, SingleR, ScType). Each method labels the cell types of a dataset. The score is the agreement with a manual expert annotation. This vignette re-runs that benchmark through beam. It reports how stable the ranking is across the modeling choices, how much of the ranking variance comes from the dataset, and whether the ranking changes between normal and cancer tissue. GPT-4 ranks first. GPT-3.5 ranks at or below the three classical tools, so the difference is specific to GPT-4 and not to language models in general.

Provenance of the bundled tables

beam does not ship the single-cell data or the model outputs. It ships two small derived tables (src/beam/data/gptcelltype2024_agreement.csv, gptcelltype2024_features.csv), reduced once from the compiled evaluation table in the reproduction repository Winnie09/GPTCelltype_Paper at a pinned commit. The article and its supplement are CC-BY-4.0; provenance and license are in src/beam/data/README.md. Cite Hou W, Ji Z. Assessing GPT-4 for cell type annotation in single-cell RNA-seq analysis. Nature Methods 2024, DOI 10.1038/s41592-024-02235-4.

To reproduce, run:

SHA=5944a41511aacd368b45448e256d9625849704df
curl -sf -o all.csv \
  "https://raw.githubusercontent.com/Winnie09/GPTCelltype_Paper/$SHA/anno/compiled/all.csv"

all.csv has one row per annotated cell type. Each method has a <method>_agreement column with the ontology-aware score (1 full match, 0.5 partial match, 0 mismatch). The reduction script is src/beam/data/reduce_gptcelltype.py:

python src/beam/data/reduce_gptcelltype.py        # fetches the pinned commit
python src/beam/data/reduce_gptcelltype.py all.csv  # or reduces a local copy

Load the tensor

load_gptcelltype reads the bundled tables into a frozen dataclass. beam treats each (source, tissue) pair as one dataset and each annotated cell type as one observation. The two metrics are the mean agreement score and the full match rate, both higher is better.

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

import beam
from beam.datasets import load_gptcelltype
from beam.cards import properties_for

g = load_gptcelltype()
print("scores shape:", g.scores.shape)
print("methods:", ", ".join(g.method_names))
print("datasets:", len(g.dataset_names))
print("metrics:", ", ".join(g.metric_ids))

scores shape: (6, 54, 2)
methods: GPT-4, GPT-4-mar2023, GPT-3.5, CellMarker2.0, SingleR, ScType
datasets: 54
metrics: cell_type_annotation_agreement, cell_type_annotation_full_match_rate

Coverage

The matrix has gaps. SingleR and ScType were not run on the Azimuth and literature datasets. The March 2023 GPT-4 endpoint covers a subset. beam reads the gaps as numpy.nan and does not impute them.

for i, m in enumerate(g.method_names):
    covered = int(np.sum(~np.isnan(g.scores[i, :, 0])))
    mean_agreement = np.nanmean(g.scores[i, :, 0])
    print(f"{m:14s} covers {covered:2d}/{len(g.dataset_names)} datasets   mean agreement {mean_agreement:.3f}")

GPT-4          covers 54/54 datasets   mean agreement 0.621
GPT-4-mar2023  covers 27/54 datasets   mean agreement 0.625
GPT-3.5        covers 54/54 datasets   mean agreement 0.452
CellMarker2.0  covers 54/54 datasets   mean agreement 0.268
SingleR        covers 36/54 datasets   mean agreement 0.450
ScType         covers 36/54 datasets   mean agreement 0.450

classical = {"CellMarker2.0", "SingleR", "ScType"}
mean_agreement = np.nanmean(g.scores[:, :, 0], axis=1)
o = np.argsort(-mean_agreement)
colors = ["tab:gray" if g.method_names[i] in classical else "tab:blue" for i in o]
fig, ax = plt.subplots(figsize=(7.0, 3.2))
ax.bar([g.method_names[i] for i in o], mean_agreement[o], color=colors)
ax.set_ylabel("mean agreement over datasets")
ax.set_xlabel("method (grey = classical, blue = LLM)")
ax.set_title("Mean annotation agreement per method")
ax.tick_params(axis="x", rotation=45)
fig.tight_layout()
plt.show()

The main comparison uses the five methods Hou and Ji compare directly (the March 2023 GPT-4 endpoint set aside), on the datasets where all five were run.

head = ["GPT-4", "GPT-3.5", "CellMarker2.0", "SingleR", "ScType"]
idx = [g.method_names.index(m) for m in head]
metrics = list(g.metric_ids)

sub = g.scores[idx]                         # 5 methods x 54 datasets x 2 metrics
complete = ~np.isnan(sub[:, :, 0]).any(axis=0)
block = sub[:, complete, :]                 # 5 x 36 x 2, no missing cells
block_datasets = [d for d, keep in zip(g.dataset_names, complete) if keep]
print(f"complete block: {block.shape[0]} methods, {block.shape[1]} datasets, {block.shape[2]} metrics")

complete block: 5 methods, 36 datasets, 2 metrics

Pull the metric cards

properties_for reads the two annotation cards. Both are higher is better, bounded in [0, 1], and use min-max normalization with the arithmetic mean across datasets.

props = properties_for(metrics)
for p in props:
    print(
        f"{p.id:36s} polarity={p.polarity:16s} scale={p.scale_type:8s} "
        f"range=({p.range_lower}, {p.range_upper}) norm={p.recommended_normalization}"
    )

cell_type_annotation_agreement       polarity=higher_is_better scale=interval range=(0, 1) norm=min_max
cell_type_annotation_full_match_rate polarity=higher_is_better scale=ratio    range=(0, 1) norm=min_max

One MCDA run on the complete cases

beam pools the two metrics across the 36 datasets and ranks the five methods. The default is equal weights with simple additive weighting (SAW). beam.rank also runs the default sensitivity battery on the same normalized matrix.

block_scores = beam.Scores(
    values=block,
    tool_names=tuple(head),
    metric_ids=tuple(metrics),
    dataset_names=tuple(block_datasets),
    layout="long",
)
run = beam.rank(block_scores, weights="equal", method="saw", sensitivity=True)

print(f"weighting={run.result.weighting}  method={run.result.method}")
order = np.argsort(run.result.ranks)
print(f"\n{'method':14s} composite  rank")
for i in order:
    print(f"{head[i]:14s} {run.result.composite[i]:.3f}      {run.result.ranks[i]}")

weighting=equal  method=saw

method         composite  rank
GPT-4          0.506      1
SingleR        0.334      2
GPT-3.5        0.323      3
ScType         0.319      4
CellMarker2.0  0.183      5

GPT-4 ranks first. SingleR, a classical tool, ranks second, and GPT-3.5 ranks third.

The card audit

card_data_consistency checks the pooled scores against the cards before any normalization: every value inside the declared [0, 1] range, every baseline and floor sane. A unit or scale error would be named here.

audit = run.card_consistency
print("pooled scores consistent with the cards:", audit.ok)
for finding in audit.findings:
    print(" ", finding.severity, finding.message)

pooled scores consistent with the cards: True

Compare across weightings and aggregation methods

Four weightings (equal, entropy, std, critic) crossed with five aggregation methods give twenty runs. A method that holds its rank down a column is stable to these choices.

from beam.mcda import run_from_registry

pooled = np.nanmean(block, axis=1)        # 5 methods x 2 metrics, complete
weightings = ["equal", "entropy", "std", "critic"]
agg_methods = ["saw", "topsis", "vikor", "promethee_ii", "comet"]

combos = [(w, m) for w in weightings for m in agg_methods]
ranks_grid = np.zeros((len(combos), len(head)), dtype=int)
row_labels = []
for i, (w, m) in enumerate(combos):
    r = run_from_registry(pooled, metrics, weights=w, method=m)
    ranks_grid[i] = r.ranks
    row_labels.append(f"{w} / {m}")

from beam import plot

plot.rank_heatmap(
    ranks_grid,
    row_names=row_labels,
    col_names=head,
    row_label="weighting / aggregation method",
    col_label="annotation method",
    title="Rank by configuration (both metrics higher is better)",
)

top_per_combo = {head[int(np.argmin(row))] for row in ranks_grid}
print("methods top-ranked in at least one of the 20 configurations:", ", ".join(sorted(top_per_combo)))

methods top-ranked in at least one of the 20 configurations: GPT-4

aggregation_agreement measures that: the Kendall tau-b agreement of the rankings the five aggregations give at equal weights, and whether they agree on the top method.

from beam.mcda import aggregation_agreement

agg = aggregation_agreement(pooled, [p.polarity for p in props], tool_names=head)
print(f"mean pairwise Kendall tau-b across aggregations: {agg.mean_pairwise_tau:.2f}")
print(f"top method unanimous across aggregations: {agg.top_is_unanimous} ({head[agg.top_tool]})")

mean pairwise Kendall tau-b across aggregations: 0.98
top method unanimous across aggregations: True (GPT-4)

Rank stability

The SMAA analysis samples weight vectors and records how often each method takes each rank. The weight-perturbation check asks the smallest single-weight change that would flip the top method.

smaa = run.smaa
top = int(np.argmin(run.result.ranks))
print(f"SMAA: {head[top]} is rank 1 in {smaa.rank_acceptability_index[top, 0] * 100:.0f} percent of weight draws")

wp = run.perturbation
print(f"top rank fragile: {wp.top_rank_is_fragile}")
if wp.top_rank_perturbation is not None:
    print(f"smallest single-weight change at the top pair: {wp.top_rank_perturbation.absolute_delta:.2f}")
else:
    print("no single-weight change within range alters the top rank")

SMAA: GPT-4 is rank 1 in 100 percent of weight draws
top rank fragile: False
no single-weight change within range alters the top rank

plot.smaa(run)

leave_one_metric_out drops each metric in turn. With two metrics this ranks on each metric alone.

loo = run.leave_one_out
print(f"top method rank stability under leave-one-metric-out: {loo.rank_stability[top]:.2f}")
print(f"most influential metric: {metrics[loo.most_influential_metric]}")

top method rank stability under leave-one-metric-out: 1.00
most influential metric: cell_type_annotation_full_match_rate

Above chance, and the noise floor

beats_random_baseline reads the chance level the card declares. The full match rate has a declared chance level of 0, since a random label almost never matches the exact manual type. Every method scores above 0. noise_floor_separation would flag method pairs closer than a declared noise floor. The annotation cards declare no floor yet, so it reports as inactive.

br = run.random_baseline
print("random-baseline check active:", br.active)
if br.active:
    for pm in br.per_metric:
        print(f"  {pm.metric:36s} {pm.n_beating}/{pm.n_observed} methods score above chance ({pm.baseline})")
    print("  methods above chance on no metric:", ", ".join(br.tools_never_beating) or "none")

nf = run.noise_floor
print("noise-floor check active:", nf.active, "(no floor declared on the annotation cards)")

random-baseline check active: True
  cell_type_annotation_full_match_rate 5/5 methods score above chance (0.0)
  methods above chance on no metric: none
noise-floor check active: False (no floor declared on the annotation cards)

Critical difference on agreement

A Demsar critical-difference diagram ranks the methods on each of the 36 datasets and asks, with a Friedman test and a Nemenyi post-hoc, whether the differences are larger than sampling noise.

from beam.mcda import critical_difference

cd = critical_difference(block[:, :, 0], "higher_is_better", tool_names=head)
print(f"Friedman p-value: {cd.friedman_pvalue:.2e}")
print(f"critical difference (alpha={cd.alpha}): {cd.critical_difference:.2f}")
print("average rank (lower is better):")
for i in np.argsort(cd.average_ranks):
    print(f"  {head[i]:14s} {cd.average_ranks[i]:.2f}")
print("cliques (methods not separable):")
for clique in cd.cliques:
    print("  " + ", ".join(head[i] for i in clique))

Friedman p-value: 4.21e-20
critical difference (alpha=0.05): 1.02
average rank (lower is better):
  GPT-4          1.17
  ScType         2.94
  SingleR        2.99
  GPT-3.5        3.14
  CellMarker2.0  4.76
cliques (methods not separable):
  ScType, SingleR, GPT-3.5

plot.critical_difference(cd)

pairwise_superiority is the effect-size companion: for each method pair, the share of the 36 datasets where one method scores higher, the readable version of the significance the diagram gives.

from beam.mcda import pairwise_superiority

ps = pairwise_superiority(block[:, :, 0], "higher_is_better", method_names=head)
print("standing (Copeland-style, higher means ahead on more pairs):")
for i in np.argsort(-ps.standing):
    print(f"  {head[i]:14s} {ps.standing[i]:.2f}")

standing (Copeland-style, higher means ahead on more pairs):
  GPT-4          0.96
  ScType         0.51
  SingleR        0.50
  GPT-3.5        0.47
  CellMarker2.0  0.06

pairwise_transitivity reads the same relation and reports whether one order is consistent with the pairwise majorities.

from beam.mcda import pairwise_transitivity

trans = pairwise_transitivity(ps)
print(f"transitive: {trans.is_transitive}; circular triads: {trans.n_circular_triads} of {trans.n_triads}")

transitive: True; circular triads: 0 of 10

The matrix below orders the methods by how many others they outperform. A transitive relation fills the upper triangle; a red cell below the diagonal marks a method that outperforms one ranked above it, which only happens inside a cycle.

plot.pairwise_majority(trans)

bayesian_sign_comparison reports, for each pair, the posterior probability that one method is practically better than another. With 36 datasets the posterior is sharper, so more pairs reach a decisive label than on the smaller benchmarks.

from beam.mcda import bayesian_sign_comparison

bayes = bayesian_sign_comparison(ps)
decisive = sum(1 for p in bayes.per_pair if p.decision != "inconclusive")
print(f"pairs with a decisive posterior at 0.95: {decisive} of {len(bayes.per_pair)}")
plot.bayesian_comparison(bayes)

pairs with a decisive posterior at 0.95: 7 of 10

Rank sensitivity: data or analyst choice

The sections above each measured part of whether the ranking holds under a different analysis. rank_sensitivity combines them: it ranks the methods on each dataset, under every weighting and aggregation, and splits each method’s rank variance into the dataset, the weighting, the aggregation and their interactions.

from beam.mcda import rank_sensitivity

ctx = run.context
rs = rank_sensitivity(
    block,
    ctx.polarity,
    normalization=list(ctx.normalization),
    bounds=list(ctx.bounds),
    baselines=list(ctx.baselines),
    targets=list(ctx.targets),
    tool_names=head,
    dataset_names=block_datasets,
)
print(f"{rs.n_combinations} combinations")
print(f"  dataset:      {rs.dataset_share:.3f} of the rank variance")
print(f"  weighting:    {rs.weighting_share:.3f}")
print(f"  aggregation:  {rs.aggregation_share:.3f}")
print(f"  interactions: {rs.interaction_share:.3f}")
print(f"  most influential factor: {rs.most_influential_factor}")

720 combinations
  dataset:      0.962 of the rank variance
  weighting:    0.001
  aggregation:  0.000
  interactions: 0.037
  most influential factor: dataset

plot.rank_sensitivity(rs)

The dataset accounts for almost all of the rank variance. The weighting and the aggregation account for less than one percent combined. Leave-one-dataset-out reports whether the pooled order changes when one dataset is dropped.

The figures so far pool the shares over all methods. The per-method breakdown attributes each method’s rank spread on its own, telling apart a method that moves with the dataset from one that moves with the analyst’s choice. The bar span is the method’s best minus worst rank.

plot.rank_sensitivity_by_tool(rs)

lodo = run.leave_one_dataset_out
print(f"top method rank stability under leave-one-dataset-out: {lodo.rank_stability[top]:.2f}")
print(f"most influential dataset: {block_datasets[lodo.most_influential_dataset]}")

top method rank stability under leave-one-dataset-out: 1.00
most influential dataset: HCA_Skeletal_muscle

specification_curve lists the rankings the grid produces. The choice-only multiverse (weighting by aggregation on the pooled matrix) shows whether GPT-4 holds the top across the modeling choices; the full grid adds the dataset.

from beam.mcda import specification_curve

curve = specification_curve(rs)
dom = curve.tool_names[curve.most_frequent_top_tool]
print(f"choices plus dataset ({curve.n_specifications} combinations): "
      f"{dom} first in {curve.most_frequent_top_fraction * 100:.0f}%, "
      f"{curve.n_distinct_top_tools} methods reach the top")

pooled = specification_curve(
    rank_sensitivity(
        run.matrix, ctx.polarity,
        normalization=list(ctx.normalization), bounds=list(ctx.bounds),
        baselines=list(ctx.baselines), targets=list(ctx.targets),
        tool_names=head,
    )
)
pdom = pooled.tool_names[pooled.most_frequent_top_tool]
print(f"choices only ({pooled.n_specifications} combinations): "
      f"{pdom} first in {pooled.most_frequent_top_fraction * 100:.0f}%")

choices plus dataset (720 combinations): GPT-4 first in 94%, 3 methods reach the top
choices only (20 combinations): GPT-4 first in 100%

Dataset concordance

The pooled order sits on top of 36 per-dataset orders. dataset_concordance recovers each of those, ranking the methods within a dataset and comparing every pair of orderings with Kendall tau-b. A high mean says the datasets broadly agree with the pooled order; a low mean says they do not. It assumes nothing about the datasets being interchangeable and needs no replicates.

conc = run.dataset_concordance
names = conc.dataset_names
print(f"mean agreement across datasets (Kendall tau-b): {conc.mean_pairwise_tau:.2f}")
print("least typical dataset:", names[conc.most_idiosyncratic_dataset])
print("number of mutually consistent dataset groups:", len(conc.concordant_groups))
print("where methods depart most from their own average rank:")
for cell in conc.notable_cells[:5]:
    side = "lower" if cell.deviation > 0 else "higher"
    print(f"  {conc.tool_names[cell.tool]} on {names[cell.dataset]}: "
          f"rank {cell.rank}, {side} than its mean {cell.mean_rank:.1f}")

mean agreement across datasets (Kendall tau-b): 0.53
least typical dataset: HCA_Lung
number of mutually consistent dataset groups: 1
where methods depart most from their own average rank:
  SingleR on tabulasapiens_Heart: rank 5, lower than its mean 3.0
  SingleR on tabulasapiens_Pancreas: rank 5, lower than its mean 3.0
  SingleR on tabulasapiens_Skin: rank 5, lower than its mean 3.0
  ScType on HCA_Heart: rank 5, lower than its mean 3.0
  ScType on HCA_Lung: rank 1, higher than its mean 3.0

plot.dataset_concordance(run)

The companion view marks where each method places higher or lower than its own typical rank, showing which methods carry the dataset disagreement rather than which ranks best overall.

plot.dataset_struggle(run)

Blind analysis

A blind analysis fixes the pipeline before the method names are known. This matters here, where one method (GPT-4) is the one under test: blinding removes the chance to tune the weighting toward it. beam.blind masks the names and shuffles the rows; beam.unblind restores them. The ranking is unchanged, and the seal fingerprint is recorded in the manifest.

from beam import blind, unblind

blinded, seal = blind(block_scores, seed=0)
blind_run = beam.rank(blinded, weights="equal", method="saw", seed=0, sensitivity=False)
restored = unblind(blind_run, seal)
named = beam.rank(block_scores, weights="equal", method="saw", seed=0, sensitivity=False)
print("ranking identical after unblinding:",
      dict(zip(named.tool_names, named.result.ranks))
      == dict(zip(restored.tool_names, restored.result.ranks)))
print("top method after unblinding:", restored.top_tool)
print("blinding fingerprint:", blind_run.manifest["blinding"]["seal_sha256"][:12])

ranking identical after unblinding: True
top method after unblinding: GPT-4
blinding fingerprint: 50b6fc613b9f

Equal versus cell-weighted pooling

beam weights each dataset equally by default. The paper’s per-cell-type accuracy weights a dataset by how many cell types it has. The dataset accounts for almost all of the variance, so the dataset weight changes the magnitudes. Here it does not change the top rank.

n_cell_types = np.array(g.features.numeric["n_cell_types"])[complete]
agreement = block[:, :, 0]
equal_pool = agreement.mean(axis=1)
weighted_pool = (agreement * n_cell_types).sum(axis=1) / n_cell_types.sum()
print(f"{'method':14s} equal-per-dataset  cell-weighted")
for i in np.argsort(-equal_pool):
    print(f"{head[i]:14s} {equal_pool[i]:.3f}             {weighted_pool[i]:.3f}")

method         equal-per-dataset  cell-weighted
GPT-4          0.603             0.564
SingleR        0.450             0.436
ScType         0.450             0.427
GPT-3.5        0.436             0.399
CellMarker2.0  0.249             0.244

Mixed-effects on agreement

rank_sensitivity decomposes the ranks. A mixed-effects model decomposes the scores: it fits agreement ~ method + (1 | dataset) in R’s lme4, separating the stable method effect from the dataset being uniformly easy or hard. The fit needs the R toolchain, so this section runs only when it is available.

from beam.heterogeneity import mixed_effects_from_matrix, r_available

if r_available():
    me = mixed_effects_from_matrix(agreement, head, block_datasets)
    print(f"model: {me.formula}")
    print(f"dataset shift (ICC): {me.icc_dataset:.2f} of the agreement variance")
    order = np.argsort(-me.method_effects)
    print("\nmethod          agreement marginal mean")
    for i in order:
        print(f"{me.method_names[i]:14s}  {me.method_effects[i]:.3f} +/- {me.method_effect_se[i]:.3f}")
else:
    me = None
    print("R with lme4 not available; skipping the mixed-effects fit.")
    print("Provision it with envs/heterogeneity.yml (conda or mamba).")

model: score ~ method + (1 | dataset)
dataset shift (ICC): 0.28 of the agreement variance

method          agreement marginal mean
GPT-4           0.603 +/- 0.021
SingleR         0.450 +/- 0.021
ScType          0.450 +/- 0.021
GPT-3.5         0.436 +/- 0.021
CellMarker2.0   0.249 +/- 0.021

Bradley-Terry tree on tissue features

A further question is whether the ranking reverses on particular kinds of data. A Bradley-Terry tree turns each dataset into pairwise method comparisons and splits the datasets by their features, so each leaf gets its own ranking and a parameter-stability test decides where a split is real. The candidate features are the source, tissue, species, sample type (normal or cancer) and the number of cell types. The fit needs the R toolchain.

from beam.heterogeneity import bradley_terry_tree, bttree_available

if bttree_available():
    numeric, categorical = g.features.aligned_to(block_datasets)
    bt = bradley_terry_tree(
        agreement,
        head,
        block_datasets,
        numeric_features=numeric,
        categorical_features=categorical,
        polarity="higher_is_better",
        minsize=6,
    )
    print(f"split found: {bt.did_split}")
    print(f"global ranking (top 3): {', '.join(bt.global_ranking()[:3])}")
    print(bt.summary())
else:
    bt = None
    print("R with psychotree not available; skipping the Bradley-Terry tree.")
    print("Provision it with envs/heterogeneity.yml (conda or mamba).")

split found: False
global ranking (top 3): GPT-4, ScType, SingleR
The Bradley-Terry tree found no dataset feature that splits the method ranking at alpha 0.05 over 36 datasets, so the ranking is reported as one Bradley-Terry model over all of them, led by GPT-4. With this many datasets the split test has few observations to work with, the same small-sample limit the critical-difference diagram shows; a benchmark with more datasets is where a split can appear.

plackett_luce fits the same preference idea as a single full-ranking model over the methods, a check on the Bradley-Terry global ranking.

from beam.heterogeneity import plackett_luce, plackett_luce_available

if plackett_luce_available():
    pl = plackett_luce(agreement, head, polarity="higher_is_better")
    print("Plackett-Luce ranking (first to last):", ", ".join(pl.ranking()))
    print("top method:", pl.top_tool())
else:
    print("R with PlackettLuce not available; skipping.")
    print("Provision it with envs/heterogeneity.yml (conda or mamba).")

Plackett-Luce ranking (first to last): GPT-4, ScType, SingleR, GPT-3.5, CellMarker2.0
top method: GPT-4

Reliability of the manual annotation

The agreement metric is defined against the manual annotation, so beam uses it as the reference and cannot check it directly. One indirect check treats the five annotators as independent raters and looks for cell types where they agree on a broad lineage that differs from the manual label. Those cell types are candidates for a fine-grained mislabel or an ambiguous reference.

import csv
import collections
from importlib import resources

text = resources.files("beam.data").joinpath("gptcelltype2024_agreement.csv").read_text("utf-8")
ann_rows = [r for r in csv.DictReader(text.splitlines()) if r["method"] in head]

by_cell = collections.defaultdict(list)
for r in ann_rows:
    if r["broadtype"]:
        by_cell[(r["dataset"], r["cell_type"], r["manual_broadtype"])].append(r["broadtype"])

candidates = []
for (dataset, cell_type, manual), predicted in by_cell.items():
    if len(predicted) < 3:
        continue
    consensus, n = collections.Counter(predicted).most_common(1)[0]
    manual_lineages = {x.strip() for x in manual.split(",")}
    if n / len(predicted) >= 0.8 and consensus not in manual_lineages:
        candidates.append((dataset, cell_type, manual, consensus, f"{n}/{len(predicted)}"))

eligible = sum(1 for v in by_cell.values() if len(v) >= 3)
print(f"cell types where the methods agree on a lineage that differs from the manual one:")
print(f"  {len(candidates)} of {eligible} cell types with at least three running methods")
print("\nexamples (dataset, cell type, manual lineage, method consensus, votes):")
for dataset, cell_type, manual, consensus, votes in candidates[:8]:
    print(f"  {dataset:22s} {cell_type:26s} {manual:18s} -> {consensus} ({votes})")

cell types where the methods agree on a lineage that differs from the manual one:
  40 of 1029 cell types with at least three running methods

examples (dataset, cell type, manual lineage, method consensus, votes):
  Azimuth_Pancreas       immune                     immune cell        -> macrophage/monocyte (3/3)
  Azimuth_Fetal_Development Stromal cells              stromal cell       -> fibroblast (3/3)
  Azimuth_Lung           CD14+ Monocyte             macrophage/monocyte -> granulocyte (3/3)
  Azimuth_Lung           Capillary                  blood vessel cell  -> endothelial cell (3/3)
  Azimuth_Lung           Mucous                     mucosal cell       -> epithelial cell (3/3)
  Azimuth_Lung           Lipofibroblast             fibroblast         -> epithelial cell (3/3)
  Azimuth_Bone_Marrow    Innate Lymphoid Cell       innate lymphoid cell -> t/nk cell (3/3)
  Azimuth_Bone_Marrow    Platelet                   platelet           -> megakaryocyte (3/3)

Most of these are a granularity difference rather than an error: the manual label is a fine type (innate lymphoid cell, platelet, capillary) and the methods agree on its parent lineage (t/nk cell, megakaryocyte, endothelial cell). In those cells the agreement score counts the method as wrong even though it named a broader correct lineage.

Partial coverage

The complete cases dropped the 18 datasets where SingleR and ScType were not run. beam can also rank on the full set without completing the matrix. Available-case SAW scores each method on the datasets where it was measured, with the weights renormalized over its support; the Skillings-Mack test is the coverage-aware Friedman test for a matrix with gaps.

from beam.mcda import skillings_mack
from beam.reporting import funky_heatmap_from_run

full_scores = beam.Scores(
    values=sub,
    tool_names=tuple(head),
    metric_ids=tuple(metrics),
    dataset_names=tuple(g.dataset_names),
    layout="long",
)
full_run = beam.rank(
    full_scores, weights="equal", method="saw", sensitivity=True, missing="available"
)
print("available-case ranking on all 54 datasets:")
for i in np.argsort(full_run.result.ranks):
    print(f"  {head[i]:14s} rank {full_run.result.ranks[i]}")

sm = skillings_mack(sub[:, :, 0], higher_is_better=True, method_names=head)
print(f"\nSkillings-Mack statistic: {sm.statistic:.1f}  p-value: {sm.p_value:.2e}")

available-case ranking on all 54 datasets:
  GPT-4          rank 1
  GPT-3.5        rank 2
  SingleR        rank 3
  ScType         rank 4
  CellMarker2.0  rank 5

Skillings-Mack statistic: 114.5  p-value: 7.84e-24

The funky heatmap shows the available-case run as a glyph table: the two metrics per method, the composite, and the rank-robustness panels, on all 54 datasets rather than the 36 complete cases.

funky_heatmap_from_run(full_run, title="Cell type annotation: available-case ranking on 54 datasets")
plt.show()

/home/runner/work/beam/beam/src/beam/reporting/figures.py:778: UserWarning: Tight layout not applied. The bottom and top margins cannot be made large enough to accommodate all Axes decorations.
  fig.tight_layout()

Cross-benchmark: LLM annotation benchmarks

Three more recent benchmarks score LLMs at cell type annotation: DeepCellSeek (Briefings in Bioinformatics 2025, 10.1093/bib/bbaf677), AnnDictionary (Nature Communications 2025, 10.1038/s41467-025-64511-x) and Single-Cell Omics Arena (SOAR, arXiv:2412.02915). They use different models, datasets and metrics, so they cannot be pooled on a shared score. The common currency is the within-benchmark rank of the methods each benchmark has in common with the others. The three classical tools (CellMarker2.0, SingleR, ScType) connect GPTCelltype, DeepCellSeek and SOAR; GPT-4o connects DeepCellSeek, AnnDictionary and SOAR. The provenance and license of each bundled table are in src/beam/data/; DeepCellSeek and SOAR scores come from a pinned source, and the SOAR and AnnDictionary folders carry a note on their license status.

import csv
from importlib import resources
from beam.datasets import load_deepcellseek

CLASSICAL = {"CellMarker2.0", "SingleR", "ScType", "Cell2Sentence"}

def ranking(scores_by_method):
    order = sorted(scores_by_method, key=lambda m: -scores_by_method[m])
    return {m: i + 1 for i, m in enumerate(order)}

benchmarks = {}
benchmarks["GPTCelltype"] = {m: float(np.nanmean(g.scores[i, :, 0]))
                             for i, m in enumerate(g.method_names)}
dcs = load_deepcellseek()
benchmarks["DeepCellSeek"] = {m: float(np.nanmean(dcs.scores[i, :, 0]))
                              for i, m in enumerate(dcs.method_names)}
ad = {r["model"]: float(r["overall_binary_celltypes_mean"]) for r in csv.DictReader(
    resources.files("beam.data").joinpath(
        "anndictionary/anndictionary_celltype_agreement.csv").read_text().splitlines())}
benchmarks["AnnDictionary"] = ad
soar = {r["model"]: float(r["R-L"]) for r in csv.DictReader(
    resources.files("beam.data").joinpath("soar/soar_readme_scores.csv").read_text().splitlines())
    if r["table_name"] == "SOAR-RNA Zero-shot" and r["R-L"]}
benchmarks["SOAR"] = soar

for name, sc in benchmarks.items():
    rk = ranking(sc)
    n = len(sc)
    classical = [m for m in sc if m in CLASSICAL]
    if classical:
        top = min(rk[m] for m in classical)
        print(f"{name:13s} {n:2d} methods; top-ranked classical ranks {top}/{n}; metric: "
              f"{'ontology agreement' if name in ('GPTCelltype','DeepCellSeek') else 'NLP overlap (ROUGE-L)'}")
    else:
        print(f"{name:13s} {n:2d} methods; LLM-only (no classical comparator)")

GPTCelltype    6 methods; top-ranked classical ranks 4/6; metric: ontology agreement
DeepCellSeek  14 methods; top-ranked classical ranks 12/14; metric: ontology agreement
AnnDictionary 14 methods; LLM-only (no classical comparator)
SOAR          11 methods; top-ranked classical ranks 6/11; metric: NLP overlap (ROUGE-L)

In the two ontology-agreement benchmarks every classical tool ranks below every LLM. SOAR scores annotation text with ROUGE-L instead, and there CellMarker2.0 ranks in the middle, above several open-source LLMs, while SingleR and ScType stay near the bottom.

shared_methods = ["GPT-4o", "GPT-4", "CellMarker2.0", "SingleR", "ScType"]
bench_names = list(benchmarks)
grid = np.full((len(shared_methods), len(bench_names)), np.nan)
for j, name in enumerate(bench_names):
    rk = ranking(benchmarks[name])
    n = len(benchmarks[name])
    for i, m in enumerate(shared_methods):
        if m in rk:
            grid[i, j] = rk[m] / n          # normalized rank, 0 first, 1 last

fig, ax = plt.subplots(figsize=(6.0, 3.2))
im = ax.imshow(grid, cmap="RdYlGn_r", aspect="auto", vmin=0, vmax=1)
ax.set_xticks(range(len(bench_names))); ax.set_xticklabels(bench_names, rotation=20, ha="right")
ax.set_yticks(range(len(shared_methods))); ax.set_yticklabels(shared_methods)
ax.set_xlabel("benchmark"); ax.set_ylabel("shared method")
ax.set_title("Normalized rank of the shared methods (0 = first, 1 = last)")
for i in range(len(shared_methods)):
    for j in range(len(bench_names)):
        if not np.isnan(grid[i, j]):
            ax.text(j, i, f"{grid[i, j]:.2f}", ha="center", va="center", fontsize=8)
fig.colorbar(im, ax=ax, fraction=0.05, label="normalized rank")
fig.tight_layout()
plt.show()

The benchmarks agree that the top-ranked LLM outranks the classical tools, but they disagree on the rest. CellMarker2.0 ranks last under ontology agreement (GPTCelltype, DeepCellSeek) and in the middle under SOAR’s ROUGE-L. GPT-4o ranks first in AnnDictionary and SOAR but in the middle in DeepCellSeek, where the newer Grok-4, GPT-5 and Claude-4.1 score higher. The ranking depends on the metric and on the model generation each benchmark tested, the disagreement beam measures. A full network meta-analysis is left out here: GPTCelltype and DeepCellSeek provide per-dataset scores, but the AnnDictionary and SOAR tables are aggregate-only, so the network lacks the per-study arms the pooled model needs.

Quantitative cross-benchmark: variance on the shared classical methods

GPTCelltype and DeepCellSeek both score the three classical tools per dataset on the same ontology-aware metric, so their scores can be pooled and the variance attributed. source_variance_decomposition fits score ~ method + (1 | benchmark) + (1 | benchmark:dataset) + (1 | method:benchmark) in lme4 over the two benchmarks and the three shared methods. It separates the benchmark level effect, the method-by-benchmark interaction (whether the method order changes between benchmarks), the dataset nested in benchmark, and the residual. SOAR and AnnDictionary stay out of this fit: SOAR scores on a different metric and AnnDictionary has no classical methods. The fit needs the R toolchain.

from beam.heterogeneity import source_variance_decomposition, r_available

shared = ["CellMarker2.0", "SingleR", "ScType"]
sv_methods, sv_datasets, sv_benchmarks, sv_scores = [], [], [], []
for bname, b in [("GPTCelltype", g), ("DeepCellSeek", dcs)]:
    for m in shared:
        i = b.method_names.index(m)
        for j, ds in enumerate(b.dataset_names):
            v = b.scores[i, j, 0]
            if not np.isnan(v):
                sv_methods.append(m)
                sv_datasets.append(ds)
                sv_benchmarks.append(bname)
                sv_scores.append(float(v))

if r_available():
    sv = source_variance_decomposition(sv_methods, sv_datasets, sv_benchmarks, sv_scores)
    total = sum(v for v in sv.variance_components.values() if not np.isnan(v))
    print(f"{len(sv_scores)} observations, 2 benchmarks, 3 shared methods")
    for k, v in sv.variance_components.items():
        if not np.isnan(v):
            print(f"  {k:18s} {v / total:.2f} of the score variance")
    print("pooled method effect:", {n: round(e, 3) for n, e in zip(sv.method_names, sv.method_effects)})
else:
    sv = None
    print("R with lme4 not available; skipping the variance decomposition.")
    print("Provision it with envs/heterogeneity.yml (conda or mamba).")

195 observations, 2 benchmarks, 3 shared methods
  benchmark:dataset  0.00 of the score variance
  method:benchmark   0.10 of the score variance
  benchmark          0.24 of the score variance
  Residual           0.66 of the score variance
pooled method effect: {'CellMarker2.0': np.float64(0.389), 'ScType': np.float64(0.483), 'SingleR': np.float64(0.511)}

if sv is not None:
    comps = {k: v for k, v in sv.variance_components.items() if not np.isnan(v)}
    total = sum(comps.values())
    fig, axes = plt.subplots(1, 2, figsize=(9.0, 3.2))
    axes[0].bar(list(comps), [comps[k] / total * 100 for k in comps], color="tab:purple")
    axes[0].set_ylabel("share of score variance (percent)")
    axes[0].set_xlabel("source")
    axes[0].set_title("Variance on the three classical methods")
    axes[0].tick_params(axis="x", rotation=30)
    x = np.arange(len(shared))
    gm = [np.nanmean(g.scores[g.method_names.index(m), :, 0]) for m in shared]
    dm = [np.nanmean(dcs.scores[dcs.method_names.index(m), :, 0]) for m in shared]
    axes[1].bar(x - 0.2, gm, 0.4, label="GPTCelltype")
    axes[1].bar(x + 0.2, dm, 0.4, label="DeepCellSeek")
    axes[1].set_xticks(x)
    axes[1].set_xticklabels(shared, rotation=20, ha="right")
    axes[1].set_ylabel("mean agreement")
    axes[1].set_title("Classical methods per benchmark")
    axes[1].legend()
    fig.tight_layout()
    plt.show()

The benchmark accounts for about a quarter of the score variance and the method-by-benchmark interaction about a tenth, so the classical-method comparison is only partly portable across the two benchmarks. DeepCellSeek scores all three classical tools higher than GPTCelltype does, a level shift, and their relative order shifts somewhat between the two. The pooled order is SingleR, then ScType, then CellMarker2.0. Most of the variance is residual, the per-dataset spread at one score per cell type.

Recommendation

On the 36 datasets where all five methods were run, GPT-4 ranks first on the mean agreement and the full match rate pooled across datasets. It keeps that position across all twenty weighting-by-aggregation configurations. The SMAA draws and the weight-perturbation check report the top rank as not fragile. SingleR, a classical reference-based tool, ranks second, and GPT-3.5 ranks with the classical tools rather than with GPT-4.

On the critical-difference diagram GPT-4 separates from the other four methods, which sit in overlapping cliques and are not separable on agreement at 36 datasets. The rank-sensitivity decomposition puts about 96 percent of the rank variance on the dataset and less than one percent on the weighting and the aggregation combined.