Metadata-Version: 2.4
Name: hillrep
Version: 0.4.0
Summary: Coverage-based Hill-number diversity estimation for immune repertoires (and any abundance data).
Project-URL: Homepage, https://github.com/KilianMaire/hillrep
Project-URL: Issues, https://github.com/KilianMaire/hillrep/issues
Author: Kilian Maire
License-Expression: MIT
License-File: LICENSE
Keywords: AIRR,BCR,TCR,diversity,extrapolation,hill-numbers,iNEXT,immune-repertoire,rarefaction,sample-coverage
Classifier: Intended Audience :: Science/Research
Classifier: License :: OSI Approved :: MIT License
Classifier: Programming Language :: Python :: 3
Classifier: Topic :: Scientific/Engineering :: Bio-Informatics
Classifier: Typing :: Typed
Requires-Python: >=3.10
Requires-Dist: numpy>=1.23
Requires-Dist: pandas>=1.5
Requires-Dist: scipy>=1.9
Provides-Extra: airr
Requires-Dist: airr>=1.4; extra == 'airr'
Provides-Extra: anndata
Requires-Dist: anndata>=0.9; extra == 'anndata'
Provides-Extra: dev
Requires-Dist: anndata>=0.9; extra == 'dev'
Requires-Dist: matplotlib>=3.6; extra == 'dev'
Requires-Dist: mypy>=1.8; extra == 'dev'
Requires-Dist: pytest>=7.4; extra == 'dev'
Requires-Dist: ruff>=0.4; extra == 'dev'
Provides-Extra: plot
Requires-Dist: matplotlib>=3.6; extra == 'plot'
Description-Content-Type: text/markdown

# hillrep

[![PyPI](https://img.shields.io/pypi/v/hillrep.svg)](https://pypi.org/project/hillrep/)
[![Python](https://img.shields.io/pypi/pyversions/hillrep.svg)](https://pypi.org/project/hillrep/)
[![License: MIT](https://img.shields.io/badge/License-MIT-yellow.svg)](LICENSE)
[![DOI](https://zenodo.org/badge/1267425013.svg)](https://zenodo.org/badge/latestdoi/1267425013)

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

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

Or from a clone of this repository:

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

## Quickstart

```python
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`.

```python
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.

```python
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

```bash
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.

```python
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:

```python
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:

```python
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

```python
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`](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](https://doi.org/10.5281/zenodo.20670550);
the badge above always resolves to the latest release. Citation metadata is in
[`CITATION.cff`](CITATION.cff).

## License

MIT
