Demo 4: STARmap to MERFISH

import CAST
import scanpy as sc
import torch,os
import numpy as np

workdir = '$demo_path' #### input the demo path

This demo showcases an example workflow aligning a STARmap sample to a MERFISH sample using CAST Mark and CAST Stack.

Load / Prepare Data

###  Load and plot the data

from CAST.visualize import plot_mid


# Output directory for the results
output_path = f'{workdir}/demo4_STARmap_to_MERFISH/demo_output'
os.makedirs(output_path,exist_ok=True)

# Loading the data
sdata = sc.read_h5ad(f'{workdir}/demo4_STARmap_to_MERFISH/data/STARmap_vs_MERFISH.h5ad')
sample_list= ['STARmap','MERFISH']
batch_key = 'batch'
coords_t = np.array(sdata.obs[['center_x', 'center_y']])

# Plotting the initial data
CAST.plot_mid(coords_t[sdata.obs[batch_key] == sample_list[0]],
         coords_t[sdata.obs[batch_key] == sample_list[1]],
         output_path=output_path,
         filename = 'Align_raw',
         title_t = [sample_list[1],
                    sample_list[0]],
         s_t = 8,scale_bar_t = None)

### Extract expression and coordinates and save the data

from CAST.utils import extract_coords_exp, sub_data_extract


# Extract the coordinates and expression data
coords_raw,exps = CAST.extract_coords_exp(sdata, batch_key = 'batch', cols = ['center_x', 'center_y'], count_layer = '.X', data_format = 'norm1e4')
coords_sub,exp_sub,sub_node_idxs = CAST.sub_data_extract(sample_list,coords_raw, exps, nodenum_t = 20000)

# Save the data
torch.save(coords_raw, f'{output_path}/coords_raw.pt')
torch.save(sub_node_idxs, f'{output_path}/sub_node_idxs.pt')
torch.save(exp_sub, f'{output_path}/exp_sub.pt')
torch.save(coords_sub, f'{output_path}/coords_sub.pt')

Preprocessing...

/home/unix/panj/miniconda3/envs/cast_dev_nonrelease_cast/lib/python3.9/site-packages/scanpy/preprocessing/_normalization.py:197: UserWarning: Some cells have zero counts
  warn(UserWarning('Some cells have zero counts'))

CAST Mark

### Run CAST Mark and plot the results

from CAST import CAST_MARK
from CAST.visualize import kmeans_plot_multiple


embed_dict = CAST_MARK(coords_sub,exp_sub,output_path,graph_strategy='delaunay')
kmeans_plot_multiple(embed_dict,sample_list,coords_sub,'demo1',output_path,k=20,dot_size = 10,minibatch=True)

Constructing delaunay graphs for 2 samples...
Training on cuda:0...

Loss: -435.040 step time=0.189s: 100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 400/400 [01:17<00:00,  5.15it/s]

Finished.
The embedding, log, model files were saved to /home/unix/panj/wanglab/jessica/CAST/demo/demo4_STARmap_to_MERFISH/demo_output
Perform KMeans clustering on 40000 cells...
Plotting the KMeans clustering results...

array([18, 11,  1, ..., 10,  2, 19], dtype=int32)

CAST Stack

### Run CAST Stack - Align the two samples

from CAST import CAST_STACK
from CAST.CAST_Stack import reg_params


# Set up parameters
query_sample = sample_list[0]
params_dist = CAST.reg_params(dataname = query_sample,
                            gpu = 0 if torch.cuda.is_available() else -1,
                            #### Affine parameters
                            iterations=150,
                            dist_penalty1=0,
                            bleeding=500,
                            d_list = [3,2,1,1/2,1/3],
                            attention_params = [None,3,1,0],
                            #### FFD parameters
                            dist_penalty2 = [0],
                            alpha_basis_bs = [500],
                            meshsize = [8],
                            iterations_bs = [40],
                            attention_params_bs = [[None,3,1,0]],
                            mesh_weight = [None])
params_dist.alpha_basis = torch.Tensor([1/1000,1/1000,1/50,5,5]).reshape(5,1).to(params_dist.device)

# Run CAST Stack
coord_final = CAST.CAST_STACK(coords_raw,embed_dict,output_path,sample_list,params_dist,sub_node_idxs = sub_node_idxs, rescale=True) # v3.0

Loss: 3801.573: 100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 150/150 [00:22<00:00,  6.67it/s]
/home/unix/panj/miniconda3/envs/cast_dev_nonrelease_cast/lib/python3.9/site-packages/torch/functional.py:568: UserWarning: torch.meshgrid: in an upcoming release, it will be required to pass the indexing argument. (Triggered internally at  /opt/conda/conda-bld/pytorch_1646756402876/work/aten/src/ATen/native/TensorShape.cpp:2228.)
  return _VF.meshgrid(tensors, **kwargs)  # type: ignore[attr-defined]
Loss: 2705.461: 100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 40/40 [00:07<00:00,  5.31it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 40/40 [00:00<00:00, 168.28it/s]

Visualize Results

### Plot the kmeans clustering of the final results

from CAST.visualize import kmeans_plot_multiple


# load the final coordinates
coord_final = torch.load(f'{output_path}/STARmap_align_to_MERFISH_coords_final.data')

# Load the embed dict and sub node indexes
sample_list= ['STARmap','MERFISH']
embed_dict = torch.load(f'{output_path}/demo_embed_dict.pt',map_location='cpu')
sub_node_idxs = torch.load(f'{output_path}/sub_node_idxs.pt')

# Extract the sub coordinates
coords_sub_new = dict()
for sample_t in sample_list:
    coords_sub_new[sample_t] = coord_final[sample_t][sub_node_idxs[sample_t],:]

# Plot the final results
kmeans_plot_multiple(embed_dict,sample_list,coords_sub_new,'demo1_new',output_path,k=20,dot_size = 10,minibatch=True,plot_strategy='stack')

Perform KMeans clustering on 40000 cells...
Plotting the KMeans clustering results...

array([18, 11,  1, ..., 10,  2, 19], dtype=int32)