hillrep 0.3.0

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

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 coverage still reachable by doubling the sample, and the extrapolation factor each sample must undergo to reach the common standardization target Cmax. A sample that must be extrapolated several-fold is one the whole comparison leans on; that is the reliability signal.

from hillrep import assess

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

Comparison assessment over 2 assemblage(s)
  verdict: NOT_COMPARABLE
  target coverage (Cmax): 0.4165
  recommendation: A fair comparison is not supported as-is: sequence the
    degenerate sample(s) deeper, or drop them.
  flagged samples:
    - tumor: degenerate (coverage 0.074, extrap x2.00)

The verdict is reliable / caution / not_comparable, driven by explicit, documented, tunable thresholds (coverage near 1 means depth bias is small; extrapolating beyond ~2-3x is unreliable, Chao et al. 2014). 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.

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.

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.

License

MIT