Multi-dataset: Differential Expression Analysis

This example demonstrates how to identify differentially expressed genes in a multi-FOV setting using spatial_query_multi.de_genes.

Two common comparisons are shown:

1. Across conditions: Compare motif+ anchor cells between two conditions (e.g., motif+ podocytes in normal vs diabetic kidney disease)

2. Within a condition: Compare motif+ vs motif− anchor cells within one condition

Dataset: CZI Kidney — normal vs autosomal dominant tubulointerstitial kidney disease (ADTKD) (spatial transcriptomics), availabel at https://cellxgene.cziscience.com/collections/8e880741-bf9a-4c8e-9227-934204631d2a

Key API: spatial_query_multi.motif_enrichement_dist/spatial_query_multi.motif_enrichment_knn and spatial_query_multi.de_genes

Setup

import warnings
warnings.filterwarnings("ignore")
import os

import anndata as ad
import numpy as np

from SpatialQuery import spatial_query_multi

Load Data & Initialize

DATA_DIR = "../data/CZI_kidney"

data_files = os.listdir(DATA_DIR)
adatas = [ad.read_h5ad(os.path.join(DATA_DIR, f)) for f in data_files]

print(f"Loaded {len(adatas)} FOVs")
adatas[0]

Loaded 84 FOVs

AnnData object with n_obs × n_vars = 12906 × 17811
    obs: 'assay_ontology_term_id', 'self_reported_ethnicity_ontology_term_id', 'is_primary_data', 'organism_ontology_term_id', 'sample', 'tissue_ontology_term_id', 'disease_state', 'sex_ontology_term_id', 'genotype', 'development_stage_ontology_term_id', 'author_cell_type', 'cell_type_ontology_term_id', 'disease_ontology_term_id', 'donor_id', 'suspension_type', 'tissue_type', 'cell_type', 'assay', 'disease', 'organism', 'sex', 'tissue', 'self_reported_ethnicity', 'development_stage', 'observation_joinid'
    var: 'feature_is_filtered', 'feature_name', 'feature_reference', 'feature_biotype', 'feature_length'
    uns: 'citation', 'schema_reference', 'schema_version', 'title'
    obsm: 'X_spatial'

# Create condition labels and filter out mitochondrial / unannotated genes
for i, adata in enumerate(adatas):
    adata.obs['disease_state_genotype'] = (
        adata.obs['disease_state'].astype(str) + '_' + adata.obs['genotype'].astype(str)
    )
    adata.var_names = adata.var['feature_name'].tolist()
    adata = adata[:, ~adata.var_names.str.startswith('mt-')]
    adata = adata[:, ~adata.var_names.str.endswith('Rik')]
    adata = adata[:, ~adata.var_names.str.startswith('Gm')]
    adatas[i] = adata

# Map verbose condition names to short labels
dataset_mapping = {
    'early diabetic kidney disease_BTBR-ob/ob': 'DKD_BTBR-ob/ob',
    'autosomal dominant tubulointersital kidney disease (ADTKD)_UMOD-KI/KI': 'ADTKD_UMOD-KI/KI',
    'control_BTBR-wt/wt': 'control_BTBR-WT/WT',
    'control_UMOD-WT/WT': 'control_UMOD-WT/WT',
}

dataset_name_col = 'disease_state_genotype'
datasets = [adata.obs[dataset_name_col].unique()[0] for adata in adatas]
dataset_names = [dataset_mapping[d] for d in datasets]

spatial_key = "X_spatial"
label_key = "cell_type"
feature_name = "feature_name"

Inspect adata.obs for disease labels and cell types. Then split by condition.

spm = spatial_query_multi(
    adatas=adatas,
    datasets=dataset_names,
    spatial_key=spatial_key,
    label_key=label_key,
    feature_name=feature_name,
    build_gene_index=False,
    if_lognorm=True,
    if_normalize_spatial_coord=True
)

Replacing _ with hyphen in DKD_BTBR-ob/ob.
Replacing _ with hyphen in DKD_BTBR-ob/ob.
Replacing _ with hyphen in control_BTBR-WT/WT.
Replacing _ with hyphen in control_UMOD-WT/WT.
Replacing _ with hyphen in ADTKD_UMOD-KI/KI.
Replacing _ with hyphen in ADTKD_UMOD-KI/KI.
Replacing _ with hyphen in ADTKD_UMOD-KI/KI.
Replacing _ with hyphen in control_BTBR-WT/WT.
Replacing _ with hyphen in control_BTBR-WT/WT.
Replacing _ with hyphen in DKD_BTBR-ob/ob.
...(492 lines omitted)...
build_gene_index is False. Using adata.X for gene expression analysis.
Log normalizing the expression data... If data is already log normalized, please set if_lognorm to False.

Define Anchor and Motif

We use macrophage as the anchor cell type. The motif includes fibroblast and thick ascending limb epithelial cell.

anchor_ct = "macrophage"
motif = ['kidney interstitial fibroblast', 'kidney loop of Henle thick ascending limb epithelial cell']
max_dist = 5
ds_control = 'control-UMOD-WT/WT'
ds_case = 'ADTKD-UMOD-KI/KI'

Retrieve Motif+ Cell IDs Per Condition

Use motif_enrichment_dist with return_cellID=True to obtain the indices of motif-positive anchor cells and their neighbors for each condition. The returned cell IDs are organized as dictionaries keyed by FOV name.

Note: In the multi-FOV setting, de_genes expects ind_group1 and ind_group2 to be dictionaries {fov_name: [cell_indices]}, not flat lists.

control_result, control_motif_id, control_center_id = spm.motif_enrichment_dist(
    ct=anchor_ct, motifs=motif, dataset=ds_control, max_dist=max_dist, return_cellID=True,
)

case_result, case_motif_id, case_center_id = spm.motif_enrichment_dist(
    ct=anchor_ct, motifs=motif, dataset=ds_case, max_dist=max_dist, return_cellID=True,
)

print(f"Control motif+ anchors: {sum(len(v) for v in control_center_id[str(sorted(motif))].values())} cells")
print(f"Case motif+ anchors: {sum(len(v) for v in case_center_id[str(sorted(motif))].values())} cells")

Control motif+ anchors: 3156 cells
Case motif+ anchors: 10698 cells

Alternative: KNN-based Cell ID Retrieval

You can also retrieve motif-positive cell IDs using KNN neighborhoods.

k = 20

control_result_knn, control_motif_id_knn, control_center_id_knn = spm.motif_enrichment_knn(
    ct=anchor_ct, motifs=motif, dataset=ds_control, k=k, return_cellID=True,
)

case_result_knn, case_motif_id_knn, case_center_id_knn = spm.motif_enrichment_knn(
    ct=anchor_ct, motifs=motif, dataset=ds_case, k=k, return_cellID=True,
)

print(f"KNN Control motif+ anchors: {sum(len(v) for v in control_center_id_knn[str(sorted(motif))].values())} cells")
print(f"KNN Case motif+ anchors: {sum(len(v) for v in case_center_id_knn[str(sorted(motif))].values())} cells")

KNN Control motif+ anchors: 1997 cells
KNN Case motif+ anchors: 8532 cells

Comparison 1: Motif+ Cells Across Conditions

Compare gene expression of motif-positive anchor cells between normal and DKD.

Key parameters (multi-FOV version):

Parameter	Description
ind_group1	Dict of {fov_name: [cell_indices]} for group 1
ind_group2	Dict of {fov_name: [cell_indices]} for group 2
method	Statistical test: 't-test', 'wilcoxon', or 'fisher'
alpha	Significance threshold

The multi-FOV de_genes automatically handles pooling across FOVs.

de_across = spm.de_genes(
    ind_group1=control_center_id[str(sorted(motif))],
    ind_group2=case_center_id[str(sorted(motif))],
    method="t-test",
    alpha=0.05,
)

print(f"DE genes (Control vs Case motif+ {anchor_ct}): {len(de_across)}")
de_across.head(10)

Testing 1669 genes ...
DE genes (Control vs Case motif+ macrophage): 1466

	gene	p_value	adj-pval	log2fc	proportion_1	proportion_2	abs_difference	de_in
0	Apoe	0.000000e+00	0.000000e+00	-2.901169	0.168885	0.623107	0.454222	group2
1	Mgp	0.000000e+00	0.000000e+00	-2.684509	0.119455	0.478314	0.358859	group2
2	Spp1	1.314239e-234	7.311551e-232	-1.984394	0.129278	0.372967	0.243689	group2
3	Lcn2	3.064570e-190	1.278692e-187	-4.911119	0.004119	0.119555	0.115436	group2
4	Igfbp7	1.892602e-178	6.317507e-176	-1.054887	0.517744	0.696766	0.179022	group2
5	Dpt	1.660872e-175	4.619992e-173	-3.769909	0.011090	0.132548	0.121458	group2
6	Lyz2	2.759899e-165	6.580389e-163	-2.049103	0.074461	0.268181	0.193720	group2
7	Tpm1	2.577798e-158	5.377930e-156	-1.991123	0.079848	0.257899	0.178051	group2
8	Sparc	1.252472e-157	2.322639e-155	-2.217716	0.058619	0.226958	0.168340	group2
9	Cfh	2.595302e-156	4.331559e-154	-1.791375	0.107414	0.299869	0.192455	group2

Comparison 2: Motif+ vs Motif− Within DKD

Compare gene expression between motif-positive and motif-negative anchor cells within the DKD condition only.

# Get non-motif anchor cell IDs for DKD FOVs
non_motif_center = {str(sorted(motif)): {}}
for sp in spm.spatial_queries:
    if sp.dataset.split("_")[0] != ds_case:
        continue
    ct_id = np.where(sp.labels == anchor_ct)[0]
    motif_ids = case_center_id[str(sorted(motif))].get(sp.dataset, [])
    non_motif_center[str(sorted(motif))][sp.dataset] = list(set(ct_id) - set(motif_ids))

de_within = spm.de_genes(
    ind_group1=case_center_id[str(sorted(motif))],
    ind_group2=non_motif_center[str(sorted(motif))],
    method="t-test",
    alpha=0.05,
)

print(f"DE genes (motif+ vs motif− in DKD): {len(de_within)}")
de_within.head(10)

Testing 1919 genes ...
DE genes (motif+ vs motif− in DKD): 1274

	gene	p_value	adj-pval	log2fc	proportion_1	proportion_2	abs_difference	de_in
0	Kap	0.000000e+00	0.000000e+00	-1.128119	0.808843	0.901811	0.092968	group2
1	Apoe	1.884793e-259	1.808458e-256	1.029396	0.623107	0.435769	0.187338	group1
2	Aldob	4.441241e-214	2.840914e-211	-1.014052	0.317162	0.520432	0.203270	group2
3	Malat1	2.169127e-185	1.040639e-182	0.648035	0.863993	0.776140	0.087854	group1
4	Napsa	2.504258e-168	9.611342e-166	-0.932830	0.326136	0.503900	0.177765	group2
5	Slc12a1	7.426424e-155	2.375218e-152	1.007723	0.388764	0.243827	0.144937	group1
6	Acsm2	2.493802e-153	6.836580e-151	-1.001984	0.218826	0.390109	0.171284	group2
7	Gpx3	4.661235e-150	1.118114e-147	-0.855314	0.335857	0.501539	0.165682	group2
8	Akr1c21	4.131642e-144	8.809578e-142	-0.837663	0.316134	0.486653	0.170519	group2
9	Acy3	3.313493e-139	6.358593e-137	-0.868884	0.278837	0.455020	0.176183	group2