Tutorial 1: Simulated Data (Spatial Patterns of Varying Complexity)

In this tutorial, we evaluate the INR module’s ability to simultaneously capture diverse spatial patterns while maintaining spatial smoothness. We use the Simulated dataset 1 (Spatial Patterns of Varying Complexity), in which we use raw expression data and cell-type labels from the Human Lymph Node Dataset, and generate four spatial distribution patterns: quadrant, stripe, arc, and layered. To ensure the reliability of the simulation experiments, we conducted three independent replications for each spatial pattern by randomly dividing the dataset.

All datasets used in this paper are available at https://doi.org/10.5281/zenodo.14854747.

Loading package

import warnings
warnings.filterwarnings('ignore')
import torch
import scanpy as sc
import numpy as np
from sklearn.metrics.cluster import adjusted_rand_score
import anndata as ad
import matplotlib.pyplot as plt
from scipy.sparse import csr_matrix

import os
## Should be replaced with the R package installation path
os.environ['R_HOME'] = 'E:/R-4.3.1'
os.environ['R_USER'] = 'E:/anaconda/lib/site-packages/rpy2'

import sys
sys.path.append(r'../..')
from Model.INR import INRModel
from Model.utils import mclust_R, reorder_categories
from Model.model import DCCAE
from Model.preprocess import fix_seed
fix_seed(2024)

Loading data

To facilitate and standardize subsequent operations, the simulated data has already undergone standard preprocessing steps during its generation. Therefore, it can be directly used as input for SpaKnit.

When selecting data, we need to specify two parameters, ‘replicate’ and ‘batch’, which have the following meanings:

• replicate: representing different spatial patterns (1-quadrant, 2-stripe, 3-arc, and 4-layered);

• batch: representing the groups into which the data are randomly divided (ranging 1, 2, and 3).

## read data
## should replace the path below according to the actual data storage situation
replicate = 1
file_fold_1 = f'../Data/Spatial_Scenario_{replicate}/simulation{replicate}_RNA'
file_fold_2 = f'../Data/Spatial_Scenario_{replicate}/simulation{replicate}_Protein'

adata_omics_1 = sc.read_h5ad(file_fold_1 + '.h5ad')
adata_omics_2 = sc.read_h5ad(file_fold_2 + '.h5ad')

batch = 1
adata_RNA = adata_omics_1[adata_omics_1.obs['batch'] == batch]
adata_ADT = adata_omics_2[adata_omics_2.obs['batch'] == batch]

print("shape of modality 1:", np.shape(adata_RNA))
print("shape of modality 2:", np.shape(adata_ADT))

shape of modality 1: (400, 3000)
shape of modality 2: (400, 31)

Training the model

Spatial Continuous Mapping(SCM) module

This module aims to capture diverse spatial patterns while ensuring spatial continuity from spatial information. In this module, spatial coordinates serve as inputs to a neural network employing three sine periodic activation functions to reconstruct each omics modality. This reconstruction process is further enhanced with kernel functions to bolster continuous representation capabilities.

By treating spatial information as prior knowledge rather than a standalone modality, the SCM module learns continuous spatial patterns of omics modalities. This design enables the SCM module to accommodate diverse combinations of modalities across varying dimensions, ranging from raw features or reduced-dimension embeddings.

## INR module training
def Spatially_Continuous_Reconstruction(adata, epoch_num):
    coords = adata.obsm['spatial'].astype(np.float32)
    node_feats = adata.X.astype(np.float32)

    coords = torch.from_numpy(coords).float()
    if isinstance(node_feats, csr_matrix):
        node_feats = node_feats.toarray()
    node_feats = torch.from_numpy(node_feats).float()

    device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
    node_feats = node_feats.to(device)
    coords = coords.to(device)

    print(f'out_dim: {node_feats.shape[1]}')

    model = INRModel(
        X=node_feats,
        spatial_coord=coords,
        device=device,
        learning_rate=1e-4,
        reg_par=0,
        epoch_num=epoch_num,
        print_train_log_info=True
    )

    reconstructed_X = model.fit()
    if 'INR' not in adata.uns:
        adata.uns['INR'] = reconstructed_X

# modality 1
Spatially_Continuous_Reconstruction(adata_omics_1, 500)

# modality 2
Spatially_Continuous_Reconstruction(adata_omics_2, 500)

## Run PCA on the spatially continuous reconstructions of 2 modalities
adata_omics_1.X = adata_omics_1.uns['INR']
adata_omics_2.X = adata_omics_2.uns['INR']

sc.tl.pca(adata_RNA, use_highly_variable=False)
sc.tl.pca(adata_ADT, use_highly_variable=False)

DCCAE module

This module focuses on capturing the nonlinear correlations among omics modalities while eliminating modality-specific noise. In the DCCAE module, the spatially reconstructed continuous omics data serve as input and the CCA Loss as well as the Reconstruction Loss are jointly optimized for model training. Finally, the feature transformation is performed through a linear CCA layer. This design enables the capture of deep nonlinear relationships among omics modalities while preserving their individual characteristics, and facilitates pairwise interaction of the integrated embeddings.

n_DCCAE = 10  # number of components for DCCAE

# Get the feature sizes of the two datasets
features1 = adata_RNA.obsm['X_pca'].shape[1]
features2 = adata_ADT.obsm['X_pca'].shape[1]

# Define the structure of the hidden layers and output size for the two modalities
layers1 = [256, 256, n_DCCAE]
layers2 = [256, 256, n_DCCAE]

# Load the PCA features from the two datasets
X = adata_RNA.obsm['X_pca'].copy()
Y = adata_ADT.obsm['X_pca'].copy()

# Specify the representations to be used
use_rep = ['DCCAE_X', "DCCAE_Y", "DCCAE"]

# Set the number of training epochs
epochs = 300

# Initialize the DCCAE model with specified parameters
dccae = DCCAE(
    input_size1=features1,  # Input size for the first modality
    input_size2=features2,  # Input size for the second modality
    n_components=n_DCCAE,    # Number of components for DCCAE
    layer_sizes1=layers1,   # Hidden layer sizes for the first modality
    layer_sizes2=layers2,   # Hidden layer sizes for the second modality
    epoch_num=epochs,       # Number of training epochs
    learning_rate=0.001     # Learning rate for training
)

# Train the DCCAE model using the two datasets X and Y
dccae.fit([X, Y])

# Transform the input data using the trained DCCAE model
Xs_transformed = dccae.transform([X, Y])

# Assign the transformed features to the AnnData objects
adata_RNA.obsm["DCCAE_X"] = Xs_transformed[0]
adata_ADT.obsm["DCCAE_Y"] = Xs_transformed[1]

# Concatenate the transformed features from both modalities
adata_RNA.obsm["DCCAE"] = np.concatenate(
    (adata_RNA.obsm["DCCAE_X"], adata_ADT.obsm["DCCAE_Y"]),
    axis=1
)

# Specify the cluster number for mclust analysis
n = 4

# Perform mclust clustering on the concatenated features (DCCAE)
mclust_R(adata_RNA, used_obsm=use_rep[2], num_cluster=n)
obs_df = adata_RNA.obs.dropna()
ARI = adjusted_rand_score(
    obs_df['clusters_mclust'], obs_df['Ground Truth']
)
print(f'n={n}, DCCAE, ARI = {ARI}')

Training Progress: 100%|██████████| 300/300 [00:21<00:00, 14.20it/s]

model training finished!
fitting ...
  |                                                                      |   0%

  |======================================================================| 100%
n=4, DCCAE, ARI = 1.0

colors_domain = [
    '#f19c79', '#2a9d8f', '#e9c46a', '#264653'
]

plt.rcParams['font.size'] = 20
plt.rcParams['font.sans-serif'] = ['Arial']

# Calculate ARI
ari_value = adjusted_rand_score(adata_RNA.obs['Ground Truth'], adata_RNA.obs['clusters_mclust'])

# Create figure
fig, ax = plt.subplots(1, 1, figsize=(4, 4))

# Plot the spatial embedding
sc.pl.embedding(
    adata_RNA,
    basis='spatial',
    color=['clusters_mclust'],
    title=None,
    s=600,
    colorbar_loc=None,
    show=False,
    ax=ax,
    palette=colors_domain,
    legend_loc=None
)

# Set title with ARI value
ax.set_title(f"SpaKnit (ARI={ari_value:.2f})")
ax.set_xlabel('')
ax.set_ylabel('')

# Hide axis borders
for spine in ax.spines.values():
    spine.set_visible(False)

# Adjust subplot parameters
plt.subplots_adjust(left=0, right=1, top=1, bottom=0)
plt.tight_layout()
plt.show()

Storing the results

adata_1 = ad.AnnData(obs=adata_RNA.obs[['clusters_mclust', 'batch']], obsm={'SpaKnit': adata_RNA.obsm['DCCAE']})
adata_2 = ad.AnnData(obs=adata_RNA.obs[['clusters_mclust', 'batch']], obsm={'SpaKnit': adata_RNA.obsm['DCCAE']})
adata_3 = ad.AnnData(obs=adata_RNA.obs[['clusters_mclust', 'batch']], obsm={'SpaKnit': adata_RNA.obsm['DCCAE']})
adata_results = adata_1.concatenate(adata_2, adata_3, batch_key='batch')

obs_df = adata_results.obs
obs_df = obs_df.rename(columns={'clusters_mclust': 'SpaKnit'})
adata_results.obs = obs_df

# results = sc.read_h5ad(f'./Results/Spatial_Scenario_{replicate}.h5ad')
# results.obs['SpaKnit'] = adata_results.obs['SpaKnit'].values
# results.obsm['SpaKnit'] = adata_results.obsm['SpaKnit']
# results.write_h5ad(f'./Results/Spatial_Scenario_{replicate}.h5ad')

AnnData object with n_obs × n_vars = 1200 × 3031
    obs: 'Ground Truth', 'batch', 'SpaGCN', 'SpatialGlue', 'STAGATE', 'MultiMAP', 'MultiVI', 'Modality1', 'Modality2', 'SpaKnit'
    obsm: 'Modality1', 'Modality2', 'MultiMAP', 'MultiVI', 'STAGATE', 'SpatialGlue', 'spatial', 'SpaKnit'

Visualizing results

# Define new_orders for each replicate
new_order_dict = {
    1: [1, 3, 4, 2],
    2: [1, 2, 3, 4],
    3: [1, 3, 2, 4],
    4: [1, 2, 4, 3]
}

colors_domain = [
    '#f19c79', '#e9c46a', '#2a9d8f', '#264653'
]

plt.rcParams['font.size'] = 20
plt.rcParams['font.sans-serif'] = ['Arial']

# Loop through each replicate
for replicate in [1, 2, 3, 4]:
    # Load the results
    adata_analysis = sc.read_h5ad(f'./Results/Spatial_Scenario_{replicate}.h5ad')

    # Select batch 0
    batch = '0'  # 0, 1, 2
    adata_analysis = adata_analysis[adata_analysis.obs['batch'] == batch]

    # Reorder categories based on the replicate's new_order
    new_order = new_order_dict[replicate]
    method = 'SpaKnit'
    reorder_categories(adata_analysis, method, new_order)

    # Calculate ARI
    ari_value = adjusted_rand_score(adata_analysis.obs['Ground Truth'], adata_analysis.obs[method])

    # Create figure
    fig, ax = plt.subplots(1, 1, figsize=(4, 4))

    # Plot the spatial embedding
    sc.pl.embedding(
        adata_analysis,
        basis='spatial',
        color=[method],
        title=None,
        s=600,
        colorbar_loc=None,
        show=False,
        ax=ax,
        palette=colors_domain,
        legend_loc=None
    )

    # Set title with ARI value
    ax.set_title(f"{method} (ARI={ari_value:.2f})")
    ax.set_xlabel('')
    ax.set_ylabel('')

    # Hide axis borders
    for spine in ax.spines.values():
        spine.set_visible(False)

    # Adjust subplot parameters
    plt.subplots_adjust(left=0, right=1, top=1, bottom=0)
    plt.tight_layout()
    plt.show()