SpatialDM Analysis of CCI
The SpatialDM documentation and package could be found here: https://spatialdm.readthedocs.io/en/latest
The paper could be found here: https://doi.org/10.1038/s41467-023-39608-w
1. Import Package
[ ]:
import os
import pandas as pd
import numpy as np
import anndata as ann
import scanpy as sc
import spatialdm as sdm
from spatialdm.datasets import dataset
import spatialdm.plottings as pl
import matplotlib.pyplot as plt
import time
from itertools import zip_longest
2. Getting Ligand-Receptor Database
Because I run the code in HPC, and that our HPC are unable to use the
sdm.extract_lr, I extract this function from the package.The L-R database csv files could be found in here: https://github.com/StatBiomed/SpatialDM/tree/main/spatialdm/datasets/LR_data
[ ]:
def extract_lr(adata, species, mean='algebra', min_cell=0, datahost='builtin'):
"""
find overlapping LRs from CellChatDB
:param adata: AnnData object
:param species: support 'human', 'mouse' and 'zebrafish'
:param mean: 'algebra' (default) or 'geometric'
:param min_cell: for each selected pair, the spots expressing ligand or receptor should be larger than the min,
respectively.
:param datahost: the host of the ligand-receptor data. 'builtin' for package built-in otherwise from figshare
:return: ligand, receptor, geneInter (containing comprehensive info from CellChatDB) dataframes \
in adata.uns
"""
if mean=='geometric':
from scipy.stats.mstats import gmean
adata.uns['mean'] = mean
if datahost == 'package':
if species in ['mouse', 'human', 'zerafish']:
datapath = './datasets/LR_data/%s-' %(species)
else:
raise ValueError("species type: {} is not supported currently. Please have a check.".format(species))
import pkg_resources
stream1 = pkg_resources.resource_stream(__name__, datapath + 'interaction_input_CellChatDB.csv.gz')
geneInter = pd.read_csv(stream1, index_col=0, compression='gzip')
stream2 = pkg_resources.resource_stream(__name__, datapath + 'complex_input_CellChatDB.csv')
comp = pd.read_csv(stream2, header=0, index_col=0)
else:
if species == 'mouse':
geneInter = pd.read_csv('/rsrch5/home/biostatistics/lku/SpatialDM/mouse_interaction_input_CellChatDB.csv', index_col=0)
comp = pd.read_csv('/rsrch5/home/biostatistics/lku/SpatialDM/mouse_complex_input_CellChatDB', header=0, index_col=0)
elif species == 'human':
geneInter = pd.read_csv('/rsrch5/home/biostatistics/lku/SpatialDM/human_interaction_input_CellChatDB.csv', header=0, index_col=0)
comp = pd.read_csv('/rsrch5/home/biostatistics/lku/SpatialDM/human_complex_input_CellChatDB.csv', header=0, index_col=0)
else:
raise ValueError("species type: {} is not supported currently. Please have a check.".format(species))
geneInter = geneInter.sort_values('annotation')
ligand = geneInter.ligand.values
receptor = geneInter.receptor.values
geneInter.pop('ligand')
geneInter.pop('receptor')
## NOTE: the following for loop needs speed up
t = []
for i in range(len(ligand)):
for n in [ligand, receptor]:
l = n[i]
if l in comp.index:
n[i] = comp.loc[l].dropna().values[pd.Series \
(comp.loc[l].dropna().values).isin(adata.var_names)]
else:
n[i] = pd.Series(l).values[pd.Series(l).isin(adata.var_names)]
if (len(ligand[i]) > 0) * (len(receptor[i]) > 0):
if mean=='geometric':
meanL = gmean(adata[:, ligand[i]].X, axis=1)
meanR = gmean(adata[:, receptor[i]].X, axis=1)
else:
meanL = adata[:, ligand[i]].X.mean(axis=1)
meanR = adata[:, receptor[i]].X.mean(axis=1)
if (sum(meanL > 0) >= min_cell) * \
(sum(meanR > 0) >= min_cell):
t.append(True)
else:
t.append(False)
else:
t.append(False)
ind = geneInter[t].index
adata.uns['ligand'] = pd.DataFrame.from_records(zip_longest(*pd.Series(ligand[t]).values)).transpose()
adata.uns['ligand'].columns = ['Ligand' + str(i) for i in range(adata.uns['ligand'].shape[1])]
adata.uns['ligand'].index = ind
adata.uns['receptor'] = pd.DataFrame.from_records(zip_longest(*pd.Series(receptor[t]).values)).transpose()
adata.uns['receptor'].columns = ['Receptor' + str(i) for i in range(adata.uns['receptor'].shape[1])]
adata.uns['receptor'].index = ind
adata.uns['num_pairs'] = len(ind)
adata.uns['geneInter'] = geneInter.loc[ind]
if adata.uns['num_pairs'] == 0:
raise ValueError("No effective RL. Please have a check on input count matrix/species.")
return
3. Run the Analysis
The below shows an example analysis on analyzing spot level data, where I have prepared the data by the following columns:
["x", "y", "celltypeA",..., "celltypeX", "cell_type(dominant cell type)", "GeneA",..., "GeneX"].
[ ]:
for num in ["slice1","slice2","slice3","slice4"]:
path_exprsn = f'/rsrch5/home/biostatistics/lku/ILIBD/data/{num}/exprsn_df.csv'
data = pd.read_csv(path_exprsn, index_col=0)
data = data.dropna(subset=['cell_type'])
data = data.set_index('cell', drop=False)
data.index.name = None
gene_expression= data.drop(columns=["x","y","cell","Excitatory_neurons","Inhibitory_neuron", "Astrocyte", "Oligodendrocyte","Oligodendrocyte_precursor_cell",
"Microglia","Pericytes", "Endothelial_cells","cell_type" ])
meta_data = data[["Excitatory_neurons","Inhibitory_neuron", "Astrocyte", "Oligodendrocyte","Oligodendrocyte_precursor_cell","Microglia","Pericytes", "Endothelial_cells" ]]
#meta_data.set_index('cell', inplace=True)
adata = sc.AnnData(X=gene_expression, obs=meta_data)
adata.obsm['spatial'] = data[['x', 'y']].to_numpy()
adata.raw = adata.copy()
sc.pp.normalize_total(adata, target_sum=1e4) # Normalizes each cell to have total count of 10,000
sc.pp.log1p(adata)
sdm.weight_matrix(adata, l=1.2, cutoff=0.2, single_cell=False) # weight_matrix by rbf kernel
# extract LR pairs
extract_lr(adata, 'human', min_cell=3)
# global selection
sdm.spatialdm_global(adata, 1000, specified_ind=None, method='both', nproc=1) # global Moran selection
sdm.sig_pairs(adata, method='permutation', fdr=True, threshold=0.1) # select significant pairs # FDR threshold
LR = adata.uns['global_res'].sort_values(by='fdr')
LR_pair = LR.index.tolist()
# local spot selection
sdm.spatialdm_local(adata, n_perm=1000, method='both', specified_ind=None, nproc=1)
sdm.sig_spots(adata, method='permutation', fdr=False, threshold=0.1)
adata.obsm['cell_type'] = adata.obs[adata.obs.columns]
sdm.write_spatialdm_h5ad(adata, filename=f'/rsrch5/home/biostatistics/lku/ILIBD/data/{num}/spatialdm_adata.h5ad')
# Create an empty list to store the Links DataFrames
all_links = []
for pair in LR_pair:
try:
# Compute local_I_R and local_I for the current pair
selected_spots_mask = adata.uns['selected_spots'].index == pair
local_I_R = adata.uns['local_stat']['local_I_R'][:, selected_spots_mask]
local_I = adata.uns['local_stat']['local_I'][:, selected_spots_mask]
# Check if either local_I_R or local_I is empty
if local_I_R.shape[1] == 0 or local_I.shape[1] == 0:
print(f"Skipping pair {pair}: local_I_R or local_I is empty.")
continue # Skip to the next iteration
# Compute ct_R and ct_L
ct_R = local_I_R * adata.obsm['cell_type']
ct_L = local_I * adata.obsm['cell_type']
# Compute sparse_ct_sum
from scipy.sparse import csc_matrix
w = adata.obsp['weight']
sparse_ct_sum = [[
(csc_matrix(w).multiply(ct_L[n1].values).T.multiply(ct_R[n2].values)).sum()
for n1 in ct_L.columns
] for n2 in ct_R.columns]
sparse_ct_sum = np.array(sparse_ct_sum)
# Create Links DataFrame
Links = pd.DataFrame({
'ligand': np.tile(ct_L.columns, ct_R.shape[1]),
'receptor': np.repeat(ct_R.columns, ct_L.shape[1]),
'value': sparse_ct_sum.reshape(1, -1)[0],
'pvalue': np.zeros(ct_L.shape[1] * ct_R.shape[1])
})
# Filter out rows with value <= 0
Links = Links.loc[Links.value > 0]
# Add LR pair name
Links['interaction_name'] = pair
# Append to the list
all_links.append(Links)
except Exception as e:
print(f"Error processing pair {pair}: {e}")
continue
# Combine all Links DataFrames into one
if all_links:
combined_links = pd.concat(all_links, ignore_index=True)
path_write_LR = f'/rsrch5/home/biostatistics/lku/ILIBD/data/{num}/spatialdm_result.csv'
combined_links.to_csv(path_write_LR , index=False)
print(f"Results saved for {num} at {path_write_LR}")
else:
print(f"No valid interactions found for {num}. Skipping CSV write.")