Stereo-seq Mouse E16.5 Embryo

Key features:

• Apply Gaussian smoothing to address data noise and sparsity
• Construct a gene co-expression network across samples using smoothed dataset
• Group genes into spatial modules with hard and soft thresholding on the co-expression network

In [1]:

import anndata as ad
import numpy as np
import pandas as pd
import scanpy as sc
import matplotlib.pyplot as plt

import smoothie
smoothie.suppress_warnings() # Suppress a few warnings for cleaner output

import anndata as ad import numpy as np import pandas as pd import scanpy as sc import matplotlib.pyplot as plt import smoothie smoothie.suppress_warnings() # Suppress a few warnings for cleaner output

Load AnnData¶

Smoothie requires data input in the AnnData format. A guide for constructing the AnnData data structure for both imaging-based and sequencing-based spatial transcriptomics outputs can be found in /guides/convert_to_anndata.md

In [ ]:

# Bash code to generate AnnData file from STOmics bin1 output: https://db.cngb.org/stomics/mosta/download/
# python -m smoothie.convert_to_anndata /path/to/E16.5_E1S1_GEM_bin1.tsv.gz \
#     --x-col x \
#     --y-col y \
#     --gene-col geneID \
#     --count-col MIDCounts \
#     --format tsv.gz \
#     --chunksize 1000000 \
#     --output mouse_E16.5_E1S1_raw_dnb.h5ad

# Bash code to generate AnnData file from STOmics bin1 output: https://db.cngb.org/stomics/mosta/download/ # python -m smoothie.convert_to_anndata /path/to/E16.5_E1S1_GEM_bin1.tsv.gz \ # --x-col x \ # --y-col y \ # --gene-col geneID \ # --count-col MIDCounts \ # --format tsv.gz \ # --chunksize 1000000 \ # --output mouse_E16.5_E1S1_raw_dnb.h5ad

In [2]:

# Load your spatial transcriptomics data
adata = ad.read_h5ad("/path/to/mouse_E16.5_E1S1_raw_dnb.h5ad")

# Load your spatial transcriptomics data adata = ad.read_h5ad("/path/to/mouse_E16.5_E1S1_raw_dnb.h5ad")

Quality Control¶

The table below shows suggested QC values given the resolution of your AnnData spatial coordinates.

• Note: If you have cell-segmented imaging / fine resolution data, then you should use the "Cell-sized" QC settings. If you have true individual transcripts or fine resolution data, then you should use the Fine / Imaging QC settings.

Resolution	min_counts (filter_cells)	min_counts (filter_genes)	min_cells (filter_genes)
"Cell-sized" (8-10 µm)	10 - 100	100 - 500	10 - 100
Binned (20-50 µm)	50 - 500	100 - 500	10 - 100
Fine / Imaging (< 1-2 µm)	1	100 - 500	10 - 100

In [3]:

# Filter spots with low UMI counts
sc.pp.filter_cells(adata, min_counts=1) # Only remove DNA nanoballs with 0 UMIs!!

# Filter genes with low counts/unique spots
sc.pp.filter_genes(adata, min_counts=100)

# Filter spots with low UMI counts sc.pp.filter_cells(adata, min_counts=1) # Only remove DNA nanoballs with 0 UMIs!! # Filter genes with low counts/unique spots sc.pp.filter_genes(adata, min_counts=100)

Normalization¶

Normalization depends on the resolution of your AnnData spatial coordinates.

• Note: If you have cell-segmented imaging / fine resolution data, then you should use the "Cell-sized" normalization method. If you have true individual transcripts or fine resolution data, then you should use the Fine / Imaging normalization method.

Resolution	Normalization Method
"Cell sized" or Binned Data (8-50 micron resolution)	CPT + Log1p normalization
Fine / Imaging (<1-2 micron resolution)	Log1p normalization only

In [4]:

# # CHOOSE: CPT + Log1p normalization (8-50 micron resolution data)
# sc.pp.normalize_total(adata, target_sum=1e3)
# sc.pp.log1p(adata)

# CHOOSE: Log1p normalization only (<1-2 micron resolution data)
sc.pp.log1p(adata)

# # CHOOSE: CPT + Log1p normalization (8-50 micron resolution data) # sc.pp.normalize_total(adata, target_sum=1e3) # sc.pp.log1p(adata) # CHOOSE: Log1p normalization only (<1-2 micron resolution data) sc.pp.log1p(adata)

Step 1: Run Gaussian Smoothing¶

Here, we run Gaussian smoothing on the spatial dataset address data sparsity and noise. Key parameters for function smoothie.run_parallelized_smoothing:

gaussian_sd: float

• The standard deviation for the Gaussian distribution, which determines smoothing degree. Generally across high-resolution S.T. platforms, a gaussian_sd corresponding to 20-30 microns is ideal.
• gaussian_sd is calculated as target_microns * micron_to_unit_conversion
  • For Slide-seq, micron_to_unit_conversion = 1.5456 (1 micron equals 1.5456 Slide-seq spatial units)
  • For Stereo-seq, micron_to_unit_conversion = 2 (1 micron equals 2 Stereo-seq spatial units)
  • A table of micron_to_unit_conversion values for spatial transcriptomics platforms is included in /guides/micron_to_unit_conversion_table.md

grid_based_or_not: boolean

• True = grid-based smoothing: smooth only at imposed hexagonal grid points. Good for subcellular to "cell-sized" resolution data (<1-10 micron).
• False = in-place smoothing: smooth at every spatial location in the dataset. Good for "cell-sized" to binned resolution data (10-50 micron).

min_spots_under_gaussian

• Minimum required number of data points within radius (3 * gaussian_sd) for smoothing to occur at a given location.
• Default is 10-100. Raise value if low density regions around the perimeter of tissue appear noisy in modules. Lower value if you have large spot size (e.g. original Visium 100 micron resolution data or binned 20-50 micron resolution).

In [5]:

# Determine gaussian_sd ---------------------------------------------------
target_microns = 20 # 20-30 microns is ideal for smoothing
micron_to_unit_conversion = 2 # IMPORTANT! For Stereo-seq, this value is 2.
# ------------------------------------------------------------------------

sm_adata = smoothie.run_parallelized_smoothing(
    adata,
    grid_based_or_not=True,  # Select True or False,
    gaussian_sd=target_microns * micron_to_unit_conversion,
    min_spots_under_gaussian=100 # Adjust as needed
)

# Determine gaussian_sd --------------------------------------------------- target_microns = 20 # 20-30 microns is ideal for smoothing micron_to_unit_conversion = 2 # IMPORTANT! For Stereo-seq, this value is 2. # ------------------------------------------------------------------------ sm_adata = smoothie.run_parallelized_smoothing( adata, grid_based_or_not=True, # Select True or False, gaussian_sd=target_microns * micron_to_unit_conversion, min_spots_under_gaussian=100 # Adjust as needed )

Checking dataset point density...
Fitting grid to tissue coordinates...
Variation in X-direction: 25549
Variation in Y-direction: 17467
Number of points in fitted grid: 218394
Gaussian smoothing will run in 36 chunks.
Smoothing chunk: 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 
Total runtime for grid-based Gaussian smoothing: 633.39 seconds.

Step 2: Measure and Cluster Gene Co-expression Network¶

First, we compute the pairwise Pearson correlation coefficient between all smoothed gene surface pairs.

Next, we evaluate optimal hyperparameter choices for network construction and clustering. The function smoothie.select_clustering_params calculates network clustering statistics across different choices for pcc_cutoff and clustering_power.

Choose values for pcc_cutoff and clustering_power that balance the following goals:

1. High gene margin scores for high gene module confidence
2. Stable number of modules (n_clusters) compared to other pcc_cutoff values
3. High gene recall for more gene-rich analyses

Refer to /docs/tutorials for examples on selecting these parameters.

In [6]:

# Compute pairwise gene correlations
pearsonR_mat, _ = smoothie.compute_correlation_matrix(sm_adata.X)

# Compute pairwise gene correlations pearsonR_mat, _ = smoothie.compute_correlation_matrix(sm_adata.X)

compute_correlation_matrix running...
Total runtime for compute_correlation_matrix: 1567.91 seconds.

In [7]:

# Clustering hyperparameter selection
smoothie.select_clustering_params(
    gene_names=sm_adata.var_names, # don't change
    pearsonR_mat=pearsonR_mat, # don't change
    output_folder=None,
    pcc_cutoffs=[0.3, 0.4, 0.5, 0.6, 0.7, 0.8], # lower pcc_cutoffs require more runtime
    clustering_powers=[1, 3, 5, 7, 9],
    min_genes_for_module=5
)

# Clustering hyperparameter selection smoothie.select_clustering_params( gene_names=sm_adata.var_names, # don't change pearsonR_mat=pearsonR_mat, # don't change output_folder=None, pcc_cutoffs=[0.3, 0.4, 0.5, 0.6, 0.7, 0.8], # lower pcc_cutoffs require more runtime clustering_powers=[1, 3, 5, 7, 9], min_genes_for_module=5 )

Running network clustering with parameters:
PCC hard cutoff: 0.3, clustering soft power: 1, 3, 5, 7, 9, 
PCC hard cutoff: 0.4, clustering soft power: 1, 3, 5, 7, 9, 
PCC hard cutoff: 0.5, clustering soft power: 1, 3, 5, 7, 9, 
PCC hard cutoff: 0.6, clustering soft power: 1, 3, 5, 7, 9, 
PCC hard cutoff: 0.7, clustering soft power: 1, 3, 5, 7, 9, 
PCC hard cutoff: 0.8, clustering soft power: 1, 3, 5, 7, 9, 

Optionally, data shuffling and smoothing can be done to help determine the lower bound for pcc_cutoff selection. Here, use the same parameters as before in smoothie.run_parallelized_smoothing.

In [ ]:

# p95, p99, p999 = smoothie.compute_shuffled_correlation_percentiles(
#     adata,
#     grid_based_or_not=True, # Choose same value as in smoothie.run_parallelized_smoothing
#     gaussian_sd=target_microns * micron_to_unit_conversion, # Choose same value as in smoothie.run_parallelized_smoothing
#     min_spots_under_gaussian=100, # Choose same value as in smoothie.run_parallelized_smoothing
#     seed=0
# )

# p95, p99, p999 = smoothie.compute_shuffled_correlation_percentiles( # adata, # grid_based_or_not=True, # Choose same value as in smoothie.run_parallelized_smoothing # gaussian_sd=target_microns * micron_to_unit_conversion, # Choose same value as in smoothie.run_parallelized_smoothing # min_spots_under_gaussian=100, # Choose same value as in smoothie.run_parallelized_smoothing # seed=0 # )

Based on the outputs above, we can select pcc_cutoff and clustering_power hyperparameters for smoothie.make_spatial_network.

pcc_cutoff: The Pearson correlation coefficient (PCC) hard threshold for network construction. Only correlations above this value are retained in the gene network. Higher cutoff values result in smaller, stronger average correlation networks.

clustering_power: A soft thresholding parameter that controls how much weaker PCC values are suppressed. Higher values favor more modules.

In [8]:

edge_list, node_label_df = smoothie.make_spatial_network(
    pearsonR_mat=pearsonR_mat, # don't change
    gene_names=sm_adata.var_names, # don't change
    pcc_cutoff=0.4,
    clustering_power=2,
    output_folder='./out'
)

# Filter node_label_df for only genes within modules of size 2 or more.
modules_df = node_label_df.groupby('module_label').filter(lambda x: len(x) >= 2)

# Examine modules
modules_df

edge_list, node_label_df = smoothie.make_spatial_network( pearsonR_mat=pearsonR_mat, # don't change gene_names=sm_adata.var_names, # don't change pcc_cutoff=0.4, clustering_power=2, output_folder='./out' ) # Filter node_label_df for only genes within modules of size 2 or more. modules_df = node_label_df.groupby('module_label').filter(lambda x: len(x) >= 2) # Examine modules modules_df

Out[8]:

	name	module_label	degree	weighted_degree	rescaled_weighted_degree	clustering_coeff	margin_score
0	Tuba1a	1	1478	812.904511	147.721621	0.247334	0.660701
1	Rtn1	1	1372	755.943484	136.671708	0.270682	0.773819
2	Lrba	1	1614	863.133571	130.268922	0.216896	0.648015
3	Tmsb10	1	1579	845.484626	128.663715	0.225472	0.686678
4	Basp1	1	1432	775.636928	126.642462	0.259535	0.769660
...	...	...	...	...	...	...	...
5939	Ndor1	269	1	0.405006	0.000070	0.000000	1.000000
5940	Pmel	270	1	0.506851	0.031714	0.000000	1.000000
5941	Dct	270	1	0.506851	0.031714	0.000000	1.000000
5942	Vkorc1	271	17	7.916375	0.511063	0.933824	0.974660
5943	Gm21974	271	8	3.707906	0.377918	1.000000	0.994388

5944 rows × 7 columns

Finally, we can import the saved edge_list.csv and node_labels.csv files into Cytoscape for interactive network visualization! More info in /guides/cytoscape_instructions.md

Visualize Smoothed Gene Expression and Modules¶

In [9]:

sm_adata = smoothie.rotate_spatial(sm_adata, angle_degrees=-90)

sm_adata = smoothie.rotate_spatial(sm_adata, angle_degrees=-90)

In [10]:

# Plot individual gene
smoothie.plot_gene(
    sm_adata=sm_adata,
    gene_name="Alb",  # Replace with your gene
    output_folder="./gene_plots",
    spot_size=50,  # Adjust for optimal visualization
)

# Plot individual gene smoothie.plot_gene( sm_adata=sm_adata, gene_name="Alb", # Replace with your gene output_folder="./gene_plots", spot_size=50, # Adjust for optimal visualization )

In [11]:

# Plot all gene modules
smoothie.plot_modules(
    sm_adata,
    node_label_df,
    output_folder='./module_plots',
    min_genes=3, # Minimum number of genes in a module to plot the module
    spot_size=50, # Adjust for optimal visualization
    plots_per_row=10,
    dpi=150 # low resolution here for github, better to choose 300-600!
)

# Plot all gene modules smoothie.plot_modules( sm_adata, node_label_df, output_folder='./module_plots', min_genes=3, # Minimum number of genes in a module to plot the module spot_size=50, # Adjust for optimal visualization plots_per_row=10, dpi=150 # low resolution here for github, better to choose 300-600! )

Analyze Specific Gene or Gene Set Correlations¶

1. Find top correlated genes to a gene of interest.
2. Build a network focused on a specific set of genes with a more permissive PCC cutoff.

In [12]:

# Find top correlated genes for a gene of interest (GOI)
GOI_correlations = smoothie.get_correlations_to_GOI(
    pearsonR_mat=pearsonR_mat, # don't change
    gene_names=sm_adata.var_names, # don't change
    GOI="Alb",
    reverse_order=False,  # False for top correlations, True for anti-correlations
    plot_histogram=True
)
    
print(f"\nTop 10 gene correlations with GOI:")
print(GOI_correlations[:10])

# Find top correlated genes for a gene of interest (GOI) GOI_correlations = smoothie.get_correlations_to_GOI( pearsonR_mat=pearsonR_mat, # don't change gene_names=sm_adata.var_names, # don't change GOI="Alb", reverse_order=False, # False for top correlations, True for anti-correlations plot_histogram=True ) print(f"\nTop 10 gene correlations with GOI:") print(GOI_correlations[:10])

Top 10 gene correlations with GOI:
[['Alb' 0.9999899864196777]
 ['Afp' 0.9688883423805237]
 ['Trf' 0.9522647261619568]
 ['Apoa2' 0.9450106024742126]
 ['Apoa1' 0.9382757544517517]
 ['Gm20425' 0.9149119853973389]
 ['Ahsg' 0.9096430540084839]
 ['Serpina1b' 0.9063860177993774]
 ['Apoc2' 0.8834957480430603]
 ['Apom' 0.8820343613624573]]

In [13]:

# Build targeted network for your gene set of interest
Slc_family_genes = [gene for gene in sm_adata.var_names if gene.startswith("Slc")]

# Build targeted network
geneset_edge_list, geneset_node_label_df = smoothie.make_geneset_spatial_network(
    pearsonR_mat=pearsonR_mat, # don't change
    gene_names=sm_adata.var_names, # don't change
    node_label_df=node_label_df, # don't change
    gene_list=Slc_family_genes,
    low_pcc_cutoff=0.2,  # More permissive pcc_cutoff than main network
    output_folder='./geneset_network_out',
    intra_geneset_edges_only=False # Only edges within gene set?
)

print(f"Gene set network: {len(geneset_edge_list)} edges")
print(f"Genes included: {geneset_node_label_df['geneset_member'].sum()}")

# Build targeted network for your gene set of interest Slc_family_genes = [gene for gene in sm_adata.var_names if gene.startswith("Slc")] # Build targeted network geneset_edge_list, geneset_node_label_df = smoothie.make_geneset_spatial_network( pearsonR_mat=pearsonR_mat, # don't change gene_names=sm_adata.var_names, # don't change node_label_df=node_label_df, # don't change gene_list=Slc_family_genes, low_pcc_cutoff=0.2, # More permissive pcc_cutoff than main network output_folder='./geneset_network_out', intra_geneset_edges_only=False # Only edges within gene set? ) print(f"Gene set network: {len(geneset_edge_list)} edges") print(f"Genes included: {geneset_node_label_df['geneset_member'].sum()}")

Gene set network: 630 edges
Genes included: 222

In [15]:

geneset_node_label_df[geneset_node_label_df['geneset_member']]

geneset_node_label_df[geneset_node_label_df['geneset_member']]

Out[15]:

	name	module_label	degree	weighted_degree	rescaled_weighted_degree	clustering_coeff	margin_score	weak_module_label	geneset_member
71	Slc25a4	1	1322	655.207787	53.767721	0.298799	0.322105	-1	True
398	Slc17a6	1	506	238.604374	11.390278	0.849521	0.477502	-1	True
488	Slc3a1	1	476	219.741037	8.116987	0.850252	-0.191400	-1	True
635	Slc1a2	1	337	153.476621	4.678490	0.902943	0.256547	-1	True
690	Slc35f1	1	341	153.698776	3.723374	0.963999	1.000000	-1	True
...	...	...	...	...	...	...	...	...	...
6088	Slc4a10	-1	-1	-1.000000	-1.000000	-1.000000	-1.000000	1	True
6089	Slc44a5	-1	-1	-1.000000	-1.000000	-1.000000	-1.000000	1	True
6090	Slc43a3	-1	-1	-1.000000	-1.000000	-1.000000	-1.000000	9	True
6091	Slc2a10	-1	-1	-1.000000	-1.000000	-1.000000	-1.000000	3	True
6092	Slc35b2	-1	-1	-1.000000	-1.000000	-1.000000	-1.000000	2	True

222 rows × 9 columns