from typing import List, Union, Optional, Dict, Tuple, Any, Literal
import numpy as np
import pandas as pd
import scipy.stats as stats
import statsmodels.stats.multitest as mt
from anndata import AnnData
from mlxtend.frequent_patterns import fpgrowth
from sklearn.preprocessing import MultiLabelBinarizer
from scipy.stats import hypergeom
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from collections import defaultdict
from sklearn.preprocessing import LabelEncoder
from statsmodels.stats.multitest import multipletests
from . import spatial_gene_covarying
from .spatial_differential_pattern import differential_analysis_motif_knn, differential_analysis_motif_dist
from .spatial_query import spatial_query
from .spatial_utils import (
find_maximal_patterns,
build_fptree_knn,
build_fptree_dist,
de_genes_scanpy,
de_genes_fisher,
)
import anndata as ad
from time import time
class spatial_query_multi:
[docs]
def __init__(
self,
adatas: List[AnnData],
datasets: List[str],
spatial_key: str,
label_key: str,
leaf_size: int = 10,
build_gene_index: bool = False,
feature_name: str = None,
if_lognorm: bool = True,
if_normalize_spatial_coord: bool = True,
):
"""
Initiate models, including setting attributes and building kd-tree for each field of view.
Parameters
----------
adatas : List[AnnData]
List of AnnData objects.
datasets : List[str]
List of dataset names.
spatial_key : str
Spatial coordination name in AnnData.obsm object.
label_key : str
Label name in AnnData.obs object.
leaf_size : int, default=10
The largest number of points stored in each leaf node.
build_gene_index : bool, default=False
Whether to build scfind index or use adata.X directly.
feature_name : str, optional
The label or key in the AnnData object's variables (var) that corresponds to the feature names.
if_lognorm : bool, default=True
Whether to log normalize the expression data.
if_normalize_spatial_coord : bool, default=True
If True, normalizes spatial coordinates so mean nearest neighbor distance equals 1.
Set to False if original spatial units should be preserved.
"""
# Each element in self.spatial_queries stores a spatial_query object
self.spatial_key = spatial_key
self.label_key = label_key
self.build_gene_index = build_gene_index
# Modify dataset names by d_0, d_2, ... for duplicates in datasets
count_dict = {}
modified_datasets = []
for dataset in datasets:
if '_' in dataset:
print(f"Replacing _ with hyphen in {dataset}.")
dataset = dataset.replace('_', '-')
if dataset in count_dict:
count_dict[dataset] += 1
else:
count_dict[dataset] = 0
mod_dataset = f"{dataset}_{count_dict[dataset]}"
modified_datasets.append(mod_dataset)
self.datasets = modified_datasets
self.spatial_queries = [spatial_query(
adata=adata,
dataset=self.datasets[i],
spatial_key=spatial_key,
label_key=label_key,
leaf_size=leaf_size,
build_gene_index=build_gene_index,
feature_name=feature_name,
if_lognorm=if_lognorm,
if_normalize_spatial_coord=if_normalize_spatial_coord,
) for i, adata in enumerate(adatas)]
[docs]
def find_fp_knn(self,
ct: str,
dataset: Union[str, List[str]] = None,
k: int = 30,
min_support: float = 0.5,
max_dist: float = 20
) -> pd.DataFrame:
"""
Find frequent patterns within the KNNs of certain cell type in multiple fields of view.
Parameters
----------
ct : str
Cell type name.
dataset : str or List[str], optional
Datasets for searching for frequent patterns.
Use all datasets if dataset=None.
k : int, default=30
Number of nearest neighbors.
min_support : float, default=0.5
Threshold of frequency to consider a pattern as a frequent pattern.
max_dist : float, default=20
The maximum distance at which points are considered neighbors.
Returns
-------
pd.DataFrame
DataFrame with frequent pattern results. Columns include:
- support: frequency of the pattern (proportion of transactions with this pattern)
- itemsets: list of cell types in the frequent pattern
"""
# Search transactions for each field of view, find the frequent patterns of integrated transactions
# start = time.time()
if_exist_label = [ct in s.labels.unique() for s in self.spatial_queries]
if not any(if_exist_label):
raise ValueError(f"Found no {self.label_key} in all datasets!")
if dataset is None:
# Use all datasets if dataset is not provided
dataset = [s.dataset.split('_')[0] for s in self.spatial_queries]
# Make sure dataset is a list
if isinstance(dataset, str):
dataset = [dataset]
# test if the input dataset name is valid
valid_ds_names = [s.dataset.split('_')[0] for s in self.spatial_queries]
for ds in dataset:
if ds not in valid_ds_names:
raise ValueError(f"Invalid input dataset name: {ds}.\n "
f"Valid dataset names are: {set(valid_ds_names)}")
# end = time.time()
# print(f"time for checking validation of inputs: {end-start} seconds")
# start = time.time()
transactions = []
for s in self.spatial_queries:
if s.dataset.split('_')[0] not in dataset:
continue
cell_pos = s.spatial_pos
labels = np.array(s.labels)
if ct not in np.unique(labels):
continue
ct_pos = cell_pos[labels == ct]
dists, idxs = s.kd_tree.query(ct_pos, k=k + 1, workers=-1)
mask = dists < max_dist
for i, idx in enumerate(idxs):
inds = idx[mask[i]]
transaction = labels[inds[1:]]
transactions.append(transaction)
# end = time.time()
# print(f"time for building {len(transactions)} transactions: {end-start} seconds.")
# start = time.time()
mlb = MultiLabelBinarizer()
encoded_data = mlb.fit_transform(transactions)
df = pd.DataFrame(encoded_data.astype(bool), columns=mlb.classes_)
# end = time.time()
# print(f"time for building df for fpgrowth: {end-start} seconds")
# start = time.time()
fp = fpgrowth(df, min_support=min_support, use_colnames=True)
if len(fp) == 0:
return pd.DataFrame(columns=['support', 'itemsets'])
# end = time.time()
# print(f"time for find fp_growth: {end-start} seconds, {len(fp)} frequent patterns.")
# start = time.time()
fp = find_maximal_patterns(fp=fp)
# end = time.time()
# print(f"time for identify maximal patterns: {end - start} seconds")
fp['itemsets'] = fp['itemsets'].apply(lambda x: list(sorted(x)))
fp.sort_values(by='support', ascending=False, inplace=True, ignore_index=True)
return fp
[docs]
def find_fp_dist(self,
ct: str,
dataset: Union[str, List[str]] = None,
max_dist: float = 20,
min_size: int = 0,
min_support: float = 0.5,
) -> pd.DataFrame:
"""
Find frequent patterns within the radius of certain cell type in multiple fields of view.
Parameters
----------
ct : str
Cell type name.
dataset : str or List[str], optional
Datasets for searching for frequent patterns.
Use all datasets if dataset=None.
max_dist : float, default=20
Maximum distance for considering a cell as a neighbor.
min_size : int, default=0
Minimum neighborhood size for each point to consider.
min_support : float, default=0.5
Threshold of frequency to consider a pattern as a frequent pattern.
Returns
-------
pd.DataFrame
DataFrame with frequent pattern results. Columns include:
- support: frequency of the pattern (proportion of transactions with this pattern)
- itemsets: list of cell types in the frequent pattern
"""
# Search transactions for each field of view, find the frequent patterns of integrated transactions
if_exist_label = [ct in s.labels.unique() for s in self.spatial_queries]
if not any(if_exist_label):
raise ValueError(f"Found no {self.label_key} in any datasets!")
if dataset is None:
dataset = [s.dataset.split('_')[0] for s in self.spatial_queries]
if isinstance(dataset, str):
dataset = [dataset]
# test if the input dataset name is valid
valid_ds_names = [s.dataset.split('_')[0] for s in self.spatial_queries]
for ds in dataset:
if ds not in valid_ds_names:
raise ValueError(f"Invalid input dataset name: {ds}.\n "
f"Valid dataset names are: {set(valid_ds_names)}")
# start = time.time()
transactions = []
for s in self.spatial_queries:
if s.dataset.split('_')[0] not in dataset:
continue
cell_pos = s.spatial_pos
labels = np.array(s.labels)
if ct not in np.unique(labels):
continue
cinds = [id for id, l in enumerate(labels) if l == ct]
ct_pos = cell_pos[cinds]
idxs = s.kd_tree.query_ball_point(ct_pos, r=max_dist, return_sorted=False, workers=-1)
for i_id, idx in zip(cinds, idxs):
transaction = [labels[i] for i in idx if i != i_id]
if len(transaction) > min_size:
transactions.append(transaction)
# end = time.time()
# print(f"time for building {len(transactions)} transactions: {end-start} seconds.")
# start = time.time()
mlb = MultiLabelBinarizer()
encoded_data = mlb.fit_transform(transactions)
df = pd.DataFrame(encoded_data.astype(bool), columns=mlb.classes_)
# end = time.time()
# print(f"time for building df for fpgrowth: {end - start} seconds")
# start = time.time()
fp = fpgrowth(df, min_support=min_support, use_colnames=True)
if len(fp) == 0:
return pd.DataFrame(columns=['support', 'itemsets'])
# end = time.time()
# print(f"time for find fp_growth: {end - start} seconds, {len(fp)} frequent patterns.")
# start = time.time()
fp = find_maximal_patterns(fp=fp)
# end = time.time()
# print(f"time for identify maximal patterns: {end - start} seconds")
fp['itemsets'] = fp['itemsets'].apply(lambda x: list(sorted(x)))
fp.sort_values(by='support', ascending=False, inplace=True, ignore_index=True)
return fp
[docs]
def motif_enrichment_knn(self,
ct: str,
motifs: Union[str, List[str], List[List[str]]] = None,
dataset: Union[str, List[str]] = None,
k: int = 30,
min_support: float = 0.5,
max_dist: float = 20,
return_cellID: bool = False,
) -> Union[pd.DataFrame, Tuple[pd.DataFrame, Dict, Dict]]:
"""
Perform motif enrichment analysis using k-nearest neighbors (KNN) in multiple fields of view.
Parameters
----------
ct : str
The cell type of the center cell.
motifs : str or List[str] or List[List[str]], optional
Specified motifs to be tested.
If motifs=None, find the frequent patterns as motifs within
the neighborhood of center cell type in each fov.
dataset : str or List[str], optional
Datasets for searching for frequent patterns and performing enrichment analysis.
Use all datasets if dataset=None.
k : int, default=30
Number of nearest neighbors to consider.
min_support : float, default=0.5
Threshold of frequency to consider a pattern as a frequent pattern.
max_dist : float, default=20
Maximum distance for neighbors.
return_cellID : bool, default=False
Indicate whether return cell IDs for each frequent pattern within the neighborhood of center cell type and center cells.
By defaults do not return cell ID.
Returns
-------
Union[pd.DataFrame, Tuple[pd.DataFrame, Dict, Dict]]
If return_cellID is False:
pd.DataFrame with motif enrichment results. Columns include:
- center: center cell type name
- motifs: list of cell types in the motif
- n_center_motif: number of center cells with motif in neighborhood
- n_center: total number of center cells
- n_motif: total number of cells with motif in neighborhood
- expectation: expected number under hypergeometric distribution
- p-values: p-value from hypergeometric test
- adj-pval: FDR-corrected p-values (when multiple motifs tested)
- if_significant: whether the enrichment is significant
If return_cellID is True:
Tuple containing:
- pd.DataFrame: enrichment results as described above
- Dict: cell IDs of motifs in center cell's neighborhood
Format: {'motif_str': {'dataset_name': [cell_ids]}}
- Dict: cell IDs of center cells with motif in neighborhood
Format: {'motif_str': {'dataset_name': [cell_ids]}}
"""
if dataset is None:
dataset = [s.dataset.split('_')[0] for s in self.spatial_queries]
if isinstance(dataset, str):
dataset = [dataset]
out = []
if_exist_label = [ct in s.labels.unique() for s in self.spatial_queries]
if not any(if_exist_label):
raise ValueError(f"Found no {self.label_key} in any datasets!")
# Check whether specify motifs. If not, search frequent patterns among specified datasets
# and use them as interested motifs
all_labels = pd.concat([s.labels for s in self.spatial_queries])
labels_unique_all = set(all_labels.unique())
if motifs is None:
fp = self.find_fp_knn(ct=ct, k=k, dataset=dataset,
min_support=min_support, max_dist=max_dist)
motifs = fp['itemsets'].tolist()
else:
if isinstance(motifs, str):
motifs = [[motifs]]
elif isinstance(motifs, list) and all(isinstance(m, str) for m in motifs):
motifs = [motifs]
# else: List[List[str]], keep as is
filtered_motifs = []
for motif in motifs:
motif_exc = [m for m in motif if m not in labels_unique_all]
if len(motif_exc) > 0:
print(f"Found no {motif_exc} in {dataset}. Ignoring them.")
valid_motif = [m for m in motif if m in labels_unique_all]
if len(valid_motif) > 0:
filtered_motifs.append(valid_motif)
if len(filtered_motifs) == 0:
raise ValueError(f"All cell types in motifs are missed in {self.label_key}.")
motifs = filtered_motifs
# Initialize dictionaries to store cell IDs if requested
motif_cell_ids = {}
center_cell_ids = {}
for motif in motifs:
n_labels = 0
n_ct = 0
n_motif_labels = 0
n_motif_ct = 0
motif = list(motif) if not isinstance(motif, list) else motif
sort_motif = sorted(motif)
if return_cellID:
motif_str = str(sort_motif)
motif_cell_ids[motif_str] = {}
center_cell_ids[motif_str] = {}
# Calculate statistics of each dataset
for fov, s in enumerate(self.spatial_queries):
if s.dataset.split('_')[0] not in dataset:
continue
cell_pos = s.spatial_pos
labels = np.array(s.labels)
labels_unique = np.unique(labels)
contain_motif = [m in labels_unique for m in motif]
if not np.all(contain_motif):
n_labels += labels.shape[0]
n_ct += np.sum(labels == ct)
continue
else:
n_labels += labels.shape[0]
label_encoder = LabelEncoder()
int_labels = label_encoder.fit_transform(labels)
int_motifs = label_encoder.transform(np.array(motif))
dists, idxs = s.kd_tree.query(cell_pos, k=k + 1, workers=-1)
num_cells = idxs.shape[0]
num_types = len(label_encoder.classes_)
valid_neighbors = dists[:, 1:] <= max_dist
filtered_idxs = np.where(valid_neighbors, idxs[:, 1:], -1)
flat_neighbors = filtered_idxs.flatten()
valid_neighbors_flat = valid_neighbors.flatten()
neighbor_labels = np.where(valid_neighbors_flat, int_labels[flat_neighbors], -1)
valid_mask = neighbor_labels != -1
neighbor_matrix = np.zeros((num_cells * k, num_types), dtype=int)
neighbor_matrix[np.arange(len(neighbor_labels))[valid_mask], neighbor_labels[valid_mask]] = 1
neighbor_counts = neighbor_matrix.reshape(num_cells, k, num_types).sum(axis=1)
# Check which cells have all motif types in their neighborhood
motif_mask = np.all(neighbor_counts[:, int_motifs] > 0, axis=1)
n_motif_labels += np.sum(motif_mask)
if ct in np.unique(labels):
int_ct = label_encoder.transform(np.array(ct, dtype=object, ndmin=1))
mask = int_labels == int_ct
center_mask = mask & motif_mask
n_motif_ct += np.sum(center_mask)
n_ct += np.sum(mask)
if return_cellID:
# Get IDs of center cells with motif in neighborhood
center_indices = np.where(center_mask)[0]
if len(center_indices) > 0:
# Get all neighbors of center cells that have the motif
center_neighbors_idxs = filtered_idxs[center_mask]
# Flatten and get unique neighbors
all_neighbors_center = center_neighbors_idxs.flatten()
all_neighbors_center = all_neighbors_center[all_neighbors_center != -1]
# Filter to only keep neighbors that are of motif cell types
motif_mask_all = np.isin(np.array(s.labels), motif)
valid_neighbors_center = all_neighbors_center[motif_mask_all[all_neighbors_center]]
id_motif_celltype = set(valid_neighbors_center)
motif_cell_ids[motif_str][s.dataset] = list(id_motif_celltype)
center_cell_ids[motif_str][s.dataset] = list(center_indices)
if ct in motif:
n_ct = round(n_ct / motif.count(ct))
hyge = hypergeom(M=n_labels, n=n_ct, N=n_motif_labels)
motif_out = {'center': ct, 'motifs': sort_motif, 'n_center_motif': n_motif_ct,
'n_center': n_ct, 'n_motif': n_motif_labels, 'expectation': hyge.mean(), 'p-values': hyge.sf(n_motif_ct)}
out.append(motif_out)
out_pd = pd.DataFrame(out)
if len(out_pd) == 1:
out_pd['if_significant'] = True if out_pd['p-values'][0] < 0.05 else False
else:
p_values = out_pd['p-values'].tolist()
if_rejected, corrected_p_values = mt.fdrcorrection(p_values,
alpha=0.05,
method='poscorr')
out_pd['adj-pval'] = corrected_p_values
out_pd['if_significant'] = if_rejected
out_pd = out_pd.sort_values(by='adj-pval', ignore_index=True)
if return_cellID:
return out_pd, motif_cell_ids, center_cell_ids
else:
return out_pd
[docs]
def motif_enrichment_dist(self,
ct: str,
motifs: Union[str, List[str], List[List[str]]] = None,
dataset: Union[str, List[str]] = None,
max_dist: float = 20,
min_size: int = 0,
min_support: float = 0.5,
return_cellID: bool = False
) -> Union[pd.DataFrame, Tuple[pd.DataFrame, Dict, Dict]]:
"""
Perform motif enrichment analysis within a specified radius-based neighborhood in multiple fields of view.
Parameters
----------
ct : str
Cell type of the center cell.
motifs : str or List[str] or List[List[str]], optional
Specified motifs to be tested.
If motifs=None, find the frequent patterns as motifs within the neighborhood of center cell type.
dataset : str or List[str], optional
Datasets for searching for frequent patterns and performing enrichment analysis.
Use all datasets if dataset=None.
max_dist : float, default=20
Maximum distance for considering a cell as a neighbor.
min_size : int, default=0
Minimum neighborhood size for each point to consider.
min_support : float, default=0.5
Threshold of frequency to consider a pattern as a frequent pattern.
return_cellID : bool, default=False
Indicate whether return cell IDs for each frequent pattern within the neighborhood of center cells.
By defaults do not return cell ID.
Returns
-------
Union[pd.DataFrame, Tuple[pd.DataFrame, Dict, Dict]]
If return_cellID is False:
pd.DataFrame with motif enrichment results. Columns include:
- center: center cell type name
- motifs: list of cell types in the motif
- n_center_motif: number of center cells with motif in neighborhood
- n_center: total number of center cells
- n_motif: total number of cells with motif in neighborhood
- expectation: expected number under hypergeometric distribution
- p-values: p-value from hypergeometric test
- adj-pval: FDR-corrected p-values (when multiple motifs tested)
- if_significant: whether the enrichment is significant
If return_cellID is True:
Tuple containing:
- pd.DataFrame: enrichment results as described above
- Dict: cell IDs of motifs in center cell's neighborhood
Format: {'motif_str': {'dataset_name': [cell_ids]}}
- Dict: cell IDs of center cells with motif in neighborhood
Format: {'motif_str': {'dataset_name': [cell_ids]}}
"""
if dataset is None:
dataset = [s.dataset.split('_')[0] for s in self.spatial_queries]
if isinstance(dataset, str):
dataset = [dataset]
out = []
if_exist_label = [ct in s.labels.unique() for s in self.spatial_queries]
if not any(if_exist_label):
raise ValueError(f"Found no {self.label_key} in any datasets!")
# Check whether specify motifs. If not, search frequent patterns among specified datasets
# and use them as interested motifs
# Properly handle different input formats for motifs
if motifs is None:
# If motifs is None, keep the existing logic to find patterns
fp = self.find_fp_dist(ct=ct, dataset=dataset, max_dist=max_dist, min_size=min_size,
min_support=min_support, )
motifs = fp['itemsets'].tolist()
else:
if isinstance(motifs, str):
motifs = [[motifs]]
elif isinstance(motifs, list) and all(isinstance(m, str) for m in motifs):
motifs = [motifs]
# At this point, motifs should be list[list[str]]
# Filter out undefined cell types from each motif
all_labels = pd.concat([s.labels for s in self.spatial_queries])
labels_unique_all = set(all_labels.unique())
filtered_motifs = []
for motif in motifs:
# Check which cell types in this motif are valid
motif_exc = [m for m in motif if m not in labels_unique_all]
if len(motif_exc) > 0:
print(f"Not found {motif_exc} in {dataset}! Ignoring them.")
# Filter the current motif to only include valid cell types
valid_motif = [m for m in motif if m in labels_unique_all]
# Only include this motif if it has at least one valid cell type
if len(valid_motif) > 0:
filtered_motifs.append(valid_motif)
# Check if we have any valid motifs left
if len(filtered_motifs) == 0:
raise ValueError(f"All cell types in motifs are missing in {self.label_key}.")
motifs = filtered_motifs
# Initialize dictionaries to store cell IDs if requested
motif_cell_ids = {}
center_cell_ids = {}
for motif in motifs:
n_labels = 0
n_ct = 0
n_motif_labels = 0
n_motif_ct = 0
motif = list(motif) if not isinstance(motif, list) else motif
sort_motif = sorted(motif)
if return_cellID:
motif_str = str(sort_motif)
motif_cell_ids[motif_str] = {}
center_cell_ids[motif_str] = {}
for s in self.spatial_queries:
if s.dataset.split('_')[0] not in dataset:
continue
cell_pos = s.spatial_pos
labels = np.array(s.labels)
labels_unique = np.unique(labels)
contain_motif = [m in labels_unique for m in motif]
if not np.all(contain_motif):
n_labels += labels.shape[0]
n_ct += np.sum(labels == ct)
continue
else:
n_labels += labels.shape[0]
# Query neighbors for all cells once (instead of using grid filtering)
idxs_all = s.kd_tree.query_ball_point(
s.spatial_pos,
r=max_dist,
return_sorted=False,
workers=-1,
)
idxs_all_filter = [np.array(ids)[np.array(ids) != i] for i, ids in enumerate(idxs_all)]
# using numpy
label_encoder = LabelEncoder()
int_labels = label_encoder.fit_transform(labels)
int_motifs = label_encoder.transform(np.array(motif))
num_cells = len(s.spatial_pos)
num_types = len(label_encoder.classes_)
# Pre-compute neighbor matrix for all cells
flat_neighbors_all = np.concatenate(idxs_all_filter)
row_indices_all = np.repeat(np.arange(num_cells), [len(neigh) for neigh in idxs_all_filter])
neighbor_labels_all = int_labels[flat_neighbors_all]
neighbor_matrix_all = np.zeros((num_cells, num_types), dtype=int)
np.add.at(neighbor_matrix_all, (row_indices_all, neighbor_labels_all), 1)
# Check which cells have all motif types in their neighborhood
motif_mask = np.all(neighbor_matrix_all[:, int_motifs] > 0, axis=1)
n_motif_labels += np.sum(motif_mask)
if ct in np.unique(labels):
int_ct = label_encoder.transform(np.array(ct, dtype=object, ndmin=1))
mask_ct = int_labels == int_ct
center_mask = mask_ct & motif_mask
n_motif_ct += np.sum(center_mask)
n_ct += np.sum(s.labels == ct)
if return_cellID:
# Get IDs of center cells with motif in neighborhood
center_indices = np.where(center_mask)[0]
if len(center_indices) > 0:
# Get motif neighbors for these center cells
motif_mask_all = np.isin(np.array(s.labels), motif)
all_neighbors_center = np.concatenate([idxs_all_filter[i] for i in center_indices])
valid_neighbors_center = all_neighbors_center[motif_mask_all[all_neighbors_center]]
id_motif_celltype = set(valid_neighbors_center)
motif_cell_ids[motif_str][s.dataset] = list(id_motif_celltype)
center_cell_ids[motif_str][s.dataset] = list(center_indices)
if ct in motif:
n_ct = round(n_ct / motif.count(ct))
hyge = hypergeom(M=n_labels, n=n_ct, N=n_motif_labels)
motif_out = {'center': ct, 'motifs': sort_motif, 'n_center_motif': n_motif_ct,
'n_center': n_ct, 'n_motif': n_motif_labels, 'expectation': hyge.mean(), 'p-values': hyge.sf(n_motif_ct)}
out.append(motif_out)
out_pd = pd.DataFrame(out)
if len(out_pd) == 1:
out_pd['if_significant'] = True if out_pd['p-values'][0] < 0.05 else False
else:
p_values = out_pd['p-values'].tolist()
if_rejected, corrected_p_values = mt.fdrcorrection(p_values,
alpha=0.05,
method='poscorr')
out_pd['adj-pval'] = corrected_p_values
out_pd['if_significant'] = if_rejected
out_pd = out_pd.sort_values(by='adj-pval', ignore_index=True)
if return_cellID:
return out_pd, motif_cell_ids, center_cell_ids
else:
return out_pd
[docs]
def find_fp_knn_fov(self,
ct: str,
dataset_i: str,
k: int = 30,
min_support: float = 0.5,
max_dist: float = 20
) -> pd.DataFrame:
"""
Find frequent patterns within the KNNs of specific cell type of interest in single field of view.
Parameters
----------
ct : str
Cell type name.
dataset_i : str
Dataset name in dataset_i format.
k : int, default=30
Number of nearest neighbors.
min_support : float, default=0.5
Threshold of frequency to consider a pattern as a frequent pattern.
max_dist : float, default=20
Maximum distance for considering a cell as a neighbor.
Returns
-------
pd.DataFrame
DataFrame with frequent pattern results. Columns include:
- support: frequency of the pattern
- itemsets: frozenset of cell types in the frequent pattern
"""
if dataset_i not in self.datasets:
raise ValueError(f"Found no {dataset_i.split('_')[0]} in any datasets.")
sp_object = self.spatial_queries[self.datasets.index(dataset_i)]
cell_pos = sp_object.spatial_pos
labels = np.array(sp_object.labels)
if ct not in np.unique(labels):
return pd.DataFrame(columns=['support', 'itemsets'])
ct_pos = cell_pos[labels == ct]
# Identify frequent patterns of cell types, including those subsets of patterns
# whose support value exceeds min_support. Focus solely on the multiplicity
# of cell types, rather than their frequency.
fp, _, _ = build_fptree_knn(
kd_tree=sp_object.kd_tree,
labels=labels,
cell_pos=ct_pos,
spatial_pos=cell_pos,
k=k,
min_support=min_support,
if_max=True,
max_dist=max_dist,
)
return fp
[docs]
def find_fp_dist_fov(self,
ct: str,
dataset_i: str,
max_dist: float = 20,
min_size: int = 0,
min_support: float = 0.5,
) -> pd.DataFrame:
"""
Find frequent patterns within the radius-based neighborhood of specific cell type of interest
in single field of view.
Parameters
----------
ct : str
Cell type name.
dataset_i : str
Dataset name in dataset_i format.
max_dist : float, default=20
Maximum distance for considering a cell as a neighbor.
min_size : int, default=0
Minimum neighborhood size for each point to consider.
min_support : float, default=0.5
Threshold of frequency to consider a pattern as a frequent pattern.
Returns
-------
pd.DataFrame
DataFrame with frequent pattern results. Columns include:
- support: frequency of the pattern
- itemsets: frozenset of cell types in the frequent pattern
"""
if dataset_i not in self.datasets:
raise ValueError(f"Found no {dataset_i.split('_')[0]} in any datasets.")
sp_object = self.spatial_queries[self.datasets.index(dataset_i)]
cell_pos = sp_object.spatial_pos
labels = sp_object.labels
if ct not in labels.unique():
return pd.DataFrame(columns=['support, itemsets'])
cinds = [id for id, l in enumerate(labels) if l == ct]
ct_pos = cell_pos[cinds]
fp, _, _ = build_fptree_dist(
kd_tree=sp_object.kd_tree,
labels=labels,
cell_pos=ct_pos,
spatial_pos=cell_pos,
max_dist=max_dist,
min_support=min_support,
if_max=True,
min_size=min_size,
cinds=cinds,
)
return fp
[docs]
def differential_analysis_knn(self,
ct: str,
datasets: List[str],
motifs: Optional[Union[str, List[str], List[List[str]]]] = None,
k: int = 30,
min_support: float = 0.5,
max_dist: float = 20,
) -> dict:
"""
Perform differential analysis of spatial motif patterns between two datasets using KNN neighborhood.
Mann-Whitney U test is used to compare the support values of each motif pattern across FOVs of two datasets,
and FDR correction is applied for multiple testing.
This function identifies motif patterns that are differentially enriched in the KNN neighborhood
of a center cell type between two conditions (e.g., disease vs control). It supports two modes:
- **Unbiased discovery** (``motifs=None``): Automatically discovers frequent patterns in each
FOV of both datasets, then tests for differential enrichment.
- **Hypothesis-driven** (``motifs`` specified): Tests user-specified motifs for differential
enrichment, allowing validation of known or hypothesized spatial patterns.
Parameters
----------
ct : str
Cell type of interest as center point.
datasets : List[str]
List of exactly 2 dataset names to compare (e.g., ['Disease', 'Control']).
motifs : str or List[str] or List[List[str]], optional
Optional user-specified motif(s) to test. Can be:
- Single cell type: 'CellTypeA'
- Single motif: ['CellTypeA', 'CellTypeB']
- Multiple motifs: [['CellTypeA'], ['CellTypeB', 'CellTypeC']]
If None, performs unbiased discovery of frequent patterns first.
k : int, default=30
Number of nearest neighbors to consider.
min_support : float, default=0.5
Threshold of frequency to consider a pattern as frequent (only used when motifs=None
for unbiased discovery mode).
max_dist : float, default=20
Maximum distance for considering a cell as a neighbor.
Returns
-------
dict
Dictionary with dataset names as keys and DataFrames as values. Each DataFrame
contains motif patterns significantly enriched in that dataset, with columns:
itemsets (motif as tuple), support_{datasets[0]}_mean, support_{datasets[1]}_mean,
adj-pval. Only patterns with adj-pval < 0.05 are included.
"""
if len(datasets) != 2:
raise ValueError("Require 2 datasets for differential analysis.")
# Check if the two datasets are valid
valid_ds_names = [s.dataset.split('_')[0] for s in self.spatial_queries]
for ds in datasets:
if ds not in valid_ds_names:
raise ValueError(f"Invalid input dataset name: {ds}.\n"
f"Valid dataset names are: {set(valid_ds_names)}")
# If motifs are specified, use the new function for user-specified motifs
if motifs is not None:
return differential_analysis_motif_knn(
spatial_queries=self.spatial_queries,
datasets_list=datasets,
ct=ct,
motifs=motifs,
k=k,
max_dist=max_dist,
)
# Otherwise, use the original flow to discover frequent patterns first
# Step 1: Collect all frequent patterns discovered across all FOVs in both datasets
all_discovered_motifs = set()
for d in datasets:
dataset_i = [ds for ds in self.datasets if ds.split('_')[0] == d]
for d_i in dataset_i:
fp_fov = self.find_fp_knn_fov(ct=ct,
dataset_i=d_i,
k=k,
min_support=min_support,
max_dist=max_dist)
if len(fp_fov) > 0:
# Collect all discovered motifs (use union instead of intersection)
for itemset in fp_fov['itemsets']:
all_discovered_motifs.add(tuple(sorted(itemset)))
# Convert to list of lists for the motif-based analysis
all_motifs = [list(motif) for motif in all_discovered_motifs]
# Step 2: Use the motif-specified pathway to compute true support values for all discovered motifs
# This ensures consistency with the specified-motif pathway
return differential_analysis_motif_knn(
spatial_queries=self.spatial_queries,
datasets_list=datasets,
ct=ct,
motifs=all_motifs,
k=k,
max_dist=max_dist,
)
[docs]
def differential_analysis_dist(self,
ct: str,
datasets: List[str],
motifs: Optional[Union[str, List[str], List[List[str]]]] = None,
max_dist: float = 20,
min_support: float = 0.5,
min_size: int = 0,
) -> dict:
"""
Perform differential analysis of spatial motif patterns between two datasets using radius-based neighborhood.
Mann-Whitney U test is used to compare the support values of each motif pattern across FOVs of two datasets,
and FDR correction is applied for multiple testing.
This function identifies motif patterns that are differentially enriched in the radius-based
neighborhood of a center cell type between two conditions (e.g., disease vs control). It supports
two modes:
- **Unbiased discovery** (``motifs=None``): Automatically discovers frequent patterns in each
FOV of both datasets, then tests for differential enrichment.
- **Hypothesis-driven** (``motifs`` specified): Tests user-specified motifs for differential
enrichment, allowing validation of known or hypothesized spatial patterns.
Parameters
----------
ct : str
Cell type of interest as center point.
datasets : List[str]
List of exactly 2 dataset names to compare (e.g., ['Disease', 'Control']).
motifs : str or List[str] or List[List[str]], optional
Optional user-specified motif(s) to test. Can be:
- Single cell type: 'CellTypeA'
- Single motif: ['CellTypeA', 'CellTypeB']
- Multiple motifs: [['CellTypeA'], ['CellTypeB', 'CellTypeC']]
If None, performs unbiased discovery of frequent patterns first.
max_dist : float, default=20
Maximum distance for considering a cell as a neighbor.
min_support : float, default=0.5
Threshold of frequency to consider a pattern as frequent (only used when motifs=None
for unbiased discovery mode).
min_size : int, default=0
Minimum neighborhood size for each center cell to be considered.
Returns
-------
dict
Dictionary with dataset names as keys and DataFrames as values. Each DataFrame
contains motif patterns significantly enriched in that dataset, with columns:
itemsets (motif as tuple), support_{datasets[0]}_mean, support_{datasets[1]}_mean,
adj-pval. Only patterns with adj-pval < 0.05 are included.
"""
if len(datasets) != 2:
raise ValueError("Require 2 datasets for differential analysis.")
# Check if the two datasets are valid
valid_ds_names = [s.dataset.split('_')[0] for s in self.spatial_queries]
for ds in datasets:
if ds not in valid_ds_names:
raise ValueError(f"Invalid input dataset name: {ds}.\n"
f"Valid dataset names are: {set(valid_ds_names)}")
# If motifs are specified, use the new function for user-specified motifs
if motifs is not None:
return differential_analysis_motif_dist(
spatial_queries=self.spatial_queries,
datasets_list=datasets,
ct=ct,
motifs=motifs,
max_dist=max_dist,
min_size=min_size,
)
# Otherwise, use the original flow to discover frequent patterns first
# Step 1: Collect all frequent patterns discovered across all FOVs in both datasets
all_discovered_motifs = set()
for d in datasets:
dataset_i = [ds for ds in self.datasets if ds.split('_')[0] == d]
for d_i in dataset_i:
fp_fov = self.find_fp_dist_fov(ct=ct,
dataset_i=d_i,
max_dist=max_dist,
min_size=min_size,
min_support=min_support,
)
if len(fp_fov) > 0:
# Collect all discovered motifs (use union instead of intersection)
for itemset in fp_fov['itemsets']:
all_discovered_motifs.add(tuple(sorted(itemset)))
# Convert to list of lists for the motif-based analysis
all_motifs = [list(motif) for motif in all_discovered_motifs]
print(f"Discovered {len(all_motifs)} motifs across the datasets for differential analysis.")
# Step 2: Use the motif-specified pathway to compute true support values for all discovered motifs
# This ensures consistency with the specified-motif pathway
return differential_analysis_motif_dist(
spatial_queries=self.spatial_queries,
datasets_list=datasets,
ct=ct,
motifs=all_motifs,
max_dist=max_dist,
min_size=min_size,
)
[docs]
def de_genes(self,
ind_group1: Dict[str, List[int]],
ind_group2: Dict[str, List[int]],
genes: Optional[Union[str, List[str]]] = None,
min_fraction: float = 0.05,
method: Literal['fisher', 't-test', 'wilcoxon'] = 'fisher',
alpha: Optional[float] = None,
) -> pd.DataFrame:
"""
Perform differential expression analysis on the given indices.
The ind_group1 and ind_group2 should be a defaultdict with keys as modified dataset names and values as
lists of indices in corresponding group.
It provides a flexible way to perform DE analysis on different datasets, e.g., across different FOVs of the
same condition, or across FOVs from different conditions.
Parameters
----------
ind_group1 : defaultdict[str, List[int]]
A defaultdict with keys as modified dataset names and values as lists of indices in corresponding group.
ind_group2 : defaultdict[str, List[int]]
A defaultdict with keys as modified dataset names and values as lists of indices in corresponding group.
genes : str or List[str], optional
Genes to be searched in the gene index.
min_fraction : float, default=0.05
The minimum fraction of cells that express a gene for it to be considered differentially expressed.
method : {'fisher', 't-test', 'wilcoxon'}, default='fisher'
The method to use for DE analysis. If build_gene_index=True, only Fisher's exact test is supported.
alpha : float, optional
Significance threshold for adjusted p-values. If None, defaults to 0.1 when using Fisher's exact test and 0.05 otherwise.
Returns
-------
pd.DataFrame
DataFrame with differential expression results. Columns include:
- gene: gene name
- proportion_1: proportion of cells expressing the gene in group 1
- proportion_2: proportion of cells expressing the gene in group 2
- p_value: p-value from statistical test
- adj-pval: FDR-corrected p-value
- de_in: which group the gene is differentially expressed in ('group1' or 'group2')
"""
if self.build_gene_index:
# Use scfind index-based method with Fisher's exact test
if method != 'fisher':
print(f"Warning: When build_gene_index=True, only Fisher's exact test is supported. Ignoring method='{method}'.")
return self._de_genes_scfind(ind_group1, ind_group2, genes, min_fraction, alpha)
else:
# Use adata.X directly with specified method
return self._de_genes_adata(ind_group1, ind_group2, genes, min_fraction, method, alpha)
def _de_genes_scfind(self,
ind_group1: Dict[str, List[int]],
ind_group2: Dict[str, List[int]],
genes: Optional[Union[str, List[str]]] = None,
min_fraction: float = 0.05,
alpha: Optional[float] = None
) -> pd.DataFrame:
"""
Perform differential expression analysis using scfind index with Fisher's exact test.
This method is used when build_gene_index=True.
"""
# For each gene, calculate the number of cells in the provided indices expressing the gene in each group
if genes is None:
genes = set.intersection(*[set(s.genes) for s in self.spatial_queries])
genes = list(genes)
print(f"Testing {len(genes)} genes with Fisher's exact test ...\n")
if alpha is None:
alpha = 0.1
n_1 = np.sum([len(ids) for ids in ind_group1.values()])
n_2 = np.sum([len(ids) for ids in ind_group2.values()])
valid_ds1 = [ds for ds in ind_group1.keys() if ds in self.datasets]
valid_ds2 = [ds for ds in ind_group2.keys() if ds in self.datasets]
# Check if there are valid datasets
if not valid_ds1:
raise ValueError("No valid datasets found in ind_group1.")
if not valid_ds2:
raise ValueError("No valid datasets found in ind_group2.")
group1_results = []
start = time()
for ds, ids in ind_group1.items():
print(f"Processing {ds} in group1...")
if ds not in valid_ds1:
print(f'Warning: {ds} is not a valid dataset name. Ignoring it.')
continue
ds_i = self.datasets.index(ds)
sp = self.spatial_queries[ds_i]
# Get counts of cells expressing each gene in this dataset
start1 = time()
genes_sp = sp.index._case_correct(genes, if_print=False)
if not genes_sp:
continue
ds_counts = sp.index.index.cell_counts_in_indices_genes(ids, genes_sp)
end1 = time()
print(f'Cell search in dataset {ds}: {end1 - start1:.2f} seconds')
genes_list = [item['gene'] for item in ds_counts]
counts_list = [item['expressed_cells'] for item in ds_counts]
# Create a DataFrame
if genes_list:
temp_df = pd.DataFrame({'gene': genes_list, 'count': counts_list})
group1_results.append(temp_df)
end = time()
print(f'Time for cell search in group1: {end - start}')
# Count cells expressing each gene in group 2
group2_results = []
start = time()
for ds, ids in ind_group2.items():
print(f"Processing {ds} in group2...")
if ds not in valid_ds2:
print(f'Warning: {ds} is not a valid dataset name. Ignoring it.')
continue
ds_i = self.datasets.index(ds)
sp = self.spatial_queries[ds_i]
# Get counts of cells expressing each gene in this dataset
genes_sp = sp.index._case_correct(genes, if_print=False)
if not genes_sp:
continue
ds_counts = sp.index.index.cell_counts_in_indices_genes(ids, genes_sp)
genes_list = [item['gene'] for item in ds_counts]
counts_list = [item['expressed_cells'] for item in ds_counts]
# Create a DataFrame in one operation
if genes_list:
temp_df = pd.DataFrame({'gene': genes_list, 'count': counts_list})
group2_results.append(temp_df)
end = time()
print(f'Time for cell search in group2: {end - start:.2f} seconds')
# Prepare data for statistical testing
# Combine all results
# start = time()
if group1_results:
group1_df = pd.concat(group1_results, ignore_index=True)
group1_agg = group1_df.groupby('gene')['count'].sum().reset_index()
group1_agg = group1_agg.rename(columns={'count': 'count_1'})
else:
group1_agg = pd.DataFrame(columns=['gene', 'count_1'])
if group2_results:
group2_df = pd.concat(group2_results, ignore_index=True)
group2_agg = group2_df.groupby('gene')['count'].sum().reset_index()
group2_agg = group2_agg.rename(columns={'count': 'count_2'})
else:
group2_agg = pd.DataFrame(columns=['gene', 'count_2'])
# Merge the two groups
merged_df = pd.merge(group1_agg, group2_agg, on='gene', how='outer').fillna(0)
# Calculate proportions
merged_df['proportion_1'] = merged_df['count_1'] / n_1
merged_df['proportion_2'] = merged_df['count_2'] / n_2
# Filter by minimum fraction
filtered_df = merged_df[(merged_df['proportion_1'] >= min_fraction) |
(merged_df['proportion_2'] >= min_fraction)].copy()
if filtered_df.empty:
print("No genes meet the minimum fraction threshold.")
return pd.DataFrame(
columns=["gene", "proportion_1", "proportion_2", "abs",
"difference", "p_value", "adj-pval", "de_in"]
)
# Calculate differences
filtered_df.loc[:, 'difference'] = filtered_df['proportion_1'] - filtered_df['proportion_2']
filtered_df.loc[:, 'abs'] = filtered_df['difference'].abs()
# For Fisher's exact test, prepare arrays for vectorized operations
count_1_array = filtered_df['count_1'].values.astype(int)
count_2_array = filtered_df['count_2'].values.astype(int)
not_count_1_array = n_1 - count_1_array
not_count_2_array = n_2 - count_2_array
# end = time()
# print(f'Time for prepare data for preparation of statistical testing: {end - start}')
# Use numpy to create all contingency tables at once
# This creates a 3D array of shape (n_rows, 2, 2)
# start = time()
contingency_tables = np.array([
[[count_1_array[i], not_count_1_array[i]],
[count_2_array[i], not_count_2_array[i]]]
for i in range(len(count_1_array))
])
# Apply Fisher's exact test - this still needs a loop but is more efficient
# Use parallelization if available (requires joblib)
def apply_fisher(table):
_, p_value = stats.fisher_exact(table, alternative='two-sided')
return p_value
# Run tests in parallel
# n_jobs=-1 uses all available cores
# p_values_parallel = Parallel(n_jobs=-1)(
# delayed(apply_fisher)(table) for table in contingency_tables
# )
# Also compute p_values sequentially (without Parallel)
p_values = [apply_fisher(table) for table in contingency_tables]
# Add p-values to DataFrame
filtered_df.loc[:, 'p_value'] = p_values
# Sort by p-value
filtered_df = filtered_df.sort_values('p_value')
# Multiple testing correction
if len(filtered_df) > 1:
adjusted_pvals = multipletests(filtered_df['p_value'], method='fdr_bh')[1]
filtered_df['adj-pval'] = adjusted_pvals
else:
filtered_df['adj-pval'] = filtered_df['p_value']
filtered_df = filtered_df[filtered_df['adj-pval']<alpha].reset_index(drop=True)
# Add information about which group shows higher expression
filtered_df['de_in'] = np.where(
(filtered_df['proportion_1'] > filtered_df['proportion_2']),
'group1',
np.where(
(filtered_df['proportion_2'] > filtered_df['proportion_1']),
'group2',
None
)
)
# end = time()
# print(f'Time for statistical testing: {end - start}')
# Return the final results
return filtered_df[["gene", "proportion_1", "proportion_2", "abs",
"difference", "p_value", "adj-pval", "de_in"]]
def _de_genes_adata(self,
ind_group1: Dict[str, List[int]],
ind_group2: Dict[str, List[int]],
genes: Optional[Union[str, List[str]]] = None,
min_fraction: float = 0.05,
method: Literal['fisher', 't-test', 'wilcoxon'] = 'fisher',
alpha: Optional[float] = None
) -> pd.DataFrame:
"""
Perform differential expression analysis using adata.X directly.
This method is used when build_gene_index=False.
For each dataset, collect the cells and concatenate them, then perform DE analysis.
"""
# Validate datasets
valid_ds1 = [ds for ds in ind_group1.keys() if ds in self.datasets]
valid_ds2 = [ds for ds in ind_group2.keys() if ds in self.datasets]
if not valid_ds1:
raise ValueError("No valid datasets found in ind_group1.")
if not valid_ds2:
raise ValueError("No valid datasets found in ind_group2.")
if alpha is None:
alpha = 0.05
# Collect all cells from group 1
all_adatas_g1 = []
for ds in valid_ds1:
if ds not in ind_group1 or len(ind_group1[ds]) == 0:
continue
ds_i = self.datasets.index(ds)
sp = self.spatial_queries[ds_i]
# Get the adata for this dataset
if sp.adata is None:
raise ValueError(f"Error: {ds} does not have adata.X. Please use use fisher's exact using indexed data.")
# Extract cells for group 1
idx_g1 = ind_group1[ds]
adata_subset = sp.adata[idx_g1].copy()
all_adatas_g1.append(adata_subset)
# Collect all cells from group 2
all_adatas_g2 = []
for ds in valid_ds2:
if ds not in ind_group2 or len(ind_group2[ds]) == 0:
continue
ds_i = self.datasets.index(ds)
sp = self.spatial_queries[ds_i]
# Get the adata for this dataset
if sp.adata is None:
raise ValueError(f"Error: {ds} does not have adata.X. Please use use fisher's exact using indexed data.")
# Extract cells for group 2
idx_g2 = ind_group2[ds]
adata_subset = sp.adata[idx_g2].copy()
all_adatas_g2.append(adata_subset)
if not all_adatas_g1:
raise ValueError("No valid adata found in group 1.")
if not all_adatas_g2:
raise ValueError("No valid adata found in group 2.")
adata_g1 = ad.concat(all_adatas_g1, join='inner')
adata_g2 = ad.concat(all_adatas_g2, join='inner')
# Combine both groups
adata_combined = ad.concat([adata_g1, adata_g2], join='inner')
# Get the overlapping genes of each data
genes_list = adata_combined.var_names.tolist()
# Create indices for combined adata
ind_combined_g1 = list(range(len(adata_g1)))
ind_combined_g2 = list(range(len(adata_g1), len(adata_combined)))
# Perform DE analysis using spatial_utils
if method == 'fisher':
results_df = de_genes_fisher(
adata_combined, genes_list, ind_combined_g1, ind_combined_g2, genes, min_fraction, alpha
)
elif method == 't-test' or method == 'wilcoxon':
results_df = de_genes_scanpy(
adata_combined, genes_list, ind_combined_g1, ind_combined_g2, genes, min_fraction, method=method, alpha=alpha
)
else:
raise ValueError(f"Invalid method: {method}. Choose from 'fisher', 't-test', or 'wilcoxon'.")
return results_df
[docs]
def plot_cell_type_distribution(
self,
dataset: Optional[Union[str, List[str]]] = None,
data_type: Literal['number', 'proportion'] = 'proportion',
colormap: str = 'tab20c',
save_path: Optional[str] = None
):
"""
Visualize the distribution of cell types across datasets using a stacked bar plot.
Parameters
----------
dataset : str or List[str], optional
Datasets for searching.
data_type : {'number', 'proportion'}, default='proportion'
Plot bar plot by number of cells or by proportions.
colormap : str, default='tab20c'
Matplotlib colormap name.
save_path : str, optional
Path to save the figure.
Returns
-------
None
Displays a stacked bar plot.
"""
if data_type not in ['number', 'proportion']:
raise ValueError("Invalid data_type. It should be one of 'number' or 'proportion'.")
if dataset is None:
dataset = [s.dataset.split('_')[0] for s in self.spatial_queries]
if isinstance(dataset, str):
dataset = [dataset]
summary = defaultdict(lambda: defaultdict(int))
valid_queries = [s for s in self.spatial_queries if s.dataset.split('_')[0] in dataset]
cell_types = set(ct for s in valid_queries for ct in s.labels.unique())
for s in valid_queries:
for cell_type in cell_types:
summary[s.dataset][cell_type] += np.sum(s.labels == cell_type)
df = pd.DataFrame(
[(ds, ct, cnt) for ds, cts in summary.items() for ct, cnt in cts.items()],
columns=['Dataset', 'Cell Type', 'Count']
)
df['dataset'] = df['Dataset'].str.split('_').str[0]
summary = (
df.groupby(['dataset', 'Cell Type'])['Count']
.sum()
.reset_index()
)
plot_data = (
summary
.pivot(index='Cell Type', columns='dataset', values='Count')
.fillna(0)
)
plot_data = plot_data.T
# sort columns by cell counts
plot_data = plot_data.loc[
plot_data.sum(axis=1).sort_values(ascending=False).index
]
if data_type != 'number':
plot_data = plot_data.div(plot_data.sum(axis=1), axis=0)
n_cols = df['Cell Type'].unique().shape[0] # bar 中 stack 的数量
cmap = cm.get_cmap(colormap)
colors = [cmap(i) for i in range(n_cols)]
n_x = plot_data.shape[0]
fig_w = max(6, n_x * 0.5)
fig_h = fig_w * 0.5
fig, ax = plt.subplots(figsize=(fig_w, fig_h))
plot_data.plot(
kind='bar',
stacked=True,
ax=ax,
edgecolor='black',
color=colors,
width=0.75,
)
ax.set_title("Distribution of cell types across datasets", fontsize=16)
ax.set_xlabel("Dataset", fontsize=12)
ax.set_ylabel(
"Number of Cells" if data_type == "number" else "Proportion of Cells",
fontsize=12,
)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.tick_params(axis='x', rotation=45, labelsize=10)
ax.legend(
title="Cell type",
loc='center left',
bbox_to_anchor=(1, 0.5),
fontsize=10,
frameon=False,
)
# plt.tight_layout(rect=[0, 0, 0.85, 1])
if save_path is not None:
plt.savefig(save_path, dpi=300, bbox_inches="tight")
plt.show()
[docs]
def plot_cell_type_distribution_fov(
self,
dataset: str,
data_type: str = 'number',
colormap: str = 'tab20c',
save_path: str = None
):
"""
Visualize the distribution of cell types across FOVs in the dataset using a stacked bar plot.
Parameters
----------
dataset : str
Dataset name.
data_type : {'number', 'proportion'}, default='number'
Plot bar plot by number of cells or by proportions.
colormap : str, default='tab20c'
Matplotlib colormap name.
save_path : str, optional
Path to save the figure.
Returns
-------
None
Displays a stacked bar plot.
"""
if data_type not in ['number', 'proportion']:
raise ValueError("Invalild data_type. It should be one of 'number' or 'proportion'.")
valid_ds_names = [s.dataset.split('_')[0] for s in self.spatial_queries]
if dataset not in valid_ds_names:
raise ValueError(f"Invalid input dataset name: {dataset}. \n"
f"Valid dataset names are: {set(valid_ds_names)}")
valid_queries = [s for s in self.spatial_queries if s.dataset.split('_')[0] == dataset]
cell_types = set([ct for s in valid_queries for ct in s.labels.unique()])
summary = defaultdict(lambda: defaultdict(int))
for s in valid_queries:
for cell_type in cell_types:
summary[s.dataset][cell_type] = np.sum(s.labels == cell_type)
df = pd.DataFrame([(dataset, cell_type, count)
for dataset, cell_types in summary.items()
for cell_type, count in cell_types.items()],
columns=['Dataset', 'Cell Type', 'Count'])
df['FOV'] = df['Dataset'].str.split('_').str[1]
summary = df.groupby(['FOV', 'Cell Type'])['Count'].sum().reset_index()
plot_data = summary.pivot(columns='Cell Type', index='FOV', values='Count').fillna(0)
# Sort the cell types by total count (descending)
row_sums = plot_data.sum(axis=1)
plot_data_sorted = plot_data.loc[row_sums.sort_values(ascending=False).index]
if data_type != 'number':
plot_data_sorted = plot_data_sorted.div(plot_data_sorted.sum(axis=1), axis=0)
n_col = df['Cell Type'].unique().shape[0] # bar 中 stack 的数量
cmap = cm.get_cmap(colormap)
colors = [cmap(i) for i in range(n_col)]
n_x = plot_data.shape[0]
fig_w = max(6, n_x * 0.5)
fig_h = fig_w * 0.4
fig, ax = plt.subplots(figsize=(fig_w, fig_h))
# Create the stacked bar plot
plot_data_sorted.plot(kind='bar',
stacked=True,
edgecolor='black',
color=colors,
ax=ax,
)
ax.set_title(f"Distribution of FOVs in {dataset} dataset", fontsize=16)
ax.set_xlabel("FOV", fontsize=12)
ax.set_ylabel(
"Number of Cells" if data_type == "number" else "Proportion of Cells",
fontsize=12,
)
ax.tick_params(axis='x', rotation=0, labelsize=10)
ax.legend(
title="Cell type",
loc='center left',
bbox_to_anchor=(1, 0.5),
fontsize=10,
)
plt.tight_layout(rect=[0, 0, 0.85, 1])
if save_path is not None:
plt.savefig(save_path, dpi=300, bbox_inches="tight")
plt.show()
[docs]
def compute_gene_gene_correlation(self,
ct: str,
motif: Union[str, List[str]],
dataset: Union[str, List[str]] = None,
genes: Optional[Union[str, List[str]]] = None,
max_dist: Optional[float] = None,
k: Optional[int] = None,
min_size: int = 0,
min_nonzero: int = 10,
alpha: Optional[float] = None
) -> pd.DataFrame:
"""
Compute gene-gene co-varying patterns between motif and center cells across multiple FOVs.
Similar to compute_gene_gene_correlation in single FOV, but:
- Aggregates center-neighbor pairs across all FOVs in specified dataset
- Uses FOV-specific cell type means for centering (NOT global means)
- Computes correlations by accumulating statistics across FOVs
This function calculates cross correlation between gene expression in motif cells that are
neighbors of the center type, motif cells that are not neighbors, and neighboring cells
without nearby motif. Aggregation is performed across all FOVs in the specified dataset.
Parameters
----------
ct : str
Cell type as the center cells.
motif : str or List[str]
Motif (names of cell types) to be analyzed.
dataset : str or List[str], optional
Datasets to include in analysis. If None, use all datasets.
genes : str or List[str], optional
List of genes to analyze. If None, uses intersection of genes across all FOVs.
max_dist : float, optional
Maximum distance for considering a cell as a neighbor. Use either max_dist or k.
k : int, optional
Number of nearest neighbors. Use either max_dist or k.
min_size : int, default=0
Minimum neighborhood size for each center cell (only used when max_dist is specified).
min_nonzero : int, default=10
Minimum number of non-zero expression values required for a gene to be included.
alpha : float, optional
Significance threshold.
Returns
-------
results_df : pd.DataFrame
DataFrame with correlation results. Columns: gene_center, gene_motif, corr_neighbor,
corr_non_neighbor, corr_center_no_motif, p_value_test1, p_value_test2,
delta_corr_test1, delta_corr_test2, combined_score, adj-pval-test1, adj-pval-test2.
"""
if not self.build_gene_index:
print('Computing covarying genes using expression data ...')
results_df = spatial_gene_covarying.compute_gene_gene_correlation_adata_multi_fov(
sq_objs=self,
ct=ct,
motif=motif,
dataset=dataset,
genes=genes,
max_dist=max_dist,
k=k,
min_size=min_size,
min_nonzero=min_nonzero,
alpha=alpha,
)
else:
print('Compting covarying genes using binary data ...')
results_df = spatial_gene_covarying.compute_gene_gene_correlation_binary_multi_fov(
sq_objs=self,
ct=ct,
motif=motif,
dataset=dataset,
genes=genes,
max_dist=max_dist,
k=k,
min_size=min_size,
min_nonzero=min_nonzero,
alpha=alpha,
)
return results_df
[docs]
def compute_gene_gene_correlation_by_type(self,
ct: str,
motif: Union[str, List[str]],
dataset: Union[str, List[str]] = None,
genes: Optional[Union[str, List[str]]] = None,
max_dist: Optional[float] = None,
k: Optional[int] = None,
min_size: int = 0,
min_nonzero: int = 10,
alpha: Optional[float] = None
) -> pd.DataFrame:
"""
Compute gene-gene cross correlation separately for each cell type in the motif across multiple FOVs.
Similar to compute_gene_gene_correlation_by_type in spatial_query.py, but aggregates across FOVs.
For each non-center cell type in the motif, compute:
- Correlation 1: Center cells with motif vs neighboring motif cells of THIS TYPE
- Correlation 2: Center cells with motif vs distant motif cells of THIS TYPE
- Correlation 3: Center cells without motif vs neighbors (same for all types)
Parameters
----------
ct : str
Cell type as the center cells.
motif : str or List[str]
Motif (names of cell types) to be analyzed. Include all cell types for neighbor finding.
dataset : str or List[str], optional
Datasets to include in analysis. If None, use all datasets.
genes : str or List[str], optional
List of genes to analyze. If None, uses intersection of genes across all FOVs.
max_dist : float, optional
Maximum distance for considering a cell as a neighbor. Use either max_dist or k.
k : int, optional
Number of nearest neighbors. Use either max_dist or k.
min_size : int, default=0
Minimum neighborhood size for each center cell (only used when max_dist is specified).
min_nonzero : int, default=10
Minimum number of non-zero expression values required for a gene to be included.
alpha : float, optional
Significance threshold.
Returns
-------
pd.DataFrame
DataFrame with correlation results for each cell type and gene pair.
Columns: cell_type, gene_center, gene_motif, corr_neighbor, corr_non_neighbor,
corr_center_no_motif, p_value_test1, p_value_test2, q_value_test1, q_value_test2,
delta_corr_test1, delta_corr_test2, reject_test1_fdr, reject_test2_fdr,
combined_score, abs_combined_score.
"""
if not self.build_gene_index:
print('Computing covarying genes using expression data ...')
results_df = spatial_gene_covarying.compute_gene_gene_correlation_by_type_adata_multi_fov(
sq_objs=self,
ct=ct,
motif=motif,
dataset=dataset,
genes=genes,
max_dist=max_dist,
k=k,
min_size=min_size,
min_nonzero=min_nonzero,
alpha=alpha
)
else:
print('Compting covarying genes using binary data ...')
results_df = spatial_gene_covarying.compute_gene_gene_correlation_by_type_binary_multi_fov(
sq_objs=self,
ct=ct,
motif=motif,
dataset=dataset,
genes=genes,
max_dist=max_dist,
k=k,
min_size=min_size,
min_nonzero=min_nonzero,
alpha=alpha
)
return results_df
[docs]
@staticmethod
def test_score_difference(
result_A: pd.DataFrame,
result_B: pd.DataFrame,
score_col: str = 'combined_score',
significance_col: str = 'if_significant',
gene_center_col: str = 'gene_center',
gene_motif_col: str = 'gene_motif',
percentile_threshold: float = 95.0,
background: Literal['Overlapping', 'Significant'] = 'Significant'
) -> pd.DataFrame:
"""
Identify gene pairs with large score differences between two covariation pattern results.
This function compares covariation scores between two groups (e.g., disease vs control,
treatment vs baseline, covarying gene pairs by distinct motif types) and identifies gene
pairs with the largest score differences using percentile-based ranking. Gene pairs in the
top percentile_threshold% and bottom (100 - percentile_threshold)% of score differences are
flagged as emperical significant ones.
Parameters
----------
result_A : pd.DataFrame
Results from compute_gene_gene_correlation or compute_gene_gene_correlation_by_type
for condition A. Must contain columns: gene_center, gene_motif, combined_score, if_significant.
result_B : pd.DataFrame
Results from compute_gene_gene_correlation or compute_gene_gene_correlation_by_type
for condition B. Must contain the same columns as result_A.
score_col : str, default='combined_score'
Name of the column containing covariation scores to compare.
significance_col : str, default='if_significant'
Name of the column indicating whether a gene pair is significant.
gene_center_col : str, default='gene_center'
Name of the column containing center gene names.
gene_motif_col : str, default='gene_motif'
Name of the column containing motif gene names.
percentile_threshold : float, default=95.0
Percentile threshold for identifying outliers. Gene pairs with score_diff in the
top percentile_threshold% (e.g., >95th percentile) or bottom (100 - percentile_threshold)%
(e.g., <5th percentile) are flagged as outliers.
background : Literal['Overlapping', 'Significant'], default='Significant'
Defines the background gene pairs for comparison:
- 'Significant': Only gene pairs significant in at least one condition are considered.
- 'Overlapping': All overlapping gene pairs between both conditions are considered.
Returns
-------
pd.DataFrame
DataFrame with gene pair comparison results. Columns include:
- gene_center: center gene name
- gene_motif: motif gene name
- score_A: covariation score in condition A
- score_B: covariation score in condition B
- score_diff: score difference (score_A - score_B)
- percentile: percentile rank of score_diff in the distribution
- is_outlier: True if percentile > percentile_threshold or < (100 - percentile_threshold)
- significant_in_A: whether the gene pair is significant in condition A
- significant_in_B: whether the gene pair is significant in condition B
- outlier_direction: 'higher_in_A' (top percentile), 'higher_in_B' (bottom percentile),
or 'not_outlier'
"""
return spatial_gene_covarying.test_score_difference(result_A, result_B, score_col, significance_col, gene_center_col, gene_motif_col, percentile_threshold, background)