1. Spatial Transcriptomics (ST) prediction tutorial

⚠️ Note:
The following is an example demonstration. Here, the ST prediction model is trained on only 5 samples, using a single inner fold and a single outer fold.
The training process runs for a maximum of 50 epochs.
In the manuscript, the model was trained for 500 epochs on 25 samples from 14 patients,
utilizing 14 outer folds and 7 inner folds, resulting in a total of 98 trained models (14 × 7).

Load packages

import numpy as np
import pandas as pd
import os
import subprocess
import sys
import pkg_resources
import time
import platform
import psutil
import pickle
from scipy.spatial import cKDTree
from pathlib import Path

# Function to get system information
import platform
import psutil
import torch

def get_system_info():
    print("\n--- System Information ---")
    print(f"Operating System: {platform.system()} {platform.release()} ({platform.version()})")
    print(f"Processor: {platform.processor()}")
    print(f"CPU Cores: {psutil.cpu_count(logical=False)} (Physical), {psutil.cpu_count(logical=True)} (Logical)")
    print(f"Total RAM: {psutil.virtual_memory().total / (1024**3):.2f} GB")
    gpu = torch.cuda.get_device_name(0) if torch.cuda.is_available() else "CPU only (no CUDA GPU detected)"
    print(f"GPU: {gpu}")
get_system_info()

1.1. Feature extraction

Load input/output folder and metadata

path2input = f'../input_data/'  # inputs folder 
path2output = f'../output_data/' # outputs folder 

# Load the meta file
meta_path = f"{path2input}metadata.csv"
metadata = pd.read_csv(meta_path)

# List all samples in meta
samples=metadata['sample_ID']
print("Samples in meta:", samples)

Feature extraction for the samples

# Measure execution time
start_time = time.time()

try:
    # Loop over the slides and run the feature extraction for each of the samples
    for i_slide in range(len(samples)):
        command = f"python ../scripts/1.ST_prediction/1.1.Feature_extraction/1main_feature_extraction.py {path2input} {path2output} {i_slide}"
        print(command)
        result = subprocess.run(command, shell=True, check=True, capture_output=True, text=True, cwd=os.getcwd())
        print("Command Output:\n", result.stdout)
except subprocess.CalledProcessError as e:
    print("Error occurred:\n", e.stderr)

# Calculate elapsed time
end_time = time.time()
elapsed_time = end_time - start_time
print(f"Execution Time: {elapsed_time:.2f} seconds")

Collect features

Create a single pickle file containing all the features

feature_path=f'{path2output}features/'
features = []

# Iterate through all files in the directory
for filename in os.listdir(feature_path):
    if filename.endswith(".pkl"):
        file_path = os.path.join(feature_path, filename)
        x=np.load(file_path, allow_pickle=True)
        features.append((filename, x))

with open(f'{path2output}/features.pkl', 'wb') as f:
    pickle.dump(features, f)

View features for one spot

features[0]

1.2. Regression: Training the regression model for one fold

Load the split file. Here, since we have 5 slides, we use a split file for 5x4 nested cross validation

# Load the .npz file
file_path = f"{path2input}train_valid_test_idx.npz"
split_file = np.load(file_path, allow_pickle = True)

# Explore the contents of each key
for key in split_file.files:
    print(f"\nKey: {key}")
    print("Shape:", split_file[key].shape)

Load the gene file

genes=pd.read_pickle(f'{path2input}gene_file.pkl')['gene']

Define the model parameters

ik_fold=0 # to fully train the model, the regression model should be trained for ik_folds=[0,1,2,3,4]
il_fold=0 # to fully train the model, the regression model should be trained for il_folds=[0,1,2,3]
max_epochs=10  # to fully train the model, the regression model should be trained with max_epochs=500

Define the command

command = f"python ../scripts/1.ST_prediction/1.2.Regression/1main_regression.py {path2input} {path2output}  {ik_fold} {il_fold} {max_epochs}"
print(command)

Execute the command

# Measure execution time
start_time = time.time()

try:
    result = subprocess.run(command, shell=True, check=True, capture_output=True, text=True, cwd=os.getcwd())
    print("Command Output:\n", result.stdout)
except subprocess.CalledProcessError as e:
    print("Error occurred:\n", e.stderr)

# Calculate elapsed time
end_time = time.time()
elapsed_time = end_time - start_time
print(f"Execution Time: {elapsed_time:.2f} seconds")

1.3 Prediction: ST prediction for a TCGA slide using a single trained model

1.3.1 Feature extraction on the test slide

We need to extract the features on the TCGA slide (like we did in step 1.1)

Define the inputs

path2input_pred = f'{path2input}slides/TCGA-OL-A5S0-01Z-00-DX1.49A7AC9D-C186-406C-BA67-2D73DE82E13B.svs'
path2output_pred = f'{path2output}features_test/'
slide_name_pred = 'OL-A5S0-01Z-00-DX1'

Define the command

command = f"python ../scripts/1.ST_prediction/1.1.Feature_extraction/1main_feature_extraction_TCGA.py {path2input_pred} {path2output_pred} {slide_name_pred}"
print(command)

Execute the command

# Measure execution time
start_time = time.time()

try:
    result = subprocess.run(command, shell=True, check=True, capture_output=True, text=True, cwd=os.getcwd())
    print("Command Output:\n", result.stdout)
except subprocess.CalledProcessError as e:
    print("Error occurred:\n", e.stderr)

# Calculate elapsed time
end_time = time.time()
elapsed_time = end_time - start_time
print(f"Execution Time: {elapsed_time:.2f} seconds")

1.3.2 ST Prediction on the test slide

Example of Spatial Transcriptomics (ST) Prediction for a TCGA Slide Using a Model with One Inner and One Outer Fold

We start by importing the necessary modules and setting the path to the prediction script.

# Define the path to the prediction script and add it to the system path
script_path = Path("../scripts/1.ST_prediction/1.3.Prediction")
sys.path.append(str(script_path))

# Import the prediction function
from TCGA_Prediction import predict_st_from_trained_model

#Define the path to the test features file
path_to_test_features = "../output_data/features_test/OL-A5S0-01Z-00-DX1.pkl"
# Load the features from the pickle file
features_list = np.load(f"{path_to_test_features}", allow_pickle=True)
# Define the outer and inner cross-validation folds
ik_folds = [0]  # Example outer fold index
il_folds = [0]  # Example inner fold index

Run the Prediction Model

Now, we use the predict_st_from_trained_model function to generate predictions for the test slide.

tcga_st_pred = predict_st_from_trained_model(ik_folds, il_folds,features_list,
                                              path2models = path2output,
                                              genes = genes)

View the Prediction Results

tcga_st_pred.head()

1.3.3. Smoothing the predictions

After obtaining thre predictiosn, we perform the smoothing step

#### Smoothing function

def smooth_genes_kdtree(slide_df, genes, radius=2, weights = 'uniform'):

    # Extract spatial coordinates and gene data
    coordinates = slide_df[['x', 'y']].values
    gene_data = slide_df[genes].values

    # Build a KDTree for fast radius neighbor searches
    tree = cKDTree(coordinates)

    # Query the tree for neighbors within the radius for each point
    indices = [tree.query_ball_point(point, r=radius) for point in coordinates]

    # Prepare an array to hold smoothed gene values
    smoothed_gene_data = np.zeros_like(gene_data)

    # Compute average values of all genes for neighbors within the radius
    for i, neighbors in enumerate(indices):
        if neighbors:
            # Calculate the mean across all selected genes for these indices
            smoothed_gene_data[i] = gene_data[neighbors].mean(axis=0)
            if weights != 'uniform':
                
                smoothed_gene_data[i] = (len(neighbors)*smoothed_gene_data[i] + (weights-1)*gene_data[i])/(len(neighbors) + weights -1)

        else:
            # Handle cases with no neighbors (could handle differently if needed)
            smoothed_gene_data[i] = gene_data[i]

    # Create a DataFrame for smoothed gene values
    smoothed_knn = pd.DataFrame(smoothed_gene_data, columns=genes, index=slide_df.index)

    # Optionally, combine original data with smoothed data
    result_df = slide_df.copy()
    result_df.update(smoothed_knn)

    return result_df

## extract x,y coordinates from spot names

def get_x_y_grid(df_pred, tile_size=640):
    slide_names=df_pred.index.str.split('_').str[0]
    
    # Split 'spot_id' to extract 'x' and 'y' coordinates    
    x = df_pred.index.str.split('_').str[-1]
    y = df_pred.index.str.split('_').str[-2]

    # Convert 'x' and 'y' from strings to floats
    x = (x.astype(float)-x.astype(float).min())//tile_size
    y = (y.astype(float)-y.astype(float).min())//tile_size
    
    out_df = pd.DataFrame({'slide_file_name':slide_names,'x':x, 'y':y}, index = df_pred.index)


    return out_df

grid_df = get_x_y_grid(tcga_st_pred)
tcga_st_pred = grid_df.join(tcga_st_pred)
tcga_st_pred = tcga_st_pred.sort_values(['slide_file_name', 'x', 'y'])

tcga_st_pred_s = smooth_genes_kdtree(tcga_st_pred, genes, radius=2, weights ='uniform')
tcga_st_pred_s.head(5)

Save the smoothed and non-smoothed predictions

tcga_st_pred_s.to_pickle(f'{path2output}smoothed_preds_tcga.pkl')
tcga_st_pred.to_pickle(f'{path2output}preds_tcga.pkl')

Python and packages versions

# Print Python version
print(f"Python version: {sys.version}")

# Print installed packages and their versions
installed_packages = {pkg.key: pkg.version for pkg in pkg_resources.working_set}
for package, version in installed_packages.items():
    print(f"{package}=={version}")