Coverage-standardized Hill diversity (tl.hill_diversity_profile)

CHANGELOG.md

@@ -12,6 +12,12 @@ and this project adheres to [Semantic Versioning][].
 
 ### Additions
 
+  - Add `tl.hill_diversity_profile` and `tl.convert_hill_table` for coverage-standardized
+    Hill-number diversity. Profiles are standardized to a common sample coverage (iNEXT
+    framework) so they are comparable across samples of different sequencing depth, and a
+    warning is raised when a fair comparison is not supported. Estimation is delegated to
+    the [hillrep](https://github.com/KilianMaire/hillrep) package, installed via the
+    `diversity` extra ([#714](https://github.com/scverse/scirpy/pull/714)).
   - Add support for TCRBLOSUM alpha/beta substitution matrices in the `tcrdist` distance metric via
     `base_matrix="tcrblosum"`, and allow configuring the substitution-to-distance cap with `distance_cap`.
 

docs/api.rst

@@ -148,6 +148,8 @@ Analyse clonal diversity
    tl.clonal_expansion
    tl.summarize_clonal_expansion
    tl.alpha_diversity
+   tl.hill_diversity_profile
+   tl.convert_hill_table
    tl.repertoire_overlap
    tl.clonotype_modularity
    tl.clonotype_imbalance

pyproject.toml

@@ -22,6 +22,7 @@ dependencies = [
   "airr>=1.6.1",
   "anndata>=0.9",
   "awkward>=2.1",
+  "hillrep>=0.3",
   # 0.10.0 and 0.10.1 have the bug described in https://github.com/igraph/python-igraph/issues/570
   "igraph!=0.10,!=0.10.1",
   "joblib>=1.3.1",
@@ -74,7 +75,7 @@ doc = [
   "myst-nb>=1.1",
   "nbconvert",
   "pycairo",
-  "scirpy[dandelion]",
+  "scirpy[dandelion]",               # hillrep for the Hill diversity tutorial section
   "scverse-misc[sphinx]>=0.1.2",
   "sphinx>=8.1",
   "sphinx-autodoc-typehints",

src/scirpy/tests/test_tools.py

@@ -1093,3 +1093,113 @@ def test_mutational_load(adata_mutation, region_vars, expected):
 def test_mutational_load_adata_not_aligned(adata_not_aligned):
     with npt.assert_raises(ValueError):
         ir.tl.mutational_load(adata_not_aligned, germline_key="germline_alignment")
+
+
+def _adata_from_counts(counts: dict[str, list[int]], mudata: bool = False):
+    """Build a minimal AnnData/MuData whose obs encodes the given per-group clonotype counts.
+
+    ``counts`` maps a group label to a list of clonotype abundances. Each abundance
+    is expanded into that many cells sharing one clonotype id, so that collapsing
+    ``obs`` back by (group, clonotype) recovers the original counts exactly.
+    """
+    records = []
+    for group, abundances in counts.items():
+        for clone_idx, n in enumerate(abundances):
+            for _ in range(n):
+                records.append([f"{group}_ct{clone_idx}", group])
+    obs = pd.DataFrame(records, columns=["clonotype_", "group"])
+    obs.index = [f"cell{i}" for i in range(len(obs))]
+    return _make_adata(obs, mudata)
+
+
+# two well-sampled groups of different depth: comparable at a common coverage < 1
+_HILL_COUNTS = {
+    "A": [30, 20, 14, 10, 7, 5, 4, 3, 2, 2, 1, 1, 1, 1, 1],
+    "B": [90, 60, 42, 30, 22, 16, 12, 9, 7, 5, 4, 3, 2, 2, 1, 1, 1, 1],
+}
+
+
+@pytest.mark.extra
+@pytest.mark.parametrize("mudata", [False, True], ids=["AnnData", "MuData"])
+def test_hill_diversity_profile(mudata):
+    import warnings
+
+    import hillrep
+
+    adata = _adata_from_counts(_HILL_COUNTS, mudata)
+
+    with warnings.catch_warnings():
+        warnings.simplefilter("error")  # the happy path must not warn
+        res = ir.tl.hill_diversity_profile(adata, groupby="group", target_col="clonotype_", q_min=0, q_max=2, q_step=1)
+
+    # the adapter must reproduce hillrep's coverage-standardized estimate exactly
+    ref = hillrep.compare(_HILL_COUNTS, level="coverage", q=[0.0, 1.0, 2.0], n_boot=0)
+    expected = ref.pivot(index="order_q", columns="assemblage", values="qD")
+
+    assert list(res.columns) == ["A", "B"]
+    for q in (0.0, 1.0, 2.0):
+        for grp in ("A", "B"):
+            npt.assert_allclose(res.loc[q, grp], expected.loc[q, grp], rtol=1e-9)
+
+
+@pytest.mark.extra
+def test_hill_diversity_profile_warns_when_not_comparable():
+    # group "A" is heavily undersampled with no doubletons: a fair comparison is not supported
+    counts = {"A": [3, 1, 1, 1, 1], "B": [120, 60, 40, 30, 20, 12, 8, 5, 3, 2, 1, 1]}
+    adata = _adata_from_counts(counts)
+    with pytest.warns(UserWarning, match="cannot be compared"):
+        ir.tl.hill_diversity_profile(adata, groupby="group", target_col="clonotype_")
+
+
+@pytest.mark.extra
+def test_convert_hill_table():
+    adata = _adata_from_counts(_HILL_COUNTS)
+    profile = ir.tl.hill_diversity_profile(adata, groupby="group", target_col="clonotype_", q_min=0, q_max=2, q_step=1)
+
+    div = ir.tl.convert_hill_table(profile, convert_to="diversity")
+    assert list(div.index) == ["Observed richness", "Shannon entropy", "Inverse Simpson", "Gini-Simpson"]
+    npt.assert_allclose(div.loc["Observed richness"].to_numpy(dtype=float), profile.loc[0].to_numpy(dtype=float))
+    npt.assert_allclose(div.loc["Shannon entropy"].to_numpy(dtype=float), np.log(profile.loc[1].to_numpy(dtype=float)))
+    npt.assert_allclose(div.loc["Inverse Simpson"].to_numpy(dtype=float), profile.loc[2].to_numpy(dtype=float))
+    npt.assert_allclose(div.loc["Gini-Simpson"].to_numpy(dtype=float), 1 - 1 / profile.loc[2].to_numpy(dtype=float))
+
+    with pytest.raises(ValueError):
+        ir.tl.convert_hill_table(profile, convert_to="not_a_mode")
+
+
+def test_convert_hill_table_modes():
+    # operates on a plain DataFrame, so it needs no optional dependency
+    profile = pd.DataFrame({"A": [10.0, 6.0, 4.0], "B": [20.0, 12.0, 8.0]}, index=[0, 1, 2])
+
+    div = ir.tl.convert_hill_table(profile, convert_to="diversity")
+    npt.assert_allclose(div.loc["Shannon entropy"].to_numpy(dtype=float), np.log(profile.loc[1].to_numpy(dtype=float)))
+    npt.assert_allclose(div.loc["Gini-Simpson"].to_numpy(dtype=float), 1 - 1 / profile.loc[2].to_numpy(dtype=float))
+
+    ef = ir.tl.convert_hill_table(profile, convert_to="evenness_factor")
+    npt.assert_allclose(ef.to_numpy(dtype=float), (profile / profile.loc[0]).to_numpy(dtype=float))
+
+    rel = ir.tl.convert_hill_table(profile, convert_to="relative_evenness")
+    npt.assert_allclose(rel.to_numpy(dtype=float), (np.log(profile) / np.log(profile.loc[0])).to_numpy(dtype=float))
+
+    with pytest.raises(ValueError, match="Invalid"):
+        ir.tl.convert_hill_table(profile, convert_to="nope")
+
+    with pytest.raises(ValueError, match="missing diversity order"):
+        ir.tl.convert_hill_table(profile.drop(index=2), convert_to="diversity")
+
+
+def test_hill_diversity_profile_requires_hillrep(monkeypatch):
+    import builtins
+
+    from scirpy.tl import _diversity
+
+    real_import = builtins.__import__
+
+    def fake_import(name, *args, **kwargs):
+        if name == "hillrep":
+            raise ImportError("simulated missing hillrep")
+        return real_import(name, *args, **kwargs)
+
+    monkeypatch.setattr(builtins, "__import__", fake_import)
+    with pytest.raises(ImportError, match="pip install hillrep"):
+        _diversity._import_hillrep()

src/scirpy/tl/__init__.py

@@ -4,7 +4,7 @@
 from ._clonotype_modularity import clonotype_modularity
 from ._clonotypes import clonotype_network, clonotype_network_igraph, define_clonotype_clusters, define_clonotypes
 from ._convergence import clonotype_convergence
-from ._diversity import alpha_diversity
+from ._diversity import alpha_diversity, convert_hill_table, hill_diversity_profile
 from ._group_abundance import group_abundance
 from ._ir_query import ir_query, ir_query_annotate, ir_query_annotate_df
 from ._mutational_load import mutational_load

src/scirpy/tl/_diversity.py

@@ -1,8 +1,10 @@
 from collections.abc import Callable
-from typing import cast
+from typing import Literal, cast
 
+import hillrep
 import numpy as np
 import pandas as pd
+from scanpy import logging
 
 from scirpy.util import DataHandler, _is_na
 
@@ -43,6 +45,149 @@ def _dxx(counts: np.ndarray, *, percentage: int):
     return i / len(freqs) * 100
 
 
+@DataHandler.inject_param_docs()
+def hill_diversity_profile(
+    adata: DataHandler.TYPE,
+    groupby: str,
+    *,
+    target_col: str = "clone_id",
+    airr_mod: str = "airr",
+    q_min: float = 0,
+    q_max: float = 2,
+    q_step: float = 1,
+) -> tuple[pd.DataFrame, "hillrep.Assessment"]:
+    """\
+    Computes a coverage-standardized Hill diversity profile for a range of diversity orders (`q`).
+
+    Hill numbers unify the common alpha diversity indices into a single family indexed
+    by an order `q` (`q=0` is observed richness, `q=1` is the exponential of Shannon
+    entropy, `q=2` is the inverse Simpson index).
+
+    Naively plugging observed frequencies into the Hill formula yields values
+    that grow with cell count, so two samples with different cell count look different
+    even when the underlying repertoire is the same.
+    This function instead standardizes all groups to a common sample coverage
+    before reporting the profile, following the iNEXT framework
+    (:cite:`Chao.2014,Hsieh.2016`),
+    so the profiles are comparable across depth.
+
+    When the groups cannot be standardized to a fully reliable shared coverage (for
+    example because one group is heavily undersampled), a warning is emitted and the
+    reason is available from the returned assessment: sometimes the honest answer is
+    that a fair comparison is not possible.
+
+    This function relies on the `hillrep <https://github.com/KilianMaire/hillrep>`__ package for the estimation.
+
+    Parameters
+    ----------
+    {adata}
+    groupby
+        Column of `obs` by which the grouping will be performed.
+    target_col
+        Column containing the clonotype annotation.
+    {airr_mod}
+    q_min
+        Lowest (start) diversity order.
+    q_max
+        Highest (end) diversity order.
+    q_step
+        Step between consecutive diversity orders.
+
+    Returns
+    -------
+    profile
+        A tidy `DataFrame` with one row per group and diversity order, with the columns
+        `assemblage` (the group), `order_q`, `m` (the standardized sample size),
+        `method`, `qD` (the Hill number), `qD_lcl`/`qD_ucl` (its confidence interval)
+        and `coverage`. This format flows directly into plotting libraries such as
+        seaborn.
+    assessment
+        A `hillrep` `Assessment` object describing whether the groups can be compared
+        fairly. Its `verdict` is one of `"reliable"`, `"caution"` or `"not_comparable"`,
+        `table` holds the per-group diagnostics (coverage, singletons/doubletons,
+        extrapolation factor) and `summary()` renders a human-readable report.
+    """
+    params = DataHandler(adata, airr_mod)
+    ir_obs = params.get_obs([target_col, groupby]).dropna(axis=0, how="any")
+    clono_counts = ir_obs.groupby([groupby, target_col], observed=True).size().reset_index(name="count")
+
+    counts_by_group = {
+        str(k): cast(pd.Series, clono_counts.loc[clono_counts[groupby] == k, "count"]).to_numpy()
+        for k in sorted(ir_obs[groupby].unique())
+    }
+
+    q_values = np.arange(q_min, q_max + q_step, q_step, dtype=float).tolist()
+
+    assessment = hillrep.assess(counts_by_group)
+    if assessment.verdict != "reliable":
+        logging.warning(
+            "The groups cannot be compared at a fully reliable common coverage. "
+            "Read the profile with caution.\n" + assessment.summary(),
+        )
+
+    tidy = hillrep.compare(counts_by_group, level="coverage", q=q_values, n_boot=0)
+
+    return tidy, assessment
+
+
+def convert_hill_table(
+    diversity_profile: pd.DataFrame,
+    convert_to: Literal["diversity", "evenness_factor", "relative_evenness"] = "diversity",
+) -> pd.DataFrame:
+    """\
+    Converts a profile from :func:`~scirpy.tl.hill_diversity_profile` into other alpha diversity indices.
+
+    See :cite:`Daly.2018` for an overview of the indices and evenness measures.
+
+    Parameters
+    ----------
+    diversity_profile
+        A `DataFrame` produced by :func:`~scirpy.tl.hill_diversity_profile`. It must
+        contain the diversity orders `0`, `1` and `2` in its index.
+    convert_to
+        Which conversion to perform:
+
+        * `"diversity"` - the classical indices (observed richness, Shannon entropy,
+          inverse Simpson, Gini-Simpson) derived from the Hill numbers.
+        * `"evenness_factor"` - each Hill number divided by the observed richness.
+        * `"relative_evenness"` - the log of each Hill number over the log of the
+          observed richness.
+
+    Returns
+    -------
+    A `DataFrame` whose rows are the requested indices (or diversity orders) and whose
+    columns are the groups.
+    """
+    for required_q in (0, 1, 2):
+        if required_q not in diversity_profile.index:
+            raise ValueError(
+                f"The profile is missing diversity order q={required_q}. "
+                "`convert_hill_table` requires the orders 0, 1 and 2."
+            )
+
+    if convert_to == "diversity":
+        richness = diversity_profile.loc[0]
+        inverse_simpson = diversity_profile.loc[2]
+        return pd.DataFrame(
+            {
+                "Observed richness": richness,
+                "Shannon entropy": np.log(diversity_profile.loc[1]),
+                "Inverse Simpson": inverse_simpson,
+                "Gini-Simpson": 1 - 1 / inverse_simpson,
+            }
+        ).T
+
+    if convert_to == "evenness_factor":
+        return diversity_profile / diversity_profile.loc[0]
+
+    if convert_to == "relative_evenness":
+        return np.log(diversity_profile) / np.log(diversity_profile.loc[0])
+
+    raise ValueError(
+        f"Invalid `convert_to` value {convert_to!r}. Choose 'diversity', 'evenness_factor' or 'relative_evenness'."
+    )
+
+
 @DataHandler.inject_param_docs()
 def alpha_diversity(
     adata: DataHandler.TYPE,