FineST enhances gene signal
Main steps
Step1: Align image patch and ST spot: using image pixels with spot coordinates
Step2: Train model on within-spots: train model on 80% - test on 20% of the data
Step3: Infer gene expression of between-spots using the trained model
Step4: Impute (smooth) on all sub-spots using measured spot expression
[1]:
import os
import time
import torch
import numpy as np
import pandas as pd
import scanpy as sc
from datetime import datetime
import json
import warnings
warnings.filterwarnings('ignore')
import logging
logging.getLogger().setLevel(logging.INFO)
[2]:
## From local
path = '/mnt/lingyu/nfs_share2/Python/'
os.chdir(str(path) + 'FineST/FineST/')
import FineST as fst
from FineST.utils import device
from FineST.datasets import dataset
import FineST.plottings as fstplt
print("FineST version: %s" %fst.__version__)
## install FineST package
# import FineST as fst
# from FineST.datasets import dataset
# import FineST.plottings as fstplt
# print("torch version: %s" %torch.__version__)
# print("FineST version: %s" %fst.__version__)
fst.setup_seed(666)
FineST version: 0.1.3
CUDA is available. GPU: NVIDIA A100-PCIE-40GB
[3]:
## Set path
path = '/mnt/lingyu/nfs_share2/Python/'
# os.chdir(f'{path}FineST/FineST/')
os.chdir(f'{path}FineST/FineST_demo/') # here I use FineST_demo
## Set logging
model_folder = 'logging/'
dir_name = model_folder + datetime.now().strftime('%Y%m%d%H%M%S%f')
print('dir_name: ', dir_name)
if not os.path.exists(dir_name):
os.makedirs(dir_name)
logger = fst.setup_logger(dir_name)
## Set parameter
parameter_file_path = 'parameter/parameters_HCC_virchow2_modify.json'
with open(parameter_file_path,"r") as json_file:
params = json.load(json_file)
logger.info("Load parameters:\n" + json.dumps(params, indent=2))
[2026-03-29 19:00:37] INFO - Load parameters:
{
"n_encoder_hidden_matrix": 256,
"n_encoder_hidden_image": 128,
"n_encoder_latent": 128,
"n_projection_hidden": 256,
"n_projection_output": 128,
"batch_size": 200,
"batch_size_pair": 64,
"n_encoder_layers": 2,
"dropout_rate": 0,
"training_epoch": 50,
"inital_learning_rate": 0.2,
"k_nearest_positives": 0,
"temperature": 0.03,
"n_input_image": 1280,
"w1": 0,
"w2": 0,
"w3": 1,
"w4": 1
}
dir_name: logging/20260329190037934881
You can set parameter_file_path as the optimal parameters obtained from your train and test model for your dataset, but here we use our trained parameter for paper results repeated (saved in parameter/parameters_HCC_virchow2_modify.json).
1. Load ST RNA-seq and Image data
Download the HCC_P7T folder from figshare. Note: run in terminal
cd /mnt/lingyu/nfs_share2/Python/FineST/FineST_dem
wget "https://figshare.com/ndownloader/articles/26763241/versions/12?folder_path=HCCP7T" -O HCCP1T.zip
unzip HCCP1T.zip
rm HCCP1T.zip
The human hepatocellular carcinoma (P1T adn P7T) from Liu, Y. et al. are publicly available in skrx2fz79n. For easier reuse, we included them in an anndata object which can be loaded directly in FineST Python package. Otherwise you can also using:
demo_data = 'HCCP1T/'adata = sc.read_visium(f'{demo_data}')[ ]:
## Convert HCC data from mtx file to AnnData object
demo_data = 'HCC_data/'
adata = fst.convert_mtx2adata(f'{demo_data}filtered_feature_bc_matrix', f'{demo_data}spatial/')
print(adata)
Warning: fail to load spatial images/scalefactors. name 'positions_path' is not defined
AnnData object with n_obs × n_vars = 3354 × 36601
obs: 'in_tissue', 'array_row', 'array_col', 'pxl_row_in_fullres', 'pxl_col_in_fullres'
var: 'gene_ids', 'feature_types'
uns: 'spatial'
obsm: 'spatial'
[ ]:
## save adata
# demo_data = 'HCC_data/'
# adata.write_h5ad(f'{demo_data}P1T_filtered_feature_bc_matrix.h5ad')
[4]:
adata = dataset.HCCP1T()
print(adata)
100%|██████████| 88.4M/88.4M [00:16<00:00, 5.76MB/s]
AnnData object with n_obs × n_vars = 3354 × 36601
obs: 'in_tissue', 'array_row', 'array_col', 'pxl_row_in_fullres', 'pxl_col_in_fullres'
var: 'gene_ids', 'feature_types'
uns: 'spatial'
obsm: 'spatial'
[ ]:
# adata = dataset.HCCP1T()
demo_data = 'HCC_data/'
adata = sc.read_h5ad(f'{demo_data}P1T_filtered_feature_bc_matrix.h5ad')
print(adata)
AnnData object with n_obs × n_vars = 3354 × 36601
obs: 'in_tissue', 'array_row', 'array_col', 'pxl_row_in_fullres', 'pxl_col_in_fullres'
var: 'gene_ids', 'feature_types'
uns: 'spatial'
obsm: 'spatial'
In FineST, we only considers 963 human genes involved in ligand-receptor pairs from CellChatDB (version.1.1.3) (mouse: 2,022 pairs, human: 1,940 pairs, zebrafish: 2,774 pairs), and extract the overlap gene’s expression profile.
[5]:
## Selected LR genes. Here gene_list can be 'LR_genes', 'HV_genes' or 'LR_HV_genes'
adata = fst.adata_LR(adata, gene_list='LR_HV_genes', n_top_genes=500)
After making an intersection with the LR genes and filtering out genes that are detected in less than 10 cells, 596 genes from the HCC datasetare used for training and prediction.
[6]:
adata_count = adata.copy()
adata = fst.adata_preprocess(adata_count, normalize=False)
print(adata)
AnnData object with n_obs × n_vars = 3354 × 1073
obs: 'in_tissue', 'array_row', 'array_col', 'pxl_row_in_fullres', 'pxl_col_in_fullres', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'pct_counts_in_top_50_genes', 'pct_counts_in_top_100_genes', 'pct_counts_in_top_200_genes', 'pct_counts_in_top_500_genes', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt'
var: 'gene_ids', 'feature_types', 'mt', 'n_cells_by_counts', 'mean_counts', 'log1p_mean_counts', 'pct_dropout_by_counts', 'total_counts', 'log1p_total_counts', 'n_cells'
uns: 'spatial'
obsm: 'spatial'
Note The spatial coordinates tissue_positions_list.csv and gene expression matrix row/barcode need to be re-ordered by the segmented image file name.
Step0: HE image feature extraction
python ./demo/Image_feature_extraction.py \
--dataset HCC \
--position_path HCC_data/spatial/tissue_positions_list.csv \
--rawimage_path HCC_data/T11-V10F06-115-A1.jpg \
--scale_image False \
--method Virchow2 \
--patch_size 112 \
--output_img HCC_data/ImgEmbeddings/pth_112_14_image \
--output_pth HCC_data/ImgEmbeddings/pth_112_14 \
--logging HCC_data/ImgEmbeddings/
[7]:
## Load ST spot position
demo_data = 'HCC_data/'
position = pd.read_csv(f'{demo_data}spatial/tissue_positions_list.csv', header=None)
position = position.rename(columns={position.columns[-2]: 'pixel_x',
position.columns[-1]: 'pixel_y'})
print("The coords of ST spot: ", position.shape)
## Order spot position by image file name
file_paths = sorted(os.listdir(f'{demo_data}/ImgEmbeddings/Bpth_112_14/'))
print("Image embedding file: ", file_paths[:3])
## Get patch position and merge position
position_image = fst.get_image_coord(file_paths, dataset_class="Visium")
position_image = fst.image_coord_merge(position_image, position, dataset_class = 'Visium')
position_order = fst.update_st_coord(position_image)
print("The coords of image patch (merged): ", position_order.shape)
print(position_order.head())
## Save the position data -- train model used
link_data = 'OrderData/'
save_dir = os.path.join(demo_data, link_data)
os.makedirs(save_dir, exist_ok=True)
# position_order.to_csv(f'{demo_data}{link_data}position_order.csv', index=False, header=False)
## Order matrix row/barcode, spatial coordinates by image coordinates
spotID_order = np.array(position_image[0])
gene_hv = np.array(adata.var_names)
matrix_order, matrix_order_df = fst.sort_matrix(adata, position_image, spotID_order, gene_hv)
## Save gene expression mateix -- train model used
# np.save(f'{demo_data}{link_data}matrix_order.npy', matrix_order_df.T)
The coords of ST spot: (4992, 6)
Image embedding file: ['HCC_10018_14318.pth', 'HCC_10020_12925.pth', 'HCC_10020_13273.pth']
The coords of image patch (merged): (3354, 4)
pixel_y pixel_x array_row array_col
0 10018 14318 81 67
1 10020 12925 81 59
2 10020 13273 81 61
3 10021 12577 81 57
4 10022 11881 81 53
(3354, 1073)
[8]:
print(position_order['pixel_y'].max(), position_order['pixel_x'].max())
14624 16065
[9]:
## Update adata with new barcodes and lications
adata = fst.update_adata_coord(adata, matrix_order_df, position_image)
print(adata)
adata_count = adata.copy()
adata_norml = fst.adata_preprocess(adata.copy(), normalize=True)
## Save adata
save_data = 'SaveData/'
save_dir = os.path.join(demo_data, save_data)
os.makedirs(save_dir, exist_ok=True)
# adata_count.write_h5ad(f'{demo_data}{save_data}adata_count.h5ad')
# adata_norml.write_h5ad(f'{demo_data}{save_data}adata_norml.h5ad')
AnnData object with n_obs × n_vars = 3354 × 1073
obs: 'in_tissue', 'array_row', 'array_col', 'pxl_row_in_fullres', 'pxl_col_in_fullres', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'pct_counts_in_top_50_genes', 'pct_counts_in_top_100_genes', 'pct_counts_in_top_200_genes', 'pct_counts_in_top_500_genes', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt'
var: 'gene_ids', 'feature_types', 'mt', 'n_cells_by_counts', 'mean_counts', 'log1p_mean_counts', 'pct_dropout_by_counts', 'total_counts', 'log1p_total_counts', 'n_cells'
uns: 'spatial'
obsm: 'spatial'
Note The above two cells only need to be run once. They aim to generate the ordered gene expression profile matrix_order.npy and the ordered ST spot coordinates position_order.csv according to image pixel coordinates, and save the original and normalized gene expression data adata_count.h5ad,adata_norml.h5ad.
[10]:
fig_path = 'Figures/'
save_dir = os.path.join(demo_data, fig_path)
os.makedirs(save_dir, exist_ok=True)
fstplt.gene_expr(adata_count, adata_count.to_df(), gene_selet='SPP1', marker='o',
s=8, figsize=(5, 4), cnt_color='Spectral_r', save_path=f'{demo_data}{fig_path}SPP1_expr_count1.pdf')
fstplt.gene_expr(adata_norml, adata_norml.to_df(), gene_selet='SPP1', marker='o',
s=8, figsize=(5, 4), cnt_color='Spectral_r', save_path=f'{demo_data}{fig_path}SPP1_expr_norml1.pdf')
# fstplt.gene_expr(adata_count, adata_count.to_df(), gene_selet='THBS1', marker='o',
# s=8, figsize=(5, 4), cnt_color='Spectral_r', save_path=f'{demo_data}{fig_path}THBS1_expr_count1.pdf')
# fstplt.gene_expr(adata_norml, adata_norml.to_df(), gene_selet='THBS1', marker='o',
# s=8, figsize=(5, 4), cnt_color='Spectral_r', save_path=f'{demo_data}{fig_path}THBS1_expr_norml1.pdf')
fstplt.gene_expr(adata_count, adata_count.to_df(), gene_selet='CD274', marker='o',
s=8, figsize=(5, 4), cnt_color='Spectral_r', save_path=f'{demo_data}{fig_path}CD274_expr_count1.pdf')
fstplt.gene_expr(adata_norml, adata_norml.to_df(), gene_selet='CD274', marker='o',
s=8, figsize=(5, 4), cnt_color='Spectral_r', save_path=f'{demo_data}{fig_path}CD274_expr_norml1.pdf')
fstplt.gene_expr(adata_count, adata_count.to_df(), gene_selet='PDCD1', marker='o',
s=8, figsize=(5, 4), cnt_color='Spectral_r', save_path=f'{demo_data}{fig_path}PDCD1_expr_count1.pdf')
fstplt.gene_expr(adata_norml, adata_norml.to_df(), gene_selet='PDCD1', marker='o',
s=8, figsize=(5, 4), cnt_color='Spectral_r', save_path=f'{demo_data}{fig_path}PDCD1_expr_norml1.pdf')
[25]:
fstplt.gene_expr(adata_count, adata_count.to_df(), gene_selet='CCR1', marker='o',
s=8, figsize=(5, 4), cnt_color='Spectral_r')
2. Train model and Internal validation
Input: three input file
image_embed_path: image embedding from Setp0Image_feature_extraction.pyspatial_pos_path: ordered ST spot coords, inf'{demo_data}{link_data}position_order.csvreduced_mtx_path: ordered gene expression, inf'{demo_data}{link_data}matrix_order.npy
[11]:
image_embed_path = f'{demo_data}ImgEmbeddings/Bpth_112_14/*.pth'
spatial_pos_path = f'{demo_data}{link_data}position_order.csv'
reduced_mtx_path = f'{demo_data}{link_data}/matrix_order.npy'
2.1 Train and test model on within spot, with plitting: 80% for train and 20% for test
[15]:
## add params
params['n_input_matrix'] = len(gene_hv)
## init the model
model = fst.FineSTModel(n_input_matrix=params['n_input_matrix'],
n_input_image=params['n_input_image'],
n_encoder_hidden_matrix=params["n_encoder_hidden_matrix"],
n_encoder_hidden_image=params["n_encoder_hidden_image"],
n_encoder_latent=params["n_encoder_latent"],
n_projection_hidden=params["n_projection_hidden"],
n_projection_output=params["n_projection_output"],
n_encoder_layers=params["n_encoder_layers"]).to(device)
## Load the data
train_loader, test_loader = fst.build_loaders(batch_size=params['batch_size'],
image_embed_path=image_embed_path,
spatial_pos_path=spatial_pos_path,
reduced_mtx_path=reduced_mtx_path,
image_clacss='Virchow2',
dataset_class='Visium64')
## Set optimizer
optimizer = torch.optim.SGD(model.parameters(), lr=params['inital_learning_rate'],
momentum=0.9, weight_decay=5e-4)
## Load loss function
l = fst.ContrastiveLoss(temperature=params['temperature'],
w1=params['w1'], w2=params['w2'], w3=params['w3'], w4=params['w4'])
## train model
(dir_name,
train_losses, test_losses,
best_epoch, best_loss) = fst.train_model_fst(params, model, train_loader, test_loader,
optimizer, l, dir_name, logger, dataset_class='Visium64')
CUDA is available. GPU: NVIDIA A100-PCIE-40GB
Finished loading all files
train/test split completed
2683 671
***** Finished building loaders *****
[2025-06-27 01:09:27] INFO - Begin Training ...
[2025-06-27 01:09:27] INFO - epoch [1/50]
Epoch: 1
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 52.62it/s]
--- 77.35094261169434 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 52.23it/s]
[2025-06-27 01:10:50] INFO - Best epoch_loss: [0: 1.6787]
[2025-06-27 01:10:50] INFO - epoch [2/50]
Epoch: 2
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 46.08it/s]
--- 80.54439234733582 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 64.99it/s]
[2025-06-27 01:12:15] INFO - Best epoch_loss: [1: 1.6349]
[2025-06-27 01:12:15] INFO - epoch [3/50]
Epoch: 3
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 70.87it/s]
--- 76.63284683227539 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 66.81it/s]
[2025-06-27 01:13:35] INFO - Best epoch_loss: [2: 1.5319]
[2025-06-27 01:13:35] INFO - epoch [4/50]
Epoch: 4
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 65.61it/s]
--- 87.64106297492981 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 50.29it/s]
[2025-06-27 01:15:07] INFO - Best epoch_loss: [3: 1.5085]
[2025-06-27 01:15:07] INFO - epoch [5/50]
Epoch: 5
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 71.46it/s]
--- 69.16080522537231 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 45.55it/s]
[2025-06-27 01:16:21] INFO - Best epoch_loss: [4: 1.4906]
[2025-06-27 01:16:21] INFO - epoch [6/50]
Epoch: 6
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 42.40it/s]
--- 70.99906039237976 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 62.97it/s]
[2025-06-27 01:17:36] INFO - Best epoch_loss: [5: 1.4803]
[2025-06-27 01:17:36] INFO - epoch [7/50]
Epoch: 7
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 69.35it/s]
--- 74.13247632980347 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 45.28it/s]
[2025-06-27 01:18:54] INFO - Best epoch_loss: [6: 1.4629]
[2025-06-27 01:18:54] INFO - epoch [8/50]
Epoch: 8
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 71.85it/s]
--- 76.93378639221191 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 67.57it/s]
[2025-06-27 01:20:14] INFO - epoch [9/50]
Epoch: 9
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 70.52it/s]
--- 70.62517213821411 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 47.11it/s]
[2025-06-27 01:21:30] INFO - Best epoch_loss: [8: 1.4567]
[2025-06-27 01:21:30] INFO - epoch [10/50]
Epoch: 10
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 46.82it/s]
--- 70.62352538108826 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 59.28it/s]
[2025-06-27 01:22:44] INFO - Best epoch_loss: [9: 1.4529]
[2025-06-27 01:22:44] INFO - epoch [11/50]
Epoch: 11
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 53.25it/s]
--- 72.42334389686584 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 56.19it/s]
[2025-06-27 01:24:00] INFO - Best epoch_loss: [10: 1.4508]
[2025-06-27 01:24:00] INFO - epoch [12/50]
Epoch: 12
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 73.60it/s]
--- 79.01900458335876 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 57.85it/s]
[2025-06-27 01:25:24] INFO - epoch [13/50]
Epoch: 13
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 63.10it/s]
--- 68.16388893127441 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 56.13it/s]
[2025-06-27 01:26:36] INFO - Best epoch_loss: [12: 1.4486]
[2025-06-27 01:26:36] INFO - epoch [14/50]
Epoch: 14
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 57.20it/s]
--- 71.6641001701355 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 58.73it/s]
[2025-06-27 01:27:52] INFO - epoch [15/50]
Epoch: 15
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 57.33it/s]
--- 74.0404303073883 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 49.19it/s]
[2025-06-27 01:29:11] INFO - Best epoch_loss: [14: 1.4457]
[2025-06-27 01:29:11] INFO - epoch [16/50]
Epoch: 16
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 63.44it/s]
--- 74.35751104354858 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 58.46it/s]
[2025-06-27 01:30:30] INFO - Best epoch_loss: [15: 1.4419]
[2025-06-27 01:30:30] INFO - epoch [17/50]
Epoch: 17
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 55.58it/s]
--- 70.4228765964508 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 47.16it/s]
[2025-06-27 01:31:45] INFO - epoch [18/50]
Epoch: 18
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 58.99it/s]
--- 77.42666721343994 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 71.13it/s]
[2025-06-27 01:33:06] INFO - Best epoch_loss: [17: 1.4393]
[2025-06-27 01:33:06] INFO - epoch [19/50]
Epoch: 19
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 70.26it/s]
--- 81.8407289981842 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 36.43it/s]
[2025-06-27 01:34:32] INFO - Best epoch_loss: [18: 1.4382]
[2025-06-27 01:34:32] INFO - epoch [20/50]
Epoch: 20
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 71.66it/s]
--- 73.52662658691406 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 43.35it/s]
[2025-06-27 01:35:51] INFO - Best epoch_loss: [19: 1.4360]
[2025-06-27 01:35:51] INFO - epoch [21/50]
Epoch: 21
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 58.64it/s]
--- 71.90704107284546 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 54.15it/s]
[2025-06-27 01:37:07] INFO - Best epoch_loss: [20: 1.4339]
[2025-06-27 01:37:07] INFO - epoch [22/50]
Epoch: 22
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 53.85it/s]
--- 78.12762641906738 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 45.05it/s]
[2025-06-27 01:38:30] INFO - Best epoch_loss: [21: 1.4308]
[2025-06-27 01:38:30] INFO - epoch [23/50]
Epoch: 23
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 47.38it/s]
--- 75.41903877258301 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 58.24it/s]
[2025-06-27 01:39:49] INFO - Best epoch_loss: [22: 1.4291]
[2025-06-27 01:39:49] INFO - epoch [24/50]
Epoch: 24
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 68.72it/s]
--- 74.24105286598206 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 41.02it/s]
[2025-06-27 01:41:08] INFO - Best epoch_loss: [23: 1.4263]
[2025-06-27 01:41:08] INFO - epoch [25/50]
Epoch: 25
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 46.79it/s]
--- 71.68878555297852 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 46.19it/s]
[2025-06-27 01:42:24] INFO - epoch [26/50]
Epoch: 26
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 47.23it/s]
--- 72.92937278747559 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 51.65it/s]
[2025-06-27 01:43:41] INFO - epoch [27/50]
Epoch: 27
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 69.41it/s]
--- 76.12817716598511 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 53.67it/s]
[2025-06-27 01:45:01] INFO - Best epoch_loss: [26: 1.4232]
[2025-06-27 01:45:01] INFO - epoch [28/50]
Epoch: 28
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 74.46it/s]
--- 79.23716235160828 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 37.92it/s]
[2025-06-27 01:46:25] INFO - epoch [29/50]
Epoch: 29
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 52.85it/s]
--- 72.1925539970398 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 53.49it/s]
[2025-06-27 01:47:42] INFO - epoch [30/50]
Epoch: 30
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 57.12it/s]
--- 76.20742964744568 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 69.32it/s]
[2025-06-27 01:49:01] INFO - Best epoch_loss: [29: 1.4222]
[2025-06-27 01:49:01] INFO - epoch [31/50]
Epoch: 31
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 76.73it/s]
--- 75.62271332740784 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 40.98it/s]
[2025-06-27 01:50:22] INFO - Best epoch_loss: [30: 1.4178]
[2025-06-27 01:50:22] INFO - epoch [32/50]
Epoch: 32
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 39.18it/s]
--- 71.54358172416687 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 47.57it/s]
[2025-06-27 01:51:39] INFO - epoch [33/50]
Epoch: 33
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 56.55it/s]
--- 74.94433379173279 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 66.82it/s]
[2025-06-27 01:52:58] INFO - Best epoch_loss: [32: 1.4178]
[2025-06-27 01:52:58] INFO - epoch [34/50]
Epoch: 34
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 71.16it/s]
--- 80.91631603240967 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 58.33it/s]
[2025-06-27 01:54:22] INFO - Best epoch_loss: [33: 1.4170]
[2025-06-27 01:54:22] INFO - epoch [35/50]
Epoch: 35
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 69.69it/s]
--- 77.2474627494812 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 47.22it/s]
[2025-06-27 01:55:45] INFO - Best epoch_loss: [34: 1.4152]
[2025-06-27 01:55:45] INFO - epoch [36/50]
Epoch: 36
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 53.90it/s]
--- 72.98612236976624 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 46.57it/s]
[2025-06-27 01:57:02] INFO - epoch [37/50]
Epoch: 37
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 53.64it/s]
--- 76.44309687614441 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 50.95it/s]
[2025-06-27 01:58:23] INFO - Best epoch_loss: [36: 1.4114]
[2025-06-27 01:58:23] INFO - epoch [38/50]
Epoch: 38
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 53.04it/s]
--- 75.97451972961426 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 45.20it/s]
[2025-06-27 01:59:44] INFO - epoch [39/50]
Epoch: 39
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 66.11it/s]
--- 80.21203470230103 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 63.27it/s]
[2025-06-27 02:01:08] INFO - epoch [40/50]
Epoch: 40
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 60.73it/s]
--- 73.20628809928894 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 53.72it/s]
[2025-06-27 02:02:26] INFO - Best epoch_loss: [39: 1.4098]
[2025-06-27 02:02:26] INFO - epoch [41/50]
Epoch: 41
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 37.40it/s]
--- 74.92529940605164 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 45.27it/s]
[2025-06-27 02:03:46] INFO - epoch [42/50]
Epoch: 42
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 57.49it/s]
--- 77.56578207015991 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 55.65it/s]
[2025-06-27 02:05:07] INFO - epoch [43/50]
Epoch: 43
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 66.97it/s]
--- 78.94461631774902 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 59.17it/s]
[2025-06-27 02:06:31] INFO - Best epoch_loss: [42: 1.4096]
[2025-06-27 02:06:31] INFO - epoch [44/50]
Epoch: 44
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 60.56it/s]
--- 71.25590491294861 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 58.89it/s]
[2025-06-27 02:07:47] INFO - Best epoch_loss: [43: 1.4081]
[2025-06-27 02:07:47] INFO - epoch [45/50]
Epoch: 45
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 53.90it/s]
--- 76.08900451660156 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 43.85it/s]
[2025-06-27 02:09:07] INFO - epoch [46/50]
Epoch: 46
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 71.57it/s]
--- 81.3759970664978 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 70.36it/s]
[2025-06-27 02:10:32] INFO - epoch [47/50]
Epoch: 47
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 64.48it/s]
--- 83.80668592453003 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 54.85it/s]
[2025-06-27 02:12:01] INFO - epoch [48/50]
Epoch: 48
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 54.96it/s]
--- 72.97136950492859 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 53.43it/s]
[2025-06-27 02:13:20] INFO - epoch [49/50]
Epoch: 49
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 53.86it/s]
--- 77.22329640388489 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 63.66it/s]
[2025-06-27 02:14:40] INFO - epoch [50/50]
Epoch: 50
train model
100%|██████████████████████████████████████████████████████████████████| 13/13 [00:00<00:00, 70.25it/s]
--- 86.04526400566101 seconds ---
test model
100%|████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00, 68.47it/s]
[2025-06-27 02:16:10] INFO - Done!, Best epoch_loss: [43: 1.4081]
[2025-06-27 02:16:10] INFO - Finished Training
Done!, Best epoch_loss: [43: 1.4081]
Training epoch time: 4003.5061 seconds
[16]:
fst.loss_curve(train_losses, test_losses, best_epoch, best_loss, max_step=5, min_step=1,
fig_size=(5, 4), format='svg', save_path=f'{demo_data}{fig_path}HCC_loss_curve.svg')
Note The above cell for section 2.1 needs to be run more than once by setting different parameters in the given .json file, if the correlation in section 3.4 is not satisfying. If one obtains a good performance, one can save the dir_name and then directly use it to load the trained model and infer gene expression, and doesn’t need to run section 2.1 again. Recommended parameter fine-tuning and its range: training_epoch: 40~60; temperature: 0.01-0.05.
3. Inference, Imputation and Evaluation on “within spot”
3.1 Infer the gene expression of within spots
You can use the dir_name obtained from section 2.1, here we use our trained logging/20250626161134717029 20250627010909403881 with inital_learning_rate: 0.2 at Goole Drive for paper results repeated.
[17]:
!pwd
/mnt/lingyu/nfs_share2/Python/FineST/FineST_demo
[17]:
dir_name
[17]:
'logging/20250626234007931504'
[12]:
dir_name = 'logging/20250627010909403881'
[13]:
## Load the trained model
model = fst.load_model(dir_name, parameter_file_path, gene_hv)
## Load test data
test_loader = fst.build_loaders_inference(batch_size=adata.shape[0],
image_embed_path=image_embed_path,
spatial_pos_path=spatial_pos_path,
reduced_mtx_path=reduced_mtx_path,
image_clacss='Virchow2',
dataset_class='Visium64')
## Inference
(matrix_profile,
reconstructed_matrix,
recon_ref_adata_image_f2,
reconstructed_matrix_reshaped,
input_coord_all) = fst.infer_model_fst(model, test_loader, logger, dataset_class='Visium64')
CUDA is available. GPU: NVIDIA A100-PCIE-40GB
[2025-11-18 17:31:35] INFO - Running inference task...
Finished loading all files
***** Finished building loaders_inference *****
device cuda:0
***** Begin perform_inference: ******
100%|██████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 3.10it/s]
***** batch_size=adata.shape[0] doen't effect *****
torch.Size([3354, 1073])
torch.Size([3354, 64, 1280])
1
1
1
***** *****
Finished extractting test data
[2025-11-18 17:31:38] INFO - Running inference task DONE!
Inference within spots: 2.8029 seconds
Reconstructed_matrix_reshaped shape: torch.Size([214656, 1073])
3.2 Get coordinates of sub-spots from within spots
We use reshape_latent_image() function to convert the super-resolved gene expression with 16X or 64X relolution reconstructed_matrix_reshaped (torch.Size([21296, 596]) or torch.Size([85184, 596])) into tensor format reconstructed_matrix_reshaped_tensor (torch.Size([1331, 16, 596]) or torch.Size([1331, 64, 596])).
[14]:
## Get the sub-spot level gene expression of all genes on all within-spots
reconstructed_matrix_reshaped_tensor, _ = fst.reshape_latent_image(reconstructed_matrix_reshaped,
dataset_class='Visium64')
(_, _, all_spot_all_variable,
C2, adata_infer) = fst.subspot_coord_expr_adata(reconstructed_matrix_reshaped_tensor,
adata, gene_hv, patch_size=112, dataset_class='Visium64')
print(adata_infer)
_, adata_infer_reshape = fst.reshape_latent_image(torch.tensor(adata_infer.X), dataset_class='Visium64')
adata_infer_spot = fst.reshape2adata(adata, adata_infer_reshape, gene_hv)
print(adata_infer_spot)
## save adata
# adata_infer.write_h5ad(f'{demo_data}{save_data}adata_infer.h5ad') # adata 85184 × 596
# adata_infer_spot.write_h5ad(f'{demo_data}{save_data}adata_infer_spot.h5ad') # adata 1331 × 596
pixel_step (half of patch_size): 7.0
AnnData object with n_obs × n_vars = 214656 × 1073
obs: 'x', 'y'
obsm: 'spatial'
AnnData object with n_obs × n_vars = 3354 × 1073
obs: 'in_tissue', 'array_row', 'array_col', 'pxl_row_in_fullres', 'pxl_col_in_fullres', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'pct_counts_in_top_50_genes', 'pct_counts_in_top_100_genes', 'pct_counts_in_top_200_genes', 'pct_counts_in_top_500_genes', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt'
var: 'gene_ids', 'feature_types', 'mt', 'n_cells_by_counts', 'mean_counts', 'log1p_mean_counts', 'pct_dropout_by_counts', 'total_counts', 'log1p_total_counts', 'n_cells'
uns: 'spatial'
obsm: 'spatial'
Note Here, adata_infer is the inferred super-resolved gene expression data with 16X or 64X solution.
[15]:
## Plot the first_spot_first_variable spiltting
(first_spot_first_variable, C,
_, _, _) = fst.subspot_coord_expr_adata(reconstructed_matrix_reshaped_tensor, adata, gene_hv,
p=0, q=0, patch_size=112, dataset_class='Visium64')
print("first_spot_first_variable shape:", first_spot_first_variable.shape)
fstplt.subspot_expr(C, first_spot_first_variable, patch_size=112, dataset_class='Visium64',
marker='s', s=250, rotation=30, fig_size=(2.5, 2.5), format='svg',
# save_path=f'{demo_data}{fig_path}spot_splitting.svg')
save_path=None)
pixel_step (half of patch_size): 7.0
first_spot_first_variable shape: (64,)
[19]:
fstplt.gene_expr_compare(adata, "CCR1", adata_infer_reshape, gene_hv, s=50, cmap='Spectral_r', save_path=None)
fstplt.gene_expr_compare(adata, "THBS1", adata_infer_reshape, gene_hv, s=50, cmap='Spectral_r', save_path=None)
[20]:
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
from scipy.stats import pearsonr
from sklearn.linear_model import LinearRegression
def plot_regression_with_stats(
x, y,
weight=None,
moran_func=None,
x_label='X',
y_label='Y',
title='Regression Plot',
figsize=(5, 5),
dpi=150,
scatter_alpha=0.5,
line_color='red',
trans=False, format='pdf', save_path=None
):
"""
Plot regression with Pearson R, p-value, (optional) Moran R and regression equation.
All annotations are in English.
"""
x = np.array(x).flatten()
y = np.array(y).flatten()
## Pearson correlation
r, p = pearsonr(x, y)
## Linear regression
X_ = x.reshape(-1, 1)
model = LinearRegression().fit(X_, y)
slope = model.coef_[0]
intercept = model.intercept_
# ## Moran R if function and weight are provided
# moran_str = ''
# if (moran_func is not None) and (weight is not None):
# try:
# moran_r = moran_func(x, y, weight)[0].item()
# moran_str = f'Moran R = {moran_r:.2f}\n'
# except Exception as e:
# moran_str = ''
# print(f"Warning: Moran R calculation failed: {e}")
## Plotting
fig, ax = plt.subplots(figsize=figsize, dpi=dpi)
sns.regplot(
x=x, y=y, ax=ax,
scatter_kws={'alpha': scatter_alpha},
line_kws={'color': line_color}
)
## Compose annotation string
annotation = (
# f"{moran_str}"
f"Pearson R = {r:.2f}\n"
f"p-value = {p:.2e}\n"
f"Regression: y = {slope:.2f}x + {intercept:.2f}"
)
ax.text(
0.95, 0.00, annotation, color='black',
ha='right', va='bottom', transform=ax.transAxes
)
ax.set_title(title)
ax.set_xlabel(x_label)
ax.set_ylabel(y_label)
plt.tight_layout()
if save_path is not None:
plt.savefig(save_path, transparent=trans, format=format,
dpi=300, bbox_inches='tight')
plt.show()
[20]:
## Show marker gene expression -- only infer from image
fstplt.gene_expr_compare(adata, "IL18", adata_infer_reshape, gene_hv, s=50, save_path=None)
fstplt.gene_expr_compare(adata, "SPP1", adata_infer_reshape, gene_hv, s=50, save_path=None)
# fstplt.gene_expr_compare(adata, "CCR1", adata_infer_reshape, gene_hv, s=50, save_path=None)
# fstplt.gene_expr_compare(adata, "THBS1", adata_infer_reshape, gene_hv, s=50, save_path=None)
# fstplt.gene_expr(adata, adata_infer_reshape, gene_selet='SPP1', marker='o',
# s=8, figsize=(5, 4), save_path=f'{demo_data}{fig_path}SPP1_expr_infer_spot_within.pdf')
# fstplt.gene_expr(adata, adata_infer_reshape, gene_selet='IL18', marker='o',
# s=8, figsize=(5, 4), save_path=f'{demo_data}{fig_path}IL18_expr_infer_spot_within.pdf')
[22]:
## save the correlation -- -- only infer from image
fstplt.mean_cor_box(adata, adata_infer_reshape, logger, save_path=None)
fstplt.mean_cor_box(adata_norml, adata_infer_reshape, logger, save_path=None)
# fstplt.mean_cor_box(adata, adata_infer_reshape, logger, save_path=f'{demo_data}{fig_path}Boxplot_infer_cor_count.pdf')
# fstplt.mean_cor_box(adata_norml, adata_infer_reshape, logger, save_path=f'{demo_data}{fig_path}Boxplot_infer_cor_norml.pdf')
[2025-06-27 09:49:10] INFO - mean correlation of spots: 0.93184632417744
mean correlation of spots: 0.93184632417744
[2025-06-27 09:49:11] INFO - mean correlation of genes: 0.22118164568459048
mean correlation of genes: 0.22118164568459048
[2025-06-27 09:49:12] INFO - mean correlation of spots: 0.7258101153916254
mean correlation of spots: 0.7258101153916254
[2025-06-27 09:49:12] INFO - mean correlation of genes: 0.14001626192891992
mean correlation of genes: 0.14001626192891992
[23]:
## Pearson correlations
fstplt.cor_hist(adata, adata_infer_spot.to_df(),
fig_size=(5, 4), trans=False, format='svg',
# save_path=None)
save_path=f'{demo_data}{fig_path}Hist_infer_cor_count.svg')
Pearson correlations: 0.22118164568459045
3.3 Impute the gene expression of sub-spots
Input: two input files:
-
adata: original spot-level gene expression data (within spots).-
adata_infer: inferred super-resolved gene expression data (within spots).Input: three input files:
-
adata_impt: = w*adata_infer + (1-w)*sudo_adata, predicted super-resolved gene expression data (here w=0.98).-
sudo_adata: Imputed data using the nearest k neighbors of within spots (here k=6).-
data_impt: predicted super-resolved gene expression data in tensor form (within spots).[16]:
# adata_smooth = fst.impute_adata(adata, adata_infer, C2, gene_hv,
# dataset_class='Visium64', weight_exponent=2)
adata_smooth = fst.impute_adata(adata_norml, adata_infer, C2, gene_hv,
dataset_class='Visium64', weight_exponent=2)
Smoothing time: 69.8896 seconds
[17]:
adata_imput, data_impt = fst.weight_adata(adata_infer, adata_smooth, gene_hv, w=0.98)
print(adata_imput)
_, data_impt_reshape = fst.reshape_latent_image(data_impt, dataset_class='Visium64')
adata_imput_spot = fst.reshape2adata(adata, data_impt_reshape, gene_hv)
print(adata_imput_spot)
## save adata
# adata_imput.write_h5ad(f'{demo_data}{save_data}_adata_imput.h5ad') # adata: 214656 × 1073
# adata_imput_spot.write_h5ad(f'{demo_data}{save_data}adata_imput_spot.h5ad') # adata: 3354 × 1073
AnnData object with n_obs × n_vars = 214656 × 1073
obs: 'x', 'y'
uns: 'spatial'
obsm: 'spatial'
AnnData object with n_obs × n_vars = 3354 × 1073
obs: 'in_tissue', 'array_row', 'array_col', 'pxl_row_in_fullres', 'pxl_col_in_fullres', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'pct_counts_in_top_50_genes', 'pct_counts_in_top_100_genes', 'pct_counts_in_top_200_genes', 'pct_counts_in_top_500_genes', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt'
var: 'gene_ids', 'feature_types', 'mt', 'n_cells_by_counts', 'mean_counts', 'log1p_mean_counts', 'pct_dropout_by_counts', 'total_counts', 'log1p_total_counts', 'n_cells'
uns: 'spatial'
obsm: 'spatial'
[18]:
fstplt.gene_expr_compare(adata, "CD274", adata_infer_reshape, gene_hv, s=40, cmap='Spectral_r',
save_path=f'{demo_data}{fig_path}CD274_expr_infer_spot_within_count.pdf')
# save_path=None)
fstplt.gene_expr_compare(adata, "PDCD1", adata_infer_reshape, gene_hv, s=40, cmap='Spectral_r',
save_path=f'{demo_data}{fig_path}PDCD1_expr_infer_spot_within_count.pdf')
# save_path=None)
[19]:
fstplt.gene_expr_compare(adata_norml, "CD274", adata_infer_reshape, gene_hv, s=40, cmap='Spectral_r',
save_path=f'{demo_data}{fig_path}CD274_expr_infer_spot_within_norml.pdf')
# save_path=None)
fstplt.gene_expr_compare(adata_norml, "PDCD1", adata_infer_reshape, gene_hv, s=40, cmap='Spectral_r',
save_path=f'{demo_data}{fig_path}PDCD1_expr_infer_spot_within_norml.pdf')
# save_path=None)
[28]:
fstplt.gene_expr_compare(adata, "CCR1", adata_infer_reshape, gene_hv, s=40, cmap='Spectral_r',
save_path=f'{demo_data}{fig_path}CCR1_expr_infer_spot_within_count.pdf')
# save_path=None)
fstplt.gene_expr_compare(adata, "CCR2", adata_infer_reshape, gene_hv, s=40, cmap='Spectral_r',
save_path=f'{demo_data}{fig_path}CCR2_expr_infer_spot_within_count.pdf')
# save_path=None)
fstplt.gene_expr_compare(adata, "THBS1", adata_infer_reshape, gene_hv, s=40, cmap='Spectral_r',
save_path=f'{demo_data}{fig_path}THBS1_expr_infer_spot_within_count.pdf')
# save_path=None)
fstplt.gene_expr_compare(adata, "SPP1", adata_infer_reshape, gene_hv, s=40, cmap='Spectral_r',
save_path=f'{demo_data}{fig_path}SPP1_expr_infer_spot_within_count.pdf')
# save_path=None)
fstplt.gene_expr_compare(adata, "IL18", adata_infer_reshape, gene_hv, s=40, cmap='Spectral_r',
save_path=f'{demo_data}{fig_path}IL18_expr_infer_spot_within_count.pdf')
# save_path=None)
[26]:
fstplt.gene_expr_compare(adata_norml, "CCR1", adata_infer_reshape, gene_hv, s=40, cmap='Spectral_r',
save_path=f'{demo_data}{fig_path}CCR1_expr_infer_spot_within_norml.pdf')
# save_path=None)
fstplt.gene_expr_compare(adata_norml, "CCR2", adata_infer_reshape, gene_hv, s=40, cmap='Spectral_r',
save_path=f'{demo_data}{fig_path}CCR2_expr_infer_spot_within_norml.pdf')
# save_path=None)
fstplt.gene_expr_compare(adata_norml, "THBS1", adata_infer_reshape, gene_hv, s=40, cmap='Spectral_r',
save_path=f'{demo_data}{fig_path}THBS1_expr_infer_spot_within_norml.pdf')
# save_path=None)
fstplt.gene_expr_compare(adata_norml, "SPP1", adata_infer_reshape, gene_hv, s=40, cmap='Spectral_r',
save_path=f'{demo_data}{fig_path}SPP1_expr_infer_spot_within_norml.pdf')
# save_path=None)
fstplt.gene_expr_compare(adata_norml, "IL18", adata_infer_reshape, gene_hv, s=40, cmap='Spectral_r',
save_path=f'{demo_data}{fig_path}IL18_expr_infer_spot_within_norml.pdf')
# save_path=None)
[29]:
fstplt.gene_expr_compare(adata_norml, "SAA1", adata_infer_reshape, gene_hv, s=40, cmap='Spectral_r',
save_path=f'{demo_data}{fig_path}SAA1_expr_infer_spot_within_norml.pdf')
# save_path=None)
# fstplt.gene_expr_compare(adata_norml, "SAA2", adata_infer_reshape, gene_hv, s=40, cmap='Spectral_r',
# save_path=f'{demo_data}{fig_path}SAA2_expr_infer_spot_within_norml.pdf')
# # save_path=None)
[29]:
# plot_regression_with_stats(
# adata_norml.to_df()['CCR1'],
# adata_norml.to_df()['THBS1'],
# # weight=adata.obsp['weight'],
# # moran_func=sdm.utils.Moran_R,
# x_label='CCR1',
# y_label='THBS1',
# title='CCR1 VS THBS1',
# save_path=f'{demo_data}{fig_path}CCR1_VS_THBS1_Original.pdf'
# )
# plot_regression_with_stats(
# adata_imput_spot.to_df()['CCR1'],
# adata_imput_spot.to_df()['THBS1'],
# # weight=adata.obsp['weight'],
# # moran_func=sdm.utils.Moran_R,
# x_label='CCR1',
# y_label='THBS1',
# title='CCR1 VS THBS1',
# save_path=f'{demo_data}{fig_path}CCR1_VS_THBS1_FineST.pdf'
# )
[ ]:
fstplt.gene_expr_compare(adata_norml, "IL18", adata_infer_reshape, gene_hv, s=40, cmap='Spectral_r',
save_path=f'{demo_data}{fig_path}IL18_expr_infer_spot_within1.pdf')
# save_path=None)
fstplt.gene_expr_compare(adata_norml, "SPP1", adata_infer_reshape, gene_hv, s=40, cmap='Spectral_r',
save_path=f'{demo_data}{fig_path}SPP1_expr_infer_spot_within1.pdf')
# save_path=None)
[27]:
# plot_regression_with_stats(
# adata_norml.to_df()['IL18'],
# adata_norml.to_df()['SPP1'],
# # weight=adata.obsp['weight'],
# # moran_func=sdm.utils.Moran_R,
# x_label='IL18',
# y_label='SPP1',
# title='IL18 VS SPP1',
# save_path=f'{demo_data}{fig_path}IL18_VS_SPP1_Original.pdf'
# )
# plot_regression_with_stats(
# adata_imput_spot.to_df()['IL18'],
# adata_imput_spot.to_df()['SPP1'],
# # weight=adata.obsp['weight'],
# # moran_func=sdm.utils.Moran_R,
# x_label='IL18',
# y_label='SPP1',
# title='IL18 VS SPP1',
# save_path=f'{demo_data}{fig_path}IL18_VS_SPP1_FineST.pdf'
# )
Note For performance evaluation, we integrate super-resolved gene expression data_impt into spot-level gene expression data_impt_reshape (i.e., from 16x sub-spot level to spot level)
3.4 Visualization: selected gene (Original vs. FineST)
[30]:
fstplt.gene_expr_compare(adata, "SPP1", data_impt_reshape, gene_hv, s=50, save_path=None)
fstplt.gene_expr_compare(adata, "IL18", data_impt_reshape, gene_hv, s=50, save_path=None)
[31]:
fstplt.gene_expr(adata, data_impt_reshape, gene_selet='SPP1', marker='o',
s=8, figsize=(5, 4), save_path=f'{demo_data}{fig_path}SPP1_expr_imput_spot_within.pdf')
fstplt.gene_expr(adata, data_impt_reshape, gene_selet='IL18', marker='o',
s=8, figsize=(5, 4), save_path=f'{demo_data}{fig_path}IL18_expr_imput_spot_within.pdf')
3.4 Correlation: Selected gene、all spots and all genes、mean correlation box
The correlation can be used as a measure to evaluate the gene expression prediction performance. Here we provide three functions:
-
sele_gene_cor: show the selected gene correlation (Reconstructed by FineST vs. Original).-
mean_cor : calculate the mean correlation of all spots or all genes.-
mean_cor_box : box plot of correlation for all spots or all genes.[32]:
fstplt.sele_gene_cor(adata, data_impt_reshape, gene_hv, gene = "SPP1",
ylabel='FineST Expression', title = "SPP1 expression", size=5,
save_path=f'{demo_data}{fig_path}SPP1_cor_imput_spot_within.pdf')
fstplt.sele_gene_cor(adata, data_impt_reshape, gene_hv, gene = "IL18",
ylabel='FineST Expression', title = "IL18 expression", size=5,
save_path=f'{demo_data}{fig_path}IL18_cor_imput_spot_within.pdf')
[33]:
logger.info("Running Gene Correlation task...")
(pearson_cor_gene,
spearman_cor_gene,
cosine_sim_gene) = fst.mean_cor(adata, data_impt_reshape, 'reconf2', sample="gene")
logger.info("Running Gene Correlation task DINE!")
[2025-06-27 08:27:46] INFO - Running Gene Correlation task...
matrix1: (3354, 1073)
matrix2: (3354, 1073)
Mean Pearson correlation coefficient--reconf2: 0.6303
[2025-06-27 08:27:47] INFO - Running Gene Correlation task DINE!
Mean Spearman correlation coefficient--reconf2: 0.4470
Mean cosine similarity--reconf2: 0.5114
[34]:
fstplt.mean_cor_box(adata, data_impt_reshape, logger, save_path=f'{demo_data}{fig_path}Boxplot_imput_cor_count.pdf')
[2025-06-27 08:27:48] INFO - mean correlation of spots: 0.9395336865190809
mean correlation of spots: 0.9395336865190809
[2025-06-27 08:27:48] INFO - mean correlation of genes: 0.6302883615430775
mean correlation of genes: 0.6302883615430775
[33]:
## Pearson correlations
fstplt.cor_hist(adata, adata_imput_spot.to_df(),
fig_size=(5, 4), trans=False, format='svg',
# save_path=None)
save_path=f'{demo_data}{fig_path}Hist_impute_cor_count.svg')
Pearson correlations: 0.6302883615430775
[35]:
#####################
# release GUP memory
#####################
fst.release_gpu_tensors(globals())
# release_gpu_tensors(locals()) # If inside a function or class method
Releasing the following GPU tensors:
reconstructed_matrix_reshaped 878.62 MB shape=(214656, 1073)
reconstructed_matrix_reshaped_tensor 878.62 MB shape=(3354, 64, 1073)
Released GPU memory occupied by the above variables.
Note Here the adata_impt is the imputed spot-level gene expression from FineST, it only contains the given 1331 within spots. The comparisons are based on the original adata patient1_adata_orignal.h5ad vs. FineST’s predicted adata patient1_adata.h5ad.
4. Infer gene expression of “within spot” and “between spot”
Spot_interpolate
cd ~/FineST_demo
conda activate FineST
python ./demo/Spot_interpolation.py --position_path HCC_data/spatial/tissue_positions_list.csv
Step0: HE image feature extraction
time python ./demo/Image_feature_extraction.py \
--dataset NEW_HCC \
--position_path HCC_data/spatial/tissue_positions_list_add.csv \
--rawimage_path HCC_data/T11-V10F06-115-A1.jpg \
--scale_image False \
--method Virchow2 \
--patch_size 112 \
--output_img HCC_data/ImgEmbeddings/NEW_pth_112_14_image \
--output_pth HCC_data/ImgEmbeddings/NEW_pth_112_14 \
--logging HCC_data/ImgEmbeddings/
4.1 Get all spot coordinates (within spot & between spot) of image embeddings
Input: two input file pathways:
-
file_paths_spot: image embedding, of within spots, from Setp0 Image_feature_extraction.py.-
file_paths_between_spot: image embedding, of between spots, from Setp0 Image_feature_extraction.py.[15]:
## Add coords for each .pth file
file_paths_spot = os.listdir(f'{demo_data}ImgEmbeddings/Bpth_112_14/')
file_paths_between_spot = os.listdir(f'{demo_data}ImgEmbeddings/NEW_pth_112_14/')
file_paths_all = file_paths_spot + file_paths_between_spot
print("file_paths_spot number: ", len(file_paths_spot))
print("file_paths_between_spot number:", len(file_paths_between_spot))
print("file_paths_all number:", len(file_paths_all))
## Merge, sort and process file paths
data_all = fst.get_image_coord_all(file_paths_all)
position_order_allspot = pd.DataFrame(data_all, columns=['pixel_y', 'pixel_x'])
## save all spots
position_order_allspot.to_csv(f'{demo_data}{link_data}position_order_all.csv', index=False, header=False)
file_paths_spot number: 3354
file_paths_between_spot number: 9625
file_paths_all number: 12979
Note The above two cells for section 4.1 only need to be run once, aiming to generate the ordered ST spot coordinates position_order_all.csv according to image pixel coordinates.
4.2 Load all spot image embeddings with their coordinates
[16]:
import glob
file_paths_spot = f'{demo_data}ImgEmbeddings/Bpth_112_14/*.pth'
file_paths_between_spot = f'{demo_data}ImgEmbeddings/NEW_pth_112_14/*.pth'
spatial_pos_path=f'{demo_data}{link_data}position_order_all.csv'
4.4 Load the trained model to infer all spots
[17]:
dir_name ='logging/20250627010909403881'
[18]:
## Load the trained model
model = fst.load_model(dir_name, parameter_file_path, gene_hv)
## load all spot (within and between) spots data
all_dataset = fst.build_loaders_inference_allimage(batch_size=len(file_paths_all),
file_paths_spot=file_paths_spot,
file_paths_between_spot=file_paths_between_spot,
spatial_pos_path=spatial_pos_path,
dataset_class='Visium64')
## inference
logger.info("Running inference tesk between spot...")
start_time = time.time()
(recon_ref_adata_image_f2,
reconstructed_matrix_reshaped,
representation_image_reshape_between_spot,
input_image_exp_between_spot,
input_coord_all) = fst.perform_inference_image_between_spot(model, all_dataset, dataset_class='Visium64')
print("--- %s seconds ---" % (time.time() - start_time))
## print
print("recon_ref_adata_image_f2:", recon_ref_adata_image_f2.shape)
logger.info("Running inference tesk between spot DONE!")
CUDA is available. GPU: NVIDIA A100-PCIE-40GB
***** Building loaders_inference between spot *****
[2025-06-27 10:42:22] INFO - Running inference tesk between spot...
Finished loading all files
***** Finished building loaders_inference *****
device cuda:0
***** Begin perform_inference: ******
100%|████████████████████████████████████████████████████████████████████| 1/1 [00:01<00:00, 1.07s/it]
***** batch_size=adata.shape[0] *****
torch.Size([12979, 64, 1280])
1
Finished extractting image_between_spot data
[2025-06-27 10:42:23] INFO - Running inference tesk between spot DONE!
--- 1.31569504737854 seconds ---
recon_ref_adata_image_f2: (12979, 1073)
4.5 Visualization all spots
[19]:
## process_and_check_duplicates of the coordinates of all spots (within and between)
spatial_loc_all = fst.get_allspot_coors(input_coord_all)
print(spatial_loc_all)
Are there any duplicate rows? : False
[[10018. 14318. ]
[10020. 12925. ]
[10020. 13273. ]
...
[ 9987.5 3612. ]
[ 9988. 3264. ]
[ 9988. 3438. ]]
[21]:
fstplt.gene_expr_allspots("SPP1", spatial_loc_all, recon_ref_adata_image_f2,
gene_hv, 'Inferred all spot', s=5, marker='o',
figsize=(8, 7), cmap='Spectral_r', save_path=None)
fstplt.gene_expr_allspots("IL18", spatial_loc_all, recon_ref_adata_image_f2,
gene_hv, 'Inferred all spot', s=5, marker='o',
figsize=(8, 7), cmap='Spectral_r', save_path=None)
SPP1 gene expression dim: (12979, 1)
SPP1 gene expression:
[[0.00843153]
[0.01055692]
[0.00843315]
...
[0.02241297]
[0.05169463]
[0.02044788]]
IL18 gene expression dim: (12979, 1)
IL18 gene expression:
[[0.00735867]
[0.00794413]
[0.00701345]
...
[0.01055571]
[0.03799286]
[0.00964013]]
[22]:
fstplt.gene_expr_allspots("SPP1", spatial_loc_all, recon_ref_adata_image_f2,
gene_hv, 'Inferred all spot', s=1.5, marker='s',
figsize=(5, 4), save_path=f'{demo_data}{fig_path}SPP1_expr_infer_allspot.pdf')
fstplt.gene_expr_allspots("IL18", spatial_loc_all, recon_ref_adata_image_f2,
gene_hv, 'Inferred all spot', s=1.5, marker='s',
figsize=(5, 4), save_path=f'{demo_data}{fig_path}IL18_expr_infer_allspot.pdf')
SPP1 gene expression dim: (12979, 1)
SPP1 gene expression:
[[0.00843153]
[0.01055692]
[0.00843315]
...
[0.02241297]
[0.05169463]
[0.02044788]]
IL18 gene expression dim: (12979, 1)
IL18 gene expression:
[[0.00735867]
[0.00794413]
[0.00701345]
...
[0.01055571]
[0.03799286]
[0.00964013]]
4.6 Visualization all sub-spots
[23]:
## Get the sub-spot level gene expression of all genes in all spots
reconstructed_matrix_reshaped_tensor, _ = fst.reshape_latent_image(reconstructed_matrix_reshaped,
dataset_class='Visium64')
print("reconstructed_matrix_reshaped_tensor size: ", reconstructed_matrix_reshaped_tensor.shape)
(_, _, all_spot_all_variable,
C2, adata_infer_all) = fst.subspot_coord_expr_adata(reconstructed_matrix_reshaped_tensor,
spatial_loc_all, gene_hv, patch_size=112, dataset_class="Visium64")
print(adata_infer_all)
reconstructed_matrix_reshaped_tensor size: torch.Size([12979, 64, 1073])
pixel_step (half of patch_size): 7.0
AnnData object with n_obs × n_vars = 830656 × 1073
obs: 'x', 'y'
obsm: 'spatial'
[24]:
fstplt.gene_expr_allspots("SPP1", C2, all_spot_all_variable, gene_hv,
'Inferred all spot', marker='s', s=0.3,
figsize=(15, 12), save_path=f'{demo_data}{fig_path}SPP1_expr_infer_allsubspot.pdf')
# fstplt.gene_expr_allspots("IL18", C2, all_spot_all_variable, gene_hv,
# 'Inferred all spot', marker='s', s=0.3,
# figsize=(15, 12), save_path=f'{demo_data}{fig_path}IL18_expr_infer_allsubspot.pdf')
SPP1 gene expression dim: (830656, 1)
SPP1 gene expression:
[[0.00870922]
[0.00938632]
[0.00817643]
...
[0.00841701]
[0.00790322]
[0.00843024]]
[25]:
#####################
# release GUP memory
#####################
fst.release_gpu_tensors(globals())
# release_gpu_tensors(locals()) # If inside a function or class method
Releasing the following GPU tensors:
input_image_exp_between_spot 4055.94 MB shape=(830656, 1280)
reconstructed_matrix_reshaped 3400.02 MB shape=(830656, 1073)
reconstructed_matrix_reshaped_tensor 3400.02 MB shape=(12979, 64, 1073)
Released GPU memory occupied by the above variables.
5. Imputate sub-spot gene expression using measured spot expression
5.1 Impute the gene expression of all sub-spots from all spots (within spots and between spots)
[26]:
# adata_smooth_all = fst.impute_adata(adata, adata_infer_all, C2, gene_hv, dataset_class='Visium64', weight_exponent=2)
adata_smooth_all = fst.impute_adata(adata_norml, adata_infer_all, C2, gene_hv, dataset_class='Visium64', weight_exponent=2)
print(adata_smooth_all)
Smoothing time: 319.0520 seconds
AnnData object with n_obs × n_vars = 830656 × 1073
obs: 'x', 'y'
var: 'gene_ids', 'feature_types', 'mt', 'n_cells_by_counts', 'mean_counts', 'log1p_mean_counts', 'pct_dropout_by_counts', 'total_counts', 'log1p_total_counts', 'n_cells'
uns: 'spatial'
obsm: 'spatial'
[27]:
adata_impt_all, data_impt_all = fst.weight_adata(adata_infer_all, adata_smooth_all, gene_hv, w=0.98)
print(adata_impt_all)
## save adata
adata_impt_all.write_h5ad(f'{demo_data}{save_data}adata_imput_all_subspot.h5ad') # aadata: 830656 × 596
AnnData object with n_obs × n_vars = 830656 × 1073
obs: 'x', 'y'
uns: 'spatial'
obsm: 'spatial'
Note Here, adata_impt_all is the imputed sub-spot-level gene expression from FineST. patient1_adata_all.h5ad contains about 1331x4x16 sub-spots. (1331 within spots, 1331x3 between spots, 16x resolution).
5.2 Convert super-resolved gene expression to spot-level and save the imputated spot-level adata
[28]:
_, adata_impt_all_reshape = fst.reshape_latent_image(data_impt_all, dataset_class='Visium64')
adata_impt_spot = fst.reshape2adata(adata, adata_impt_all_reshape, gene_hv, spatial_loc_all)
## save adata
adata_impt_spot.write_h5ad(f'{demo_data}{save_data}adata_imput_all_spot.h5ad') # adata: 5039 × 596
Note Here, adata_impt_spot is the imputed sub-spot-level gene expression from FineST. patient1_adata_all_spot.h5ad contains about 1331x4 spots. (1331 within spots, 1331x3 between spots), it is integrated from patient1_adata_all.h5ad using the reshape_latent_image function.
5.3 Visualization: gene at all spot
[29]:
fstplt.gene_expr_allspots("SPP1", spatial_loc_all, adata_impt_all_reshape,
gene_hv, 'FineST all spot', marker='o', s=2.5,
figsize=(5, 4), cmap='Spectral_r', save_path=f'{demo_data}{fig_path}SPP1_expr_imput_allspot1.pdf')
fstplt.gene_expr_allspots("IL18", spatial_loc_all, adata_impt_all_reshape,
gene_hv, 'FineST all spot', marker='o', s=2.5,
figsize=(5, 4), cmap='Spectral_r', save_path=f'{demo_data}{fig_path}IL18_expr_imput_allspot1.pdf')
SPP1 gene expression dim: (12979, 1)
SPP1 gene expression:
[[0.00877749]
[0.01101463]
[0.00826448]
...
[0.03007508]
[0.06530798]
[0.03069134]]
IL18 gene expression dim: (12979, 1)
IL18 gene expression:
[[0.0074036 ]
[0.02174616]
[0.00687318]
...
[0.01149232]
[0.04141673]
[0.01114439]]
5.4 Visualization: gene at all sub-spot
[ ]:
fstplt.gene_expr_allspots("SPP1", C2, adata_impt_all.X, gene_hv,
'FineST all sub-spot', marker='o', s=0.3, figsize=(15, 12),
save_path=f'{demo_data}{fig_path}SPP1_expr_imput_allsubspot.pdf')
fstplt.gene_expr_allspots("IL18", C2, adata_impt_all.X, gene_hv,
'FineST all sub-spot', marker='o', s=0.3, figsize=(15, 12),
save_path=f'{demo_data}{fig_path}IL18_expr_imput_allsubspot.pdf')
[31]:
#####################
# release GUP memory
#####################
fst.release_gpu_tensors(globals())
# release_gpu_tensors(locals()) # If inside a function or class method
No GPU tensors detected.
6. Infer and impute gene expression at “single-cell” level
Note Section 6 directly infers super-resolved gene expression based on image features from nuclei segmentation. The image embeddings sc_pth_14_14 is based on adata_imput_all_spot.adata from Nuclei segmentation . The results NucleiSegments can be downloaded Goole Drive.
Nuclei segmentation
time python ./demo/StarDist_nuclei_segmente.py \
--tissue HCC_allspot \
--out_dir HCC_data/NucleiSegments \
--adata_path HCC_data/SaveData/adata_imput_all_spot.h5ad \
--img_path HCC_data/T11-V10F06-115-A1.jpg
Later
time python ./demo/StarDist_nuclei_segmente.py \
--tissue HCC_allspot_p075 \
--out_dir HCC_data/NucleiSegments \
--adata_path HCC_data/SaveData/adata_imput_all_spot.h5ad \
--img_path HCC_data/T11-V10F06-115-A1.jpg \
--prob_thresh 0.75
Step0: HE image feature extraction
time python ./demo/Image_feature_extraction.py \
--dataset sc_HCC \
--position_path HCC_data/NucleiSegments/HCC_allspot/position_all_tissue_sc.csv \
--rawimage_path HCC_data/T11-V10F06-115-A1.jpg \
--scale_image False \
--method Virchow2 \
--patch_size 14 \
--output_img HCC_data/ImgEmbeddings/sc_pth_14_14_image \
--output_pth HCC_data/ImgEmbeddings/sc_pth_14_14 \
--logging HCC_data/ImgEmbeddings/
use threshold 0.75 to shorten the nuclei number
time python ./demo/Image_feature_extraction.py \
--dataset sc_HCC_p075 \
--position_path HCC_data/NucleiSegments/HCC_allspot_p075/position_all_tissue_sc.csv \
--rawimage_path HCC_data/T11-V10F06-115-A1.jpg \
--scale_image False \
--method Virchow2 \
--patch_size 14 \
--output_img HCC_data/ImgEmbeddings/scp075_pth_14_14_image \
--output_pth HCC_data/ImgEmbeddings/scp075_pth_14_14 \
--logging HCC_data/ImgEmbeddings/
6.1 Infer the gene expression of at single-cell resolution
You can use the dir_name just obtained from Section 2 (Model training), but just like Section 3, here we use our trained logging/20250622170940952737 for paper results repeated.
[16]:
dir_name ='logging/20250627010909403881'
[34]:
# file_paths_sc = sorted(os.listdir(f'{demo_data}/ImgEmbeddings/sc_pth_14_14/'))
file_paths_sc = sorted(os.listdir(f'{demo_data}/ImgEmbeddings/scp075_pth_14_14/'))
print("Image embedding file: ", file_paths_sc[:2])
## Image patch position
data_all_sc = fst.get_image_coord_all(file_paths_sc)
spatial_loc_sc = pd.DataFrame(data_all_sc, columns=['pixel_y', 'pixel_x'])
## save .csv file
spatial_loc_sc.to_csv(f'{demo_data}{link_data}/position_order_sc.csv', index=False, header=False)
Image embedding file: ['sc_HCC_p075_10000.0_5201.0.pth', 'sc_HCC_p075_10000.0_5506.0.pth']
[35]:
## opt - 66029; p075 - 14647
len(file_paths_sc)
[35]:
14647
Note The above two cells for section 6.1 only need to be run once, aiming to generate the ordered ST spot coordinates position_order_sc.csv according to image pixel coordinates.
[36]:
import glob
image_embed_path_sc = f'{demo_data}/ImgEmbeddings/scp075_pth_14_14/*.pth'
spatial_pos_path_sc = f'{demo_data}{link_data}/position_order_sc.csv'
[19]:
# image_paths = glob.glob(image_embed_path_sc)
# image_paths.sort()
# image_paths
[20]:
# spatial_pos_csv = pd.read_csv(spatial_pos_path_sc, sep=",", header=None)
# spatial_pos_csv
14624 16065
[37]:
## load models
model = fst.load_model(dir_name, parameter_file_path, gene_hv)
## load all data
all_dataset_sc = fst.build_loaders_inference_allimage(batch_size=len(file_paths_sc),
file_paths_spot=image_embed_path_sc,
spatial_pos_path=spatial_pos_path_sc,
dataset_class='VisiumSC')
## inference
logger.info("Running inference tesk single-nuclei...")
start_time = time.time()
(recon_ref_adata_image_f2,
reconstructed_matrix_reshaped,
representation_image_reshape_between_spot,
input_image_exp_between_spot,
input_coord_all) = fst.perform_inference_image_between_spot(model, all_dataset_sc, dataset_class='VisiumSC')
print("--- %s seconds ---" % (time.time() - start_time))
## print
print("recon_ref_adata_image_f2 shape:", recon_ref_adata_image_f2.shape)
logger.info("Running inference tesk single-nuclei DONE!")
CUDA is available. GPU: NVIDIA A100-PCIE-40GB
***** Building loaders_inference sc image *****
[2025-06-29 22:58:23] INFO - Running inference tesk single-nuclei...
Finished loading all files
***** Finished building loaders_inference *****
device cuda:0
***** Begin perform_inference: ******
100%|████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 1.78it/s]
[2025-06-29 22:58:23] INFO - Running inference tesk single-nuclei DONE!
***** batch_size=adata.shape[0] *****
torch.Size([14647, 1, 1280])
1
Finished extractting image_between_spot data
--- 0.7193112373352051 seconds ---
recon_ref_adata_image_f2 shape: (14647, 1073)
[38]:
spatial_loc_sc = fst.get_allspot_coors(input_coord_all)
reconstructed_matrix_reshaped_tensor, _ = fst.reshape_latent_image(reconstructed_matrix_reshaped,
dataset_class='VisiumSC')
(_, _, all_spot_all_variable,
C2, adata_infer_sc) = fst.subspot_coord_expr_adata(reconstructed_matrix_reshaped_tensor,
spatial_loc_sc, gene_hv, patch_size=14, dataset_class="VisiumSC")
adata_smooth_sc = fst.impute_adata(adata_norml, adata_infer_sc, C2, gene_hv, dataset_class='VisiumSC', weight_exponent=2)
# adata_smooth_sc = fst.impute_adata(adata, adata_infer_sc, C2, gene_hv, dataset_class='VisiumSC', weight_exponent=2)
print(adata_smooth_sc)
Are there any duplicate rows? : False
pixel_step (half of patch_size): 7.0
Smoothing time: 6.3607 seconds
AnnData object with n_obs × n_vars = 14647 × 1073
obs: 'x', 'y'
var: 'gene_ids', 'feature_types', 'mt', 'n_cells_by_counts', 'mean_counts', 'log1p_mean_counts', 'pct_dropout_by_counts', 'total_counts', 'log1p_total_counts', 'n_cells'
uns: 'spatial'
obsm: 'spatial'
[23]:
# print(spatial_loc_sc)
# print(spatial_loc_sc[:,0].max(), spatial_loc_sc[:,1].max())
[[10000. 5250.49012146]
[10000. 5456.03825281]
[10000.14259663 3906.76939118]
...
[ 9999.69174785 4531.17247899]
[ 9999.74650563 3600.82598116]
[ 9999.84602469 6070.81414429]]
14690.0 16131.0
[24]:
# input_coord_all
[24]:
[[tensor([10000.0000, 10000.0000, 10000.1426, ..., 9999.6917, 9999.7465,
9999.8460], dtype=torch.float64),
tensor([5250.4901, 5456.0383, 3906.7694, ..., 4531.1725, 3600.8260,
6070.8141], dtype=torch.float64)]]
[25]:
# print(spatial_loc_sc)
# print(spatial_loc_sc[:,0].max(), spatial_loc_sc[:,1].max())
[[10000. 5250.49012146]
[10000. 5456.03825281]
[10000.14259663 3906.76939118]
...
[ 9999.69174785 4531.17247899]
[ 9999.74650563 3600.82598116]
[ 9999.84602469 6070.81414429]]
14690.0 16131.0
[27]:
# print(adata_impt_sc.obsm['spatial'])
# print(adata_impt_sc.obsm['spatial'][:,0].max(), adata_impt_sc.obsm['spatial'][:,1].max())
[45]:
adata_impt_sc, data_impt_sc = fst.weight_adata(adata_infer_sc, adata_smooth_sc, gene_hv, w=0.5)
print(adata_impt_sc)
## save adata
adata_impt_sc.write_h5ad(f'{demo_data}{save_data}adata_imput_all_sc.h5ad') # adata: 66029 × 596
AnnData object with n_obs × n_vars = 14647 × 1073
obs: 'x', 'y'
uns: 'spatial'
obsm: 'spatial'
[41]:
fstplt.gene_expr_allspots("IL18", spatial_loc_sc, recon_ref_adata_image_f2,
gene_hv, 'Inferred single-cell', marker='o', s=0.5,
figsize=(5, 4), cmap='Spectral_r', save_path=f'{demo_data}{fig_path}SPP1_expr_infer_allsc1.pdf')
IL18 gene expression dim: (14647, 1)
IL18 gene expression:
[[0.07839042]
[0.00958498]
[0.00892444]
...
[0.00970013]
[0.00887754]
[0.00888959]]
[42]:
# Save inferred expression
fstplt.gene_expr_allspots("SPP1", spatial_loc_sc, recon_ref_adata_image_f2, gene_hv,
'Inferred single-cell', s=0.6,
figsize=(5, 4), save_path=f'{demo_data}{fig_path}SPP1_expr_infer_allsc.pdf')
fstplt.gene_expr_allspots("IL18", spatial_loc_sc, recon_ref_adata_image_f2, gene_hv,
'Inferred single-cell', s=0.6,
figsize=(5, 4), save_path=f'{demo_data}{fig_path}IL18_expr_infer_allsc.pdf')
SPP1 gene expression dim: (14647, 1)
SPP1 gene expression:
[[0.0389441 ]
[0.03739255]
[0.00954874]
...
[0.01013606]
[0.00885651]
[0.00959768]]
IL18 gene expression dim: (14647, 1)
IL18 gene expression:
[[0.07839042]
[0.00958498]
[0.00892444]
...
[0.00970013]
[0.00887754]
[0.00888959]]
[47]:
## Save imputed expression
_, adata_impt_sc_reshape = fst.reshape_latent_image(data_impt_sc, dataset_class='VisiumSC')
fstplt.gene_expr_allspots("SPP1", spatial_loc_sc, adata_impt_sc_reshape, gene_hv,
'FineST single-cell', s=0.6,
figsize=(5, 4),
# save_path=f'{demo_data}{fig_path}SPP1_expr_imput_allsc.pdf'
)
fstplt.gene_expr_allspots("IL18", spatial_loc_sc, adata_impt_sc_reshape, gene_hv,
'FineST single-cell', s=0.6,
figsize=(5, 4),
# save_path=f'{demo_data}{fig_path}IL18_expr_imput_allsc.pdf'
)
SPP1 gene expression dim: (14647, 1)
SPP1 gene expression:
[[0.06748754]
[0.08309467]
[0.14246341]
...
[0.00506803]
[0.00442825]
[0.00479884]]
IL18 gene expression dim: (14647, 1)
IL18 gene expression:
[[0.06113877]
[0.0173814 ]
[0.00446222]
...
[0.00485007]
[0.00443877]
[0.00444479]]
Note The above all .h5ad files SaveData can be downloaded from Goole Drive.
[44]:
fst.release_gpu_tensors(globals())
Releasing the following GPU tensors:
input_image_exp_between_spot 71.52 MB shape=(14647, 1280)
reconstructed_matrix_reshaped 59.95 MB shape=(14647, 1073)
reconstructed_matrix_reshaped_tensor 59.95 MB shape=(14647, 1, 1073)
Released GPU memory occupied by the above variables.
map subspot to nuclei
[44]:
import os
path = '/mnt/lingyu/nfs_share2/Python/'
os.chdir(f'{path}FineST/FineST_demo/')
import numpy as np
import scanpy as sc
import pandas as pd
demo_data = 'HCC_data/'
save_data = 'SaveData/'
adata_impt_all_subspot = sc.read_h5ad(f'{demo_data}{save_data}adata_imput_all_subspot.h5ad')
adata_impt_all_sc = sc.read_h5ad(f'{demo_data}{save_data}adata_imput_all_sc.h5ad')
print(adata_impt_all_subspot)
print(adata_impt_all_sc)
AnnData object with n_obs × n_vars = 830656 × 1073
obs: 'x', 'y'
uns: 'spatial'
obsm: 'spatial'
AnnData object with n_obs × n_vars = 66029 × 1073
obs: 'x', 'y'
uns: 'spatial'
obsm: 'spatial'
[45]:
import scanpy as sc
from scipy.spatial import cKDTree
def map_subspot_to_nuclei(adata_subspot, adata_nuclei,
inherit_uns_from="subspot", spatial_key="spatial"):
"""
Map each subspot to its nearest nuclei and keep unique mappings.
Parameters
----------
adata_subspot : AnnData
AnnData object containing subspot data (to be mapped).
adata_nuclei : AnnData
AnnData object containing nuclei data (reference).
spatial_key : str
Key in obsm for spatial coordinates.
Returns
-------
adata_map : AnnData
AnnData object with only the uniquely mapped subspots,
including spatial coordinates.
"""
## Build tree for nuclei coordinates, find nearest nuclei index for each subspot
tree = cKDTree(adata_nuclei.obsm[spatial_key])
_, closest_indices = tree.query(adata_subspot.obsm[spatial_key], k=1)
## Get the corresponding coordinates of mapped nuclei, Keep only unique mappings
mapped_coords = adata_nuclei.obsm[spatial_key][closest_indices]
_, unique_idx = np.unique(mapped_coords, axis=0, return_index=True)
## Construct new AnnData object with unique mappings
adata_map = sc.AnnData(
adata_subspot.X[unique_idx],
var=adata_subspot.var.copy(),
obs=adata_subspot.obs.iloc[unique_idx].copy()
)
adata_map.obsm[spatial_key] = mapped_coords[unique_idx]
adata_map.var_names = adata_subspot.var_names
# Copy .uns["spatial"]
if inherit_uns_from == "subspot":
adata_map.uns["spatial"] = adata_subspot.uns["spatial"].copy()
elif inherit_uns_from == "nuclei":
adata_map.uns["spatial"] = adata_nuclei.uns["spatial"].copy()
else:
raise ValueError("inherit_uns_from must be 'subspot' or 'nuclei'.")
return adata_map
[46]:
adata_map = map_subspot_to_nuclei(adata_impt_all_subspot, adata_impt_all_sc)
adata_map.write_h5ad(f'{demo_data}{save_data}adata_imput_all_sc_mapped.h5ad')
print(adata_map)
AnnData object with n_obs × n_vars = 65675 × 1073
obs: 'x', 'y'
uns: 'spatial'
obsm: 'spatial'
[47]:
fstplt.gene_expr_allspots("IL18", adata_map.obsm['spatial'], adata_map.to_df(), gene_hv,
'FineST single-cell map', s=0.6, figsize=(5, 4),
# save_path=None)
save_path=f'{demo_data}{fig_path}/IL18_expr_imput_allsc_mapped.pdf')
fstplt.gene_expr_allspots("SPP1", adata_map.obsm['spatial'], adata_map.to_df(), gene_hv,
'FineST single-cell map', s=0.6, figsize=(5, 4),
# save_path=None)
save_path=f'{demo_data}{fig_path}/SPP1_expr_imput_allsc_mapped.pdf')
IL18 gene expression dim: (65675, 1)
IL18 gene expression:
[[0.04592751]
[0.01063765]
[0.00969013]
...
[0.009544 ]
[0.00829374]
[0.00671693]]
SPP1 gene expression dim: (65675, 1)
SPP1 gene expression:
[[0.16072396]
[0.10252475]
[0.02213156]
...
[0.01143425]
[0.01013806]
[0.00764641]]
[50]:
fstplt.gene_expr_allspots("SPP1", adata_map.obsm['spatial'], adata_map.to_df(),
gene_hv, 'FineST single-cell map', marker='o', s=0.5,
figsize=(5, 4), cmap='Spectral_r', save_path=f'{demo_data}{fig_path}/SPP1_expr_imput_allsc_mapped.pdf')
fstplt.gene_expr_allspots("IL18", adata_map.obsm['spatial'], adata_map.to_df(),
gene_hv, 'FineST single-cell map', marker='o', s=0.5,
figsize=(5, 4), cmap='Spectral_r', save_path=f'{demo_data}{fig_path}/IL18_expr_imput_allsc_mapped.pdf')
SPP1 gene expression dim: (65675, 1)
SPP1 gene expression:
[[0.16072396]
[0.10252475]
[0.02213156]
...
[0.01143425]
[0.01013806]
[0.00764641]]
IL18 gene expression dim: (65675, 1)
IL18 gene expression:
[[0.04592751]
[0.01063765]
[0.00969013]
...
[0.009544 ]
[0.00829374]
[0.00671693]]
[52]:
fstplt.gene_expr_allspots("CCR1", adata_map.obsm['spatial'], adata_map.to_df(),
gene_hv, 'FineST single-cell map', marker='o', s=0.5,
figsize=(5, 4), cmap='Spectral_r', save_path=f'{demo_data}{fig_path}/CCR1_expr_imput_allsc_mapped.pdf')
fstplt.gene_expr_allspots("THBS1", adata_map.obsm['spatial'], adata_map.to_df(),
gene_hv, 'FineST single-cell map', marker='o', s=0.5,
figsize=(5, 4), cmap='Spectral_r', save_path=f'{demo_data}{fig_path}/THBS1_expr_imput_allsc_mapped.pdf')
CCR1 gene expression dim: (65675, 1)
CCR1 gene expression:
[[0.03009469]
[0.00929491]
[0.00862425]
...
[0.01149619]
[0.00875341]
[0.00685632]]
THBS1 gene expression dim: (65675, 1)
THBS1 gene expression:
[[0.05908929]
[0.0349021 ]
[0.0125149 ]
...
[0.04234706]
[0.01004844]
[0.02690993]]
