Coverage-standardized Hill diversity (tl.hill_diversity_profile)

[KilianMaire]

Close #535

Summary

Adds coverage-standardized Hill-number diversity to scirpy.tl:

• tl.hill_diversity_profile computes a Hill diversity profile over a range of orders q, standardized to a common sample coverage so that profiles are comparable across samples of different sequencing depth.
• tl.convert_hill_table converts a profile into the classical alpha diversity indices (observed richness, Shannon entropy, inverse Simpson, Gini-Simpson) and into evenness measures.

This builds on #535 and addresses the open question raised there about sequencing-depth correction. Credit to @MKanetscheider for the original hill_diversity_profile / convert_hill_table design and for the convert_hill_table conversions requested by @FFinotello; this PR keeps those public signatures and replaces the estimation engine.

Motivation

The plug-in Hill estimator grows with sampling depth at every order of q (richness most strongly, but inverse Simpson too), so two samples sequenced to different depth show different profiles even when the underlying repertoire is identical. This is the confounder discussed in #535, and it is more acute for scRNA-seq where the number of cells varies a lot between samples.

The established fix is the iNEXT framework (Chao et al. 2014; Hsieh, Ma & Chao 2016): estimate Hill numbers and standardize them to a common sample coverage before comparing. tl.hill_diversity_profile standardizes all groups to a shared coverage (iNEXT's Cmax rule) and returns the standardized profile.

What is new compared to #535

• The naive plug-in engine is replaced with the coverage-standardized estimator.
• When the groups cannot be standardized to a fully reliable shared coverage (for example one group is heavily undersampled), a warning is raised rather than silently returning a number.
• Unit tests are added for both functions, including the warning path. The numeric output is checked against the underlying estimator to a relative tolerance of 1e-9, for both AnnData and MuData.

Implementation notes

• Estimation is delegated to the hillrep package, whose kernels are validated against R iNEXT 3.0.2 to a relative tolerance of 1e-6. hillrep is pure numpy/scipy/pandas and is added as an optional dependency under the existing diversity extra, alongside scikit-bio.
• The backend is isolated behind a single private helper (_coverage_hill_profile), so swapping the dependency for a vendored estimator later would be a localized change.
• Group counts are extracted via DataHandler, so both AnnData and MuData are supported (tests cover both).

Open question for maintainers

One design decision I would like your call on:

• (A) scirpy takes hillrep as an optional dependency (as in this PR), so the validated kernels live in one place and stay in sync with the R reference; or
• (B) the coverage estimator is vendored into _diversity.py, keeping scirpy dependency-free at the cost of duplicating the math.

This PR implements (A) because hillrep is dependency-light, but it is your dependency policy and your call. Switching to (B) would be localized to _coverage_hill_profile.

API

import scirpy as ir

profile = ir.tl.hill_diversity_profile(
    adata, groupby="sample", target_col="clone_id", q_min=0, q_max=2, q_step=1
)
indices = ir.tl.convert_hill_table(profile, convert_to="diversity")

hill_diversity_profile returns a DataFrame with one row per diversity order q and one column per group, which plots directly with seaborn and flows into convert_hill_table.

Checklist

• CHANGELOG.md updated
• Tests added (for the new functions, including the not-comparable warning)
• Tutorial updated (Hill diversity section added to the 5k BCR tutorial)

[KilianMaire]

Note on the red CI: all the tests added in this PR pass on every environment (test_hill_diversity_profile[AnnData], [MuData], test_convert_hill_table, test_hill_diversity_profile_warns_when_not_comparable).

The two failing items are pre-existing on main and unrelated to this PR:

• test_io.py::test_convert_dandelion fails with ImportError: Please install dandelion (sc-dandelion not available in the test env).
• test_plotting.py fails at collection with duplicate parametrization of 'adata_clonotype'.

Both also fail on the latest main run and on other open PRs (e.g. #711) with the same messages, so they are not introduced here. Happy to rebase once CI on main is green again.

[grst]

Sorry for being slow with reviews, I still need to take a closer look. I'll also take care of CI.

One thing I was wondering at first glance: Do you think it makes sense to also provide a corresponding plotting function (scirpy.pl.hill_diversity_profile)? But if it's just a seaborn oneliner, I would also be fine with just adding that to the tutorial.

[review-notebook-app]

Check out this pull request on 

See visual diffs & provide feedback on Jupyter Notebooks.

---

Powered by ReviewNB

[KilianMaire]

Thanks! I went with the tutorial route for now, since the output of hill_diversity_profile is a tidy DataFrame (one row per q, one column per group) and plotting it is essentially a seaborn lineplot oneliner. A dedicated pl.hill_diversity_profile would mostly be a thin wrapper around that, so I would rather not add the maintenance surface unless you prefer it for consistency with pl.alpha_diversity. Happy to add it as a follow-up if you do.

I filled in the previously empty "Clonotype Diversity" section of the 5k BCR tutorial: it computes the coverage-standardized profile across patient status groups, plots it with seaborn, and shows the conversion to the classical indices via convert_hill_table. I ran the full tutorial up to that section locally to confirm the new cells execute and the outputs are committed.

The only open design point from my side is the dependency question in the description (hillrep as an optional dep vs vendoring the estimator), whenever you get to it. No rush.

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 78.47%. Comparing base (120dd83) to head (fb37ffe).
⚠️ Report is 1 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #714      +/-   ##
==========================================
+ Coverage   78.31%   78.47%   +0.16%     
==========================================
  Files          51       51              
  Lines        4607     4646      +39     
==========================================
+ Hits         3608     3646      +38     
- Misses        999     1000       +1     

Files with missing lines	Coverage Δ
src/scirpy/tl/__init__.py	100.00% <100.00%> (ø)
src/scirpy/tl/_diversity.py	94.31% <100.00%> (+4.52%)	⬆️

... and 1 file with indirect coverage changes

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[grst]

The only open design point from my side is the dependency question in the description (hillrep as an optional dep vs vendoring the estimator), whenever you get to it. No rush.

It of course also depends a bit on your commitment to maintain the package in the future! But given that hillrep adds no additional transitive dependencies and doesn't have compiled code, I'm fine with adding it as a dependency.

[KilianMaire]

Great, thanks. And yes, I am committed to maintaining hillrep: it is my package, under active development, and the estimators are validated against R iNEXT to a relative tolerance of 1e-6 with the golden values committed, so regressions are caught.

Two follow-ups I just pushed:

• The tutorial-execution job was failing because the docs environment did not have hillrep (it only pulls the doc group, not test). Fixed by adding scirpy[diversity] to the doc dependency group, so the new tutorial section runs on CI.
• Added tests for the remaining convert_hill_table branches (the evenness modes and the missing-order error) and for the missing-hillrep import error. These need no optional dependency, so they run in every environment.

Let me know if you would like anything changed in the tutorial section or the API.

[KilianMaire]

I merged current main into the branch, so the PR is mergeable again. Two notes on CI:

hatch-test.py3.14-pre is green again. The earlier BrokenProcessPool errors were transient pre-release breakage and are gone after the merge.

Test tutorials is now red, and I don't think it comes from this PR. Both tutorial_5k_bcr and tutorial_3k_tcr fail with the same error, raised inside muon:

AttributeError: 'AnnData' object has no attribute '_X'
  muon/_core/preproc.py:758, in _filter_attr -> if data._X is not None:

anndata 0.13.0 (released 2026-07-07) removed the private _X attribute (I checked: hasattr(AnnData(...), "_X") is True on 0.12.19, False on 0.13.2), and muon 0.1.7 requires anndata with no upper bound, so the docs env resolves 0.13.2.

Things that point away from this PR:

• tutorial_3k_tcr is not touched by this branch and fails identically.
• Compiling the doc group with and without my scirpy[diversity] line resolves the same anndata==0.13.2 / muon==0.1.7, so the diversity extra isn't moving the resolution. Neither hillrep nor scikit-bio constrains anndata.

This should affect main as well — the last tutorials run there is from 2026-07-01, before anndata 0.13 was released, so it still looks green.

I've left it alone since it's your CI and outside the scope of this PR, but happy to add a pin here if you'd prefer that.

[grst]

I again started looking at this. I simplified the implementation a bit, I don't think the additional abstraction into _coverage_hill_profile is necessary. I also make hillrep a hard dependency, since it adds no additional transitive dependencies. I also think that your original output (a 'tidy' table) makes more sense than the pivoted dataframe, and that it is useful for users to have access to the assessment object.

I have mixed feelings about the convert_hill_table function (I know the design is not yours, but taken from #535). It's certainly useful to convert to the standard indices (Shannon etc.), but I don't like that it outputs different things based on the convert_to argument.

What do you think about the evenness_factor and relative_evenness metrics? Are they useful? If yes, they could easily become additional columns in the "tidy" dataframe. In that case, do you already want to add them on the hillrep side?

So much on the implementation side.

---

On the results side, I'm confused:

• The results for Mild, Ciritcal, Moderate, Severe are constant across all q
• Mild has a much higher Hill number than the other groups. There is absolutely no biological justification for this. I would expect them all to be similar, with maybe Asymptomatic, or Severe/Critical differing from the rest.

If we look at the normalized shannon entropy (which should be somewhat similar to the depth-adjusted hill number at q = 1, it's basically the same for all groups, just moderate is a little lower.

Do you have any insights on this?
The sample size is not large, but >450 cells per group is also not nothing, and small, real datasets could be in the same ballpark. With cells as the sampling unit, we are simply not in the same range as with reads in a bulk sequencing experiment.

>>> {k: len(x) for k, x in counts_by_group.items()}
{'Asymptomatic': 450,
 'Critical': 935,
 'Mild': 1147,
 'Moderate': 1039,
 'Severe': 1203}

[KilianMaire]

Thanks for picking this up again, and for the simplification, I agree with essentially all of it. Taking the implementation points first, then the results, because the results question turned out to be the interesting one and the answer changed my mind about what assess should do.

Implementation

Tidy output and exposing the Assessment, both better than what I had. Dropping _coverage_hill_profile is right; the indirection wasn't earning anything.

convert_hill_table - I share your discomfort, and there's now a concrete problem: it still does diversity_profile.loc[0] / .loc[1] / .loc[2], which assumed the pivoted frame. With hill_diversity_profile returning tidy it will raise (or silently index a row label) for every caller. So it has to change anyway, which is a good moment to drop the convert_to switch. My preference: keep "diversity" as its own small function, and make evenness columns rather than a conversion.

evenness_factor / relative_evenness - yes, useful, and yes, they belong on the hillrep side; there's no evenness there at all today. They're also a good guard rail: on the data below they'd sit at ~0.99 for every group, which is itself the diagnostic that something is wrong. One caveat before I add them, qD/⁰D and ln qD/ln ⁰D are the Daly 2018 forms, but the normalised evenness profiles in Chao & Ricotta 2019 (10.1002/ecy.2852) are the coherent counterpart to Hill numbers and are what iNEXT.4steps reports. I'd rather implement those and expose the simple ratios alongside. Happy to go either way if you have a preference.

The results

Your instinct is right, but not for the reason I first assumed, and one part of your reading needs correcting.

First, those are clonotypes, not cells. len(x) is the length of the abundance vector. The actual counts:

group	cells (n)	clonotypes (S_obs)	f1	f2	largest clone
Asymptomatic	456	450	444	6	2
Critical	949	941	937	2	5
Mild	1152	1147	1143	3	3
Moderate	1155	1096	1069	18	11
Severe	1288	1241	1221	15	17

450 clonotypes for 456 cells. There is essentially no clonal expansion in this 5k subsample. Sample coverage is 0.008–0.074, and f0_hat estimates 16 000–219 000 undetected clonotypes.

Why the profile is flat. Cmax = min_i Ĉ(2n_i) = 0.0130. Standardising to 1.3 % coverage gives wildly different m:

group	m	method	qD(q=0)	qD(q=2)	qD(q=0)/m
Moderate	57.7	Rarefaction	57.4	56.9	0.993
Severe	73.0	Rarefaction	72.5	71.8	0.993
Asymptomatic	225.0	Rarefaction	223.6	222.2	0.994
Critical	1031.5	Extrapolation	1022.5	997.3	0.991
Mild	2304.0	Extrapolation	2287.0	2257.0	0.993

At those depths every cell drawn is a distinct clonotype, so qD ≈ m for every q. The curve can't decrease - there's no abundance structure left to measure. What's plotted is the standardisation depth, not diversity.

Why Mild is on top, this is where I was wrong. I assumed noise. The direction is actually defensible: Mild is the least clonally expanded group (0.78 % of cells in clones ≥ 2; asymptotic 1/λ̂ = 110 496) and Moderate the most (7.45 %; 1/λ̂ = 3 830). Ranked by expansion, Mild really is the most diverse. Observed ²D/⁰D says the same thing: 0.994 for Mild against 0.810 / 0.794 for Moderate / Severe. So a reviewer could defend the ordering.

What is not defensible is the magnitude, because it isn't a property of the repertoires. Ratio of Mild to Moderate at q=0, as a function of the target coverage, across the whole admissible range:

target C	0.002	0.005	0.0130	0.030	0.050	0.0745
Mild / Moderate	29.8	31.9	39.9	35.3	23.7	14.4

The "40×" is a property of where Cmax happened to land. And Cmax here is decided by three doubletons in Mild — merging two of Mild's singletons into one doubleton moves Cmax from 0.01302 to 0.01648 and Critical's m from 1032 to 1812. Chao & Jost's replication principle requires the standardised ratio to be roughly independent of the standardisation level; here it varies by a factor of 2.8 inside the admissible window.

For completeness I also checked the other rule: Chao et al. 2014 Box 1 recommends Cbase = max(max_i Ĉ(n_i), min_i Ĉ(2n_i)) = 0.0745, not the min_i Ĉ(2n_i) that iNEXT::estimateD (and therefore hillrep) uses. That gives 14× instead of 40× - but requires extrapolating Critical and Mild to 16× and 14× their size. No choice of standardisation rule rescues this dataset.

What the honest comparison looks like. Rarefying to equal size (m = 456) instead of equal coverage:

order_q   Asymptomatic  Critical   Mild  Moderate  Severe
0.0              450.0     453.5  455.1     441.3   445.4
1.0              447.8     452.1  454.8     431.2   436.6
2.0              444.3     449.2  454.1     407.6   411.3

All but identical, Moderate/Severe slightly lower, exactly your normalised Shannon entropy (0.9992 / 0.9991 / 0.9997 / 0.9929 / 0.9940). The two methods agree; it's coverage standardisation that manufactures the discrepancy.

This is not an implementation issue. I reproduced the whole thing in R with iNEXT 3.0.2 via estimateD(base="coverage"): identical Cmax (0.013016177898) and identical qD to five decimals for all five groups.

Two upstream causes, and the second one is the real one

1. The clonotype definition. Switching define_clonotypes from receptor_arms="all" to "any" on the same cells moves Cmax from 0.013 to 0.243, makes the profile decrease in q again, and reverses the ranking (Mild goes from first to last). I am not suggesting "any" is the right call, merging on a single chain manufactures expansion. The point is that this diagnostic is dominated by a decision taken upstream of any diversity estimator. In 10x scTCR a substantial fraction of cells recover only one chain; under "all" those become distinct clonotypes and inflate singletons mechanically.

2. The groups pool five donors each. This is the one I'd emphasise. Each severity group aggregates 5 patients, and 0 of 4875 clonotypes is shared between any two patients. So each "assemblage" is a union of five disjoint communities, which is not what Chao's estimators assume - they assume one community sampled multinomially. Pooling structurally inflates f1 and suppresses f2: a clonotype private to one donor cannot become a doubleton by aggregation. It shows in the per-patient breakdown, 8 of the 25 patients have f2 = 0 outright, and the group-level doubleton counts are concentrated in single patients (11 of Moderate's 18 come from one patient, 8 of Severe's 15 from another). The entire comparison rests on whether one patient per group happened to show some expansion.

So the sample size isn't really 450–1288. It's five samples of 30–437 cells each, from five disjoint repertoires.

What I changed in hillrep

You've found a real gap, and it's in the part I claimed as the package's main value. assess returned caution here, which is far too soft: at 1 % coverage there is nothing to be cautious about, the comparison simply isn't defined. I've implemented:

• depth_ratio = standardised richness / standardisation depth. When it approaches 1, the returned diversity is the depth. It separates the regimes cleanly and needs no threshold folklore: this data sits at 0.991–0.994, the 12 bulk TCR-β repertoires I test against sit at 0.940–0.947, and the receptor_arms="any" version at 0.837–0.861.
• coverage_gain_2n = 1 - exp(-2 f2/f1), the largest relative reduction of the coverage deficit a doubled sample can buy (Chao & Jost 2012, p. 2540). Here it is 0.4–3.3 %.
• A depth_dominated flag yielding not_comparable, and a recommendation that no longer says "sequence deeper" — because it wouldn't help, but points at the clonotype definition and at grouping that pools individuals whose classes cannot co-occur.

On your data it now returns:

verdict: NOT_COMPARABLE
target coverage (Cmax): 0.0130
reasons:
  - Asymptomatic, Critical, Mild, Moderate, Severe: at Cmax=0.0130 the standardized
    richness is 99.4% of the standardization depth, so the reported diversity is the
    standardization depth rather than a property of the assemblage; the profile is
    flat in q and the between-sample ranking is a ranking of m.

Also fixed along the way: a silent NaN when the asymptotic estimator is infinite (reachable at q ≥ 3 on repertoires whose largest clone has two members, including this dataset), and the completeness column is now documented as an upper bound, since Chao1 estimates a lower bound on the undetected-class count. And one thing worth knowing for the x2 extrapolation flag: Chao et al. 2014 state that rule for q = 0 only, and say extrapolation to the asymptote is statistically safe for q ≥ 1. hillrep applies the conservative q = 0 rule to all orders; that's why Mild was flagged. I've documented it rather than changed it, but it's arguable.

This is on a branch, not released yet, it'll go out as 0.4.0 before this merges, so the pinned dependency can require it.

Suggestion for the docs

I think this strengthens rather than weakens the case for exposing the Assessment object, which was your call. But the tutorial dataset can't carry the method: at 5k cells pooled across donors it is squarely in the regime where the honest answer is "not with this data". Three options, in my order of preference:

1. Use this dataset as the worked example of the assessment saying no, and demonstrate compare(level="size") as what to do instead. It's an unusual thing for a package to lead with and it's the most useful thing here.
2. Group per sample/patient rather than pooling by severity, so each assemblage is one repertoire.
3. Use a dataset with real clonal expansion (TILs, antigen-sorted, chronic infection), where f2/f1 is substantial and the framework is in its working range.

Happy to write whichever you prefer.