Differential Abundance

Which cell populations change in frequency between two conditions?

There are two ways to answer that, and this tutorial runs both:

1. Count cells per cluster. The classical answer, and the question most people actually arrive with — “are naive B cells rarer in old mice?” You annotate cells, count them per sample, and compare the proportions across conditions. Simple, interpretable, and entirely dependent on the annotation you chose.

2. Score every cell state (Kompot). kompot.da estimates the cell-state density of each condition over a continuous representation and reports a per-cell log fold change with calibrated significance. It says whether any cell states are more or less common in one condition — regardless of the annotation. No clusters required.

The two are not interchangeable, and Section 3 comes back to why: counting tests enrichment for an annotated group of cells, one proportion per cluster, while kompot.da asks, for every cell state in the data, whether it is more or less common under one condition than the other. The answers usually correlate, but they are different quantities — a cell type that is depleted overall can still contain sub-states that occur only in that condition, which the count cannot see once it aggregates to the type.

The method is described in the paper. For gene-level changes, see Getting Started (DE).

%matplotlib inline

import anndata as ad
import matplotlib.pyplot as plt
import numpy as np
import palantir
import pandas as pd
import scanpy as sc
from scipy import stats

import kompot

# Consistent, clean plotting style across the tutorial set
plt.rcParams["axes.spines.right"] = False
plt.rcParams["axes.spines.top"] = False
plt.rcParams["image.cmap"] = "Spectral_r"

DATA_PATH = "../data/murine_bone_marrow_aging.h5ad"
GROUPING_COLUMN = "Age"                       # condition column in adata.obs
CONDITIONS = ["Young", "Old"]                 # first entry is the reference
CELL_TYPE_COLUMN = "highres_celltype"         # fine annotation
SAMPLE_COLUMN = "Replicate"                   # biological replicate (mouse)
LIBRARY_COLUMN = "Sample"                     # one sorted library = compartment x mouse
COMPARTMENT_COLUMN = "Compartment"            # FACS gate: HSPC or Mature
DIMENSIONALITY_REDUCTION = "DM_EigenVectors"  # cell-state representation

MIN_CELLS = 50  # ignore cell types too rare to count reliably

Load Data

import os
from pathlib import Path

import requests
from tqdm.auto import tqdm

Path(DATA_PATH).parent.mkdir(parents=True, exist_ok=True)

if not os.path.exists(DATA_PATH):
    print("Downloading dataset (~2.4 GB) from Zenodo ...")
    url = "https://zenodo.org/records/15587768/files/murine_bone_marrow_aging.h5ad?download=1"
    response = requests.get(url, stream=True)
    total = int(response.headers.get("content-length", 0))
    with open(DATA_PATH, "wb") as file, tqdm(total=total, unit="B", unit_scale=True) as bar:
        for chunk in response.iter_content(chunk_size=1 << 20):
            file.write(chunk)
            bar.update(len(chunk))

adata = ad.read_h5ad(DATA_PATH)
palantir.utils.run_diffusion_maps(adata, pca_key="X_pca_harmony", n_components=40)
adata

AnnData object with n_obs × n_vars = 8090 × 16285
    obs: 'Compartment', 'Replicate', 'Age', 'Sample', 'Info', 'batch', 'doublet_score', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt', 'total_counts_hb', 'pct_counts_hb', 'S_score', 'G2M_score', 'phase', 'leiden', 'phenograph', 'highres_celltype', 'midres_celltype'
    var: 'gene_ids', 'feature_types', 'genome', 'mt', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts', 'hb', 'highly_variable', 'means', 'dispersions', 'dispersions_norm', 'highly_variable_nbatches', 'highly_variable_intersection'
    uns: 'Age_colors', 'Compartment_colors', 'DMEigenValues', 'DM_EigenValues', 'Info_colors', 'README', 'Replicate_colors', 'Sample_colors', 'batch_colors', 'draw_graph', 'highres_celltype_colors', 'hvg', 'leiden', 'leiden_colors', 'midres_celltype_colors', 'neighbors', 'pca', 'phase_colors', 'umap'
    obsm: 'AbCapture', 'DM_EigenVectors', 'HTO', 'X_draw_graph_fa', 'X_pca', 'X_pca_harmony', 'X_pca_noregression', 'X_umap'
    varm: 'PCs'
    layers: 'MAGIC_imputed_data', 'cc_counts', 'logged_counts', 'normalized_counts', 'raw_counts'
    obsp: 'DM_Kernel', 'DM_Similarity', 'connectivities', 'distances'

sc.pl.umap(adata, color=[CELL_TYPE_COLUMN, GROUPING_COLUMN], frameon=False, wspace=1.3)

1. Counting cells per cluster

The familiar analysis: count how many cells of each annotated type each animal contributes, turn the counts into within-sample proportions, and compare Young against Old across replicates.

One detail decides whether the numbers mean anything. Cells here were FACS-sorted into two compartments before sequencing (Compartment: HSPC and Mature), so each mouse yields two libraries and the HSPC : Mature ratio is set by the sorting gate, not by the animal. Proportions are therefore only comparable within a compartment. Every cell type in this dataset sits ≥91% inside one compartment, so we score each type against the libraries of its own compartment — three Young mice versus three Old mice.

obs = adata.obs[adata.obs[GROUPING_COLUMN].isin(CONDITIONS)]

# Which compartment does each cell type live in, and how many cells do we have?
by_compartment = pd.crosstab(obs[CELL_TYPE_COLUMN], obs[COMPARTMENT_COLUMN])
home = by_compartment.idxmax(axis=1)                     # majority compartment per type
keep = by_compartment.sum(axis=1) >= MIN_CELLS

# Composition of every library: what fraction of it is cell type X?
counts = pd.crosstab(obs[LIBRARY_COLUMN], obs[CELL_TYPE_COLUMN])
props = counts.div(counts.sum(axis=1), axis=0) * 100     # percent of library

library = obs.groupby(LIBRARY_COLUMN, observed=True)[
    [GROUPING_COLUMN, COMPARTMENT_COLUMN]
].agg(lambda s: s.iloc[0])

props.head()

highres_celltype	HSC	LMPP	GMP	Myelo P	Neutrophil	Mono P	cMoP	Monocyte	Non Classical Monocyte	cDC	...	Mature Naive B cell	Memory B cell	Plasma cell	NK	T cell	CD8 TEM	Naive CD8 T cell	CD4 T cell	Treg	ILC
Sample
HSPC_Old_1	23.310023	23.076923	6.293706	8.857809	0.0	8.857809	0.0	0.0	0.0	0.0	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
HSPC_Old_2	15.887850	26.635514	6.074766	8.411215	0.0	8.411215	0.0	0.0	0.0	0.0	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
HSPC_Old_3	16.222760	26.392252	7.021792	8.232446	0.0	6.779661	0.0	0.0	0.0	0.0	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
HSPC_Young_1	0.931677	22.360248	7.763975	8.695652	0.0	18.012422	0.0	0.0	0.0	0.0	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
HSPC_Young_2	0.266667	28.533333	6.933333	10.400000	0.0	7.466667	0.0	0.0	0.0	0.0	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0

5 rows × 31 columns

def count_test(cell_types, home, props, library):
    """Young vs. Old proportion per cell type, within that type's own compartment."""
    ref, alt = CONDITIONS
    rows = {}
    for ct in cell_types:
        same_gate = library[COMPARTMENT_COLUMN] == home[ct]
        y = props.loc[same_gate & (library[GROUPING_COLUMN] == ref), ct].values
        o = props.loc[same_gate & (library[GROUPING_COLUMN] == alt), ct].values
        rows[ct] = {
            f"{ref}_%": y.mean(),
            f"{alt}_%": o.mean(),
            "log2FC": np.log2((o.mean() + 1e-3) / (y.mean() + 1e-3)),
            "p": stats.ttest_ind(o, y, equal_var=False).pvalue,   # Welch, n = 3 vs 3
        }
    return pd.DataFrame(rows).T.sort_values("log2FC")


counted = count_test(by_compartment.index[keep], home, props, library)
counted.round(3)

	Young_%	Old_%	log2FC	p
Naive CD8 T cell	4.962	0.053	-6.520	0.082
Mature Naive B cell	14.758	3.626	-2.025	0.039
CLP	13.053	7.192	-0.860	0.031
cMoP	2.921	1.623	-0.847	0.230
Ery P	23.337	15.522	-0.588	0.012
Monocyte	9.096	6.296	-0.531	0.241
Mono P	11.406	8.016	-0.509	0.416
GMP	7.555	6.463	-0.225	0.063
Myelo P	9.621	8.500	-0.179	0.141
T cell	1.870	1.736	-0.107	0.737
pDC	1.963	1.945	-0.013	0.972
MKP	9.777	9.838	0.009	0.954
LMPP	24.420	25.368	0.055	0.713
Neutrophil	42.094	47.101	0.162	0.216
Pre-B cell	3.093	4.049	0.389	0.295
Immature B cell	5.112	7.595	0.571	0.150
Pre-B II cell	4.689	9.733	1.054	0.093
Memory B cell	0.430	2.941	2.773	0.015
CD8 TEM	0.106	2.625	4.613	0.070
HSC	0.571	18.474	5.014	0.017

from matplotlib.colors import to_rgba

sig = counted["p"] < 0.05
# solid = significant, faded = not; red = enriched in Old, blue = depleted
bar_colors = [to_rgba("#C44E52" if f > 0 else "#4C72B0", 1.0 if s else 0.35)
              for f, s in zip(counted["log2FC"], sig)]

fig, ax = plt.subplots(figsize=(6, 7))
ax.barh(counted.index, counted["log2FC"], color=bar_colors)
for y, (ct, row) in enumerate(counted.iterrows()):
    if sig[ct]:
        ax.text(row["log2FC"] + np.sign(row["log2FC"]) * 0.15, y, "*",
                va="center", ha="center", fontsize=13)
ax.axvline(0, color="k", lw=1)
ax.set(xlabel=f"log2 fold change in proportion ({CONDITIONS[1]} / {CONDITIONS[0]})",
       title="Per-cluster cell counting\n(* = p < 0.05, Welch t-test, n = 3 vs 3)")
plt.tight_layout()

print(f"significant cell types: {int(sig.sum())} of {len(counted)}")

significant cell types: 5 of 20

The classic aging phenotype falls straight out: HSCs expand sharply in old marrow (log2FC +5.0) and memory B cells rise, while the naive lymphoid compartment shrinks (mature naive B cells, CLPs) — as do erythroid progenitors.

Now read the p column, because that is where the count starts costing you:

• Only 5 of 20 cell types clear p < 0.05, and the largest effect in the table is not one of them. Naive CD8 T cells fall from 5.0% of the mature compartment to 0.05% — a ~90-fold depletion — at p = 0.08, and CD8 TEM cells rise ~25-fold at p = 0.07. (The second-largest effect, HSC, does clear it, at p = 0.02.) A per-cluster proportion crushes every cell into a single number per animal, so with three mice per condition even enormous effects can miss.

• The answer is only as good as the annotation. Change the clustering resolution and you change the result — as we are about to demonstrate.

• It says nothing about where inside a cluster the change happened.

2. Differential abundance without clusters

kompot.da never sees the labels. It estimates the cell-state density of each condition over DM_EigenVectors and scores every cell. See kompot.da.

da_results = kompot.da(
    adata,
    groupby=GROUPING_COLUMN,
    condition1=CONDITIONS[0],
    condition2=CONDITIONS[1],
    obsm_key=DIMENSIONALITY_REDUCTION,
)

[2026-07-11 23:12:20,388] [INFO    ] Condition 1 (Young): 2,917 cells
[2026-07-11 23:12:20,390] [INFO    ] Condition 2 (Old): 3,116 cells
[2026-07-11 23:12:20,399] [INFO    ] Fitting density estimator for condition 1...
[2026-07-11 23:13:17,059] [INFO    ] Fitting density estimator for condition 2...
[2026-07-11 23:14:03,696] [INFO    ] This run will have `run_id=0`.

Per-cell results land in adata.obs:

• kompot_da_Young_to_Old_lfc — log fold change (positive = enriched in Old)

• kompot_da_Young_to_Old_lfc_direction — up / down / unchanged

c1, c2 = CONDITIONS
direction = f"kompot_da_{c1}_to_{c2}_lfc_direction"
lfc = f"kompot_da_{c1}_to_{c2}_lfc"

sc.pl.embedding(
    adata, "umap",
    color=[direction, lfc],
    color_map="RdBu_r", vcenter=0, frameon=False,
)

Volcano plot

Effect size vs. significance, one point per cell (volcano_da):

kompot.plot.volcano_da(adata, color=CELL_TYPE_COLUMN)

[2026-07-11 23:14:04,055] [INFO    ] Found DA run info for run_id=-1
[2026-07-11 23:14:04,056] [INFO    ] Found lfc_key='kompot_da_Young_to_Old_lfc' from run info
[2026-07-11 23:14:04,056] [INFO    ] Found ptp_key='kompot_da_Young_to_Old_neg_log10_lfc_ptp' from run info
[2026-07-11 23:14:04,057] [INFO    ] Successfully inferred fields: {'lfc_key': 'kompot_da_Young_to_Old_lfc', 'ptp_key': 'kompot_da_Young_to_Old_neg_log10_lfc_ptp'}

Show only the significantly changed cells, split into enriched vs. depleted with mgroups (kompot.plot.embedding):

kompot.plot.embedding(
    adata, "umap",
    color=CELL_TYPE_COLUMN,
    mgroups=[{direction: [d]} for d in ["up", "down"]],
    frameon=False, wspace=.5,
)

[2026-07-11 23:14:05,095] [INFO    ] Selected 585 cells out of 8,090 total cells.
[2026-07-11 23:14:05,190] [INFO    ] Selected 490 cells out of 8,090 total cells.

Summary by cell type

The per-cell calls can of course still be summarized over an annotation — that is a read-out, not an input (direction_barplot):

kompot.plot.direction_barplot(adata, category_column=CELL_TYPE_COLUMN)

[2026-07-11 23:14:05,838] [INFO    ] Found DA run info for run_id=-1
[2026-07-11 23:14:05,839] [INFO    ] Found direction_key='kompot_da_Young_to_Old_lfc_direction' from run info
[2026-07-11 23:14:05,840] [INFO    ] Successfully inferred fields: {'direction_key': 'kompot_da_Young_to_Old_lfc_direction'}
[2026-07-11 23:14:05,840] [INFO    ] Using DA run 0 for direction_barplot: comparing Young to Old
[2026-07-11 23:14:05,841] [INFO    ] Creating direction barplot: comparing Young to Old
[2026-07-11 23:14:05,841] [INFO    ] Using fields - category_column: 'highres_celltype', direction_column: 'kompot_da_Young_to_Old_lfc_direction'

3. A different question

It is tempting to read Sections 1 and 2 as two routes to the same number and to check whether they agree. They are not. The per-cluster count asks how many cells fall in this cluster — a property of the annotation you supplied. kompot.da asks how common is this cell state — a property of where cells sit in expression space, estimated per cell and without reference to any labels.

The two usually correlate: a cluster that genuinely expands contains cell states that became more common, so both move together. That correlation is why the results often look interchangeable. But they are not measuring the same quantity, and the gap is easy to miss:

• A cluster can hold steady in size while the states inside it shift, if one subpopulation expands as another contracts. Reduced to a single proportion, those cancel; per cell, both are visible.

• The two can disagree simply because a cluster boundary need not follow the boundary along which abundance actually changed.

Neither is wrong, and neither supersedes the other. Reach for the per-cluster count when you have a trusted annotation and want a familiar per-cell-type number; reach for kompot.da when you want the shift located at whatever resolution it happens to occupy, without committing to labels first.

4. Accounting for sample variance

With biological replicates, add sample_col so the test reflects between-animal variation — the same concern that motivates the replicate-level t-test in the counting section, now applied per cell. The test becomes more conservative, so we relax the threshold slightly. This adds a _sample_var set of results as run_id=1:

da_sample_var = kompot.da(
    adata,
    groupby=GROUPING_COLUMN,
    condition1=c1, condition2=c2,
    obsm_key=DIMENSIONALITY_REDUCTION,
    sample_col=SAMPLE_COLUMN,
    threshold=kompot.DAThresholdSettings(ptp_threshold=0.1),  # relaxed from 0.05
)

[2026-07-11 23:14:07,676] [INFO    ] Results with result_key='kompot_da' already exist in the dataset. Previous run was at 2026-07-11T23:14:03.636291 comparing Young to Old. Fields that will be overwritten: obs.kompot_da_Young_to_Old_lfc, obs.kompot_da_Young_log_density, obs.kompot_da_Old_log_density. Note: Only fields NOT affected by sample variance (like lfc, log_density) will be overwritten since they don't use the sample variance suffix. These results will likely be identical if other parameters haven't changed. This is a partial rerun with sample variance added to a previous analysis with matching parameters. Set overwrite=False to prevent overwriting or overwrite=True to silence this message.
[2026-07-11 23:14:07,678] [INFO    ] Condition 1 (Young): 2,917 cells
[2026-07-11 23:14:07,678] [INFO    ] Condition 2 (Old): 3,116 cells
[2026-07-11 23:14:07,682] [INFO    ] Using sample column 'Replicate' for sample variance estimation
[2026-07-11 23:14:07,684] [INFO    ] Found 3 unique sample(s) in condition 1
[2026-07-11 23:14:07,685] [INFO    ] Found 3 unique sample(s) in condition 2
[2026-07-11 23:14:07,686] [INFO    ] Fitting density estimator for condition 1...
[2026-07-11 23:14:27,714] [INFO    ] Fitting density estimator for condition 2...
[2026-07-11 23:14:53,472] [INFO    ] Sample variance estimation automatically enabled due to provided sample indices
[2026-07-11 23:14:53,474] [INFO    ] Setting up sample variance estimation...
[2026-07-11 23:14:53,477] [INFO    ] Found 3 unique groups for variance estimation
[2026-07-11 23:14:53,478] [INFO    ] 3 groups have sufficient cells (>= 2) for variance estimation
[2026-07-11 23:14:53,479] [INFO    ] Training estimator for group 1 with 944 cells
[2026-07-11 23:15:03,574] [INFO    ] Training estimator for group 2 with 1,007 cells
[2026-07-11 23:15:12,572] [INFO    ] Training estimator for group 3 with 966 cells
[2026-07-11 23:15:19,902] [INFO    ] Found 3 unique groups for variance estimation
[2026-07-11 23:15:19,904] [INFO    ] 3 groups have sufficient cells (>= 2) for variance estimation
[2026-07-11 23:15:19,904] [INFO    ] Training estimator for group 1 with 1,057 cells
[2026-07-11 23:15:30,214] [INFO    ] Training estimator for group 2 with 1,019 cells
[2026-07-11 23:15:39,994] [INFO    ] Training estimator for group 3 with 1,040 cells
[2026-07-11 23:15:56,065] [INFO    ] Computing sample-specific variance for condition 1...
[2026-07-11 23:16:00,929] [INFO    ] Computing sample-specific variance for condition 2...
[2026-07-11 23:16:05,447] [INFO    ] This run will have `run_id=1`.

Fold change is unchanged; significance (the z-score) becomes more conservative:

sc.pl.embedding(
    adata, "umap",
    color=[
        f"kompot_da_{c1}_to_{c2}_lfc",
        f"kompot_da_{c1}_to_{c2}_lfc_zscore",
        f"kompot_da_{c1}_to_{c2}_lfc_zscore_sample_var",
    ],
    color_map="RdBu_r", vcenter=0, ncols=3, frameon=False,
)

Volcano plots, standard vs. sample-variance:

kompot.plot.volcano_da(adata, color=CELL_TYPE_COLUMN, run_id=0, title="Standard (pooled samples)")

[2026-07-11 23:16:06,034] [INFO    ] Found DA run info for run_id=0
[2026-07-11 23:16:06,034] [INFO    ] Found lfc_key='kompot_da_Young_to_Old_lfc' from run info
[2026-07-11 23:16:06,035] [INFO    ] Found ptp_key='kompot_da_Young_to_Old_neg_log10_lfc_ptp' from run info
[2026-07-11 23:16:06,036] [WARNING ] Field 'kompot_da_Young_to_Old_lfc' has been written by 2 different runs, indicating potential overwrites
[2026-07-11 23:16:06,036] [WARNING ] Field 'kompot_da_Young_to_Old_neg_log10_lfc_ptp' was last written by run 0, but current context expects run 1. The field may have been overwritten.
[2026-07-11 23:16:06,037] [INFO    ] Successfully inferred fields: {'lfc_key': 'kompot_da_Young_to_Old_lfc', 'ptp_key': 'kompot_da_Young_to_Old_neg_log10_lfc_ptp'}
[2026-07-11 23:16:06,038] [WARNING ] Field inference completed with 2 warnings

kompot.plot.volcano_da(adata, color=CELL_TYPE_COLUMN, run_id=1, title="With sample variance")

[2026-07-11 23:16:07,018] [INFO    ] Found DA run info for run_id=1
[2026-07-11 23:16:07,019] [INFO    ] Found lfc_key='kompot_da_Young_to_Old_lfc' from run info
[2026-07-11 23:16:07,019] [INFO    ] Found ptp_key='kompot_da_Young_to_Old_neg_log10_lfc_ptp_sample_var' from run info
[2026-07-11 23:16:07,020] [WARNING ] Field 'kompot_da_Young_to_Old_lfc' has been written by 2 different runs, indicating potential overwrites
[2026-07-11 23:16:07,021] [INFO    ] Successfully inferred fields: {'lfc_key': 'kompot_da_Young_to_Old_lfc', 'ptp_key': 'kompot_da_Young_to_Old_neg_log10_lfc_ptp_sample_var'}
[2026-07-11 23:16:07,022] [WARNING ] Field inference completed with 1 warnings

Fewer cells clear significance under sample variance, but effect sizes and the overall cell-type pattern are preserved — the surviving calls are the robust ones.

kompot.cleanup(adata)
adata.write_h5ad("../data/murine_bone_marrow_aging_da.h5ad")

[2026-07-11 23:16:07,999] [WARNING ] No run history found for de analysis.

Where to go deeper

• Getting Started — differential expression.

• DE with sample variance — the gene-level counterpart to the section above.

• Smooth expression — the GP primitive that underlies both DA and DE.

• Paper · Documentation