- Baidu Research, Sunnyvale, CA, United States

Human and animal tissues consist of heterogeneous cell types that organize and interact in highly structured manners. Bulk and single-cell sequencing technologies remove cells from their original microenvironments, resulting in a loss of spatial information. Spatial transcriptomics is a recent technological innovation that measures transcriptomic information while preserving spatial information. Spatial transcriptomic data can be generated in several ways. RNA molecules are measured by *in situ* sequencing, *in situ* hybridization, or spatial barcoding to recover original spatial coordinates. The inclusion of spatial information expands the range of possibilities for analysis and visualization, and spurred the development of numerous novel methods. In this review, we summarize the core concepts of spatial genomics technology and provide a comprehensive review of current analysis and visualization methods for spatial transcriptomics.

## 1 Introduction

Quantification of gene expression has important applications across various aspects of biology. Understanding the spatial distribution of gene expression has helped to answer fundamental questions in developmental biology (Asp et al., 2019; Rödelsperger et al., 2021), pathology (Maniatis et al., 2019; Chen et al., 2020), cancer microenvironment (Berglund et al., 2018; Thrane et al., 2018; Ji et al., 2020; Moncada et al., 2020), and neuroscience (Shah et al., 2016; Moffitt et al., 2018; Close et al., 2021). Two widely used methods for gene expression quantification are fluorescent *in situ* hybridization (FISH) and next-generation sequencing. With FISH, fluorescently-labeled RNA sequences are used as probes to identify its naturally occurring complementary sequence in cells while preserving the spatial location of the target sequences (Schwarzacher and Heslop-Harrison, 2000). Traditionally, the number of target sequences simultaneously identified by *in situ* hybridization is restricted by the number of fluorescent channels, making this method suitable for targeted gene detection. On the other hand, next-generation sequencing methods use a shotgun approach to quantify RNA molecules across the entire transcriptome (Metzker, 2010). To achieve transcriptome-wide quantification, RNA must be first isolated and purified, which removes RNA molecules from their native microenvironment. Even with single-cell sequencing, where the cellular origin of RNA molecules is preserved, spatial information of cells can only be inferred but not directly measured (Shapiro et al., 2013; Gawad et al., 2016).

Various approaches have been made to measure gene expression while preserving spatial information. Tomo-seq applied the principle of tomography to measure spatial transcriptomic information in 3D. In tomo-seq, tissue samples are sliced by cryosection and measured with RNA-seq. Each measurement corresponds to the average gene expression within a slice. Measurements are taken along multiple axes to reconstruct pixel-wise 3D gene expression information. (Junker et al., 2014). LCM-seq isolates single cells with laser capture microscopy (LCM) and measures captured cells with single-cell RNA-sequencing. LCM can capture cells of desired types and with specific spatial locations of the tissue specimen (Nichterwitz et al., 2016). While these methods retain the spatial location of RNA-seq measurements, they suffer from high labor costs and incomplete spatial coverage. In this review, we cover recent advances in spatial transcriptomic methods that attempt to address these challenges. In addition, we provide a comprehensive review of analysis and visualization techniques for spatial transcriptomic datasets.

The following sections are organized as follows (Figure 1). Section 2 discusses the latest developments in experimental spatial transcriptomic technologies. Section 3 discusses preprocessing of spatial transcriptomic data, an essential step prior to any analysis or visualization. Section 4 dissects methods whose inputs are gene expression without spatial coordinates. This includes dimensionality reduction, clustering, and cell-type identification. Section 5 describes methods whose inputs are gene expression combined with spatial coordinates. This includes identification of spatially coherent gene expression patterns and identification of spatial domains. Section 6 describes methods that analyze the interaction between cells or genes. All methods reviewed are listed in Table 1. This includes the identification of cell-to-cell communication and gene interaction. We note that other reviews on spatial transcriptomic technology (Dries et al., 2021a) have been published during the peer review of this article.

**FIGURE 1**. An overview of spatial transcriptomic tasks. **(A)** Spatial transcriptomic datasets map gene expression measurements to their respective locations. **(B)** A spatial transcriptomic dataset can be analyzed in gene expression space, irrespective of spatial locations. Tasks such as clustering and cell-type identification fall into this category. **(C)** Spatial information can be used jointly with gene expression to detect spatial expression patterns and spatial domains. **(D)** These two sources of information can also be used to detect cell-cell and gene-gene interactions.

**TABLE 1**. Current analysis and visualization tools for spatial transcriptomic datasets (accession date: 12/22/2021).

## 2 Spatial Transcriptomic Technologies

Integration of spatial information with transcriptome-wide quantification has given rise to the emerging field of spatial transcriptomics. Currently, spatial transcriptome quantification falls into three broad categories (Table 2). First, spatial barcoding methods ligate oligonucleotide barcodes with known spatial locations to RNA molecules prior to sequencing (Ståhl et al., 2016; Rodriques et al., 2019; Vickovic et al., 2019; Liu et al., 2020; Chen et al., 2021; Cho et al., 2021; Stickels et al., 2021). Both barcodes and RNA molecules are jointly sequenced, and spatial information of sequenced RNA molecules can be recovered from associated barcodes. Second, *in situ* hybridization methods coupled with combinatorial indexing can vastly increase the number of RNA species identified (Lubeck et al., 2014; Chen et al., 2015; Moffitt et al., 2016; Eng et al., 2019). The latest *in situ* hybridization methods can detect around 10,000 RNA species from a given sample (Eng et al., 2019). Third, *in situ* sequencing method uses fluorescent-based direct sequencing to read out base pair information from RNA molecules in their original spatial location (Lee et al., 2014; Wang et al., 2018).

Several metrics need to be considered when selecting a method for a specific application (Table 2). Methods employing *in situ* hybridization provide subcellular resolution. Leveraging super-resolution microscopy, *in situ* hybridization methods can achieve a resolution of ∼10nm, sufficient to distinguish single RNA molecules (Schermelleh et al., 2019). In addition, *in situ* methods require no PCR amplification of cDNA, thus avoiding amplification bias. However, the number of RNA species detected by *in situ* methods is limited by the indexing scheme. The current detection limit is ∼10,000 genes but will likely improve in the future. Furthermore, the area examined by *in situ* methods is limited by the field-of-view of the microscope objective lens. In contrast, spatial barcoding followed by shotgun sequencing can in principle sample the whole transcriptome. This is ideal if the target molecules are unknown *a priori*. Spatial barcoding can also examine larger tissue areas, making it ideal for larger samples such as tissue slices from the brain. However, the density of measurement spots limits the spatial resolution of current spatial barcoding methods, ranging from multicellular to subcellular. In addition, shotgun sequencing inevitably suffers from PCR amplification bias (Aird et al., 2011), as well as “dropout” when sequencing read depth is insufficient (Kim et al., 2020). Thus far, we have provided an overall picture of different spatial transcriptomic methods and their characteristics. Because this review focuses on analysis and visualization of spatial transcriptomics, readers who wish to understand the experimental details can refer to comprehensive reviews elsewhere (Crosetto et al., 2015).

## 3 Preprocessing

Spatial transcriptomic datasets add a new dimension to transcriptomic analyses. Spatial coordinates of cells enable novel analyses such as spatial differential expression (Svensson et al., 2018) and cell-cell interaction (Cang and Nie, 2020). Similar to single-cell RNA-seq datasets, a spatial transcriptomic dataset can be represented by a gene-by-cell count matrix. A second matrix of coordinates is attached to the cell dimension of the count matrix to represent spatial information. Comprehensive toolkits such as Space Ranger can process raw sequence reads into count matrices. Taking a microscope image and FASTQ files as input, Space Ranger can perform alignment, tissue and fiducial detection, barcode/UMI counting, and feature-spot matrices generation.

Various preprocessing steps may be performed prior to any analysis. First, genes and cells may be filtered based on a threshold specific to the dataset. For example, a cell may be removed if it has 1) less than 1,000 expressed genes or 2) a high proportion of mitochondria RNA. A gene may be removed if it is detected in less than ten cells (Wolf et al., 2018; Lun et al., 1000). Transformation of count data may be performed according to downstream modeling assumptions. Methods modeling raw counts do not require any transformation (Sun et al., 2020). Otherwise, gene expression per cell may be normalized to have the same total library size such that expression levels are comparable across cells. The gene expression matrix may then be log-transformed and be regressed against confounders such as batch effect, percentage of mitochondria genes, and other technical variations. Although preprocessing steps mentioned above are widely adopted, the exact configuration should follow input data modality and modeling assumptions, and there is no one-size-fits-all strategy.

### 3.1 Gene Expression Normalization

Current spatial transcriptomic techniques introduce unwanted technical artifacts. Raw data commonly exhibit spot-to-spot variation and high dropout rates, which may impact downstream analyses. Several normalization strategies have been created to address these challenges. Due to the similarity between spatial transcriptomics and scRNA-seq, many normalization methods for spatial transcriptomics data are inspired by scRNA-seq studies.

A widely-used normalization tool is scran, a method based on the summation of expression values and deconvolution of pooled size factors (Lun et al., 2016). In the first step, expression values of all cells in the data set are averaged to serve as a reference. The cells are then partitioned into different pools, where the summation of expression values in each pool is normalized against the reference to generate a pool-based size factor. A linear system can be constructed by repeating the above normalization over multiple pools. Finally, the normalized cell-based counterparts can be calculated by solving the linear system with standard least-squares methods, i.e., deconvolving the pool-based size factor to individual cells. By representing the individual cells with multiple pools of cells, scran is capable of avoiding estimation inaccuracy in the presence of stochastic zeroes and is robust to differentially expressed genes. Similar to scran, a number of methods adopt the global scale factor strategy, where one normalization factor is applied to each cell, and all genes in this cell share the same factor. When the relationship between transcript-specific expression and sequencing depth is not shared across genes, such strategy will likely lead to overcorrection for weakly and moderately expressed genes. To address the problem, Bacher et al. proposed SCnorm, a quantile-regression based method that can estimate the dependence of expression on sequencing depth for each gene (Bacher et al., 2017). Then genes are grouped based on the similarity of dependence, and a second quantile regression is used to estimate a shared scale factor within each group.

Lytal et al. conducted an empirical survey to evaluate the effectiveness of seven single-cell normalization methods. Based on the experimental results over several real and simulated data sets, the study concludes that there is no “one-size-fits-all” normalization technique for every data set (Lytal et al., 2020). Further, Saiselet et al. investigated whether normalization is warranted for spatial transcriptomic datasets. They discovered that variation of total read counts is related to morphology and local cell density. Therefore, total counts per spot are biologically informative and do not necessarily need to be normalized out (Saiselet et al., 2020).

## 4 Analysis and Visualization in the Expression Domain

A first step in the spatial transcriptomic analysis is to identify the cell type (for datasets of single-cell resolution) or cell mixture (for datasets of multicellular resolution) of each spatial unit or spot. Cell type identification usually starts with the dimensionality reduction technique to reduce time and space complexity for downstream analysis. The reduced representations are used to cluster cells based on the assumption that cells of the same type fall into the same cluster.

### 4.1 Clustering

The selection of clustering techniques is critical for obtaining good clustering results. Certain methods with assumptions about cluster shapes may not be suitable for spatial genomic data. For example, K-means clustering assumes that the shapes of clusters are spherical and that clusters are of similar size (Kanungo et al., 2002), and Gaussian mixture models assume that points with each cluster follow a Gaussian distribution (Reynolds, 2009). These assumptions are rarely satisfied by spatial transcriptomic data.

Agglomerative clustering methods are a class of methods that iteratively aggregate data points into clusters. These methods do not carry assumptions about the shape and size of clusters. At each iteration, data points are aggregated to optimize a pre-defined metric. Popular agglomerative clustering methods include hierarchical agglomerative clustering (Johnson, 1967) and community detection methods such as Louvain (Blondel et al., 2008) and Leiden (Traag et al., 2019) algorithms. Hierarchical agglomerative clustering is initialized by treating each point as its own cluster. Each iteration aggregates two clusters with the closest distance to form a new cluster until no clusters can be merged. Community detection methods, i.e., Louvain (Blondel et al., 2008) and Leiden (Traag et al., 2019) algorithms, have seen wide adoption in the single-cell and the spatial transcriptomics community. Both algorithms try to iteratively maximize the modularity, which can be understood as the difference between the number of observed and expected edges. Intuitively, a tightly connected community or cluster should have a large number of observed edges relative to the expected number of edges. The Louvain algorithm is initialized by assigning each node to its own community. At each iteration, each node moves from its own community to all neighboring communities, and changes in modularity are calculated. The node is moved to the community, which results in the largest increase in *H*. At the end of each iteration, a new network is built by aggregating all nodes within the same community, and a new iteration begins. The procedure will terminate when the increase in *H* can no longer be achieved.

These general-purpose methods can be combined into more sophisticated pipelines tailored towards single-cell clustering. SC3 is an ensemble clustering method in which multiple clustering outcomes are merged into a consensus. SC3 first calculates distance matrices using the Euclidean distance, as well as Pearson and Spearman correlations. Spectral clustering is performed on these distance matrices with a varying number of eigenvectors. These results were combined to assign a consensus cluster membership to each point (Kiselev et al., 2017). Seurat uses a smart local moving (SLM) algorithm (Waltman and Van Eck, 2013) to perform modularity-based clustering. Seurat first constructs a distance matrix based on canonical correlation vectors and a shared nearest neighbor (SNN) graph based on the distance matrix. The SNN graph is used as an input to the SLM algorithm to find clusters (Butler et al., 2018). SIMLR calculates a distance matrix as a weighted sum of multiple distance kernels and solves for a similarity matrix to minimize the product between the distance and similarity matrices. To ensure a fixed number of connected components, SIMLR uses constrained optimization to encourage a block diagonal structure in the similarity matrix (Wang et al., 2017).

### 4.2 Identification of Cell Types

Identification of cell types starts by defining cell-type specific genes or marker genes. A straightforward approach is to perform differential expression analysis (McCarthy et al., 2012; Love et al., 2014) between all pairs of clusters. Genes that are consistently over-expressed in one cluster are considered the cluster’s marker genes. This is the approach implemented in scran (Lun et al., 1000) and Mast (Finak et al., 2015).

Another method, scGeneFit, uses a label-aware compression method to find marker genes (Dumitrascu et al., 2021). Given cell-by-gene expression matrix and corresponding cell labels inferred from clustering results, scGeneFit finds a projection onto a lower-dimensional space, in which cells with the same labels are closer in the lower-dimensional space than cells with different labels. The projection is constrained such that the axes in the lower-dimensional space align with a single gene. Therefore, the marker genes will be the set of axes in the lower-dimensional space that best conserves label structures. The marker genes can then be matched with an expert-curated list of cell-type specific genes to infer cell types (Kim and Volsky, 2005; Subramanian et al., 2005). Other methods directly map unknown cell types onto a reference dataset, bypassing the target gene identification step. Scmap projects the query cells onto the reference cell types from other experiments and datasets (Kiselev et al., 2018). The known reference cluster is represented by its centroid, and the projection is carried out by a fast approximate k-nearest-neighbor (KNN) search by cluster using product quantization (Jégou et al., 2010), where a similarity matrix between the query cell and reference clusters is used as the distance in KNN search. Another reference-based method is SingleR (Aran et al., 2019). The method proceeds by first identifying variable genes among cell types in the reference set. Next, SingleR calculates the Spearman correlation between each single cell and the reference variable genes. Multiple correlation coefficients within each cell type are aggregated to form one correlation per reference cell type per single cell. Only the top 80% of correlation values are selected to remove random noise. In the fine-tuning step, the correlation analysis is iteratively re-run but only for the top cell types from the previous step, and the lowly correlated cell types are removed. Eventually, the cell type with the top correlation is assigned to the query single-cell. Using SingleR, the authors identified a novel disease-associated macrophage subgroup between monocyte-derived and alveolar macrophages. Cortal et al. proposed a clustering-free multivariate statistical method named Cell-ID for gene signature extraction and cell identification (Cortal et al., 2021). Cell-ID first performs a dimensionality reduction on the cell-by-gene expression matrix using the multiple correspondence analysis (MCA). Both cells and genes are simultaneously projected in a common low-dimensional space, where the distance between a gene and a cell represents the specific degree between them. According to the distance, Cell-ID can build up a gene-rank for each cell, and the top-ranked genes are defined as the cell’s gene signature. With the gene signature of the query cell, Cell-ID can perform automatic cell type and functional annotation via the hypergeometric tests against reference marker gene lists and/or gene signatures of reference single-cell datasets. The authors demonstrated the consistently reproducible gene signatures across diverse benchmarks, which helps to improve biological interpretation at the individual cell level. Unlike the above approaches, JSTA uses deep learning for cell-type identification and incorporates three distinct and interactive components: a segmentation map and two deep neural network-based cell type classifiers for pixel-level and cell-level classification (Littman et al., 2021). JSTA first trains a taxonomy-based cell-level classifier with the external data from the Neocortical Cell Type Taxonomy (NCTT) set (Yuste et al., 2020). Then the segmentation map and pixel-level classifier are iteratively refined with an expectation-maximization (EM) algorithm. Specifically, the segmentation map is initialized by a classical image segmentation algorithm watershed (Roerdink and Meijster, 2000) and paired with the trained cell-level type classifier to predict the current cell (sub)types. Given the local mRNA density at each pixel as the input, the pixel-level classifier is optimized to closely match each pixel’s current cell type assignment. Next, the updated pixel-level classifier reclassifies the cell types of all border pixels, and the resulting segmentation map requires an update of the cell-level classification, which further triggers an update of pixel-level classifier training. This learning process is repeated until convergence. The eventual segmentation map tends to maximize consistency between local RNA density and cell-type expression priors. Abdellaal et al. benchmarked 22 broadly used cell identification methods on 27 publicly available single-cell RNA data. Interested readers are referred to (Abdelaal et al., 2019).

### 4.3 Visualization of Gene Expression in Low Dimensions

The identified clusters can be visualized to ensure cells assigned to the same cluster are close in expression space. Dimensionality reduction techniques are necessary to project the high dimensional data into 2D or 3D. Principal component analysis (PCA) is widely adopted in the single-cell and spatial transcriptomic literature (Wold et al., 1987). This method identifies linear combinations of the original dimensions, or principal components (PC), that maximize the projection variance from data points onto the principal components (Figure 2A). The principal components can be computed in an iterative way: the first PC can point in any direction to maximize the variance of projections, and each subsequent PC is orthogonal to previous PCs (Tsuyuzaki et al., 2020).

**FIGURE 2**. Comparison between principal component analysis and t-SNE. **(A)** Principal component analysis iteratively identifies vectors that minimize the sum of squared distances to the direction of the vector. Each vector is orthogonal to all previously selected vectors. **(B)** t-SNE calculates a pairwise similarity based on the probability density function of the Gaussian distribution in the original high dimension space. The points are randomly projected to a low dimensional space and iteratively refined so that the similarity in low dimension matches that in high dimension. At each iteration, similar pairs attract, and dissimilar pairs repel each other.

In contrast to PCA, manifold learning is a class of non-linear dimensionality reduction techniques that aims to project the data to a lower dimension while maintaining the distance relations in the original high-dimension space; points close to each other in the original space will be close in the low-dimensional space (Figure 2B). Uniform manifold approximate and projection (UMAP) and t-distributed stochastic neighbor embedding (t-SNE) are two manifold learning methods widely adopted in single-cell and spatial transcriptomic literature (Van der Maaten and Hinton, 2008; McInnes et al., 2018). Both methods follow a two-step procedure. In the first step, a similarity matrix is computed based on a pre-defined distance metric. In the second step, all data points are placed in a low-dimensional Euclidean space such that the structure of the similarity matrix is preserved. This step is initialized by randomly placing data points in the low-dimensional space. At each iteration, data points are moved according to the similarity matrix from the high-dimensional space; points with high similarity in the high-dimensional space will attract, and those with low similarity will repel. Because optimization is done iteratively, UMAP and t-SNE results are stochastic and vary between runs. Random seeds are needed for reproducibility. The two methods differ in their construction of similarity matrix. In t-SNE, a distance matrix is calculated according to probability density functions (PDF) of the Gaussian distribution in the high-dimension space and PDFs of the t-distribution in the low-dimension embedding. In UMAP, an adjacency matrix is constructed by extending a sphere whose radius depends on the local density of nearby points; two points are connected if their spheres overlap. In practice, UMAP is faster than t-SNE and tends to preserve the high-dimensional structure better.

## 5 Analysis and Visualization in the Spatial Domain

An important question in spatial transcriptomic data analysis is to identify genes whose expression follow coherent spatial patterns. Genes with spatial expression patterns are critical determinants of polarity and anatomical structures. For example, the gene *wingless* is a member of the *wnt* family that plays a central role in anterior-posterior pattern generation during the embryonic development of *Drosophila melanogaster*. It is expressed in alternating stripes across the entire embryo (van den Heuvel et al., 1989). Another example is the neocortex of mammalian brains, which contain six distinct layers. Each layer consists of different types of neurons and glial cells that express cell-type specific marker genes (Lui et al., 2011). Spatial transcriptomic data enables unbiased transcriptome-wide identification of spatially expressed genes, but it is excessively labor-intensive to visually examine all genes. This prompted the development methods including SpatialDE (Svensson et al., 2018), trendsceek (Edsgärd et al., 2018), and Spark (Sun et al., 2020).

### 5.1 Identification of Genes with Spatial Expression Patterns

SpatialDE (Svensson et al., 2018) uses a Gaussian process to model gene expression levels. Intuitively, a Gaussian process model treats all data points as observations from a random variable that follows a multivariate Gaussian (MVN) distribution (Wang, 2020). To test whether expression levels follow a spatial pattern, the authors specify a null model, in which the covariance matrix is diagonal, and an alternative model, in which the covariance matrix follows a radial basis function kernel:

where *K*(*x*_{i}*, x*_{j}) is the covariance between *i*th and *j*th measurement; *x*_{i} and *x*_{j} represent the spatial coordinates of the *i*th and the *j*th measurement;

Trendsceek (Edsgärd et al., 2018) uses a marked point process model in which each point of measurement, or a spot, is treated as a point process, and each point is marked with a gene expression value. To decide whether a gene whose expression follows a spatial pattern, trendsceek test whether the probability of finding two marks given the distance between two points deviates from what would be expected if the marks were randomly distributed over points. To calculate the null distribution given no spatial pattern, trendsceek implements a sampling procedure in which marks are permuted with the location of points fixed. In practice, such sampling procedure is computationally expensive and makes trendsceek only suitable for small datasets.

Spark (Sun et al., 2020) uses a generalized linear spatial model (GLSM) to directly model count data (McCullagh and Nelder, 1983; Gotway and Stroup, 1997), which results in better power than SpatialDE. A simplified model is presented below:

Where *y(s)* is the gene expression of sample *s*. *x(s)* represents covariates such as batch effect and library size for sample *s*. The second term *b(s)* is the spatial correlation pattern modeled as a Gaussian process. The last term *Monte Carlo* methods are the gold standard for parameter estimation for GLSM but are computationally expensive. Instead, the authors developed a penalized quasi-likelihood (PQL) estimation procedure to make computation tractable for large datasets (Wedderburn, 1974; McCullagh and Nelder, 1983; Breslow and Clayton, 1993). Spark produces well-calibrated *p*-values and is more powerful than trendsceek and SpatialDE through a series of simulation experiments.

### 5.2 Identification of Spatial Domains

Spatially coherent domains often underly important anatomical regions (Figure 3A). A motivating example is the histological staining of cancer tissue slides. Cancer regions and normal tissues can be visually distinguished due to differential affinities to staining agents. This enables pathologists to grade and stage individual cancer tissue slides according to the location and size of the cancer regions (Fletcher, 2007). Spatial transcriptomics enables histology-like identification of spatial domains. Regular histology slides can be visualized conveniently with RGB pixels. In contrast, spatial transcriptomic data cannot be directly visualized because each spot (i.e., pixel) in spatial transcriptomic data has a dimension equal to the number of genes. This prompts the development of methods to detect spatial domains, including BayesSpace (Zhao et al., 2021), SpatialDE (Svensson et al., 2018), and a hidden Markov random field (HMRF) method (Zhu et al., 2018).

**FIGURE 3**. Visualization of gene expression in the Euclidean space. **(A)** Spatially coherent genes and spatial domains can be visualized as 2D images. **(B)** Spot deconvolution methods estimate the proportion of each cell type within each spot. Pie charts are routinely used to represent cell type proportions within each spot. **(C)** Spot super-resolution methods estimate the cell type of sub-pixels based on correlation with neighbor spots. In this case, each spot of the original dataset is divided into nine spots in the super-resolved dataset.

The three methods share a common assumption that hidden spatial domains can be described by latent variables, which are not directly observed but can be inferred from observed gene expression values. However, these methods use different modeling assumptions to infer latent variables. Zhu et al. (Zhu et al., 2018) developed an HMRF-based method, a widely adopted model in the image processing community to identify patterns in 2D images (Li, 2000; Blake et al., 2011), to identify spatial domains. An HMRF has two components: it uses a Markov random field to describe the joint distribution of latent variables and a set of observed examples that depends on them. The latent variables are assumed to satisfy the Markov property, in which any node in the network is conditionally independent of other nodes given its neighbors. Following this assumption, a Markov random field of latent variables can be decomposed into a set of subgraphs, called cliques, which gives rise to the observed gene expression. The parameters of the model by Zhu et al. are estimated with an EM algorithm (Dempster et al., 1977; Moon, 1996).

Both SpatialDE (Svensson et al., 2018) and BayesSpace (Zhao et al., 2021) model observed gene expression values as a mixture of Gaussian random variables. The means of the Gaussian random variables are determined by the spatial domain membership. In SpatialDE, the mean expression value of each spatial domain is described by a Gaussian process, whose covariance follows a radial basis function kernel. The observed expression follows a Gaussian distribution centered around the mean expression value of a given spatial domain. The posterior distribution of parameters and the latent spatial domain membership is estimated by variational inference. Different from SpatialDE, BayesSpace uses a diagonal matrix to model the covariance of the mean expression of each spatial domain. The observed gene expression is modeled as a Gaussian random variable centered around the mean expression and has a diagonal covariance matrix modeled as a Wishart random variable. BayesSpace uses a Markov chain Monte Carlo (MCMC) method to estimate model parameters (Geyer, 1992).

While the above methods consider spatial genes and spatial domain as two separate tasks, SpaGCN proposed a graph convolutional network-based (GCN) approach to address these two tasks jointly (Hu et al., 2021b). With the integration of gene expression, spatial location, and histology information, SpaGCN models spatial dependency of gene expression for clustering analysis of spatial domains and identification of domain enriched spatial variable genes (SVG) or meta genes. SpaGCN first converts the spatial transcriptomics data into an undirected weighted graph of spots, and the graph structure represents the spatial dependency of the data. Next, a GCN (Kipf and Welling, 2016) is utilized to aggregate gene expression information from the neighboring spots and update every spot’s representation. Then, SpaGCN adopts an unsupervised clustering algorithm (Xie et al., 2016) to cluster the spots iteratively, and each identified cluster will be considered as a spatial domain. The resulting domains guide the differential expression analysis to detect the SVG or meta genes with enriched expression patterns in the identified domains.

### 5.3 Spot Deconvolution and Super-resolution

Because spots in the spatial transcriptomic dataset may not correspond to cell boundaries, several additional features can be included when plotting on the spatial domain. When the spatial transcriptomic measurement technology has a multicellular resolution, spots can be decomposed into constituent cell types. A 2D array of pie charts can be used to represent the cell types’ percentages of spots, as demonstrated in DSTG (Figure 3B). To enable the investigation of cellular architecture at higher resolution, DSTG uses a GCN to uncover the cellular compositions within each spot (Song and Su, 2021). DSTG first leverages single-cell RNA-seq data to construct pseudo spatial transcriptomic (pseudo-ST) data by selecting two to eight single cells from the same tissue and combining their transcriptomic profiles. This pseudo-ST data is designed to mimic the cell mixture in the real spatial transcriptomic data and provide the basis for model training. Via canonical correlation analysis, DSTG identifies a link graph of spots with the integration of the pseudo-ST data and the real spatial transcriptomic data. A GCN (Kipf and Welling, 2016) iteratively updates the representation of each spot by aggregating its neighborhoods’ information. The GCN model is trained in a semi-supervised manner, where the known cell compositions of the pseudo-ST nodes are served as the labeled data, and the real spatial transcriptomic nodes are the prediction targets. The resulting cell type proportions can be displayed as a pie chart at each spot (Figure 3B).

While cell type deconvolution provides an estimation of cellular constituents, it does not directly increase the resolution of the dataset. BayesSpace uses a Bayesian model to increase the resolution to the subspot level, which approaches single-cell resolution with the Visium platform (Figure 3C). The model specification is similar to the spatial domain detection model described above, except that unit of analysis is the subspot rather than the spot. Since gene expression is not observed at the subspot level, BayesSpace models it as another latent variable and estimates it using MCMC. The increase in resolution is different across measurement technology. For square spots (Ståhl et al., 2016), BayesSpace by default divides each spot into nine subspots. For hexagonal spots like Visium, they are divided into six subspots by default. The subspots can be visualized in Euclidean space similar to regular spots.

### 5.4 Visualization in Euclidean Space

After obtaining spatial genes and domains, visualization in the Euclidean space is relatively straightforward. Spatial genes can be visualized by plotting their log-transformed expression values. Spatial domains can be colored by mean expression values or by their identities. Several packages such as Giotto (Dries et al., 2021b), Scanpy (Wolf et al., 2018), Seurat (Hao et al., 2021), and Squidpy (Palla et al., 2021) provide functionalities to plot spatial transcriptomic data in Euclidean space.

## 6 Analysis and Visualization in the Interaction Domain

Cell signaling describes the process in which cells send, receive, process, and transmit signals within the environment and with themselves. Based on the signaling distance and the sender-receiver identities, cell signaling can be classified into autocrine, paracrine, endocrine, intracrine, and juxtacrine (Bradshaw and Dennis, 2009). It serves critical functions in development (Wei et al., 2004), immunity (Dustin and Chan, 2000), and homeostasis (Taguchi and White, 2008) across all organisms. For example, the Hedgehog signaling pathway is involved in tissue patterning and orientation, and aberrant activations of hedgehog signaling lead to several types of cancers (Taipale and Beachy, 2001). Single-cell datasets enable correlation analysis to unravel cell-to-cell interaction (Krishnaswamy et al., 2014; Friedman et al., 2018; Wirka et al., 2019). Due to the lack of spatial information, single-cell analysis cannot distinguish short-distance (juxtacrine and paracrine) and long-distance (endocrine) signaling. Spatial transcriptomic datasets provide the spatial coordinate of each cell or spot and enable spatial dissection of cell signaling.

### 6.1 Cell-to-Cell Interaction

Cell signaling frequently occurs between cells in spatial proximity. Giotto takes spatial proximity into consideration to identify cell-to-cell interaction. It first constructs a spatial neighborhood network to identify cell types that occur in spatial proximity. Each node of the network represents a cell, and pair of neighboring cells are connected through an edge. The neighbors of each cell can be determined by extending a circle of a predefined radius, selecting the k-nearest neighbors, or constructing a Delaunay network (Chen and Xu, 2004). Cell types connected in the network more than expected are considered interacting. Giotto permutes the cell type labels without changing the topology of the network and calculates the expected frequencies between every pair of cell types. *p*-values are derived based on where the observed frequency falls on the distribution of expected frequencies.

Another method, SpaOTsc, leverages both single-cell and spatial transcriptomic data for a comprehensive profile of spatial interaction (Cang and Nie, 2020). It uses an optimal transport algorithm to map single-cell to spatial transcriptomic data. An optimal transport is a function that maps a source distribution to a target distribution while minimizing the amount of effort with respect to a predefined cost function (Villani, 2009). SpaOTsc generates a cost function based on the expression profile dissimilarity of shared genes across the single-cell and the spatial transcriptomic datasets. The optimal transport plan maps single cells onto spatial locations. SpaOTsc then formulates cell-to-cell communication as a second optimal transport problem between sender and receiver cells. The expression of ligand and receptor genes are used to estimate sender and receiver cells, and the spatial distance is the cost function. The resultant optimal transport plan represents the likelihood of cell-to-cell communication.

### 6.2 Ligand-Receptor Pairing

Another aspect of cell signaling is the pairings between ligands and receptors. Giotto identifies ligand-receptor pairs whose mean expression is higher than expected. To obtain the observed expression of ligand-receptor pairs for a pair of cell types, Giotto averages the expression of ligand in all sender cells and the expression of receptors in all receiver cells in proximity of the sender cells. Giotto then permutes the location of cells to obtain an expected expression of the ligand-receptor pair. A *p*-value can be obtained by mapping the observed expression onto the distribution of expected expression. Different from Giotto, SpaOTsc uses a partial information decomposition (PID) approach to determine gene-to-gene interaction. Intuitively, PID decomposes the mutual information between multiple source variables and a target variable into unique information contributed by each source variable, redundant information shared by many source variables, and synergistic information contributed by the combination of source variables (Kunert-Graf et al., 2020). SpaOTsc estimates the unique information from a source gene to a target gene that is within a predefined spatial distance, taking into consideration all other genes. Yuan *et al.* proposed a method called GCNG (Yuan and Bar-Joseph, 2020) to infer the extracellular gene relationship using Graph convolutional neural networks (GCN). Single-cell spatial expression data is represented as a graph of cells. Cell locations are encoded as a binary cell adjacency matrix with a selected distance threshold, and expression of gene pairs within each cell is encoded as corresponding node features. A GCN is used to combine the graph structure and node information as input and predict whether the studied gene pair can interact. The deep learning model is trained in a supervised manner, where positive samples are built from known ligand-receptor pairs, and negative samples are randomly selected from non-interacting genes. In addition to the methods specifically designed to utilize the spatial expression information for cell-to-cell interaction, many other tools developed for expression data without spatial information can also be applied to the spatial transcriptomic data. Interested readers can refer to a recent review of these methods (Armingol et al., 2021).

### 6.3 Visualization of Interactions Between Cells and Genes

Cell-to-cell and gene-to-gene interactions are naturally represented as networks and correlation matrices (Figure 4). Integrative packages such as Giotto (Dries et al., 2021b) provide functions to visualize cell-to-cell and gene-to-gene interactions as heatmaps, dot plots, or networks. A heatmap is a visual depiction of a matrix whose values are represented as colored boxes on a grid. With heatmaps, large blocks of highly connected cells or genes can be visually identified. A dot plot is similar to a heatmap, except that the boxes are replaced by dots of varying sizes. A dot plot can use both the size and the color of each dot to represent values in each interaction. Different from heatmaps and dot plots, networks use nodes to represent cells or genes and edges to represent their interactions. The widths and colors of edges can be used to describe the strength of interactions. Besides Giotto, the igraph package is widely adopted for network visualization and provides programming interfaces in R, python, C/C++, and Mathematica (Csardi and Nepusz, 2006). Cytoscape is another widely used package to visualize complex network interaction. Its graphical user interface makes it easy to manipulate and examine nodes and edges in the network (Shannon et al., 2003).

**FIGURE 4**. Identification of cell-cell and receptor-ligand interaction. **(A)** Cell-to-cell interaction can be identified by the correlation of gene expression values between cell pairs. **(B)** Receptor-ligand interaction can be identified by the correlation between genes in interacting cell types.

## 7 Discussions

Spatial transcriptomic technologies have made tremendous progress in recent years. Although earlier technologies are restricted by the number of profiled genes (Chen et al., 2015; Moffitt et al., 2016; Shah et al., 2016) or the spatial resolution (Ståhl et al., 2016), current methods can profile the whole transcriptome at single-cell or subcellular resolution (Liu et al., 2020; Chen et al., 2021; Cho et al., 2021). While available commercialized methods (Visium) cannot achieve cellular resolution, we believe newer technologies will soon be production-ready. As commercial platforms become more affordable, we believe the speed at which spatial transcriptomic datasets become publicly available will only accelerate. For example, phase two of the Brain Initiative Cell Census Network (BICCN) will map the spatial organization of more than 300,000 cells from the mouse’s primary motor cortex (Marx, 2021). Large-scale projects to comprehensively profile spatial gene expression are currently limited, but we envision that these projects will expand in three directions. First, more model organisms will be profiled, enabling comparative analysis of cell types and their spatial organizations across evolution. Second, more organ and tissue types will be profiled for a comprehensive understanding of spatial expression architecture. Third, cell states (e.g., stimulated vs resting) and disease states (cancer vs normal) will be profiled to understand cellular activation and disease pathology.

As spatial transcriptomic datasets become more abundant, meta-analysis across published datasets will become commonplace. Methods to remove batch effects are needed to account for technical confounders across datasets. Unlike bulk and single-cell sequencing, batch effects in spatial transcriptomic data must account for correlation across space. Further, the batch effect may also occur on companion histology images, and methods to jointly analyze histology image and spatial transcriptomic data are required. Although several methods have been developed for batch effect removal in bulk (Leek et al., 2012; Stegle et al., 2012) and single-cell (Korsunsky et al., 2019; Li et al., 2020) sequencing, it is still an under-explored area for spatial transcriptomics.

Histopathology is widely adopted across various domains of medicine and is considered the gold standard for certain diagnoses such as cancer staging (Edge et al., 2010). However, histology is limited by the type and number of cellular features delineated by staining agents. Spatial transcriptomics extends histology to test for both imaging and molecular features and may enable testing for oncogenic driver mutations critical for determining cancer subtypes. A recent method named SpaCell integrates both histology and spatial transcriptomic information to predict cancer staging (Tan et al., 2020). In this method, histological images are tiled into patches, where each patch corresponds to a spatial transcriptomic spot in a tissue. A convolutional neural network is used to extract image features from each patch, and combine the features with the spot gene count. A subsequent deep network is applied to predict the disease stages. We envision that spatial transcriptomics will become a diagnostic routine as it becomes more affordable and the clinical interpretation becomes more streamlined.

In this review, we surveyed state-of-the-art methods for spatial transcriptomic data analysis and visualization, and categorized them into three main categories according to the way their output is visualized. It is unlikely that we covered all available methods for spatial transcriptomics, but we hope this review will serve as a stepping stone and attract more researchers to this field.

## Author Contributions

BL conceived the project. BL and YL performed literature review. BL and YL wrote the paper. LZ made the illustrations.

## Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

## References

Abdelaal, T., Michielsen, L., Cats, D., Hoogduin, D., Mei, H., Reinders, M. J. T., et al. (2019). A Comparison of Automatic Cell Identification Methods for Single-Cell RNA Sequencing Data. *Genome Biol.* 20 (1), 194–219. doi:10.1186/s13059-019-1795-z

Aird, D., Ross, M. G., Chen, W.-S., Danielsson, M., Fennell, T., Russ, C., et al. (2011). Analyzing and Minimizing PCR Amplification Bias in Illumina Sequencing Libraries. *Genome Biol.* 12 (2), R18. doi:10.1186/gb-2011-12-2-r18

Aran, D., Looney, A. P., Liu, L., Wu, E., Fong, V., Hsu, A., et al. (2019). Reference-based Analysis of Lung Single-Cell Sequencing Reveals a Transitional Profibrotic Macrophage. *Nat. Immunol.* 20 (2), 163–172. doi:10.1038/s41590-018-0276-y

Armingol, E., Officer, A., Harismendy, O., and Lewis, N. E. (2021). Deciphering Cell-Cell Interactions and Communication from Gene Expression. *Nat. Rev. Genet.* 22 (2), 71–88. doi:10.1038/s41576-020-00292-x

Asp, M., Giacomello, S., Larsson, L., Wu, C., Fürth, D., Qian, X., et al. (2019). A Spatiotemporal Organ-wide Gene Expression and Cell Atlas of the Developing Human Heart. *Cell* 179 (7), 1647–1660. doi:10.1016/j.cell.2019.11.025

Bacher, R., Chu, L.-F., Leng, N., Gasch, A. P., Thomson, J. A., Stewart, R. M., et al. (2017). SCnorm: Robust Normalization of Single-Cell RNA-Seq Data. *Nat. Methods* 14 (6), 584–586. doi:10.1038/nmeth.4263

Berglund, E., Maaskola, J., Schultz, N., Friedrich, S., Marklund, M., Bergenstråhle, J., et al. (2018). Spatial Maps of Prostate Cancer Transcriptomes Reveal an Unexplored Landscape of Heterogeneity. *Nat. Commun.* 9 (1), 2419. doi:10.1038/s41467-018-04724-5

Blake, A., Kohli, P., and Rother, C. (2011). *Markov Random fields for Vision and Image Processing*. MIT press.

Blondel, V. D., Guillaume, J.-L., Lambiotte, R., and Lefebvre, E. (2008). Fast Unfolding of Communities in Large Networks. *J. Stat. Mech.* 2008 (10), P10008. doi:10.1088/1742-5468/2008/10/p10008

Breslow, N. E., and Clayton, D. G. (1993). Approximate Inference in Generalized Linear Mixed Models. *J. Am. Stat. Assoc.* 88 (421), 9–25. doi:10.1080/01621459.1993.10594284

Butler, A., Hoffman, P., Smibert, P., Papalexi, E., and Satija, R. (2018). Integrating Single-Cell Transcriptomic Data across Different Conditions, Technologies, and Species. *Nat. Biotechnol.* 36 (5), 411–420. doi:10.1038/nbt.4096

Cang, Z., and Nie, Q. (2020). Inferring Spatial and Signaling Relationships between Cells from Single Cell Transcriptomic Data. *Nat. Commun.* 11 (1), 2084–2097. doi:10.1038/s41467-020-15968-5

Chen, A., Liao, S., Cheng, M., Ma, K., Wu, L., Lai, Y., et al. (2021). Large Field of View-Spatially Resolved Transcriptomics at Nanoscale Resolution. *bioRxiv*, 2021. doi:10.1101/2021.01.17.427004

Chen, K. H., Boettiger, A. N., Moffitt, J. R., Wang, S., and Zhuang, X. (2015). RNA Imaging. Spatially Resolved, Highly Multiplexed RNA Profiling in Single Cells. *Science* 348 (6233), aaa6090. doi:10.1126/science.aaa6090

Chen, W.-T., Lu, A., Craessaerts, K., Pavie, B., Sala Frigerio, C., Corthout, N., et al. (2020). Spatial Transcriptomics and *In Situ* Sequencing to Study Alzheimer's Disease. *Cell* 182 (4), 976–991. doi:10.1016/j.cell.2020.06.038

Cho, C.-S., Xi, J., Si, Y., Park, S.-R., Hsu, J.-E., Kim, M., et al. (2021). Microscopic Examination of Spatial Transcriptome Using Seq-Scope. *Cell* 184 (13), 3559–3572. doi:10.1016/j.cell.2021.05.010

Close, J. L., Long, B. R., and Zeng, H. (2021). Spatially Resolved Transcriptomics in Neuroscience. *Nat. Methods* 18 (1), 23–25. doi:10.1038/s41592-020-01040-z

Cortal, A., Martignetti, L., Six, E., and Rausell, A. (2021). Gene Signature Extraction and Cell Identity Recognition at the Single-Cell Level with Cell-ID. *Nat. Biotechnol.* 39 (9), 1095–1124. doi:10.1038/s41587-021-00896-6

Crosetto, N., Bienko, M., and van Oudenaarden, A. (2015). Spatially Resolved Transcriptomics and beyond. *Nat. Rev. Genet.* 16 (1), 57–66. doi:10.1038/nrg3832

Csardi, G., and Nepusz, T. (2006). The Igraph Software Package for Complex Network Research. *InterJournal, complex Syst.* 1695 (5), 1–9.

Dempster, A. P., Laird, N. M., and Rubin, D. B. (1977). Maximum Likelihood from Incomplete Data via theEMAlgorithm. *J. R. Stat. Soc. Ser. B (Methodological)* 39 (1), 1–22. doi:10.1111/j.2517-6161.1977.tb01600.x

Dries, R., Zhu, Q., Dong, R., Eng, C-H. L., Li, H., Liu, K., et al. (2021). Giotto: a Toolbox for Integrative Analysis and Visualization of Spatial Expression Data. *Genome Biol.* 22 (1), 1–31. doi:10.1186/s13059-021-02286-2

Dries, R., Chen, J., Del Rossi, N., Khan, M. M., Sistig, A., and Yuan, G.-C. (2021). Advances in Spatial Transcriptomic Data Analysis. *Genome Res.* 31 (10), 1706–1718. doi:10.1101/gr.275224.121

Dumitrascu, B., Villar, S., Mixon, D. G., and Engelhardt, B. E. (2021). Optimal Marker Gene Selection for Cell Type Discrimination in Single Cell Analyses. *Nat. Commun.* 12 (1), 1186. doi:10.1038/s41467-021-21453-4

Dustin, M. L., and Chan, A. C. (2000). Signaling Takes Shape in the Immune System. *Cell* 103 (2), 283–294. doi:10.1016/s0092-8674(00)00120-3

Edge, S. B., Byrd, D. R., Carducci, M. A., Compton, C. C., Fritz, A., and Greene, F. (2010). *AJCC Cancer Staging Manual*. New York: Springer.

Edsgärd, D., Johnsson, P., and Sandberg, R. (2018). Identification of Spatial Expression Trends in Single-Cell Gene Expression Data. *Nat. Methods* 15 (5), 339–342. doi:10.1038/nmeth.4634

Eng, C.-H. L., Lawson, M., Zhu, Q., Dries, R., Koulena, N., Takei, Y., et al. (2019). Transcriptome-scale Super-resolved Imaging in Tissues by RNA seqFISH+. *Nature* 568 (7751), 235–239. doi:10.1038/s41586-019-1049-y

Finak, G., McDavid, A., Yajima, M., Deng, J., Gersuk, V., Shalek, A. K., et al. (2015). MAST: a Flexible Statistical Framework for Assessing Transcriptional Changes and Characterizing Heterogeneity in Single-Cell RNA Sequencing Data. *Genome Biol.* 16 (1), 278–291. doi:10.1186/s13059-015-0844-5

Fletcher, C. D. (2007). *Diagnostic Histopathology of Tumors: 2-volume Set with CD-ROMs*. Elsevier Health Sciences.

Friedman, C. E., Nguyen, Q., Lukowski, S. W., Helfer, A., Chiu, H. S., Miklas, J., et al. (2018). Single-Cell Transcriptomic Analysis of Cardiac Differentiation from Human PSCs Reveals HOPX-dependent Cardiomyocyte Maturation. *Cell Stem Cell* 23 (4), 586–598. doi:10.1016/j.stem.2018.09.009

Gawad, C., Koh, W., and Quake, S. R. (2016). Single-cell Genome Sequencing: Current State of the Science. *Nat. Rev. Genet.* 17 (3), 175–188. doi:10.1038/nrg.2015.16

Geyer, C. J. (1992). Practical Markov Chain Monte Carlo. *Stat. Sci.*, 473–483. doi:10.1214/ss/1177011137

Gotway, C. A., and Stroup, W. W. (1997). A Generalized Linear Model Approach to Spatial Data Analysis and Prediction. *J. Agric. Biol. Environ. Stat.* 2, 157–178. doi:10.2307/1400401

Hao, Y., Hao, S., Andersen-Nissen, E., Mauck, W. M., Zheng, S., Butler, A., et al. (2021). Integrated Analysis of Multimodal Single-Cell Data. *Cell*. doi:10.1016/j.cell.2021.04.048

Hu, J., Li, X., Coleman, K., Schroeder, A., Ma, N., Irwin, D. J., et al. (2021). SpaGCN: Integrating Gene Expression, Spatial Location and Histology to Identify Spatial Domains and Spatially Variable Genes by Graph Convolutional Network. *Nat. Methods* 18 (11), 1342–1351. doi:10.1038/s41592-021-01255-8

Hu, J., Schroeder, A., Coleman, K., Chen, C., Auerbach, B. J., and Li, M. (2021). Statistical and Machine Learning Methods for Spatially Resolved Transcriptomics with Histology. *Comput. Struct. Biotechnol. J.* 19, 3829–3841. doi:10.1016/j.csbj.2021.06.052

J. Xie, R. Girshick, and A. Farhadi (Editors) (2016). “Unsupervised Deep Embedding for Clustering Analysis,” *International Conference on Machine Learning* (PMLR).

Jégou, H., Douze, M., and Schmid, C. (2010). Improving Bag-Of-Features for Large Scale Image Search. *Int. J. Comput. Vis.* 87 (3), 316–336.

Ji, A. L., Rubin, A. J., Thrane, K., Jiang, S., Reynolds, D. L., Meyers, R. M., et al. (2020). Multimodal Analysis of Composition and Spatial Architecture in Human Squamous Cell Carcinoma. *Cell* 182 (6), 1661–1662. doi:10.1016/j.cell.2020.08.043

Johnson, S. C. (1967). Hierarchical Clustering Schemes. *Psychometrika* 32 (3), 241–254. doi:10.1007/bf02289588

Junker, J. P., Noël, E. S., Guryev, V., Peterson, K. A., Shah, G., Huisken, J., et al. (2014). Genome-wide RNA Tomography in the Zebrafish Embryo. *Cell* 159 (3), 662–675. doi:10.1016/j.cell.2014.09.038

Kanungo, T., Mount, D. M., Netanyahu, N. S., Piatko, C. D., Silverman, R., and Wu, A. Y. (2002). An Efficient K-Means Clustering Algorithm: Analysis and Implementation. *IEEE Trans. Pattern Anal. Machine Intell.* 24 (7), 881–892. doi:10.1109/tpami.2002.1017616

Kim, S. Y., and Volsky, D. J. (2005). PAGE: Parametric Analysis of Gene Set Enrichment. *BMC bioinformatics* 6 (1), 144–156. doi:10.1186/1471-2105-6-144

Kim, T. H., Zhou, X., and Chen, M. (2020). Demystifying "Drop-Outs" in Single-Cell UMI Data. *Genome Biol.* 21 (1), 196. doi:10.1186/s13059-020-02096-y

Kipf, T. N., and Welling, M. (2016). *Semi-supervised Classification with Graph Convolutional Networks* arXiv preprint arXiv:160902907.

Kiselev, V. Y., Kirschner, K., Schaub, M. T., Andrews, T., Yiu, A., Chandra, T., et al. (2017). SC3: Consensus Clustering of Single-Cell RNA-Seq Data. *Nat. Methods* 14 (5), 483–486. doi:10.1038/nmeth.4236

Kiselev, V. Y., Yiu, A., and Hemberg, M. (2018). Scmap: Projection of Single-Cell RNA-Seq Data across Data Sets. *Nat. Methods* 15 (5), 359–362. doi:10.1038/nmeth.4644

Korsunsky, I., Millard, N., Fan, J., Slowikowski, K., Zhang, F., Wei, K., et al. (2019). Fast, Sensitive and Accurate Integration of Single-Cell Data with Harmony. *Nat. Methods* 16 (12), 1289–1296. doi:10.1038/s41592-019-0619-0

Krishnaswamy, S., Spitzer, M. H., Mingueneau, M., Bendall, S. C., Litvin, O., Stone, E., et al. (2014). Conditional Density-Based Analysis of T Cell Signaling in Single-Cell Data. *Science* 346 (6213), 1250689. doi:10.1126/science.1250689

Kunert-Graf, J., Sakhanenko, N., and Galas, D. (2020). Partial Information Decomposition and the Information Delta: A Geometric Unification Disentangling Non-pairwise Information. *Entropy* 22 (12), 1333. doi:10.3390/e22121333

Lee, J. H., Daugharthy, E. R., Scheiman, J., Kalhor, R., Yang, J. L., Ferrante, T. C., et al. (2014). Highly Multiplexed Subcellular RNA Sequencing *In Situ*. *Science* 343 (6177), 1360–1363. doi:10.1126/science.1250212

Leek, J. T., Johnson, W. E., Parker, H. S., Jaffe, A. E., and Storey, J. D. (2012). The Sva Package for Removing Batch Effects and Other Unwanted Variation in High-Throughput Experiments. *Bioinformatics* 28 (6), 882–883. doi:10.1093/bioinformatics/bts034

Li, X., Wang, K., Lyu, Y., Pan, H., Zhang, J., Stambolian, D., et al. (2020). Deep Learning Enables Accurate Clustering with Batch Effect Removal in Single-Cell RNA-Seq Analysis. *Nat. Commun.* 11 (1), 2338–2352. doi:10.1038/s41467-020-15851-3

Littman, R., Hemminger, Z., Foreman, R., Arneson, D., Zhang, G., Gómez-Pinilla, F., et al. (2021). Joint Cell Segmentation and Cell Type Annotation for Spatial Transcriptomics. *Mol. Syst. Biol.* 17 (6), e10108. doi:10.15252/msb.202010108

Liu, Y., Yang, M., Deng, Y., Su, G., Enninful, A., Guo, C. C., et al. (2020). High-Spatial-Resolution Multi-Omics Sequencing via Deterministic Barcoding in Tissue. *Cell* 183 (6), 1665–1681. doi:10.1016/j.cell.2020.10.026

Love, M. I., Huber, W., and Anders, S. (2014). Moderated Estimation of Fold Change and Dispersion for RNA-Seq Data with DESeq2. *Genome Biol.* 15 (12), 550. doi:10.1186/s13059-014-0550-8

Lubeck, E., Coskun, A. F., Zhiyentayev, T., Ahmad, M., and Cai, L. (2014). Single-cell *In Situ* RNA Profiling by Sequential Hybridization. *Nat. Methods* 11 (4), 360–361. doi:10.1038/nmeth.2892

Lui, J. H., Hansen, D. V., and Kriegstein, A. R. (2011). Development and Evolution of the Human Neocortex. *Cell* 146 (1), 18–36. doi:10.1016/j.cell.2011.06.030

Lun, A. T., Bach, K., and Marioni, J. C. (2016). Pooling across Cells to Normalize Single-Cell RNA Sequencing Data with many Zero Counts. *Genome Biol.* 17 (1), 75–14. doi:10.1186/s13059-016-0947-7

Lun, A. T., McCarthy, D. J., and Marioni, J. C. (10002016). A Step-by-step Workflow for Low-Level Analysis of Single-Cell RNA-Seq Data with Bioconductor. *F1000Res* 5, 2122. doi:10.12688/f1000research.9501.2

Lytal, N., Ran, D., and An, L. (2020). Normalization Methods on Single-Cell RNA-Seq Data: an Empirical Survey. *Front. Genet.* 11, 41. doi:10.3389/fgene.2020.00041

Maniatis, S., Äijö, T., Vickovic, S., Braine, C., Kang, K., Mollbrink, A., et al. (2019). Spatiotemporal Dynamics of Molecular Pathology in Amyotrophic Lateral Sclerosis. *Science* 364 (6435), 89–93. doi:10.1126/science.aav9776

Marx, V. (2021). Method of the Year: Spatially Resolved Transcriptomics. *Nat. Methods* 18 (1), 9–14. doi:10.1038/s41592-020-01033-y

McCarthy, D. J., Chen, Y., and Smyth, G. K. (2012). Differential Expression Analysis of Multifactor RNA-Seq Experiments with Respect to Biological Variation. *Nucleic Acids Res.* 40 (10), 4288–4297. doi:10.1093/nar/gks042

McCullagh, P., and Nelder, J. A. (1983). *Generalized Linear Models*. Thames, Oxfordshire, England: Routledge.

McInnes, L., Healy, J., and Melville, J. (2018). *Umap: Uniform Manifold Approximation and Projection for Dimension Reduction*. arXiv preprint arXiv:180203426.

Metzker, M. L. (2010). Sequencing Technologies - the Next Generation. *Nat. Rev. Genet.* 11 (1), 31–46. doi:10.1038/nrg2626

Moffitt, J. R., Bambah-Mukku, D., Eichhorn, S. W., Vaughn, E., Shekhar, K., Perez, J. D., et al. (2018). Molecular, Spatial, and Functional Single-Cell Profiling of the Hypothalamic Preoptic Region. *Science* 362 (6416). doi:10.1126/science.aau5324

Moffitt, J. R., Hao, J., Wang, G., Chen, K. H., Babcock, H. P., and Zhuang, X. (2016). High-throughput Single-Cell Gene-Expression Profiling with Multiplexed Error-Robust Fluorescence *In Situ* Hybridization. *Proc. Natl. Acad. Sci. USA* 113 (39), 11046–11051. doi:10.1073/pnas.1612826113

Moncada, R., Barkley, D., Wagner, F., Chiodin, M., Devlin, J. C., Baron, M., et al. (2020). Integrating Microarray-Based Spatial Transcriptomics and Single-Cell RNA-Seq Reveals Tissue Architecture in Pancreatic Ductal Adenocarcinomas. *Nat. Biotechnol.* 38 (3), 333–342. doi:10.1038/s41587-019-0392-8

Moon, T. K. (1996). The Expectation-Maximization Algorithm. Available from: http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=& arnumber=543975.

Nichterwitz, S., Chen, G., Aguila Benitez, J., Yilmaz, M., Storvall, H., Cao, M., et al. (2016). Laser Capture Microscopy Coupled with Smart-Seq2 for Precise Spatial Transcriptomic Profiling. *Nat. Commun.* 7 (1), 12139–12150. doi:10.1038/ncomms12139

Palla, G., Spitzer, H., Klein, M., Fischer, D., Schaar, A. C., Kuemmerle, L. B., et al. (2021). *Squidpy: A Scalable Framework for Spatial Single Cell Analysis*. Laurel Hollow, New York: bioRxiv.

Rao, A., Barkley, D., França, G. S., and Yanai, I. (2021). Exploring Tissue Architecture Using Spatial Transcriptomics. *Nature* 596 (7871), 211–220. doi:10.1038/s41586-021-03634-9

Reynolds, D. (2009). Gaussian Mixture Models. *Encyclopedia of biometrics* 741, 659–663. doi:10.1007/978-0-387-73003-5_196

Rödelsperger, C., Ebbing, A., Sharma, D. R., Okumura, M., Sommer, R. J., and Korswagen, H. C. (2021). Spatial Transcriptomics of Nematodes Identifies Sperm Cells as a Source of Genomic novelty and Rapid Evolution. *Mol. Biol. Evol.* 38 (1), 229–243. doi:10.1093/molbev/msaa207

Rodriques, S. G., Stickels, R. R., Goeva, A., Martin, C. A., Murray, E., Vanderburg, C. R., et al. (2019). Slide-seq: A Scalable Technology for Measuring Genome-wide Expression at High Spatial Resolution. *Science* 363 (6434), 1463–1467. doi:10.1126/science.aaw1219

Roerdink, J. B., and Meijster, A. (2000). The Watershed Transform: Definitions, Algorithms and Parallelization Strategies. *Fundamenta informaticae* 41 (12), 187–228. doi:10.3233/fi-2000-411207

Saiselet, M., Rodrigues-Vitória, J., Tourneur, A., Craciun, L., Spinette, A., Larsimont, D., et al. (2020). Transcriptional Output, Cell-type Densities, and Normalization in Spatial Transcriptomics. *J. Mol. Cel. Biol.* 12 (11), 906–908. doi:10.1093/jmcb/mjaa028

Schermelleh, L., Ferrand, A., Huser, T., Eggeling, C., Sauer, M., Biehlmaier, O., et al. (2019). Super-resolution Microscopy Demystified. *Nat. Cel Biol* 21 (1), 72–84. doi:10.1038/s41556-018-0251-8

Schwarzacher, T., and Heslop-Harrison, P. (2000). *Practical in Situ Hybridization*. BIOS Scientific Publishers Ltd.

Shah, S., Lubeck, E., Zhou, W., and Cai, L. (2016). *In Situ* Transcription Profiling of Single Cells Reveals Spatial Organization of Cells in the Mouse Hippocampus. *Neuron* 92 (2), 342–357. doi:10.1016/j.neuron.2016.10.001

Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: a Software Environment for Integrated Models of Biomolecular Interaction Networks. *Genome Res.* 13 (11), 2498–2504. doi:10.1101/gr.1239303

Shapiro, E., Biezuner, T., and Linnarsson, S. (2013). Single-cell Sequencing-Based Technologies Will Revolutionize Whole-Organism Science. *Nat. Rev. Genet.* 14 (9), 618–630. doi:10.1038/nrg3542

Song, Q., and Su, J. (2021). DSTG: Deconvoluting Spatial Transcriptomics Data through Graph-Based Artificial Intelligence. *Brief. Bioinform.* 22 (5), bbaa414. doi:10.1093/bib/bbaa414

Ståhl, P. L., Salmén, F., Vickovic, S., Lundmark, A., Navarro, J. F., Magnusson, J., et al. (2016). Visualization and Analysis of Gene Expression in Tissue Sections by Spatial Transcriptomics. *Science* 353 (6294), 78–82. doi:10.1126/science.aaf2403

Stegle, O., Parts, L., Piipari, M., Winn, J., and Durbin, R. (2012). Using Probabilistic Estimation of Expression Residuals (PEER) to Obtain Increased Power and Interpretability of Gene Expression Analyses. *Nat. Protoc.* 7 (3), 500–507. doi:10.1038/nprot.2011.457

Stickels, R. R., Murray, E., Kumar, P., Li, J., Marshall, J. L., Di Bella, D. J., et al. (2021). Highly Sensitive Spatial Transcriptomics at Near-Cellular Resolution with Slide-seqV2. *Nat. Biotechnol.* 39 (3), 313–319. doi:10.1038/s41587-020-0739-1

Subramanian, A., Tamayo, P., Mootha, V. K., Mukherjee, S., Ebert, B. L., Gillette, M. A., et al. (2005). Gene Set Enrichment Analysis: A Knowledge-Based Approach for Interpreting Genome-wide Expression Profiles. *Proc. Natl. Acad. Sci.* 102 (43), 15545–15550. doi:10.1073/pnas.0506580102

Sun, S., Zhu, J., and Zhou, X. (2020). Statistical Analysis of Spatial Expression Patterns for Spatially Resolved Transcriptomic Studies. *Nat. Methods* 17 (2), 193–200. doi:10.1038/s41592-019-0701-7

Svensson, V., Teichmann, S. A., and Stegle, O. (2018). SpatialDE: Identification of Spatially Variable Genes. *Nat. Methods* 15 (5), 343–346. doi:10.1038/nmeth.4636

Taguchi, A., and White, M. F. (2008). Insulin-like Signaling, Nutrient Homeostasis, and Life Span. *Annu. Rev. Physiol.* 70, 191–212. doi:10.1146/annurev.physiol.70.113006.100533

Taipale, J., and Beachy, P. A. (2001). The Hedgehog and Wnt Signalling Pathways in Cancer. *Nature* 411 (6835), 349–354. doi:10.1038/35077219

Tan, X., Su, A., Tran, M., and Nguyen, Q. (2020). SpaCell: Integrating Tissue Morphology and Spatial Gene Expression to Predict Disease Cells. *Bioinformatics* 36 (7), 2293–2294. doi:10.1093/bioinformatics/btz914

Thrane, K., Eriksson, H., Maaskola, J., Hansson, J., and Lundeberg, J. (2018). Spatially Resolved Transcriptomics Enables Dissection of Genetic Heterogeneity in Stage III Cutaneous Malignant Melanoma. *Cancer Res.* 78 (20), 5970–5979. doi:10.1158/0008-5472.CAN-18-0747

Traag, V. A., Waltman, L., and van Eck, N. J. (2019). From Louvain to Leiden: Guaranteeing Well-Connected Communities. *Sci. Rep.* 9 (1), 5233. doi:10.1038/s41598-019-41695-z

Tsuyuzaki, K., Sato, H., Sato, K., and Nikaido, I. (2020). Benchmarking Principal Component Analysis for Large-Scale Single-Cell RNA-Sequencing. *Genome Biol.* 21 (1), 9. doi:10.1186/s13059-019-1900-3

van den Heuvel, M., Nusse, R., Johnston, P., and Lawrence, P. A. (1989). Distribution of the Wingless Gene Product in Drosophila Embryos: a Protein Involved in Cell-Cell Communication. *Cell* 59 (4), 739–749. doi:10.1016/0092-8674(89)90020-2

Van der Maaten, L., and Hinton, G. (2008). Visualizing Data Using T-SNE. *J. machine Learn. Res.* 9 (11).

Vickovic, S., Eraslan, G., Salmén, F., Klughammer, J., Stenbeck, L., Schapiro, D., et al. (2019). High-definition Spatial Transcriptomics for *In Situ* Tissue Profiling. *Nat. Methods* 16 (10), 987–990. doi:10.1038/s41592-019-0548-y

Waltman, L., and Van Eck, N. J. (2013). A Smart Local Moving Algorithm for Large-Scale Modularity-Based Community Detection. *The Eur. Phys. J. B* 86 (11), 1–14. doi:10.1140/epjb/e2013-40829-0

Wang, B., Zhu, J., Pierson, E., Ramazzotti, D., and Batzoglou, S. (2017). Visualization and Analysis of Single-Cell RNA-Seq Data by Kernel-Based Similarity Learning. *Nat. Methods* 14 (4), 414–416. doi:10.1038/nmeth.4207

Wang, J. (2020). *An Intuitive Tutorial to Gaussian Processes Regression*. arXiv preprint arXiv:200910862.

Wang, X., Allen, W. E., Wright, M. A., Sylwestrak, E. L., Samusik, N., Vesuna, S., et al. (2018). Three-dimensional Intact-Tissue Sequencing of Single-Cell Transcriptional States. *Science* 361 (6400). doi:10.1126/science.aat5691

Wedderburn, R. W. M. (1974). Quasi-Likelihood Functions, Generalized Linear Models, and the Gauss-Newton Method. *Biometrika* 61 (3), 439–447. doi:10.2307/2334725

Wei, C.-J., Xu, X., and Lo, C. W. (2004). Connexins and Cell Signaling in Development and Disease. *Annu. Rev. Cel Dev. Biol.* 20, 811–838. doi:10.1146/annurev.cellbio.19.111301.144309

Wirka, R. C., Wagh, D., Paik, D. T., Pjanic, M., Nguyen, T., Miller, C. L., et al. (2019). Atheroprotective Roles of Smooth Muscle Cell Phenotypic Modulation and the TCF21 Disease Gene as Revealed by Single-Cell Analysis. *Nat. Med.* 25 (8), 1280–1289. doi:10.1038/s41591-019-0512-5

Wold, S., Esbensen, K., and Geladi, P. (1987). Principal Component Analysis. *Chemometrics Intell. Lab. Syst.* 2 (1-3), 37–52. doi:10.1016/0169-7439(87)80084-9

Wolf, F. A., Angerer, P., and Theis, F. J. (2018). SCANPY: Large-Scale Single-Cell Gene Expression Data Analysis. *Genome Biol.* 19 (1), 15. doi:10.1186/s13059-017-1382-0

Yuan, Y., and Bar-Joseph, Z. (2020). GCNG: Graph Convolutional Networks for Inferring Gene Interaction from Spatial Transcriptomics Data. *Genome Biol.* 21 (1), 300–316. doi:10.1186/s13059-020-02214-w

Yuste, R., Hawrylycz, M., Aalling, N., Aguilar-Valles, A., Arendt, D., Armañanzas, R., et al. (2020). A Community-Based Transcriptomics Classification and Nomenclature of Neocortical Cell Types. *Nat. Neurosci.* 23 (12), 1456–1468. doi:10.1038/s41593-020-0685-8

Zhao, E., Stone, M. R., Ren, X., Guenthoer, J., Smythe, K. S., Pulliam, T., et al. (2021). Spatial Transcriptomics at Subspot Resolution with BayesSpace. *Nat. Biotechnol.* 39 (11), 1375–1384. doi:10.1038/s41587-021-00935-2

Keywords: spatial transcriptomics, single-cell RNA-seq (scRNA-seq), clustering, cell-type identification, dimensionality reduction, spatial expression pattern, spatial interaction, visualization

Citation: Liu B, Li Y and Zhang L (2022) Analysis and Visualization of Spatial Transcriptomic Data. *Front. Genet.* 12:785290. doi: 10.3389/fgene.2021.785290

Received: 29 September 2021; Accepted: 24 December 2021;

Published: 27 January 2022.

Edited by:

Guangchuang Yu, Southern Medical University, ChinaReviewed by:

Xin Li, Shanghai Institute of Nutrition and Health (CAS), ChinaLu Zhang, Hong Kong Baptist University, Hong Kong SAR, China

Copyright © 2022 Liu, Li and Zhang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Boxiang Liu, jollier.liu@gmail.com

^{†}These authors have contributed equally to this work