Metadata-Version: 2.4
Name: hillrep
Version: 0.2.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: 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

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):

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

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

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

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

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