Differential Expression with Sample Variance

The natural next step after Getting Started.

With biological replicates, variation between samples of the same condition is often larger than technical noise. Ignoring it inflates significance and hurts reproducibility. Kompot builds a per-sample estimator, measures sample-to-sample covariance, and folds it into the Mahalanobis distance — yielding conservative, replicate-aware calls. The derivation is in the paper.

This dataset has 3 replicates per age group.

Cost warning — sample variance is priced per gene, twice over. Without it, every gene is scored against a single shared posterior covariance of shape (n_landmarks, n_landmarks). With it, every gene gets its own such matrix per condition, and two costs follow that respond to different levers:

  • memory, 2 × n_landmarks² × n_genes × 8 bytes — about 0.37 GiB per gene at the default n_landmarks=5000. StorageSettings(store_arrays_on_disk=True) keeps this out of RAM entirely, and n_landmarks shrinks it quadratically;

  • compute, one Cholesky factorisation per gene instead of one in total. store_arrays_on_disk does not touch this, but lowering n_landmarks does (~0.016 s/gene at 500 against ~2.1 s at 5 000, single-threaded), as does analysing fewer genes.

So the prescription is a cheap first pass over all genes, then a sample-variance pass restricted to the top ones, which is what this notebook does. The levers and measured plans at realistic sizes: Planning Memory and Disk.

[1]:

%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

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"

[2]:

DATA_PATH = "../data/murine_bone_marrow_aging.h5ad"
SAMPLE_COLUMN = "Replicate"             # biological-replicate identifiers
GROUPING_COLUMN = "Age"
CONDITIONS = ["Young", "Old"]
CELL_TYPE_COLUMN = "highres_celltype"
DIMENSIONALITY_REDUCTION = "DM_EigenVectors"
LAYER_FOR_EXPRESSION = "logged_counts"

Load Data

[3]:

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

[3]:

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'

Replicates per condition:

[4]:

adata.obs.groupby([GROUPING_COLUMN, SAMPLE_COLUMN]).size().unstack(fill_value=0)

/tmp/ipykernel_651/3107917919.py:1: FutureWarning: The default of observed=False is deprecated and will be changed to True in a future version of pandas. Pass observed=False to retain current behavior or observed=True to adopt the future default and silence this warning.
  adata.obs.groupby([GROUPING_COLUMN, SAMPLE_COLUMN]).size().unstack(fill_value=0)

[4]:
Replicate 1 2 3
Age
Mid 1035 1022 0
Old 1057 1019 1040
Young 944 1007 966

Pass 1 — Baseline DE Without Sample Variance

Run standard DE over all genes first. This is the cheap pass, and it produces the Mahalanobis ranking that pass 2 is restricted to (run_id=0):

[5]:

de_standard = kompot.de(
    adata,
    groupby=GROUPING_COLUMN,
    condition1=CONDITIONS[0],
    condition2=CONDITIONS[1],
    layer=LAYER_FOR_EXPRESSION,
    obsm_key=DIMENSIONALITY_REDUCTION,
    gp=kompot.GPSettings(batch_size=0),
)

[2026-07-11 01:30:39,288] [INFO    ] Condition 1 (Young): 2,917 cells
[2026-07-11 01:30:39,289] [INFO    ] Condition 2 (Old): 3,116 cells
[2026-07-11 01:30:43,908] [INFO    ] Fitting expression estimator for condition 1...
[2026-07-11 01:31:12,552] [INFO    ] Fitting expression estimator for condition 2...
[2026-07-11 01:32:15,129] [INFO    ] FDR analysis complete: 192/16285 genes significantly DE at FDR < 0.05
[2026-07-11 01:32:15,131] [INFO    ] Mahalanobis distance threshold for FDR < 0.05: 5.8215
[2026-07-11 01:32:21,822] [INFO    ] This run will have `run_id=0`.

Pass 2 — DE With Sample Variance, Restricted to the Top Genes

Add sample_col. Kompot now builds one covariance matrix per gene per condition, so the run is priced per gene and we pull every lever available:

  1. subset to top genes from pass 1 — this is the only lever that touches both the memory and the per-gene factorisation, which is why it is the one that decides whether the run is possible at all. Fifty genes keeps this tutorial fast; on a real study, on the order of 200 is the working recommendation;

  2. store_arrays_on_disk=True (StorageSettings) — Kompot then assembles one gene’s covariance at a time instead of holding the tensors, which takes memory almost back to the no-sample-variance level. With dask installed nothing is even written to disk;

  3. lower n_landmarks if the run is still too slow — memory is quadratic in it, and it cuts the per-gene factorisation too, so it buys on both axes;

  4. skip FDR (null_genes=0) — not yet calibrated for sample variance, and null genes would be charged at the same per-gene rate as real ones. With sample_col set this is already the default.

Always check the plan first with dry_run=True, which reports memory, disk, and the fields that would be written, without running anything:

The stored outputs of the next two cells were cleared, deliberately.

They were captured on Kompot 0.8.0 and they assert a resource plan the fixed estimator does not produce: the dry run showed an 18.63 GB disk requirement that is now zero on the Dask path, and the run logged a per-gene Dask path that LazyGeneCovariance deliberately no longer takes. A blank cell that says why is better than a number that lies. Re-run the notebook to regenerate them for your own data.

Every result below is unaffected. The 0.9.0 fix is numerically identical — bitwise so on the in-memory and Dask paths — so the comparison table, the volcano plots and the heatmap are still exactly the ones this notebook was written around.

[ ]:

c1, c2 = CONDITIONS
mahal = f"kompot_de_{c1}_to_{c2}_mahalanobis"
top_genes = adata.var.sort_values(mahal, ascending=False).head(50).index

plan = kompot.de(
    adata,
    groupby=GROUPING_COLUMN,
    condition1=c1, condition2=c2,
    layer=LAYER_FOR_EXPRESSION,
    obsm_key=DIMENSIONALITY_REDUCTION,
    sample_col=SAMPLE_COLUMN,
    genes=top_genes,
    gp=kompot.GPSettings(batch_size=0),
    storage=kompot.StorageSettings(store_arrays_on_disk=True),
    fdr=kompot.FDRSettings(null_genes=0),
    dry_run=True,
)

If the plan fits, run it (run_id=1):

[ ]:

de_sample_var = kompot.de(
    adata,
    groupby=GROUPING_COLUMN,
    condition1=c1, condition2=c2,
    layer=LAYER_FOR_EXPRESSION,
    obsm_key=DIMENSIONALITY_REDUCTION,
    sample_col=SAMPLE_COLUMN,
    genes=top_genes,
    gp=kompot.GPSettings(batch_size=0),
    storage=kompot.StorageSettings(store_arrays_on_disk=True),
    fdr=kompot.FDRSettings(null_genes=0),
)

Compare

Sample-variance results carry a _sample_var suffix. Fold changes are unchanged; only significance shifts:

[8]:

cols = [
    f"kompot_de_{c1}_to_{c2}_mean_lfc",
    f"kompot_de_{c1}_to_{c2}_mahalanobis",
    f"kompot_de_{c1}_to_{c2}_mahalanobis_sample_var",
]
comparison = adata.var.loc[top_genes, cols].sort_values(cols[1], ascending=False)
comparison.head(20)

[8]:
kompot_de_Young_to_Old_mean_lfc kompot_de_Young_to_Old_mahalanobis kompot_de_Young_to_Old_mahalanobis_sample_var
H2-Q7 0.371810 15.224920 10.043567
Cd74 0.105968 12.781286 10.791307
H2-Aa 0.142574 12.641082 8.861583
H2-Ab1 0.153822 12.307942 8.316963
Igkc 0.016244 11.416240 8.666152
H2-Eb1 0.123908 11.371004 7.737303
AW112010 0.238251 11.366238 7.150069
S100a9 -0.051261 10.385744 9.975386
Ifitm3 0.110643 10.134062 8.776124
S100a8 -0.049699 10.089391 9.740630
H2-Q6 0.235504 9.790485 6.996295
Cd52 -0.038509 9.649612 9.086571
Aldh1a1 0.147759 9.493316 7.973998
Ifitm1 0.128764 9.438788 8.316066
Ifitm2 0.086038 9.357804 8.315513
Ighm -0.052318 9.204495 7.853406
Cd79a -0.001314 9.182240 7.399574
Apoe -0.039501 9.142740 8.440726
Fos 0.159999 8.999856 8.108269
Gm47283 0.188797 8.930669 5.335491

Genes that lose significance are the ones that vary most across replicates. Volcano plots, standard vs. sample-variance:

[9]:

kompot.plot.volcano_de(adata, run_id=0, title="Standard (pooled samples)")

[2026-07-11 01:36:30,400] [INFO    ] Found DE run info for run_id=0
[2026-07-11 01:36:30,401] [INFO    ] Found mean_lfc_key='kompot_de_Young_to_Old_mean_lfc' from run info
[2026-07-11 01:36:30,402] [INFO    ] Found mahalanobis_key='kompot_de_Young_to_Old_mahalanobis' from run info
[2026-07-11 01:36:30,404] [WARNING ] Field 'kompot_de_Young_to_Old_mean_lfc' has been written by 2 different runs, indicating potential overwrites
[2026-07-11 01:36:30,405] [WARNING ] Field 'kompot_de_Young_to_Old_mahalanobis' was last written by run 0, but current context expects run 1. The field may have been overwritten.
[2026-07-11 01:36:30,406] [INFO    ] Successfully inferred fields: {'mean_lfc_key': 'kompot_de_Young_to_Old_mean_lfc', 'mahalanobis_key': 'kompot_de_Young_to_Old_mahalanobis'}
[2026-07-11 01:36:30,407] [WARNING ] Field inference completed with 2 warnings
[2026-07-11 01:36:30,409] [INFO    ] Using DE run 0: comparing Young to Old
[2026-07-11 01:36:30,423] [INFO    ] Using data columns from var - lfc: 'kompot_de_Young_to_Old_mean_lfc', score: 'kompot_de_Young_to_Old_mahalanobis'
[2026-07-11 01:36:30,428] [INFO    ] Highlighting 192 genes marked as DE (109 up, 83 down)
[2026-07-11 01:36:30,441] [INFO    ] Labeling top 10 genes by score
../_images/notebooks_03_sample_variance_17_1.png 
[10]:

kompot.plot.volcano_de(
    adata, run_id=1,
    n_top_genes=40,          # no FDR threshold here (null_genes=0)
    gene_labels=20,
    title="With sample variance (top 50 genes)",
)

[2026-07-11 01:36:30,748] [INFO    ] Found DE run info for run_id=1
[2026-07-11 01:36:30,749] [INFO    ] Found mean_lfc_key='kompot_de_Young_to_Old_mean_lfc' from run info
[2026-07-11 01:36:30,749] [INFO    ] Found mahalanobis_key='kompot_de_Young_to_Old_mahalanobis_sample_var' from run info
[2026-07-11 01:36:30,752] [WARNING ] Field 'kompot_de_Young_to_Old_mean_lfc' has been written by 2 different runs, indicating potential overwrites
[2026-07-11 01:36:30,754] [INFO    ] Successfully inferred fields: {'mean_lfc_key': 'kompot_de_Young_to_Old_mean_lfc', 'mahalanobis_key': 'kompot_de_Young_to_Old_mahalanobis_sample_var'}
[2026-07-11 01:36:30,754] [WARNING ] Field inference completed with 1 warnings
[2026-07-11 01:36:30,755] [INFO    ] Using DE run 1: comparing Young to Old
[2026-07-11 01:36:30,767] [INFO    ] Using data columns from var - lfc: 'kompot_de_Young_to_Old_mean_lfc', score: 'kompot_de_Young_to_Old_mahalanobis_sample_var'
[2026-07-11 01:36:30,770] [INFO    ] Highlighting top 40 genes by kompot_de_Young_to_Old_mahalanobis_sample_var
[2026-07-11 01:36:30,779] [INFO    ] Labeling top 20 genes by score
../_images/notebooks_03_sample_variance_18_1.png

Sample-to-Sample Variability

Rank genes by how much significance they lose, then compare their expression between two replicates of the same condition:

[11]:

comparison["decrease"] = comparison[cols[1]] - comparison[cols[2]]
variable_genes = comparison.nlargest(20, "decrease").index
list(variable_genes)

[11]:

['H2-Q7',
 'AW112010',
 'H2-Ab1',
 'H2-Aa',
 'H2-Eb1',
 'Gm47283',
 'Ifi27',
 'Cirbp',
 'H2-Q6',
 'Igkc',
 'H2-K1',
 'Trbc2',
 'Cd74',
 'Lilr4b',
 'Cd3g',
 'Ms4a4b',
 'Ebf1',
 'Cd79a',
 'Ly6d',
 'Cd3d']

[12]:

kompot.plot.heatmap(
    adata,
    genes=variable_genes,
    groupby=CELL_TYPE_COLUMN,
    condition_column=SAMPLE_COLUMN,
    condition1="1", condition1_name="Replicate 1",
    condition2="3", condition2_name="Replicate 3",
    vmin="p2", vmax="p98",
    cluster_rows=False,
)

[2026-07-11 01:36:31,110] [INFO    ] Inferred layer='logged_counts' from run information
[2026-07-11 01:36:31,111] [INFO    ] Creating split heatmap with 20 genes/features
[2026-07-11 01:36:31,112] [INFO    ] Using display names: 'Replicate 1' for condition1, 'Replicate 3' for condition2
[2026-07-11 01:36:31,113] [INFO    ] Using expression data from layer: 'logged_counts'
[2026-07-11 01:36:31,236] [INFO    ] Applying gene-wise z-scoring (standard_scale='var')
../_images/notebooks_03_sample_variance_21_1.png

These genes differ substantially between replicates of the same age — exactly why their significance drops under sample variance.

Best Practices

Use sample variance when you have ≥3 replicates per condition and want conservative, reproducible calls. Effect sizes stay identical; only confidence changes, and the genes that survive the correction are the robust, replicate-consistent signals.

Always run it as a second pass. The levers, and what each one actually buys:

Lever

Memory

Per-gene compute

Practical setting

genes=

linear

linear

top ~200 from pass 1

StorageSettings(store_arrays_on_disk=True)

removes the tensor term

—

set disk_storage_dir to real scratch

GPSettings(n_landmarks=…)

quadratic

reduces it; measure rather than extrapolate

1 000–2 000 when 5 000 is too slow

FDRSettings(null_genes=0)

removes 2 000 genes’ worth

likewise

already the default once sample_col is set

Note what the table says about genes: it is the only row that appears in both cost columns. Offloading solves the memory axis outright — a whole-transcriptome run goes from 7 552 GiB to 101 GiB — and leaves the per-gene factorisation untouched, which is a long job rather than an impossible one (~11.6 h at 5 000 landmarks, single-threaded). So restricting the gene list is about spending your compute on the genes you will actually read, not about clearing a wall.

And price it before you run it: dry_run=True on the exact call you intend to make takes a second. Full treatment, with measured plans: Planning Memory and Disk.

[13]:

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

[2026-07-11 01:36:33,094] [INFO    ] Cleaning up all 2 run(s)
[2026-07-11 01:36:33,248] [INFO    ] Cleaned up 3 field(s) from run 0:
[2026-07-11 01:36:33,249] [INFO    ]   layers (3 field(s)):
[2026-07-11 01:36:33,250] [INFO    ]     - kompot_de_Young_smoothed
[2026-07-11 01:36:33,250] [INFO    ]     - kompot_de_Old_smoothed
[2026-07-11 01:36:33,251] [INFO    ]     - kompot_de_Young_to_Old_fold_change
[2026-07-11 01:36:33,255] [INFO    ] Cleaned up 2 field(s) from run 1:
[2026-07-11 01:36:33,256] [INFO    ]   layers (2 field(s)):
[2026-07-11 01:36:33,256] [INFO    ]     - kompot_de_Young_std
[2026-07-11 01:36:33,257] [INFO    ]     - kompot_de_Old_std
[2026-07-11 01:36:33,258] [INFO    ] Total: Cleaned up 5 field(s) across 2 run(s)

Where to go deeper