Starfysh tutorial on real dataset
Azizi Lab
Siyu He, Yinuo Jin
02-02-2023
This is a tutorial on an example real Spatial Transcriptomics (ST) data (CID44971_TNBC) from Wu et al., 2021.
Overview
Starfysh performs cell-type deconvolution followed by various downstream analyses to discover spatial interactions in tumor microenvironment. Specifically, Starfysh looks for anchor spots (presumably with the highest compositions of one given cell type) informed by user-provided gene signatures (see example) as priors to guide the deconvolution inference, which further enables downstream analyses such as sample integration, spatial hub characterization, cell-cell interactions, etc. This tutorial focuses on the deconvolution task. Overall, Starfysh provides the following options:
Base feature:
Auxiliary Variational AutoEncoder (AVAE): Spot-level deconvolution with expected cell types and corresponding annotated signature gene sets (default)
Optional:
Archetypal Analysis (AA):
If gene signatures are not provided
Unsupervised cell type annotation
If gene signatures are provided but require refinement:
Novel cell type / cell state discovery (complementary to known cell types from the signatures)
Refine known marker genes by appending archetype-specific differentially expressed genes, and update anchor spots accordingly
Product-of-Experts (PoE) integration: Example updated soon!
Multi-modal integrative predictions with expression & histology image by leverging additional side information (e.g. cell density) from H&E image.
[1]:
import sys
IN_COLAB = "google.colab" in sys.modules
if IN_COLAB:
!pip3 install scanpy
!pip install git+https://github.com/azizilab/starfysh.git
!pip install scikit-image --upgrade
from google.colab import drive
drive.mount('/content/drive')
import sys
# Please specify the colab notebook directory
sys.path.append('/content/drive/MyDrive/Starfysh')
[2]:
import os
import numpy as np
import pandas as pd
import torch
[3]:
import matplotlib.pyplot as plt
import matplotlib.font_manager
from matplotlib import rcParams
import seaborn as sns
sns.set_style('white')
font_list = []
fpaths = matplotlib.font_manager.findSystemFonts()
for i in fpaths:
try:
f = matplotlib.font_manager.get_font(i)
font_list.append(f.family_name)
except RuntimeError:
pass
font_list = set(font_list)
plot_font = 'Helvetica' if 'Helvetica' in font_list else 'FreeSans'
rcParams['font.family'] = plot_font
rcParams.update({'font.size': 10})
rcParams.update({'figure.dpi': 300})
rcParams.update({'figure.figsize': (3,3)})
rcParams.update({'savefig.dpi': 500})
import warnings
warnings.filterwarnings('ignore')
Load starfysh
[4]:
from starfysh import (AA, utils, plot_utils, post_analysis)
from starfysh import starfysh as sf_model
(1). load data and marker genes
File Input: - Spatial transcriptomics - Count matrix: adata - (Optional): Paired histology & spot coordinates: img, map_info
Annotated signatures (marker genes) for potential cell types:
gene_sig
Starfysh is built upon scanpy and Anndata. The common ST/Visium data sample folder consists a expression count file (usually filtered_featyur_bc_matrix.h5), and a subdirectory with corresponding H&E image and spatial information, as provided by Visium platform.
For example, our example real ST data has the following structure:
├── ../data
bc_signatures_version_1013.csv
├── P1A_ER:
\__ filtered_feature_bc_mactrix.h5
├── spatial:
\__ aligned_fiducials.jpg
detected_tissue_image.jpg
scalefactors_json.json
tissue_hires_image.png
tissue_lowres_image.png
tissue_positions_list.csv
For data that doesn’t follow the common visium data structure (e.g. missing filtered_featyur_bc_matrix.h5 or the given .h5ad count matrix file lacks spatial metadata), please construct the data as Anndata synthesizing information as the example simulated data shows:
[Note]: If you’re running this tutorial locally, please download the sample data and signature gene sets, and save it in the relative path ../data (otherwise please modify the data_path defined in the cell below):
[5]:
# Specify data paths
data_path = '../data/'
sample_id = 'CID44971_TNBC'
sig_name = 'bc_signatures_version_1013.csv'
[6]:
# Load expression counts and signature gene sets
adata, adata_normed = utils.load_adata(data_folder=data_path,
sample_id=sample_id, # sample id
n_genes=2000 # number of highly variable genes to keep
)
gene_sig = pd.read_csv(os.path.join(data_path, sig_name))
gene_sig = utils.filter_gene_sig(gene_sig, adata.to_df())
gene_sig.head()
[2023-02-26 14:34:03] Preprocessing1: delete the mt and rp
[2023-02-26 14:34:03] Preprocessing2: Normalize
[2023-02-26 14:34:03] Preprocessing3: Logarithm
[2023-02-26 14:34:03] Preprocessing4: Find the variable genes
[6]:
| Basal | LumA | LumB | MBC | Normal epithelial | Tcm | Tem | Tfh | Treg | Activated CD8 | ... | Plasmablasts | MDSC | Monocytes | cDC | pDC | CAFs MSC iCAF-like | CAFs myCAF-like | PVL differentiated | PVL immature | Endothelial | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | EMP1 | SH3BGRL | UGCG | COL11A2 | KRT14 | CCR7 | IL7R | CXCL13 | TNFRSF4 | CD69 | ... | IGKV3-15 | ITGAM | LYZ | CD80 | IL3RA | APOD | COL1A1 | ACTA2 | CCL19 | ACKR1 |
| 1 | TAGLN | HSPB1 | ARMT1 | SDC1 | KRT17 | LTB | ANXA1 | NMB | LTB | CCR7 | ... | IGHG1 | CD33 | IL1B | CD86 | LILRA4 | DCN | COL1A2 | TAGLN | RGS5 | FABP4 |
| 2 | TTYH1 | PHGR1 | ISOC1 | FBN2 | LTF | IL7R | CXCR4 | NR3C1 | IL32 | CD27 | ... | IGKV1-5 | ARG1 | G0S2 | CCR7 | CD123 | PTGDS | COL3A1 | MYL9 | IGFBP7 | PLVAP |
| 3 | RTN4 | SOX9 | GDF15 | MMP1 | KRT15 | SARAF | KLRB1 | DUSP4 | BATF | BTLA | ... | IGKV3-20 | NOS2 | TYROBP | CD1A | TCF4 | CFD | LUM | TPM2 | NDUFA4L2 | RAMP2 |
| 4 | TK1 | CEBPD | ZFP36 | FABP5 | PTN | SELL | TNFAIP3 | TNFRSF18 | FOXP3 | CD40LG | ... | IGKV3-11 | CD68 | FCN1 | CD1C | IRF7 | LUM | SFRP2 | NDUFA4L2 | CCL2 | VWF |
5 rows × 29 columns
If there’s no input signature gene sets, Starfysh defines “archetypal marker genes” as signatures. Please refer to the following code snippet and see details in section (3).
aa_model = AA.ArchetypalAnalysis(adata_orig=adata_normed)
archetype, arche_dict, major_idx, evs = aa_model.compute_archetypes(r=40)
gene_sig = aa_model.find_markers(n_markers=30, display=False)
gene_sig = utils.filter_gene_sig(gene_sig, adata.to_df())
gene_sig.head()
[7]:
adata
[7]:
AnnData object with n_obs × n_vars = 1162 × 18858
obs: 'orig.ident', 'nCount_RNA', 'nFeature_RNA', 'sample', '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: 'features', 'highly_variable'
[8]:
# Load spatial information
img_metadata = utils.preprocess_img(data_path,
sample_id,
adata_index=adata.obs.index,
hchannel=False
)
img, map_info, scalefactor = img_metadata['img'], img_metadata['map_info'], img_metadata['scalefactor']
umap_df = utils.get_umap(adata, display=True)
... storing 'sample' as categorical
[9]:
plt.figure(figsize=(6, 6), dpi=200)
plt.imshow(img)
plt.show()
[10]:
map_info.head()
[10]:
| array_row | array_col | imagerow | imagecol | sample | |
|---|---|---|---|---|---|
| AACATTGGTCAGCCGT-1 | 3 | 17 | 1361 | 2511 | CID44971_TNBC |
| CATCGAATGGATCTCT-1 | 3 | 19 | 1361 | 2620 | CID44971_TNBC |
| CGGGTTGTAGCTTTGG-1 | 3 | 21 | 1361 | 2730 | CID44971_TNBC |
| CCTAAGTGTCTAACCG-1 | 2 | 22 | 1266 | 2784 | CID44971_TNBC |
| TCTGTGACTGACCGTT-1 | 3 | 23 | 1361 | 2839 | CID44971_TNBC |
(2). Preprocessing (finding anchor spots)
Identify anchor spot locations.
Instantiate parameters for Starfysh model training: - Raw & normalized counts after taking highly variable genes - filtered signature genes - library size & spatial smoothed library size (log-transformed) - Anchor spot indices (anchors_df) for each cell type & their signature means (sig_means)
[11]:
# Parameters for training
visium_args = utils.VisiumArguments(adata,
adata_normed,
gene_sig,
img_metadata,
n_anchors=60,
window_size=3,
sample_id=sample_id
)
adata, adata_normed = visium_args.get_adata()
anchors_df = visium_args.get_anchors()
[2023-02-26 14:34:12] Subsetting highly variable & signature genes ...
[2023-02-26 14:34:12] Smoothing library size by taking averaging with neighbor spots...
[2023-02-26 14:34:12] Retrieving & normalizing signature gene expressions...
[2023-02-26 14:34:12] Identifying anchor spots (highly expression of specific cell-type signatures)...
Check the dimensions of preprocessed data after selecting the union of highly variable genes & user defined signatures:
[12]:
adata.shape
[12]:
(1162, 2608)
Visualize the spatial data
[13]:
plot_utils.plot_spatial_feature(adata,
map_info,
visium_args.log_lib,
label='log library size'
)
[15]:
plot_utils.plot_spatial_feature(adata,
map_info,
visium_args.win_loglib,
label='windowed log library size',
)
plot raw gene expression:
[14]:
plot_utils.plot_spatial_gene(adata,
map_info,
gene_name='IL7R')
Visualize anchor spots
[15]:
plot_utils.plot_anchor_spots(umap_df,
visium_args.pure_spots,
visium_args.sig_mean,
bbox_x=2
)
(3). Optional: Archetypal Analysis
Overview: If users don’t provide annotated gene signature sets with cell types, Starfysh identifies candidates for cell types via archetypal analysis (AA). The underlying assumption is that the geometric “extremes” are identified as the purest cell types, whereas all other spots are mixture of the “archetypes”. If the users provide the gene signature sets, they can still optionally apply AA to refine marker genes and update anchor spots for known cell types. In addition, AA can identify & assign potential novel cell types / states. Here are the features provided by the optional archetypal analysis: - Finding archetypal spots & assign 1-1 mapping to their closest anchor spot neighbors - Finding archetypal marker genes & append them to marker genes of annotated cell types - Assigning novel cell type / cell states as the most distant archetypes
Overall, Starfysh provides the archetypal analysis as a complementary toolkit for automatic cell-type annotation & signature gene completion:
If signature genes aren’t provided: Archetypal analysis defines the geometric extrema of the data as major cell types with corresponding marker genes.
If complete signature genes are known: Users can skip this section and use only the signature priors
If signature genes are incomplete or need refinement: Archetypal analysis can be applied to
Refine signatures of certain cell types
Find novel cell types / states that haven’t been provided from the input signature
If signature genes aren’t provided
Note: - Intrinsic Dimension (ID) estimator is implemented to estimate the lower-bound for the number of archetypes \(k\), followed by elbow method with iterations to identify the optimal \(k\). By default, a conditional number is set as 30; if you believe there are potentially more / fewer cell types, please increase / decrease cn accordingly.
Major cell types & corresponding markers are represented by the inferred archetypes:
aa_model = AA.ArchetypalAnalysis(adata_orig=adata_normed)
archetype, arche_dict, major_idx, evs = aa_model.compute_archetypes(r=40)
# (1). Find archetypal spots & archetypal clusters
arche_df = aa_model.find_archetypal_spots(major=True)
# (2). Define "signature genes" as marker genes associated with each archetypal cluster
gene_sig = aa_model.find_markers(n_markers=30, display=False)
gene_sig.head()
If complete signature genes are known
Users can skip th section & run Starfysh
If signature genes are incomplete or require refinement
In this tutorial, we’ll show an example of applying best-aligned archetypes to existing ``anchors`` of given cell type(s) to append signature genes.
[16]:
aa_model = AA.ArchetypalAnalysis(adata_orig=adata_normed)
archetype, arche_dict, major_idx, evs = aa_model.compute_archetypes(r=40)
# (1). Find archetypal spots & archetypal clusters
arche_df = aa_model.find_archetypal_spots(major=True)
# (2). Find marker genes associated with each archetypal cluster
markers_df = aa_model.find_markers(n_markers=30, display=False)
# (3). Map archetypes to closest anchors (1-1 per cell type)
map_df, map_dict = aa_model.assign_archetypes(anchors_df)
# (4). Optional: Find the most distant archetypes that are not assigned to any annotated cell types
distant_arches = aa_model.find_distant_archetypes(anchors_df, n=3)
[2023-02-26 14:34:21] Computing intrinsic dimension to estimate k...
30 components are retained using conditional_number=30.00
[2023-02-26 14:34:22] Estimating lower bound of # archetype as 12...
[2023-02-26 14:34:36] Calculating UMAPs for counts + Archetypes...
[2023-02-26 14:34:41] Calculating UMAPs for counts + Archetypes...
[2023-02-26 14:34:43] 0.6859 variance explained by raw archetypes.
Merging raw archetypes within 40 NNs to get major archetypes
[2023-02-26 14:34:43] Finding 20 nearest neighbors for each archetype...
[2023-02-26 14:34:43] Finding 30 top marker genes for each archetype...
[17]:
plot_utils.plot_evs(evs, kmin=aa_model.kmin)
Visualize archetypes
[18]:
aa_model.plot_archetypes(do_3d=False, major=True, disp_cluster=False)
[2023-02-26 14:34:45] No handles with labels found to put in legend.
[18]:
(<Figure size 1800x1200 with 1 Axes>, <AxesSubplot:>)
Visualize overlapping ratio between archetypal & anchor spots:
[19]:
aa_model.plot_mapping(map_df)
[19]:
<seaborn.matrix.ClusterGrid at 0x7f9d632667c0>
Application: appending marker genes Append archetypal marker genes with the best-aligned anchors:
[20]:
visium_args = utils.refine_anchors(
visium_args,
aa_model,
thld=0.7, # alignment threshold
n_genes=5,
n_iters=1
)
# Get updated adata & signatures
adata, adata_normed = visium_args.get_adata()
gene_sig = visium_args.gene_sig
cell_types = gene_sig.columns
[2023-02-26 14:34:46] Finding 50 top marker genes for each archetype...
[2023-02-26 14:34:48] Recalculating anchor spots (highly expression of specific cell-type signatures)...
Refining round 1...
appending 5 genes in arch_0 to Treg...
appending 5 genes in arch_10 to CAFs myCAF-like...
Run starfysh
We perform n_repeat random restarts and select the best model with lowest loss for parameter c (inferred cell-type proportions):
(1). Model parameters
[21]:
n_repeats = 3
epochs = 200
patience = 50
device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
(2). Model training
Users can choose to run the one of the following Starfysh model without/with histology integration: - Base Starfysh deconvolution with ST expression matrix - Starfysh deconvolution with ST expression matrix & histology integration
Base Starfysh deconvolution with ST expression matrix
[23]:
# Run models
model, loss = utils.run_starfysh(visium_args,
n_repeats=n_repeats,
epochs=epochs,
patience=patience,
device=device
)
[2023-02-23 15:36:15] Running Starfysh with 3 restarts, choose the model with best parameters...
[2023-02-23 15:36:15] === Restart Starfysh 1 ===
[2023-02-23 15:36:17] Initializing model parameters...
[2023-02-23 15:36:26] Epoch[10/200], train_loss: 2757.1400, train_reconst: 1549.9088, train_u: 4.7961,train_z: 18.4716,train_c: 1154.8984,train_n: 33.8611
[2023-02-23 15:36:35] Epoch[20/200], train_loss: 2473.3621, train_reconst: 1495.9336, train_u: 5.5051,train_z: 17.9046,train_c: 928.7905,train_n: 30.7333
[2023-02-23 15:36:44] Epoch[30/200], train_loss: 2330.8574, train_reconst: 1471.7674, train_u: 6.2065,train_z: 18.1975,train_c: 810.6913,train_n: 30.2012
[2023-02-23 15:36:53] Epoch[40/200], train_loss: 2247.9081, train_reconst: 1440.6176, train_u: 6.9305,train_z: 18.5305,train_c: 759.2430,train_n: 29.5171
[2023-02-23 15:37:02] Epoch[50/200], train_loss: 2324.9906, train_reconst: 1419.3684, train_u: 7.8082,train_z: 18.8234,train_c: 857.7535,train_n: 29.0453
[2023-02-23 15:37:11] Epoch[60/200], train_loss: 2117.2935, train_reconst: 1400.2355, train_u: 8.9160,train_z: 19.2183,train_c: 668.7168,train_n: 29.1230
[2023-02-23 15:37:20] Epoch[70/200], train_loss: 2056.6418, train_reconst: 1385.9130, train_u: 10.2061,train_z: 19.9288,train_c: 621.5048,train_n: 29.2951
[2023-02-23 15:37:28] Epoch[80/200], train_loss: 2078.7415, train_reconst: 1375.3640, train_u: 11.6821,train_z: 20.2772,train_c: 653.6566,train_n: 29.4436
[2023-02-23 15:37:37] Epoch[90/200], train_loss: 1981.2158, train_reconst: 1358.2318, train_u: 13.3211,train_z: 20.6815,train_c: 573.3214,train_n: 28.9811
[2023-02-23 15:37:46] Epoch[100/200], train_loss: 1955.6103, train_reconst: 1342.9327, train_u: 15.1606,train_z: 20.9966,train_c: 562.3097,train_n: 29.3714
[2023-02-23 15:37:55] Epoch[110/200], train_loss: 1948.3726, train_reconst: 1335.1401, train_u: 17.1632,train_z: 21.1413,train_c: 562.8060,train_n: 29.2851
[2023-02-23 15:38:04] Epoch[120/200], train_loss: 1922.7624, train_reconst: 1331.7687, train_u: 19.3093,train_z: 21.8316,train_c: 540.2918,train_n: 28.8703
[2023-02-23 15:38:13] Epoch[130/200], train_loss: 1908.0791, train_reconst: 1320.3909, train_u: 21.6107,train_z: 22.1016,train_c: 536.2219,train_n: 29.3647
[2023-02-23 15:38:22] Epoch[140/200], train_loss: 1917.7096, train_reconst: 1316.0702, train_u: 24.0703,train_z: 22.2168,train_c: 550.2750,train_n: 29.1476
[2023-02-23 15:38:31] Epoch[150/200], train_loss: 1886.2947, train_reconst: 1303.7384, train_u: 26.6950,train_z: 22.2046,train_c: 531.8323,train_n: 28.5195
[2023-02-23 15:38:40] Epoch[160/200], train_loss: 1873.9099, train_reconst: 1294.2157, train_u: 29.4463,train_z: 22.8038,train_c: 528.0124,train_n: 28.8780
[2023-02-23 15:38:49] Epoch[170/200], train_loss: 1852.1266, train_reconst: 1296.0656, train_u: 32.2922,train_z: 23.0411,train_c: 504.3316,train_n: 28.6883
[2023-02-23 15:38:58] Epoch[180/200], train_loss: 1844.4900, train_reconst: 1285.0703, train_u: 35.2073,train_z: 23.1973,train_c: 507.3918,train_n: 28.8306
[2023-02-23 15:39:07] Epoch[190/200], train_loss: 1825.8041, train_reconst: 1277.0771, train_u: 38.2100,train_z: 23.5896,train_c: 496.7026,train_n: 28.4348
[2023-02-23 15:39:16] Epoch[200/200], train_loss: 1804.8194, train_reconst: 1267.3563, train_u: 41.2891,train_z: 23.6020,train_c: 485.0773,train_n: 28.7838
[2023-02-23 15:39:16] === Finished training ===
[2023-02-23 15:39:16] === Restart Starfysh 2 ===
[2023-02-23 15:39:16] Initializing model parameters...
[2023-02-23 15:39:25] Epoch[10/200], train_loss: 2771.4032, train_reconst: 1546.1694, train_u: 5.3345,train_z: 17.7110,train_c: 1174.0513,train_n: 33.4715
[2023-02-23 15:39:34] Epoch[20/200], train_loss: 2568.9384, train_reconst: 1514.5348, train_u: 5.7845,train_z: 17.2977,train_c: 1006.7286,train_n: 30.3773
[2023-02-23 15:39:42] Epoch[30/200], train_loss: 2353.1475, train_reconst: 1475.0976, train_u: 6.3714,train_z: 17.4276,train_c: 830.6116,train_n: 30.0107
[2023-02-23 15:39:51] Epoch[40/200], train_loss: 2266.3132, train_reconst: 1450.3741, train_u: 7.1521,train_z: 17.7185,train_c: 768.6599,train_n: 29.5607
[2023-02-23 15:40:00] Epoch[50/200], train_loss: 2215.4436, train_reconst: 1426.8405, train_u: 8.1002,train_z: 18.4010,train_c: 740.8346,train_n: 29.3675
[2023-02-23 15:40:09] Epoch[60/200], train_loss: 2145.0040, train_reconst: 1411.8011, train_u: 9.1891,train_z: 18.5509,train_c: 685.1804,train_n: 29.4716
[2023-02-23 15:40:18] Epoch[70/200], train_loss: 2071.2608, train_reconst: 1389.6598, train_u: 10.4682,train_z: 19.2782,train_c: 632.4492,train_n: 29.8736
[2023-02-23 15:40:27] Epoch[80/200], train_loss: 2066.3306, train_reconst: 1373.4050, train_u: 11.9056,train_z: 19.6237,train_c: 644.3054,train_n: 28.9965
[2023-02-23 15:40:36] Epoch[90/200], train_loss: 2027.0461, train_reconst: 1364.6049, train_u: 13.4306,train_z: 20.1261,train_c: 613.1041,train_n: 29.2110
[2023-02-23 15:40:45] Epoch[100/200], train_loss: 2039.8596, train_reconst: 1350.0161, train_u: 15.1033,train_z: 20.5216,train_c: 640.4962,train_n: 28.8257
[2023-02-23 15:40:53] Epoch[110/200], train_loss: 2021.6711, train_reconst: 1340.9502, train_u: 16.8720,train_z: 21.0727,train_c: 630.4961,train_n: 29.1521
[2023-02-23 15:41:02] Epoch[120/200], train_loss: 1955.6004, train_reconst: 1327.4103, train_u: 18.7711,train_z: 21.1771,train_c: 578.3904,train_n: 28.6227
[2023-02-23 15:41:11] Epoch[130/200], train_loss: 2008.8703, train_reconst: 1321.9498, train_u: 20.7826,train_z: 21.6331,train_c: 636.7317,train_n: 28.5556
[2023-02-23 15:41:20] Epoch[140/200], train_loss: 1922.1904, train_reconst: 1313.9041, train_u: 22.9024,train_z: 22.0497,train_c: 557.6942,train_n: 28.5424
[2023-02-23 15:41:34] Epoch[150/200], train_loss: 1942.7256, train_reconst: 1308.1408, train_u: 25.0997,train_z: 22.1656,train_c: 583.6355,train_n: 28.7836
[2023-02-23 15:41:43] Epoch[160/200], train_loss: 1885.2464, train_reconst: 1298.1472, train_u: 27.3562,train_z: 22.4174,train_c: 536.1287,train_n: 28.5531
[2023-02-23 15:41:52] Epoch[170/200], train_loss: 1837.8245, train_reconst: 1290.5161, train_u: 29.7023,train_z: 22.7446,train_c: 496.2513,train_n: 28.3124
[2023-02-23 15:42:01] Epoch[180/200], train_loss: 1845.8878, train_reconst: 1280.7779, train_u: 32.2259,train_z: 22.8609,train_c: 513.8417,train_n: 28.4073
[2023-02-23 15:42:10] Epoch[190/200], train_loss: 1820.2334, train_reconst: 1281.3175, train_u: 34.7438,train_z: 23.0705,train_c: 487.6215,train_n: 28.2238
[2023-02-23 15:42:19] Epoch[200/200], train_loss: 1824.2379, train_reconst: 1276.1209, train_u: 37.3222,train_z: 23.5247,train_c: 496.0574,train_n: 28.5349
[2023-02-23 15:42:19] === Finished training ===
[2023-02-23 15:42:19] === Restart Starfysh 3 ===
[2023-02-23 15:42:19] Initializing model parameters...
[2023-02-23 15:42:28] Epoch[10/200], train_loss: 2828.1074, train_reconst: 1544.8992, train_u: 5.7779,train_z: 18.6729,train_c: 1229.5659,train_n: 34.9694
[2023-02-23 15:42:37] Epoch[20/200], train_loss: 2530.7295, train_reconst: 1504.1783, train_u: 6.4677,train_z: 17.9246,train_c: 977.1099,train_n: 31.5168
[2023-02-23 15:42:45] Epoch[30/200], train_loss: 2347.7388, train_reconst: 1470.5720, train_u: 7.2472,train_z: 17.8295,train_c: 829.4506,train_n: 29.8868
[2023-02-23 15:42:55] Epoch[40/200], train_loss: 2259.7625, train_reconst: 1445.8083, train_u: 8.1376,train_z: 18.2212,train_c: 766.5617,train_n: 29.1713
[2023-02-23 15:43:04] Epoch[50/200], train_loss: 2184.3394, train_reconst: 1425.3061, train_u: 9.1511,train_z: 18.4206,train_c: 711.4393,train_n: 29.1735
[2023-02-23 15:43:13] Epoch[60/200], train_loss: 2139.7678, train_reconst: 1408.0426, train_u: 10.3143,train_z: 18.8218,train_c: 683.8678,train_n: 29.0356
[2023-02-23 15:43:22] Epoch[70/200], train_loss: 2097.4261, train_reconst: 1391.2586, train_u: 11.5474,train_z: 19.0560,train_c: 658.0438,train_n: 29.0678
[2023-02-23 15:43:31] Epoch[80/200], train_loss: 2077.8687, train_reconst: 1377.3639, train_u: 12.9338,train_z: 19.6265,train_c: 652.1278,train_n: 28.7505
[2023-02-23 15:43:40] Epoch[90/200], train_loss: 2026.1179, train_reconst: 1366.2854, train_u: 14.4203,train_z: 19.8512,train_c: 611.1232,train_n: 28.8581
[2023-02-23 15:43:49] Epoch[100/200], train_loss: 2015.5042, train_reconst: 1350.7212, train_u: 16.0561,train_z: 20.3470,train_c: 615.3903,train_n: 29.0457
[2023-02-23 15:43:58] Epoch[110/200], train_loss: 1980.4622, train_reconst: 1345.4858, train_u: 17.7693,train_z: 20.5726,train_c: 585.3474,train_n: 29.0565
[2023-02-23 15:44:07] Epoch[120/200], train_loss: 1956.4255, train_reconst: 1336.0091, train_u: 19.6240,train_z: 20.9087,train_c: 570.9680,train_n: 28.5396
[2023-02-23 15:44:16] Epoch[130/200], train_loss: 1918.6394, train_reconst: 1327.1537, train_u: 21.5526,train_z: 21.0591,train_c: 541.8596,train_n: 28.5669
[2023-02-23 15:44:25] Epoch[140/200], train_loss: 1936.7608, train_reconst: 1317.8962, train_u: 23.5975,train_z: 21.6194,train_c: 568.5682,train_n: 28.6770
[2023-02-23 15:44:34] Epoch[150/200], train_loss: 1930.8683, train_reconst: 1310.6734, train_u: 25.7242,train_z: 21.4405,train_c: 570.2975,train_n: 28.4570
[2023-02-23 15:44:43] Epoch[160/200], train_loss: 1915.0496, train_reconst: 1305.1493, train_u: 27.9648,train_z: 21.8778,train_c: 559.4349,train_n: 28.5877
[2023-02-23 15:44:52] Epoch[170/200], train_loss: 1948.3021, train_reconst: 1303.8845, train_u: 30.2412,train_z: 22.4456,train_c: 593.3461,train_n: 28.6260
[2023-02-23 15:45:01] Epoch[180/200], train_loss: 1842.3880, train_reconst: 1289.2384, train_u: 32.6075,train_z: 22.4516,train_c: 502.2100,train_n: 28.4880
[2023-02-23 15:45:10] Epoch[190/200], train_loss: 1868.3176, train_reconst: 1287.8066, train_u: 35.0768,train_z: 22.6185,train_c: 529.0843,train_n: 28.8082
[2023-02-23 15:45:19] Epoch[200/200], train_loss: 1832.6485, train_reconst: 1283.3994, train_u: 37.6051,train_z: 22.9518,train_c: 497.8525,train_n: 28.4448
[2023-02-23 15:45:19] === Finished training ===
Starfysh deconvolution with histology integration
We perform integrative inference of expression and the paired histology image using the Products of Expert (PoE) model (Lee and van der Schaar, 2021), which projects different data modalities into the joint latent space to aid deconvolution with actual spatial information.
[22]:
model, loss_poe = utils.run_starfysh(visium_args,
n_repeats=n_repeats,
epochs=epochs,
patience=patience,
poe=True,
device=device)
[2023-02-26 14:34:49] Running Starfysh with 3 restarts, choose the model with best parameters...
[2023-02-26 14:34:49] === Restart Starfysh 1 ===
[2023-02-26 14:34:51] Initializing model parameters...
[2023-02-26 14:35:07] Epoch[10/200], train_loss: 5106204.3294, train_reconst: 3491.4993, train_u: 0.0000,train_z: 17.4829,train_c: 1205.3407,train_n: 32.3482
[2023-02-26 14:35:21] Epoch[20/200], train_loss: 1734606.5557, train_reconst: 3285.0179, train_u: 0.0000,train_z: 19.6014,train_c: 1021.8931,train_n: 30.5615
[2023-02-26 14:35:34] Epoch[30/200], train_loss: 661569.3970, train_reconst: 3175.6771, train_u: 0.0000,train_z: 17.6588,train_c: 897.7205,train_n: 29.8804
[2023-02-26 14:35:47] Epoch[40/200], train_loss: 371853.4155, train_reconst: 3200.9416, train_u: 0.0000,train_z: 24.1707,train_c: 871.4137,train_n: 29.4798
[2023-02-26 14:36:00] Epoch[50/200], train_loss: 328906.4683, train_reconst: 3120.3758, train_u: 0.0000,train_z: 19.6875,train_c: 709.0122,train_n: 28.9512
[2023-02-26 14:36:14] Epoch[60/200], train_loss: 303952.7880, train_reconst: 3242.5707, train_u: 0.0000,train_z: 26.6282,train_c: 731.7766,train_n: 29.4965
[2023-02-26 14:36:27] Epoch[70/200], train_loss: 298386.0866, train_reconst: 3083.1203, train_u: 0.0000,train_z: 22.3257,train_c: 731.4536,train_n: 28.6894
[2023-02-26 14:36:40] Epoch[80/200], train_loss: 289565.0165, train_reconst: 3192.4708, train_u: 0.0000,train_z: 24.7036,train_c: 916.7221,train_n: 28.7957
[2023-02-26 14:36:54] Epoch[90/200], train_loss: 275127.5141, train_reconst: 2986.4965, train_u: 0.0000,train_z: 20.2450,train_c: 638.0327,train_n: 28.4446
[2023-02-26 14:37:07] Epoch[100/200], train_loss: 289832.7872, train_reconst: 3259.8671, train_u: 0.0000,train_z: 23.5122,train_c: 941.4290,train_n: 28.6065
[2023-02-26 14:37:20] Epoch[110/200], train_loss: 267550.6119, train_reconst: 3008.3280, train_u: 0.0000,train_z: 22.0322,train_c: 662.3952,train_n: 28.6267
[2023-02-26 14:37:34] Epoch[120/200], train_loss: 277357.8822, train_reconst: 2930.8920, train_u: 0.0000,train_z: 20.7796,train_c: 621.1819,train_n: 28.4497
[2023-02-26 14:37:47] Epoch[130/200], train_loss: 271227.2901, train_reconst: 2900.8132, train_u: 0.0000,train_z: 21.6138,train_c: 612.1154,train_n: 28.4958
[2023-02-26 14:38:00] Epoch[140/200], train_loss: 258140.4742, train_reconst: 2911.9901, train_u: 0.0000,train_z: 21.2212,train_c: 556.3544,train_n: 28.2061
[2023-02-26 14:38:13] Epoch[150/200], train_loss: 258732.2734, train_reconst: 2877.1855, train_u: 0.0000,train_z: 20.9135,train_c: 546.2261,train_n: 28.4047
[2023-02-26 14:38:27] Epoch[160/200], train_loss: 258729.8321, train_reconst: 2895.3583, train_u: 0.0000,train_z: 21.4752,train_c: 639.5163,train_n: 28.2653
[2023-02-26 14:38:40] Epoch[170/200], train_loss: 253291.8625, train_reconst: 2832.7214, train_u: 0.0000,train_z: 21.6450,train_c: 611.1145,train_n: 27.9201
[2023-02-26 14:38:53] Epoch[180/200], train_loss: 246630.5414, train_reconst: 2787.8257, train_u: 0.0000,train_z: 22.6931,train_c: 528.0627,train_n: 27.5186
[2023-02-26 14:39:06] Epoch[190/200], train_loss: 246964.7171, train_reconst: 2748.4639, train_u: 0.0000,train_z: 22.4038,train_c: 495.1112,train_n: 27.4699
[2023-02-26 14:39:19] Epoch[200/200], train_loss: 252978.1453, train_reconst: 2759.4455, train_u: 0.0000,train_z: 22.3224,train_c: 505.0142,train_n: 27.6357
[2023-02-26 14:39:19] === Finished training ===
[2023-02-26 14:39:19] === Restart Starfysh 2 ===
[2023-02-26 14:39:20] Initializing model parameters...
[2023-02-26 14:39:34] Epoch[10/200], train_loss: 5144781.3530, train_reconst: 3436.6599, train_u: 0.0000,train_z: 15.1901,train_c: 1274.1064,train_n: 32.0365
[2023-02-26 14:39:47] Epoch[20/200], train_loss: 1768550.6166, train_reconst: 3314.3621, train_u: 0.0000,train_z: 20.4146,train_c: 991.7330,train_n: 30.2856
[2023-02-26 14:40:00] Epoch[30/200], train_loss: 656172.5570, train_reconst: 3172.1109, train_u: 0.0000,train_z: 18.1547,train_c: 942.2212,train_n: 29.1228
[2023-02-26 14:40:14] Epoch[40/200], train_loss: 393922.9626, train_reconst: 3131.2998, train_u: 0.0000,train_z: 20.5968,train_c: 849.3791,train_n: 28.7167
[2023-02-26 14:40:27] Epoch[50/200], train_loss: 336856.9571, train_reconst: 3034.6232, train_u: 0.0000,train_z: 21.0180,train_c: 748.5228,train_n: 28.3635
[2023-02-26 14:40:40] Epoch[60/200], train_loss: 317561.7829, train_reconst: 3065.4551, train_u: 0.0000,train_z: 22.1590,train_c: 737.5154,train_n: 28.6295
[2023-02-26 14:40:54] Epoch[70/200], train_loss: 303843.1577, train_reconst: 3011.6274, train_u: 0.0000,train_z: 20.6004,train_c: 665.7587,train_n: 28.3140
[2023-02-26 14:41:07] Epoch[80/200], train_loss: 295229.1970, train_reconst: 2960.5557, train_u: 0.0000,train_z: 21.7357,train_c: 679.9993,train_n: 28.2011
[2023-02-26 14:41:20] Epoch[90/200], train_loss: 285783.2625, train_reconst: 3002.2135, train_u: 0.0000,train_z: 21.2462,train_c: 645.0632,train_n: 28.3669
[2023-02-26 14:41:34] Epoch[100/200], train_loss: 285386.6601, train_reconst: 2900.2783, train_u: 0.0000,train_z: 21.2149,train_c: 615.8687,train_n: 27.8645
[2023-02-26 14:41:47] Epoch[110/200], train_loss: 277198.2175, train_reconst: 2874.3669, train_u: 0.0000,train_z: 21.5431,train_c: 595.8658,train_n: 27.8273
[2023-02-26 14:42:00] Epoch[120/200], train_loss: 293643.0207, train_reconst: 3007.1657, train_u: 0.0000,train_z: 30.4903,train_c: 1043.7038,train_n: 29.0997
[2023-02-26 14:42:14] Epoch[130/200], train_loss: 272539.2865, train_reconst: 2923.8001, train_u: 0.0000,train_z: 26.5804,train_c: 639.8240,train_n: 28.5748
[2023-02-26 14:42:27] Epoch[140/200], train_loss: 293134.9835, train_reconst: 2887.3094, train_u: 0.0000,train_z: 25.2206,train_c: 610.0453,train_n: 27.7821
[2023-02-26 14:42:40] Epoch[150/200], train_loss: 278049.2724, train_reconst: 2877.5660, train_u: 0.0000,train_z: 23.9600,train_c: 619.1370,train_n: 27.8579
[2023-02-26 14:42:54] Epoch[160/200], train_loss: 276827.8915, train_reconst: 2916.7042, train_u: 0.0000,train_z: 26.0134,train_c: 580.2669,train_n: 27.5959
[2023-02-26 14:43:07] Epoch[170/200], train_loss: 258571.3330, train_reconst: 2784.7745, train_u: 0.0000,train_z: 23.0786,train_c: 531.4982,train_n: 27.2842
[2023-02-26 14:43:21] Epoch[180/200], train_loss: 259716.3813, train_reconst: 2715.3769, train_u: 0.0000,train_z: 23.3675,train_c: 547.0606,train_n: 27.3004
[2023-02-26 14:43:34] Epoch[190/200], train_loss: 257581.0188, train_reconst: 2758.3383, train_u: 0.0000,train_z: 23.7173,train_c: 565.6867,train_n: 27.7095
[2023-02-26 14:43:48] Epoch[200/200], train_loss: 255758.3611, train_reconst: 2693.5002, train_u: 0.0000,train_z: 22.9820,train_c: 529.0302,train_n: 27.2861
[2023-02-26 14:43:48] === Finished training ===
[2023-02-26 14:43:48] === Restart Starfysh 3 ===
[2023-02-26 14:43:49] Initializing model parameters...
[2023-02-26 14:44:02] Epoch[10/200], train_loss: 5148713.9679, train_reconst: 3279.5194, train_u: 0.0000,train_z: 12.6271,train_c: 1471.2258,train_n: 33.9123
[2023-02-26 14:44:15] Epoch[20/200], train_loss: 1760294.2466, train_reconst: 3107.1012, train_u: 0.0000,train_z: 16.9424,train_c: 1117.6885,train_n: 30.8703
[2023-02-26 14:44:29] Epoch[30/200], train_loss: 655864.3851, train_reconst: 3036.9442, train_u: 0.0000,train_z: 17.8658,train_c: 986.6379,train_n: 29.9048
[2023-02-26 14:44:42] Epoch[40/200], train_loss: 386880.2046, train_reconst: 2959.8202, train_u: 0.0000,train_z: 17.8974,train_c: 878.5079,train_n: 29.2460
[2023-02-26 14:44:55] Epoch[50/200], train_loss: 327343.0536, train_reconst: 2903.2014, train_u: 0.0000,train_z: 19.6593,train_c: 832.2568,train_n: 28.9550
[2023-02-26 14:45:09] Epoch[60/200], train_loss: 307915.7929, train_reconst: 2859.3862, train_u: 0.0000,train_z: 19.7203,train_c: 757.9050,train_n: 28.6089
[2023-02-26 14:45:22] Epoch[70/200], train_loss: 312378.1334, train_reconst: 2861.3384, train_u: 0.0000,train_z: 21.2654,train_c: 721.5658,train_n: 28.4705
[2023-02-26 14:45:36] Epoch[80/200], train_loss: 289900.4369, train_reconst: 2998.6637, train_u: 0.0000,train_z: 23.7672,train_c: 738.7895,train_n: 28.8920
[2023-02-26 14:45:49] Epoch[90/200], train_loss: 298646.3017, train_reconst: 2830.5458, train_u: 0.0000,train_z: 21.2444,train_c: 697.9101,train_n: 28.6755
[2023-02-26 14:46:02] Epoch[100/200], train_loss: 291199.0684, train_reconst: 2778.5722, train_u: 0.0000,train_z: 22.3973,train_c: 703.6338,train_n: 28.4261
[2023-02-26 14:46:16] Epoch[110/200], train_loss: 286329.9253, train_reconst: 2746.4538, train_u: 0.0000,train_z: 22.0374,train_c: 690.9767,train_n: 28.3981
[2023-02-26 14:46:29] Epoch[120/200], train_loss: 284016.9132, train_reconst: 2760.8564, train_u: 0.0000,train_z: 21.7751,train_c: 665.7985,train_n: 28.2625
[2023-02-26 14:46:43] Epoch[130/200], train_loss: 267556.9261, train_reconst: 2714.6238, train_u: 0.0000,train_z: 20.8677,train_c: 647.5391,train_n: 28.1156
[2023-02-26 14:46:56] Epoch[140/200], train_loss: 280491.7154, train_reconst: 2798.0905, train_u: 0.0000,train_z: 28.5670,train_c: 708.6395,train_n: 28.5049
[2023-02-26 14:47:09] Epoch[150/200], train_loss: 271827.7954, train_reconst: 2779.3461, train_u: 0.0000,train_z: 24.5762,train_c: 644.0507,train_n: 28.4178
[2023-02-26 14:47:23] Epoch[160/200], train_loss: 270827.5840, train_reconst: 2727.5231, train_u: 0.0000,train_z: 22.0569,train_c: 632.6442,train_n: 28.0516
[2023-02-26 14:47:36] Epoch[170/200], train_loss: 256638.2929, train_reconst: 2643.7381, train_u: 0.0000,train_z: 21.6659,train_c: 585.5864,train_n: 27.9296
[2023-02-26 14:47:49] Epoch[180/200], train_loss: 260109.5633, train_reconst: 2582.2564, train_u: 0.0000,train_z: 22.1102,train_c: 604.8425,train_n: 27.5478
[2023-02-26 14:48:05] Epoch[190/200], train_loss: 250769.0629, train_reconst: 2573.4146, train_u: 0.0000,train_z: 22.6120,train_c: 555.2672,train_n: 27.6903
[2023-02-26 14:48:18] Epoch[200/200], train_loss: 252123.1843, train_reconst: 2547.7461, train_u: 0.0000,train_z: 22.5638,train_c: 532.3200,train_n: 27.5801
[2023-02-26 14:48:18] === Finished training ===
(3). Downstream analysis
Base Starfysh Deconvolution
Parse Starfysh inference output
[ ]:
adata, adata_normed = visium_args.get_adata()
inference_outputs, generative_outputs = sf_model.model_eval(model,
adata,
visium_args,
poe=False,
device=device)
Visualize starfysh deconvolution results
Gene sig mean vs. inferred prop
[25]:
n_cell_types = gene_sig.shape[1]
idx = np.random.randint(0, n_cell_types)
post_analysis.gene_mean_vs_inferred_prop(inference_outputs,
visium_args,
idx=idx
)
Spatial visualizations:
Inferred density on Spatial map:
[26]:
plot_utils.pl_spatial_inf_feature(adata, feature='ql_m', cmap='Blues')
Inferred cell-type proportions (spatial map):
[27]:
plot_utils.pl_spatial_inf_feature(adata,
feature='qc_m',
# To display for specific cell types:
# factor = cell type or factor = [cell type1, ...]
vmax=0.1)
Inferred cell-type proportions on Z-space (UMAP):
[ ]:
plot_utils.pl_spatial_inf_feature(adata,
feature='qz_m',
# To display for specific cell types:
# factor = Cell type or factor = [Cell type1, ...]
factor=['Basal', 'LumA', 'MBC', 'Normal epithelial'],
spot_size=3,
vmax=0.2)
Infer cell-type specific expressions from each spot
[30]:
pred_exprs = sf_model.model_ct_exp(model,
adata,
visium_args,
device=device)
Plot spot-level expression (e.g. IL7R from Effector Memory T cells (Tem)):
[ ]:
sample_gene = 'IL7R'
sample_cell_type = 'Tem'
plot_utils.pl_spatial_inf_gene(adata,
factor=sample_cell_type,
feature=sample_gene)
Starfysh Deconvolution with histology integration
Parse Starfysh inference outputs
[23]:
inference_outputs_poe, generative_outputs_poe = sf_model.model_eval(model,
adata,
visium_args,
poe=True,
device=device)
Visualize starfysh deconvolution results
Gene sig mean vs. inferred prop
[24]:
n_cell_types = gene_sig.shape[1]
idx = np.random.randint(0, n_cell_types)
post_analysis.gene_mean_vs_inferred_prop(inference_outputs_poe,
visium_args,
idx=idx
)
Spatial visualizations:
[25]:
plot_utils.pl_spatial_inf_feature(adata, feature='ql_m', cmap='Blues')
Inferred cell-type proportions (spatial map):
[26]:
plot_utils.pl_spatial_inf_feature(adata,
feature='qc_m',
# To display for specific cell types:
# factor = cell type or factor = [cell type1, ...]
vmax=0.1)
Inferred cell-type proportions on Z-space (UMAP):
[27]:
plot_utils.pl_spatial_inf_feature(adata,
feature='qz_m',
# To display for specific cell types:
# factor = Cell type or factor = [Cell type1, ...]
factor=['Basal', 'LumA', 'MBC', 'Normal epithelial'],
vmax=0.2)
Infer cell-type specific expressions from each spot
[28]:
pred_exprs = sf_model.model_ct_exp(model,
adata,
visium_args,
poe=True,
device=device)
Plot spot-level expression (e.g. IL7R from Effector Memory T cells (Tem)):
[29]:
sample_gene = 'IL7R'
sample_cell_type = 'Tem'
plot_utils.pl_spatial_inf_gene(adata,
factor=sample_cell_type,
feature=sample_gene)
Save model & inferred parameters
[50]:
# Specify output directory
outdir = './results/'
if not os.path.exists(outdir):
os.mkdir(outdir)
# save the model
torch.save(model.state_dict(), os.path.join(outdir, 'starfysh_model.pt'))
# save `adata` object with inferred parameters
adata.write(os.path.join(outdir, 'st.h5ad'))