A Practical Workflow for Spatial Transcriptomics Data Analysis: From Data Acquisition to Advanced Analyses

Summary

This protocol presents a reproducible workflow for analyzing spatial transcriptomics data, guiding users from public data acquisition and Seurat-based quality control through integration, spatial feature detection, cell-type deconvolution, region-of-interest annotation, and cell–cell communication analysis, with practical checkpoints that support transparent execution.

Abstract

Spatial transcriptomics (ST) profiles genome-wide gene expression while preserving the two-dimensional spatial context of mRNA molecules within tissue sections, enabling studies of tissue architecture and microenvironment-associated biology. However, ST analysis remains challenging because data import, quality control, integration, deconvolution, spatial statistics, and visualization often require multiple software environments and reproducible parameter choices. This protocol presents a practical computational workflow for public ST datasets in R, beginning with data acquisition and software setup and proceeding through Seurat-based data loading, quality control, normalization, multi-sample integration, clustering, and spatially variable gene analysis. The workflow then applies complementary deconvolution strategies, including reference-guided SPOTlight analysis and unsupervised STdeconvolve topic modeling, followed by Giotto-based spatial cell-cell communication analysis and interactive region-of-interest (ROI) selection using a custom Python Dash application. By emphasizing script-based execution, explicit parameter rationales, expected outputs, and troubleshooting checkpoints, the protocol provides an adaptable framework for standard array-based ST datasets and related platforms after dataset- and platform-specific parameter evaluation.

Introduction

Spatial transcriptomics (ST) is a transformative family of technologies that measures genome-wide gene expression while retaining the spatial coordinates of messenger RNA (mRNA) molecules within tissue sections. ST methods include sequencing-based approaches that use position-barcoded arrays and in situ imaging approaches that map transcriptional signals within intact tissue microenvironments^1,2. By preserving spatial context, ST enables analysis of tissue architecture, cellular neighborhood organization, cell-cell communication, and microenvironment-associated biological processes that cannot be fully resolved after tissue dissociation³.

The rapid growth of public ST data repositories has created unprecedented opportunities for secondary analysis and method development³. Resources such as the CROST database curate hundreds of spatially resolved transcriptomic datasets across multiple species and technology platforms, while specialized collections like STOmicsDB focus on particular methodologies such as Stereo-seq^4,5. Despite this data abundance, computational analysis remains challenging due to the complexity of spatial data structures, diversity of analytical tools, and technical hurdles in implementing reproducible workflows^{6,7,8,9,10,11}.

To address the limitations of relying on a single software environment, an integrated computational workflow utilizing complementary analytical tools is presented herein. Existing comprehensive ST analysis ecosystems primarily include Seurat, Giotto, and Python-based frameworks such as Squidpy^6,7,12. Although Python-based tools like Squidpy offer extensive functionalities for spatial graph analysis, consolidating the primary computational pipeline within a single programming language environment minimizes cross-language technical hurdles. Consequently, the core pipeline is implemented mainly in R to reduce cross-language technical hurdles. Within this R-based workflow, Seurat is used for data loading, quality control, normalization, dimensionality reduction, visualization, and multi-sample integration, reflecting its common use in single-cell and spatial transcriptomics workflows. Giotto is then used for spatial network construction and ligand-receptor-based cell-cell communication analysis. Therefore, this pipeline links Seurat-based pre-processing and integration with Giotto-based spatial analysis while keeping the data transfer between the two toolsets explicit and reproducible.

Within this framework, two complementary deconvolution strategies are implemented: SPOTlight, a reference-guided method that uses scRNA-seq data to estimate cell-type proportions, and STdeconvolve, an unsupervised topic-modeling method that identifies latent transcriptional patterns^8,11. The outputs provide complementary views of spatial cellular heterogeneity but are not treated as a quantitative cross-validation unless users perform the optional concordance analysis described in the protocol. A custom Python Dash application, Select Spatial Spots, is integrated for interactive region-of-interest (ROI) annotation and exports standard coordinate-based annotation files that can be used in subsequent downstream analyses.

Regarding practical applicability, this workflow is intended primarily for standard array-based ST data (e.g., 55 µm resolution Visium) and may be adaptable to other tissue types after parameter evaluation. Key limitations should be considered before analysis. First, the reference-guided deconvolution module depends on a high-quality, tissue-matched scRNA-seq reference. Second, subcellular or near-single-cell platforms may require modified pre-processing, spatial bin aggregation, or image-based cell segmentation before integration². The representative mouse colon dataset is used as a demonstration case to show how the workflow can evaluate spatial domains and marker-defined tissue organization rather than as evidence of universal platform compatibility.

Protocol

All biological datasets analyzed in this protocol are publicly available and used strictly for demonstration. Specific data accessions and source repositories are provided in the relevant steps. The original datasets were generated by the original investigators in compliance with institutional ethical guidelines applicable to each source study. See the Table of Materials to verify all required software and R-package versions.

Hardware requirements: The computational memory required for this workflow scales with the number of samples and spots analyzed. For a typical spatial transcriptomics dataset (e.g., approximately 3,000 spots per sample across up to three samples), a standard workstation with a minimum of 16 GB of RAM is sufficient to execute the pipeline. However, 32 GB of RAM or higher is highly recommended to ensure optimal performance and stability, particularly during memory-intensive computational steps such as SCTransform normalization and matrix factorization during deconvolution.

1. Data acquisition and directory structure preparation

1. Obtain public spatial transcriptomics datasets (e.g., GEO accession GSE169749, with a submission date of March 26, 2021, and last update date of March 6, 2022) and single-cell RNA-seq reference data (e.g., GEO accession GSE264408, with a submission date of April 19, 2024, and last update date of December 10, 2024).
2. Download the raw data archive (e.g., GSE169749_RAW.tar) from the repository and extract the archive.
   NOTE: For the representative dataset used in this protocol, the extracted directory contains multiple files for each sample, with filenames typically incorporating the GEO sample accession (e.g., GSM5213483).
3. Organize the files into a standardized directory structure compatible with the Seurat Load10X_Spatial function.
   1. Create a main data directory (e.g., ./data/).
   2. Create a dedicated subdirectory for each sample (e.g., ./data/sample_A1/).
   3. Transfer (copy or move) the following essential files for each sample into the respective subdirectory.
      1. Place the filtered_feature_bc_matrix.h5 file into the main subdirectory of the sample.
      2. Create a spatial/ subfolder within each sample directory.
      3. Place the following files into the spatial/ subfolder: tissue_positions_list.csv.gz, scalefactors_json.json.gz, and tissue_hires_image.png.gz.
   4. Decompress the .gz files in the spatial/ folder. Ensure the original filenames remain exactly as required by the Load10X_Spatial function (e.g., tissue_positions_list.csv).
      NOTE: The final, organized structure for seamless loading must follow this pattern, using sample A1 as an example:
      ./data/sample_A1/filtered_feature_bc_matrix.h5
      ./data/sample_A1/spatial/tissue_positions_list.csv
      ./data/sample_A1/spatial/scalefactors_json.json
      ./data/sample_A1/spatial/tissue_hires_image.png

2. Software environment setup

1. Install R (version 4.4.3 or higher).
2. Obtain the analysis scripts (1_ReadSpatialData.R through 8_Giotto_Communication.R) from the GitHub repository (https://github.com/LeafLight/SpatialTranscriptomicsWorkFlow, at commit hash 2d85e18 for reproducibility).
3. Install the required R packages from CRAN and Bioconductor by executing the script setup.R.
   NOTE: A comprehensive list of the specific package versions utilized in this workflow is provided in the Table of Materials and is also available in the session_info.txt file within the associated GitHub repository (https://github.com/LeafLight/SpatialTranscriptomicsWorkFlow).
4. Install Giotto Suite by executing the installation commands provided in the official documentation.
   1. Navigate to the official installation URL to retrieve setup scripts: https://giottosuite.com/articles/installation.html
   2. Initialize the required Python environment and system dependencies according to the step-by-step instructions.
5. Install additional packages for enhanced visualization.
   1. Execute: install.packages(c("ggprism", "colorBlindness"))
6. Install the custom Select Spatial Spots tool. Ensure the operating system is Windows, macOS, or Linux, and Python version 3.8 or higher is installed. This step is optional.
   1. Obtain the tool by navigating to the GitHub repository (https://github.com/LeafLight/SelectSpatialSpots, at commit hash d20946e for reproducibility) and downloading the source code.
   2. Navigate into the tool's directory and install the Python dependencies: pip install -r requirements.txt
7. Install the glmGamPoi package by executing BiocManager::install("glmGamPoi") to accelerate the SCTransform normalization.
   NOTE: The core workflow requires only the packages listed in steps 2.1–2.3 and 2.6. The package Rfast2 is used to speed up the computation of Moran's I statistic. The optional packages in step 2.5 are for generating publication-ready plots with a prism-style theme (ggprism) and for accessing a colorblind-friendly palette (colorBlindness). The tool in step 2.6 enables the interactive spot selection functionality. Execute the sessionInfo() command and save the complete console output to a text file. This documents the exact software versions and package dependencies, which are critical for ensuring long-term reproducibility.

3. Spatial data loading and quality control (1_ReadSpatialData.R, 2_SpatialDataQC.R)

1. Read the spatial data into Seurat objects.
   1. Use Read10X_Image to manually load the high-resolution tissue image, specifying the image.dir and image.name.
   2. Use Load10X_Spatial with the image parameter set to the image object created in step 3.1.1 to create the Seurat object.
2. Calculate quality control metrics.
   1. Compute the percentage of mitochondrial reads using PercentageFeatureSet with the pattern ^mt-.
3. Visualize and interpret the data based on QC metrics.
   1. Generate violin plots of nCount_Spatial, nFeature_Spatial, and percent.mt using VlnPlot.
   2. Create spatial feature plots of these metrics using SpatialFeaturePlot to identify spots outside of the tissue area.
      NOTE: Visual checkpoint: Spots outside the tissue area typically display low UMI counts (nCount_Spatial < 500) and low gene detection (nFeature_Spatial < 200)
   3. For demonstration purposes, apply filters to remove low-quality spots (e.g., subset(seurat_obj, subset = nFeature_Spatial > 200 & nCount_Spatial > 500)). This is an optional step.
      NOTE: The primary goal of QC in spatial transcriptomics is to identify and annotate technical artifacts, such as spots outside the tissue. For the demonstration Visium colon data, nFeature_Spatial > 200 and nCount_Spatial > 500 remove low-complexity/background spots. For tissues with low RNA content or degraded sections, inspect violin and spatial feature plots before raising thresholds; for dense or high-UMI tissues, stricter thresholds may be appropriate. Avoid filtering solely on high expression values unless doublets, tissue folds, or clear artifacts are visible. Filtering is often not recommended as it physically removes spatial locations, which may disrupt the continuity of the tissue structure for downstream spatial analysis.

4. Data pre-processing, integration, and clustering (3_IntegrationAndClustering.R)

1. Normalize and pre-process individual samples.
   1. Apply the SCTransform normalization to each sample separately with assay = Spatial.
2. Integrate multiple samples.
   1. Prepare the list of SCT-normalized objects for integration. Ensure each object has an "RNA" assay by copying the "Spatial" assay: spatial_list[[1]][["RNA"]] <- spatial_list[[1]][["Spatial"]]
      ​NOTE: Copying the "Spatial" assay into a standard "RNA" assay slot serves as a necessary workaround to ensure compatibility with downstream integration functions originally designed for single-cell Seurat objects.
   2. Use SelectIntegrationFeatures and PrepSCTIntegration to identify shared variable features.
   3. Find integration anchors using FindIntegrationAnchors with normalization.method = "SCT".
   4. Integrate the data using IntegrateData with normalization.method = "SCT".
3. Perform dimensionality reduction and clustering on the integrated assay.
   1. Run PCA on the integrated data using RunPCA.
   2. Determine the optimal number of principal components (PCs) for downstream analysis by calculating the cumulative variance explained. Identify the elbow point programmatically (e.g., the point where cumulative variance exceeds 90% and marginal gain drops below 0.1%), using the code in 3_IntegrationAndClustering.R, lines 36-38. The resulting number of PCs is hereafter referred to as pc.use.
   3. Run RunUMAP with dims = 1:pc.use.
   4. Cluster the cells using FindNeighbors with dims = 1:pc.use and FindClusters with resolution = 0.5. Adjust resolution only after inspecting cluster stability and marker gene coherence.
   5. Perform differential expression analysis between target groups (e.g., "B1_colon_d14" versus "A1_colon_d0") using the FindMarkers function.
      ​NOTE: Visual checkpoint: A successful integration will display appropriate mixing of samples in the UMAP plot while maintaining biologically distinct clusters.
4. Identify spatially variable genes.
   1. For each original sample, run FindSpatiallyVariableFeatures using the "moransi" method on the "SCT" assay to compute spatial autocorrelation.

5. Single-cell reference data pre-processing (4_scDataPreProcessing.R)

1. Read the single-cell RNA-seq count matrix using Read10X and create a Seurat object.
2. Perform standard QC, normalization, and clustering.
   1. Calculate the percentage of mitochondrial reads (percent.mt) and filter cells (e.g., nFeature_RNA > 200 & nFeature_RNA < 7500 & percent.mt < 25).
   2. Normalize data using SCTransform, setting vars.to.regress = "percent.mt".
   3. Run PCA, UMAP, and cluster the cells using the dynamic PC selection method described in step 4.3.2.
3. Annotate cell types.
   1. Calculate module scores for canonical cell-type marker genes using AddModuleScore.
   2. Annotate clusters based on the module scores and known biology. Alternatively, import pre-computed annotations from metadata.

6. Reference-guided Deconvolution with SPOTlight (5_SPOTlight_Deconv.R)

1. Prepare data for SPOTlight.
   1. Convert the annotated single-cell Seurat object and the spatial Seurat object to SingleCellExperiment objects.
   2. Log-normalize the single-cell data using logNormCounts.
2. Run SPOTlight deconvolution.
   1. Identify hypervariable genes (HVGs) on the single-cell data using modelGeneVar and getTopHVGs.
   2. Compute cell-type marker genes using scoreMarkers and filter for high-quality markers (e.g., mean.AUC > 0.8).
   3. Downsample the single-cell reference for each cell type to a manageable number (e.g., 50 cells) to reduce computational time.
   4. Execute deconvolution using the SPOTlight function with weight_id = "mean.AUC", group_id = "cluster", and gene_id = "gene".
3. Visualize and export the results.
   1. Extract the deconvolution result matrix (cell-type proportions per spot).
   2. Use plotSpatialScatterpie to visualize the cellular composition on the spatial coordinates.
   3. Add the deconvolution results to the spatial Seurat object's metadata using AddMetaData. Expected output: a deconvolution matrix with cell-type proportion columns for each spatial spot, spatial scatterpie plots showing local cellular composition, and a Seurat object containing the deconvolution proportions as metadata.

7. Reference-free Deconvolution with STdeconvolve (7_STdeconvolve.R)

1. Prepare the spatial data.
   1. Extract the raw count matrix from the spatial Seurat object using GetAssayData with slot = "counts".
   2. Remove low-quality spots and genes using cleanCounts from STdeconvolve (e.g., min.lib.size = 100).
2. Identify the latent cell types.
   1. Further filter the corpus for genes expressed in a minimum fraction of spots using restrictCorpus (e.g., removeAbove=1.0, removeBelow = 0.05).
   2. Fit a Latent Dirichlet Allocation (LDA) model across a range of potential topic numbers (K) (e.g., Ks = seq(2, 9, by = 1)) using fitLDA.
   3. Select the optimal model based on minimum perplexity using optimalModel with opt = "min".
3. Analyze and visualize the results.
   1. Extract the cell-type (topic) proportions (theta) and gene profiles (beta) from the optimal model using getBetaTheta.
   2. To aid in the biological interpretation of the deconvoluted topics, import region-of-interest (ROI) annotations generated by the "Select Spatial Spots" tool (see step 9.3.3). Use these annotations as the groups parameter in the vizAllTopics function to project the deconvoluted cell-type proportions back onto the spatial coordinates and color-code the spots by their ROI. This is an optional step.
      NOTE: Expected output: STdeconvolve should return theta values representing topic proportions per spot and beta values representing gene profiles for each topic. If users wish to quantitatively compare the unsupervised STdeconvolve results with the reference-guided SPOTlight outputs (from Step 6), they can export both proportion matrices (theta from STdeconvolve and the SPOTlight matrix) and compute per-spot correlation or concordance metrics (e.g., Pearson or Spearman correlation) using standard R functions. Quantitative comparison between SPOTlight and STdeconvolve is optional and is not required for executing the core workflow; however, the relevant output matrices are available from Steps 6.3.1 and 7.3.1 for users who wish to perform such an analysis.

8. Spatial cell-cell communication using Giotto(8_Giotto_Communication.R)

1. Convert the Seurat object to a Giotto object.
   1. Use the createGiottoObject function, providing the raw count matrix and spatial coordinates.
2. Pre-process the Giotto object and add deconvolution results.
   1. Normalize the data using normalizeGiotto.
   2. Add the cell-type annotations (e.g., the major cell type from SPOTlight deconvolution) to the cell metadata using addCellMetadata.
3. Infer spatially-informed cell-cell communication.
   1. Create a spatial network using createSpatialNetwork with method = "Delaunay" and name = "Delaunay_network".
   2. Load a ligand-receptor database into the R environment. Ensure the loaded object is formatted as a data frame with two distinct columns representing the ligand and receptor gene symbols (e.g., a pre-compiled mouse network from Zenodo: https://zenodo.org/api/records/15168114/files/lr_network_mouse.csv/content, version v6, published on April 7, 2025).
   3. Run exprCellCellcom with cluster_column = "celltype_major" to identify significant ligand-receptor interactions between cell types that are in spatial proximity. Expected output: a result table containing ligand-receptor pairs, source and target cell-type combinations, log2 fold-change values, and adjusted P values; retain significant interactions (e.g., p.adj < 0.05) for visualization.

9. Interactive spot selection with select spatial spots (6_SelectSpatialSpots.R)

1. Prepare data for the interactive tool.
   1. Extract the spatial coordinates from the Seurat object using GetTissueCoordinates.
2. Format and export the data.
   1. Format the coordinate data frame to contain the exact column headers: CELL_ID, X, and Y.
   2. Export the formatted data frame to a CSV file.
3. Perform region-of-interest (ROI) analysis.
   1. Launch the custom Select Spatial Spots Dash application and load the CSV file.
   2. Interactively select spots based on spatial location.
   3. Export the list of selected spots and their assigned group/ROI labels as a new CSV file. Expected output: a CSV file with spot identifiers, X/Y coordinates, and assigned group/ROI labels that can be matched back to the Seurat object by CELL_ID.

Results

Workflow implementation and data integration illustrate major tissue features

The computational workflow was applied to mouse colon spatial transcriptomics data to illustrate the expected outputs across the analytical stages. As depicted in the workflow schematic (Figure 1), the pipeline began with data acquisition and quality control, where spatial feature plots delineated tissue boundaries (Figure 2A,B

Discussion

This protocol provides a comprehensive computational workflow for spatial transcriptomics data analysis that balances analytical depth with practical accessibility. The step-by-step approach guides researchers through the complete analytical pipeline, from initial data acquisition to advanced spatial analyses, while emphasizing critical decision points and potential pitfalls.

Several steps in the protocol warrant particular attention due to their impact on downstream results. The quality contr...

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
ggplot2	Posit Software, PBC	v4.0.0(CRAN)	Advanced data visualization
Giotto	Dries Lab	v4.2.2 (GitHub)	Spatial network and cell-cell communication analysis
patchwork	Thomas Lin Pedersen	v1.3.2 (CRAN)	Plot composition and arrangement
R software	R Foundation for Statistical Computing	v4.4.3	Core execution environment (macOS aarch64)
scater	Davis McCarthy et al.	v1.34.1 (Bioconductor)	Single-cell quality control and visualization
scran	Aaron Lun et al.	v1.34.0 (Bioconductor)	Single-cell variance modeling and marker detection
Select Spatial Spots (Custom Python Tool)	LeafLight	v1.0.0 (GitHub)	Interactive spatial region-of-interest (ROI) selection (https://github.com/LeafLight/SelectSpatialSpots)
Seurat	Satija Lab	v5.3.0 (CRAN)	Spatial data preprocessing, integration, and clustering
SeuratObject	Satija Lab	v5.2.0 (CRAN)	Data structures for single-cell and spatial data
SingleCellExperiment	Bioconductor Core Team	v1.28.1 (Bioconductor)	Standardized data container for scRNA-seq
SPOTlight	Marc Elosua-Bayes et al.	v1.10.0 (Bioconductor)	Reference-guided spatial deconvolution
Stdeconvolve	Jean Fan Lab	v1.3.2 (Bioconductor)	Unsupervised latent topic modeling
tidyverse	Posit Software, PBC	v2.0.0 (CRAN)	Core data manipulation and formatting suite