hillrep

Coverage-based Hill-number diversity for immune repertoires (and any abundance data), in Python.

hillrep brings the iNEXT estimation framework (Chao et al. 2014; Hsieh, Ma & Chao 2016) to Python with first-class support for AIRR-seq data. It computes Hill-number diversity profiles with sample-coverage-based rarefaction and extrapolation and bootstrap confidence intervals, so you can compare the clonal diversity of repertoires sequenced to different depths without the bias that wrecks naive index comparisons.

The problem it solves

Richness, Shannon, and Simpson indices are all biased by sequencing depth: a repertoire sequenced deeper looks more diverse simply because more rare clones were observed. Comparing raw indices across samples of unequal depth is one of the most common statistical mistakes in repertoire analysis.

The correct fix is to standardize all samples to a common sample coverage (the fraction of the assemblage represented in the sample) and report Hill numbers, the effective number of equally-abundant clones. That machinery has lived almost entirely in R (iNEXT, immunarch, alakazam). In Python you previously had only fragments: scikit-bio has a point Hill estimator and Chao1, pyrepseq has Chao1/overlap, but neither does coverage-based rarefaction/extrapolation with confidence intervals, and neither is AIRR-aware.

Install

pip install hillrep                 # core
pip install "hillrep[airr,plot]"    # + AIRR schema reader and plotting

Or from a clone of this repository:

pip install -e ".[airr,plot]"

Quickstart

from hillrep import AbundanceCounts, estimate, compare

# an abundance vector is the per-clonotype read/UMI count
counts = AbundanceCounts([95, 40, 21, 13, 8, 5, 3, 3, 2, 2, 1, 1, 1, 1])

# size-based rarefaction/extrapolation curve with 95% CIs (q = 0, 1, 2)
curve = estimate(counts, q=(0, 1, 2))

# compare several repertoires at a common coverage (the fair comparison)
result = compare(
    {"patient_A": counts_a, "patient_B": counts_b},
    level="coverage",   # standardize to the largest coverage common to all
)

compare standardizes every assemblage to a shared coverage and returns one diversity estimate per (assemblage, q) with a confidence interval:

assemblage       m        method  order_q     qD  qD_lcl  qD_ucl  coverage
 patient_A 365.618 Extrapolation        0 20.367   5.299  35.435     0.987
 patient_A 365.618 Extrapolation        1  5.853   4.768   6.938     0.987
 patient_B 110.000 Extrapolation        0  8.248   3.878  12.617     0.987
 patient_B 110.000 Extrapolation        1  3.970   2.925   5.015     0.987

Is the comparison even valid? assess

Before reporting any diversity number, ask whether a fair comparison is possible at all. assess is the one thing hillrep knows that a naive index never tells you, turned into a verdict. It reports, per sample, what iNEXT's DataInfo shows (size, richness, singletons/doubletons, sample coverage) plus what it does not: the extrapolation factor each sample must undergo to reach the common standardization target Cmax, how much coverage a doubled sample could actually buy (coverage_gain_2n), and depth_ratio — the standardized richness divided by the standardization depth.

That last one is the sharpest signal. When a sample holds almost no repeated classes, every individual drawn at Cmax is a distinct class, so the standardized Hill number equals m for every order q: the profile is flat and the ranking between samples is a ranking of sampling depth, not of diversity. depth_ratio near 1 says exactly that, and it is the primary trigger for not_comparable.

from hillrep import assess

report = assess({"blood": counts_blood, "tumor": counts_tumor})
print(report.summary())
report.table   # per sample: coverage, coverage_gain_2n, extrap_factor, depth_ratio, flag

Comparison assessment over 2 assemblage(s)
  verdict: NOT_COMPARABLE
  target coverage (Cmax): 0.9707
  recommendation: A fair comparison is not supported as-is: sequence the
    degenerate sample(s) deeper, or drop them. Standardizing anyway rests on an
    unstable undetected-class correction.
  reasons:
    - tumor: too few repeated classes for a stable undetected-class estimate
      (no doubletons or a single class).
  flagged samples:
    - tumor: degenerate (coverage 0.962, extrap x2.00, qD0/m 0.048, gain at 2n 0.0%)

The verdict is reliable / caution / not_comparable, driven by explicit, documented, tunable thresholds. Only one of the quantities behind them comes from a published cutoff — coverage_gain_2n = 1 - exp(-2 f2/f1), the doubling bound of Chao & Jost 2012 — and the module docstring says plainly which numbers are hillrep heuristics and how they were calibrated. Note also that the x2 extrapolation rule of Chao et al. 2014 is stated there for q = 0 only; hillrep applies it to every order, which is conservative. Every number composes the same kernels validated against iNEXT; the value added is the decision layer. On the command line: hillrep assess rearrangements.tsv (verdict to stderr, the table to stdout so it can be piped).

AIRR-seq

Point it at an AIRR Rearrangement TSV (the format from IgBLAST, MiXCR, 10x conversions). hillrep defines a clonotype, groups by repertoire_id, and estimates diversity per repertoire.

from hillrep.airr import read_rearrangement, clonotype_counts
from hillrep import compare

df = read_rearrangement("rearrangements.tsv")
# clonotype = (junction_aa, V-gene, J-gene); abundance = duplicate_count
reps = clonotype_counts(df, by="repertoire_id")
result = compare(reps, level="coverage")

The clonotype definition and abundance weighting are explicit choices you can override (clone_key=, abundance=), not hidden defaults. Gene/allele calls are reduced to gene level (IGHV3-23*01 becomes IGHV3-23), and ambiguous comma-separated calls keep their first element.

Command line

hillrep estimate rearrangements.tsv --by repertoire_id --q 0 1 2
hillrep compare  rearrangements.tsv --level coverage --format json
hillrep estimate counts.txt --q 0          # one abundance per line

Single-cell (AnnData / scirpy)

If your TCR/BCR data lives in an AnnData (the scanpy/scirpy stack), point hillrep straight at it. scirpy.tl.alpha_diversity only offers depth-biased indices; from_anndata brings the coverage-standardized estimators to the same object.

from hillrep import from_anndata, compare

# the single-cell sampling unit is the cell: a clonotype's abundance is its cell count
reps = from_anndata(adata, groupby="sample")   # reads adata.obs["clone_id"]
compare(reps, level="coverage")                # fair across samples of unequal size

Only adata.obs is read, so this works across scirpy versions and needs no hard anndata dependency (install the extra with pip install "hillrep[anndata]" if you do not already have it). The clonotype column is auto-resolved (clone_id, then the legacy clonotype); override with target_col=. Each cell counts as one unless you pass a numeric abundance= column (e.g. UMI counts). For overlap between single-cell samples, clonotype_matrix_from_anndata(adata, groupby="sample") builds the clonotype-by-sample table that overlap_matrix consumes.

Diversity profiles and overlap

The whole Hill profile (qD versus q) is more informative than any single order:

from hillrep import hill_profile
prof = hill_profile(reps, at="coverage", coverage=0.9)   # tidy qD-vs-q, comparable

Repertoire overlap (beta diversity) answers "how similar are two repertoires?" It needs the shared clonotypes, so it works on a clonotype-by-repertoire table:

from hillrep.airr import clonotype_matrix
from hillrep import overlap_matrix

mat = clonotype_matrix(df, by="repertoire_id")     # clonotypes x repertoires
overlap_matrix(mat, method="morisita-horn")        # symmetric similarity matrix

Methods: morisita-horn, bray-curtis, jaccard, sorensen (matched to vegan), and the Chao bias-corrected chao-sorensen / chao-jaccard (matched to fossil), which estimate the similarity of the complete assemblages, correcting for shared clones missed by undersampling.

Plotting

import matplotlib.pyplot as plt
from hillrep import estimate
from hillrep.plotting import plot_rarefaction

curve = estimate(reps, q=0)
plot_rarefaction(curve, order_q=0)   # solid = rarefaction, dashed = extrapolation, band = CI
plt.show()

How this maps to iNEXT

iNEXT	hillrep
iNEXT(x, q, datatype="abundance")$iNextEst$size_based	estimate(x, q)
estimateD(x, base="coverage")	compare(x, level="coverage")
estimateD(x, base="size")	compare(x, level="size")
ChaoRichness / ChaoShannon / ChaoSimpson	asymptotic_hill(x, q)
Chat.Ind (sample coverage)	sample_coverage(x, m)
diversity profile (qD vs q)	hill_profile(x, at=...)
DataInfo (data summary) + a comparison-reliability verdict	assess(x)

Every estimator is a direct port of the corresponding iNEXT kernel and is unit-tested against iNEXT 3.0.2 output to a relative tolerance of 1e-6 on the deterministic point estimates: observed/asymptotic richness, Shannon and Simpson; the rarefaction, extrapolation (tested over the 1 to 2n range) and coverage curves; and the estimateD coverage- and size-standardized comparison. The ground-truth values are generated by scripts/gen_golden.R and committed under tests/golden/, so the test suite runs without R.

In addition, tests/test_differential_inext.py is a differential test that generates a dozen fresh random assemblages from several clone-size distributions (n up to several thousand) and checks hillrep against R iNEXT on each. It requires R + iNEXT and is skipped when they are absent (so it runs locally, not in CI); agreement is asserted to relative 1e-5 on the deterministic estimators.

Bootstrap confidence intervals use the same construction as iNEXT but a different random stream, so they are not bit-identical; only their width is checked, to agree with iNEXT's within a factor. They are implemented and exercised, not validated to the 1e-6 tolerance the deterministic kernels are.

Validation report

docs/hillrep-validation.pdf is a reproducible report (python scripts/make_report.py) that validates hillrep against iNEXT and walks through six use cases: depth-bias correction, fair multi-sample comparison across unequal sequencing depths, an AIRR Rearrangement pipeline, real public TCR-beta repertoires from the AIRR Data Commons, the comparison-reliability assessment (assess) on that real data, and robustness on extreme inputs. Highlights:

• Reproduces iNEXT 3.0.2 to a maximum relative error of ~3e-11 on the canonical ecology data.
• Naive richness varies by ~50% across sequencing depths of the same repertoire; the coverage-standardized estimate varies by ~10%.
• The same repertoire sequenced deeper looks ~50% more diverse by a naive index (a pure artifact); coverage standardization removes the confound.
• On 12 real public TCR-beta repertoires (AIRR Data Commons, depths 1k-19k), naive richness tracks sequencing depth almost perfectly; hillrep standardizes to a common coverage (compressing the spread) and flags that the repertoires are coverage-limited (5-16%), so even the standardized estimate is uncertain there. The data are fetched by scripts/fetch_real_airr.py and cached under data/ with full provenance.
• On those same real repertoires, assess returns the verdict caution and flags all 12 as coverage-limited, rather than reporting them as if directly comparable. Their depth_ratio is 0.94: undersampled, but repeated clonotypes still carry information. Single-cell repertoires pooled across donors reach 0.99 — every cell a distinct clonotype — and assess returns not_comparable there instead. The two regimes look identical to a naive index; that separation is the point.

Scope and honest limitations

What is implemented and verified:

• Hill numbers of any real order q >= 0 (0 = richness, 1 = exp-Shannon, 2 = inverse-Simpson, and everything in between or beyond), validated against iNEXT across a fine q-grid including non-integer orders.
• Continuous Hill diversity profiles (hill_profile, qD versus q) at the observed, a fixed size, or a common coverage.
• A comparison-reliability assessment (assess): a per-sample coverage table plus a reliable / caution / not_comparable verdict with documented thresholds, so you know whether a fair comparison is supported before you report one.
• Size-based and coverage-based rarefaction and extrapolation.
• Bootstrap confidence intervals (normal approximation, the iNEXT construction).
• Pairwise repertoire overlap: Morisita-Horn, Bray-Curtis, Jaccard, Sorensen (matched to vegan) and the Chao bias-corrected estimators (matched to fossil).
• AIRR Rearrangement ingestion, single-cell AnnData ingestion (scanpy/scirpy, via from_anndata), a CLI, and plotting helpers.

What is not in this version:

• Incidence (presence/absence) data; phylogenetic or functional Hill numbers.
• Full multi-assemblage Hill-number beta diversity with coverage standardization (the iNEXT.beta3D framework); overlap here is the classical pairwise indices plus the Chao bias-corrected estimators.

Caveats worth stating plainly:

• Extrapolation beyond roughly 2-3x the observed sample size becomes unreliable; the confidence intervals widen accordingly but the point estimate should be treated with caution. The committed golden tests cover extrapolation up to 2x n.
• The bootstrap intervals are asymptotic; with very few singletons/doubletons the undetected-class estimates (and therefore the intervals) are unstable. This is a property of the method, not of the implementation.
• Degenerate inputs behave like iNEXT: an assemblage with no repeated clones (all singletons) has infinite asymptotic inverse-Simpson diversity, and its q=2 extrapolation grows with the target size rather than converging. These are properties of the estimator on data that carries no abundance information.
• compare(level="coverage") floors any assemblage that cannot reach the target coverage to its closest attainable sample size and emits a warning; always read the coverage column, not just qD.
• Rarefaction for q=1 costs roughly O(largest clone abundance) per evaluated point. The kernel is vectorized (a single point on a million-read repertoire takes a fraction of a second), but a full 40-knot bootstrapped curve on a very deep repertoire with a dominant clone can still take a while; reduce n_points or n_boot if needed.

Design and architecture

The decisions below are the substance of hillrep; the equations come from the iNEXT papers, the engineering and the abstractions are what this package adds.

A port, not an R wrapper. Wrapping iNEXT through rpy2 would have been less code, but it forces an R installation on every user and keeps the data in R's world. The Python single-cell stack (scanpy, scirpy, the scverse) holds repertoire data in AnnData and pandas; an rpy2 bridge would mean serializing those across the language boundary on every call. A native port is pip-installable, has no R runtime dependency, and lets the estimators read the data structures Python users already hold. The cost (re-deriving and re-validating every kernel) is paid once; the validation suite (below) is how we keep it honest.

One numeric core, several adapters. Everything funnels through a single validated type, AbundanceCounts (a checked integer abundance vector with memoized frequency-of-frequencies). The estimators know only that type. The AIRR reader (hillrep.airr) and the AnnData reader (hillrep.anndata) are thin adapters that each collapse their own input into AbundanceCounts, so a new data source is a new adapter, never a change to the math. This is why single-cell support was ~150 lines: the cell-as-sampling-unit decision lives entirely in the adapter.

Numerical stability over literal transcription. Several published iNEXT expressions are numerically dangerous as written. The asymptotic Shannon term, taken literally, suffers catastrophic cancellation in a bracket and overflow in a (1-A)^-(n-1) prefactor at large n; hillrep rewrites it as a bounded integral with a closed form that is finite for every n and A. The general-order correction is re-indexed so the overflowing prefactor cancels and its binomial is advanced by recurrence (no gamma function at a large argument). These rewrites are invisible in the API and are exactly what a naive reimplementation gets wrong; the differential test against live R iNEXT is what caught them.

Validate against the reference, then cut the dependency. Every deterministic estimator is checked against R iNEXT 3.0.2 to 1e-6. The reference values are generated by committed R scripts and stored as golden JSON, so the test suite runs without R (fast CI, no R toolchain for contributors) while the ground truth is still R. A separate, R-gated differential test exercises fresh random assemblages to catch regressions the fixed goldens would miss.

Two questions, two layers. "What is the diversity?" (estimate, compare, hill_profile) and "can these even be compared?" (assess) are deliberately separate. The estimators stay pure numeric functions; assess is a decision layer on top that composes them into a verdict. Keeping the verdict out of the estimators means the thresholds are explicit and tunable, and the numeric kernels stay free of policy.

A scope boundary, on purpose. hillrep tracks the validated iNEXT kernels rather than inventing estimators. Coverage-standardized multi-assemblage beta diversity (the iNEXT.beta3D framework), incidence data, and phylogenetic diversity are left out rather than shipped half-validated, because the whole value proposition is "trust this to standardize correctly". The honest-limitations section above states where the method itself (not the implementation) runs out of road.

References

• Chao, A., Gotelli, N. J., Hsieh, T. C., et al. (2014). Rarefaction and extrapolation with Hill numbers. Ecological Monographs 84(1), 45-67.
• Hsieh, T. C., Ma, K. H., & Chao, A. (2016). iNEXT: an R package for rarefaction and extrapolation of species diversity. Methods in Ecology and Evolution 7, 1451-1456.
• Chao, A., & Jost, L. (2012). Coverage-based rarefaction and extrapolation. Ecology 93(12), 2533-2547.
• Chao, A., Wang, Y. T., & Jost, L. (2013). Entropy and the species accumulation curve. Methods in Ecology and Evolution 4(11), 1091-1100. (the asymptotic Shannon estimator)

hillrep is an independent reimplementation and is not affiliated with the iNEXT authors. If you use it, please also cite the papers above.

Citation

If you use hillrep, please cite the archived software (and the iNEXT papers above). The all-versions DOI is 10.5281/zenodo.20670550; the badge above always resolves to the latest release. Citation metadata is in CITATION.cff.

License

MIT