Coverage-standardized Hill diversity (tl.hill_diversity_profile)

CHANGELOG.md

@@ -10,6 +10,15 @@ and this project adheres to [Semantic Versioning][].
 
 ## Unreleased
 
+### 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)).
+
 ### Performance improvements
 
   - Speed up identity distance metric computation for comparisons between two different sequence arrays.

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

@@ -43,7 +43,7 @@ dependencies = [
 ]
 optional-dependencies.cupy = [ "cupy-cuda12x" ]
 optional-dependencies.dandelion = [ "sc-dandelion>=0.3.5" ]
-optional-dependencies.diversity = [ "scikit-bio>=0.5.7" ]
+optional-dependencies.diversity = [ "hillrep>=0.3", "scikit-bio>=0.5.7" ]
 optional-dependencies.parasail = [
   # parasail 1.2.1 fails to be installd on MacOS
   "parasail!=1.2.1",
@@ -79,6 +79,7 @@ doc = [
   "nbconvert",
   "pycairo",
   "sc-dandelion>=0.5",
+  "scirpy[diversity]",               # hillrep for the Hill diversity tutorial section
   "sphinx>=7,<7.4",                  # https://github.com/sphinx-doc/sphinx/issues/12589
   "sphinx-autodoc-typehints",
   "sphinx-book-theme>=1",

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,5 +1,6 @@
+import warnings
 from collections.abc import Callable
-from typing import cast
+from typing import Literal, cast
 
 import numpy as np
 import pandas as pd
@@ -43,6 +44,175 @@ def _dxx(counts: np.ndarray, *, percentage: int):
     return i / len(freqs) * 100
 
 
+def _import_hillrep():
+    """Lazily import :mod:`hillrep`, with an actionable error if it is missing."""
+    try:
+        import hillrep
+    except ImportError:
+        raise ImportError(
+            "Coverage-based Hill diversity requires the `hillrep` package. "
+            "You can install it with `pip install hillrep`."
+        ) from None
+    return hillrep
+
+
+def _coverage_hill_profile(
+    counts_by_group: dict[str, np.ndarray],
+    q_values: list[float],
+) -> pd.DataFrame:
+    """Coverage-standardized Hill profile for a set of groups.
+
+    This is the only place that talks to the estimation backend. It standardizes
+    all groups to a common sample coverage (iNEXT's ``Cmax`` rule) and returns the
+    point estimates of the Hill numbers, and it warns when the groups cannot be
+    compared fairly at a shared coverage.
+
+    Returns a ``DataFrame`` indexed by diversity order ``q`` with one column per group.
+    """
+    hillrep = _import_hillrep()
+
+    assessment = hillrep.assess(counts_by_group)
+    if assessment.verdict != "reliable":
+        warnings.warn(
+            "The groups cannot be compared at a fully reliable common coverage. "
+            "Read the profile with caution.\n" + assessment.summary(),
+            stacklevel=3,
+        )
+
+    tidy = hillrep.compare(counts_by_group, level="coverage", q=q_values, n_boot=0)
+    profile = tidy.pivot(index="order_q", columns="assemblage", values="qD")
+    profile.index.name = None
+    profile.columns.name = None
+    # preserve the input group order rather than the alphabetical pivot order
+    return profile[list(counts_by_group)]
+
+
+@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,
+) -> pd.DataFrame:
+    """\
+    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 sequencing depth, so two samples
+    sequenced to different depth 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
+    `(Chao et al. 2014; Hsieh, Ma & Chao 2016) <https://doi.org/10.1890/13-0133.1>`__,
+    so the profiles are comparable across depth. The estimation is delegated to the
+    `hillrep <https://github.com/KilianMaire/hillrep>`__ package.
+
+    When the groups cannot be standardized to a fully reliable shared coverage (for
+    example because one group is heavily undersampled), a warning is emitted: sometimes
+    the honest answer is that a fair comparison is not possible.
+
+    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
+    -------
+    A `DataFrame` with one row per diversity order `q` and one column per group. The
+    output flows directly into :func:`~scirpy.tl.convert_hill_table` and into plotting
+    libraries such as seaborn.
+    """
+    params = DataHandler(adata, airr_mod)
+    ir_obs = params.get_obs([target_col, groupby])
+    ir_obs = ir_obs.loc[~_is_na(ir_obs[target_col]), :]
+    clono_counts = ir_obs.groupby([groupby, target_col], observed=True).size().reset_index(name="count")
+
+    counts_by_group = {
+        str(k): cast(
+            np.ndarray,
+            cast(pd.Series, clono_counts.loc[clono_counts[groupby] == k, "count"]).to_numpy(),
+        )
+        for k in sorted(ir_obs[groupby].dropna().unique())
+    }
+
+    q_values = [float(q) for q in np.arange(q_min, q_max + q_step, q_step)]
+    return _coverage_hill_profile(counts_by_group, q_values)
+
+
+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 `Daly et al. 2018 <https://doi.org/10.1093/bib/bbx019>`__ 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,