Quickstart#
To follow along with this tutorial, you can download the data from here.
Installation#
You can install SegTraQ via pip:
pip install segtraq
Let’s start by importing the required packages and setting the number of CPU cores to use for parallel processing.
[1]:
import warnings
import anndata as ad
import spatialdata as sd
import segtraq
# Suppress specific warnings for cleaner output
warnings.simplefilter(action="ignore", category=FutureWarning)
# Use all available CPU cores
segtraq.settings.n_jobs = -1
To get started with SegTraQ, we first need to load our data into a spatialdata object. Spatialdata-io provides readers for commonly used technologies, such as 10x Genomics Xenium, Vizgen MERSCOPE, and Nanostring CosMx. Segmentation tools such as Proseg also provide the option to export their results in spatialdata format directly. For more information on this, please refer to the io notebook.
[2]:
sdata = sd.read_zarr("../../data/xenium_5K_data/proseg3.zarr")
sdata
[2]:
SpatialData object, with associated Zarr store: /g/huber/projects/CODEX/segtraq/data/20260204_SegTraQ_sdata/xenium_5K_data/proseg3.zarr
├── Images
│ └── 'image': DataArray[cyx] (1, 1795, 1519)
├── Points
│ └── 'transcripts': DataFrame with shape: (2189794, 10) (3D points)
├── Shapes
│ ├── 'cell_boundaries': GeoDataFrame shape: (17956, 1) (2D shapes)
│ ├── 'cell_boundaries_z0': GeoDataFrame shape: (17602, 2) (2D shapes)
│ ├── 'cell_boundaries_z1': GeoDataFrame shape: (17824, 2) (2D shapes)
│ ├── 'cell_boundaries_z2': GeoDataFrame shape: (17818, 2) (2D shapes)
│ ├── 'cell_boundaries_z3': GeoDataFrame shape: (16381, 2) (2D shapes)
│ └── 'nucleus_boundaries': GeoDataFrame shape: (18608, 2) (2D shapes)
└── Tables
└── 'table': AnnData (17956, 5095)
with coordinate systems:
▸ 'global', with elements:
image (Images), transcripts (Points), cell_boundaries (Shapes), cell_boundaries_z0 (Shapes), cell_boundaries_z1 (Shapes), cell_boundaries_z2 (Shapes), cell_boundaries_z3 (Shapes), nucleus_boundaries (Shapes)
Next, we can create a SegTraQ object by passing the spatialdata object to the SegTraQ() constructor. This has the advantage that we only need to set keywords like cell IDs or transcript IDs once, and they will be used for all subsequent analyses. If you do not know how to set the keywords, you can just use segtraq.SegTraQ(sdata) and the constructor will tell you which keywords need to be set. More information on this is available in the io notebook.
[3]:
st = segtraq.SegTraQ(
sdata,
images_key="image", # where the image is stored
points_cell_id_key="assignment", # what cells are called in the transcript table
points_background_id=None, # what unassigned transcripts are called
points_gene_key="gene", # what genes are called in the transcript table
points_qv_key=None, # what transcript quality is called in the transcript table
tables_area_key=None, # what the cell areas are called in the cell table
tables_cell_id_key="cell", # what cells are called in the expression matrix
shapes_cell_id_key="cell", # what cells are called in the segmentation shapes
tables_centroid_x_key="centroid_x", # what the cell centroids are called
tables_centroid_y_key="centroid_y", # what the cell centroids are called
)
/g/huber/users/lazic/src/SegTraQ/src/segtraq/SegTraQ.py:160: RuntimeWarning: No area column specified for tables. Area will be automatically computed from shapes.
validate_spatialdata(
/g/huber/users/lazic/src/SegTraQ/src/segtraq/utils.py:2139: UserWarning: Duplicate IDs detected in column 'cell_id' for shapes 'nucleus_boundaries'. Instead of using cell_id as index, resetting the current index and renaming it to `segtraq_id`.
nucleus_shapes = _ensure_index(
/g/huber/users/lazic/src/SegTraQ/src/segtraq/SegTraQ.py:188: UserWarning: Filtering control probes matching prefixes ('NegControlProbe_', 'antisense_', 'NegControlCodeword', 'BLANK_', 'Blank-', 'NegPrb', 'DeprecatedCodeword_', 'UnassignedCodeword_', 'Intergenic_Region_') and transcripts with qv < 20 from transcript points in-place; control probes are also removed from the expression table if present. Configure `min_qv`, `control_prefixes`, and `inplace` via `filter_kwargs`.
sdata_new = _filter_control_and_low_quality_transcripts(
/g/huber/users/lazic/src/SegTraQ/src/segtraq/SegTraQ.py:188: RuntimeWarning: Gene sets differ between assigned transcript points and the expression table. 0 genes occur only in points (e.g. []), and 3 occur only in the table (e.g. ['GABBR2', 'TSPY1', 'HPV16-E7']). Genes only occuring in the table can be introduced by cropping (`SpatialData`), but differences may also be due to the filtering used to generate the transcript data and expression matrix. Please check that `filter_kwargs`, particularly `min_qv` and `control_prefixes`, are consistent with the filtering used to generate the expression table.
sdata_new = _filter_control_and_low_quality_transcripts(
/g/huber/users/lazic/src/SegTraQ/src/segtraq/SegTraQ.py:188: RuntimeWarning: Cell IDs differ between assigned transcript points and the expression table. 0 cells occur only in points (e.g. []), and 191 occur only in the table (e.g. [np.uint32(17921), np.uint32(2055), np.uint32(9229), np.uint32(15375), np.uint32(16916)]). Cells only occuring in the table can reflect zero-count cells, but differences may also be due to the filtering used to generate the transcript data and expression matrix. Please check that `filter_kwargs`, particularly `min_qv` and `control_prefixes`, are consistent with the filtering used to generate the expression table.
sdata_new = _filter_control_and_low_quality_transcripts(
Now that we have a SegTraQ object, we can call run_all() to compute all SegTraQ metrics in one go.
[4]:
st.run_all()
/tmp/ipykernel_461365/3888125565.py:1: UserWarning: Skipping `run_supervised`: metric(s) could not be computed (`adata_ref` is required when running supervised module.).
st.run_all()
As you can see, SegTraQ raises a warning that some metrics could not be computed. This is because some metrics require additional information, such as an scRNA-seq reference dataset or a 3D segmentation. These can be passed into the run_all() method as keyword arguments.
[5]:
# reference scRNA-seq dataset for computing supervised metrics
adata_ref = ad.read_h5ad("../../data/xenium_5K_data/BC_scRNAseq_Janesick.h5ad")
# additional arguments to compute volume metrics
volume_kwargs = {
"run_ovrlpy": True,
"heterotypic_overlap_kwargs": {
"shapes_key_list": [
"cell_boundaries_z0",
"cell_boundaries_z1",
"cell_boundaries_z2",
"cell_boundaries_z3",
]
},
}
# run_all() with additional arguments
st.run_all(adata_ref=adata_ref, ref_cell_type="celltype_major", ref_raw_counts_layer="raw", volume_kwargs=volume_kwargs)
WARNING: adata.X seems to be already log-transformed.
Running vertical adjustment
Creating gene expression embeddings for visualization
determining pseudocells
found 2850 pseudocells
sampling expression:
0%| | 0/4 [00:00<?, ?it/s]
25%|███████████████████████████████████████████████████████████████▌ | 1/4 [00:02<00:07, 2.45s/it]
50%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████ | 2/4 [00:04<00:04, 2.36s/it]
75%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▌ | 3/4 [00:07<00:02, 2.38s/it]
100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 4/4 [00:09<00:00, 2.38s/it]
Modeling 12 pseudo-celltype clusters;
Creating signal integrity map
0%| | 0/16 [00:00<?, ?it/s]
6%|███████████████▊ | 1/16 [00:05<01:18, 5.24s/it]
12%|███████████████████████████████▋ | 2/16 [00:10<01:13, 5.26s/it]
19%|███████████████████████████████████████████████▍ | 3/16 [00:15<01:06, 5.13s/it]
25%|███████████████████████████████████████████████████████████████▎ | 4/16 [00:17<00:45, 3.75s/it]
31%|███████████████████████████████████████████████████████████████████████████████ | 5/16 [00:21<00:44, 4.01s/it]
38%|██████████████████████████████████████████████████████████████████████████████████████████████▉ | 6/16 [00:27<00:46, 4.67s/it]
44%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████▋ | 7/16 [00:33<00:45, 5.08s/it]
50%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▌ | 8/16 [00:38<00:41, 5.15s/it]
56%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▎ | 9/16 [00:44<00:36, 5.28s/it]
62%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▌ | 10/16 [00:49<00:31, 5.24s/it]
69%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▎ | 11/16 [00:55<00:27, 5.43s/it]
75%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████ | 12/16 [01:00<00:21, 5.45s/it]
81%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▊ | 13/16 [01:06<00:16, 5.59s/it]
88%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▌ | 14/16 [01:12<00:11, 5.56s/it]
94%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▎ | 15/16 [01:16<00:05, 5.27s/it]
100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 16/16 [01:20<00:00, 5.06s/it]
WARNING: adata.X seems to be already log-transformed.
INFO Creating graph using `None` transform and `1` libraries.
/g/huber/users/lazic/src/SegTraQ/src/segtraq/SegTraQ.py:1892: RuntimeWarning: neighbors_key='spatial_connectivities' not found in adata.obsp. A neighborhood graph based on Delaunay will be computed.
return sp.marker_balanced_accuracy(
All metrics are written to the spatialdata object.
[6]:
st.sdata.tables["table"]
[6]:
AnnData object with n_obs × n_vars = 17956 × 5095
obs: 'cell', 'original_cell_id', 'centroid_x', 'centroid_y', 'centroid_z', 'component', 'volume', 'surface_area', 'scale', 'region', 'leiden_subset_cells17765_res0.2_seed42', 'leiden_subset_cells11191_res0.2_seed0', 'leiden_subset_cells11191_res0.2_seed1', 'leiden_subset_cells11191_res0.2_seed2', 'leiden_subset_cells11191_res0.2_seed3', 'leiden_subset_cells11191_res0.2_seed4', 'transferred_cell_type', 'pearson_score', 'gene_count', 'mean_transcripts_per_gene', 'transcript_count', 'transcript_density', 'num_polygons', 'cell_area', 'perimeter', 'circularity', 'solidity', 'convexity', 'elongation', 'eccentricity', 'compactness', 'num_nuclei', 'similarity_nucleus_cell', 'similarity_nucleus_cell_p_value', 'nucleus_id', 'iou', 'nucleus_fraction', 'similarity_nucleus_cytoplasm', 'similarity_nucleus_cytoplasm_p_value', 'border_admixture_score', 'border_admixture_p_value', 'similarity_top_bottom', 'similarity_top_bottom_p_value', 'vertical_signal_integrity', 'heterotypic_overlap_area', 'heterotypic_overlap_fraction', 'marker_balanced_accuracy', 'n_evaluated_positive_markers', 'n_evaluated_negative_markers', 'contamination_strength', 'num_outside_cell_all_genes', 'num_in_nucleus_overlap_all_genes', 'num_in_cytoplasm_all_genes', 'perc_outside_cell_all_genes', 'perc_nucleus_all_genes', 'perc_cytoplasm_all_genes', 'distance_to_cell_centroid_norm_all_genes', 'distance_to_cell_membrane_norm_all_genes', 'skew_dist_to_cell_membrane_all_genes'
var: 'gene', 'total_count', 'lambda_bg_0', 'lambda_bg_1', 'lambda_bg_2', 'lambda_bg_3', 'total', 'unassigned', 'perc_unassigned'
uns: 'spatialdata_attrs', 'proseg_run', 'num_cells', 'num_transcripts', 'num_genes', 'perc_unassigned_transcripts', 'cluster_connectedness', 'silhouette_score', 'mean_purity', 'mean_ari', 'segtraq_markers', 'spatial_neighbors', 'contamination_strength_matrix', 'contamination_evaluable_cells_matrix', 'mutually_exclusive_coexpression_rate'
obsm: 'spatial'
obsp: 'transition_counts', 'spatial_connectivities', 'spatial_distances'
You can now use these metrics to perform quality control and filtering of your data, for example by filtering out spurious samples, cells, or comparing segmentation methods. More details on the available metrics and how to use them can be found in the module-specific notebooks.
Session Info#
[7]:
print(sd.__version__) # spatialdata
0.8.0