hillrep 0.2.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

From a clone of this repository (a PyPI release is planned):

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

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

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

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=...)

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 five 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, 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.

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.
• 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, 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