Given any spatial reference atlas consisting of binary or continuous gene expression levels for a biological system on a set of locations with known spatial coordinates, and a scRNA-seq dataset consisting of binary or continuous gene expression levels for the same biological system, we introduce a Deep-learning based Environment for the Extraction of Positional information from scRNA-seq data (DEEPsc) to impute the spatial information onto the scRNA-seq data.

In DEEPsc, we first select a common set of genes from the reference atlas and scRNA-seq data, then perform dimensionality reduction via principal component analysis (PCA) on the reduced reference atlas to shorten training time (Figure 1A). The scRNA-seq data is then projected into the same PCA space on which we learn a metric for comparison between cells and spatial positions. The DEEPsc network accepts a concatenated feature vector for a single cell and a single position and returns a likelihood the input cell originated from the input position. The network contains two fully connected hidden layers with N nodes each, where N is the number of principal components kept from PCA, and a single node in the output layer. Sigmoid activation functions are applied to each node, including the output node, so that the resulting output is in [0,1] and can be interpreted as a likelihood that the input cell originated from the input spatial position. To train the DEEPsc network, we use the spatial position feature vectors as simulated scRNA-seq data for comparison (Figure 1B). Each simulated cell is compared pairwise with every position in the spatial reference atlas; if the simulated cell is an exact match to a given position, the target output is 1 (a high likelihood of origin), and if the simulated cell and chosen position are not an exact match, the target output is 0 (a low likelihood of origin). Training is terminated when the error on a randomly chosen validation set is no longer improving.

After training the DEEPsc network, a feature vector associated with an actual cell from the scRNA-seq data is fed in as input and compared to each position in the reference atlas individually. We display the results as a heatmap on the schematic diagram of the biological system, choosing the spatial position with the largest likelihood of origin according to DEEPsc as the determined origin of the cell. This process is repeated for each cell in the scRNA-seq dataset to assign spatial origins of all cells (Figure 1C).

Each of the highlighted methods to impute spatial data onto scRNA-seq data, including DEEPsc, can be essentially boiled down to the following: For some tissue with a well-defined standard spatial structure, given known binary or continuous expression levels of G genes at each of P spatial locations (the reference atlas), calculate a correspondence score, S, of how similar each of C cells in an scRNA-seq dataset is to each of the P positions in the atlas. That is, define a function, S:[0,1]^G×[0,1]^G→[0,1], such that S(c_i,p_j);i = 1,2,…,C;j = 1,2,…,P; which describes the likelihood that cell c_i originated from position p_j, based on the similarity of the expression vectors of the cell and position.

To quantify how well a given method performs for a given spatial reference atlas, we use the reference atlas itself as simulated single cell data; that is, we generate a simulated scRNA-seq dataset with C=P cells, each an exact copy of a reference atlas position. This allows us to treat the simulated scRNA-seq data as having a known spatial origin, against which we can compare the output of each method. We define a system-adaptive, comprehensive performance score, consisting of three penalty terms: accuracy, which determines whether or not the known spatial origin was given a high likelihood of origin; precision, which determines whether or not other locations were given low likelihoods of origin; and robustness, which determines how sensitive a mapping method is to random noise in the input data. Each penalty term is represented by a number in [0,1], with 0 being no penalty and 1 being a worst-case scenario. The performance score is defined as E=1P⁢∑i=1Ei, where

S_i,j = S(c_i,p_j) is the correspondence score of cell c_i to position p_j, and E_i is interpreted as the error in the mapping of cell c_i. The quantity σ^∗ in the robustness term is calculated by determining the accuracy and precision penalty terms with no Gaussian noise added to the input data, then calculating the same two penalties with various levels of Gaussian noise with standard deviation σ ∈ [0,1]. The quantity σ^∗ is defined to be the level of Gaussian noise required to raise the mean of the accuracy and precision penalties by 0.1 from their values with no added noise, or σ^∗ = 1, whichever is smallest. The exponent of four in the robustness term was chosen empirically such that the robustness term does not dominate the performance score, keeping in mind that expression levels are normalized to [0,1] before calculating the correspondence scores, so e.g., σ^∗ = 0.5 means a method requires noise on the order of half of the expression levels to raise the precision and accuracy penalties by 0.1. The performance score has a range of [0,1], where a performance score of E=1 represents an ideal mapping that maps a cell to its known location with high confidence, to all other locations with low confidence, and does so in a manner robust to noise. An illustration of each term in the performance score is shown in Figure 2.

This performance score is limited by the fact that it relies on ground truth knowledge of the spatial origin of a single cell/spot to determine the performance of a given mapping method. However, this ground truth knowledge is not available for actual scRNA-seq data. To more directly quantify the mapping performance on actual scRNA-seq datasets, we use a measure of predictive reproducibility, obtained from a k-fold cross validation scheme, in which we randomly split the common genes in the reference atlas and scRNA-seq data into k folds and calculate the correspondence score for each method using all but one fold. The correspondence scores are then used as coefficients in a weighted sum to predict the value of the dropped-out genes in each fold for each cell (scRNA-seq predictive reproducibility) or each spatial position (atlas predictive reproducibility) and determine the error in the predicted expression level. The predicted expression of gene k in cell c_i is computed as c^i(k)=∑j=1PSi,j(k)⁢pj(k)/∑j=1PSi,j(k) and the predicted expression of gene k in position p_j is computed as p^j(k)=∑i=1CSi,j(k)⁢ci(k)/∑i=1CSi,j(k) where Si,j(k) is the correspondence score between cell c_i and position p_j with genes in folds not containing gene k and ci(k) and pj(k) are the known expression values of gene k from the scRNA-seq and the spatial atlas data, respectively. To accommodate the sparsity of data, we compute the predictive reproducibility scores separately for cells or positions with zero expression values and with positive expression values. For example, we measure the predictive reproducibility for the task of reproducing gene k in scRNA-seq data on cells with zero expression using Rs⁢c⁢_⁢z⁢e⁢r⁢o(k)=1-∑i∈Is⁢c⁢_⁢z⁢e⁢r⁢o(k)|c^i(k)-ci(k)|/|Is⁢c⁢_⁢z⁢e⁢r⁢o(k)| where Is⁢c⁢_⁢z⁢e⁢r⁢o(k)={i:ci(k)=0}. Taking the average over all common genes results in a single score R_{sc_ zero}, and in the same manner, we define R_{sc_ nonzero}, R_{atlas_ zero}, and R_{atlas_ nonzero}. When producing predictive reproducibility scores, we use the same k-fold split across all methods to ensure a fair comparison.

Using the performance score, we benchmarked the methods developed by Achim et al. (2015) and Satija et al. (2015) (Seurat v1), Karaiskos et al. (2017) (DistMap), and Peng et al. (2016) together with our DEEPsc method and applied them to four different biological systems: the zebrafish embryo (Satija et al., 2015), the Drosophila embryo (Karaiskos et al., 2017), the murine hair follicle (Joost et al., 2016), and the murine frontal cortex, downloaded from the 10x Genomics Spatial Gene Expression Datasets. The reference atlas for the zebrafish embryo contains the binarized expression of 47 genes on 64 spatial bins that assemble half of the hemisphere of the 6hpf embryo (Satija et al., 2015). The Drosophila embryo reference atlas contains 84 genes on 3,039 spatial positions (Karaiskos et al., 2017). The spatial reference atlas generated with the Visium technology (Ståhl et al., 2016) for the murine frontal cortex contains 32,285 genes on 961 spatial positions (a subset presenting the frontal cortex from the original data), from which we kept 2755 genes from the 3,000 most variable genes in spatial data that are also present in scRNA-seq data. Segmenting a standard diagram of the follicle into 233 spatial positions and using FISH imaging of eight genes identified as spatially localized (Joost et al., 2016), we manually defined a continuous reference atlas for the follicle (section “Materials and Methods”). For mapping methods requiring a binary reference atlas, we defined a cutoff expression of 0.2 to be considered on in this follicle reference atlas of follicle. We further implemented several baseline methods for benchmark comparisons, including several methods using predefined metrics where the correspondence score S is defined to be the 2-norm, infinity norm, or mean percent difference in the space of common genes between the input cell and spatial position. We also implemented a large margin nearest neighbor (LMNN) method that learns a linear metric (section “Materials and Methods”). Figure 3 shows a scatter plot of the penalty terms constituting the performance score of each implemented method on each of the four biological systems, as well as the average for each method across all four systems. Table 1 includes the numerical values for each penalty term, as well as the calculated performance score for each method. Figure 4 includes example heatmaps of simulated cells for each of the biological systems. The penalty terms for the individual locations are shown in Figure 5.

The majority of methods were able to project the simulated scRNA-seq cells to their known spatial origins with high accuracy. Specifically, Seurat v1 and DistMap achieve high performance scores in the zebrafish embryo and Drosophila embryo datasets that they were originally applied to, respectively. Designed to be a system-adaptive method, DEEPsc has the best average performance score across the four datasets (Table 1). Moreover, while some methods are stronger in terms robustness or precision, DEEPsc attains a balance between robustness and precision (Figure 3). This balance is especially important when investigating the impact of cellular spatial neighborhood on cell fate acquisition. It is desired to narrow down the inferred spatial neighborhood (precision) and at the same time reduce the sensitivity to noise (robustness). The high precision and robustness of DEEPsc is consistently observed across all locations in the dataset (Figure 5). Finally, it is worth mentioning that DEEPsc has a significantly higher robustness in the follicle dataset which has the smallest number of genes and is the noisiest among the four datasets.

We now map actual scRNA-seq data for each system and calculate the predictive reproducibility for each method (Table 2 and Figure 6). For the follicle, the scRNA-seq data contains 1,422 cells with 26,024 genes measured containing the eight genes in the spatial atlas (Joost et al., 2016). For the Drosophila embryo, we used the scRNA-seq dataset with 1,297 cells and 8,924 genes among which all the 84 spatial genes are present (Karaiskos et al., 2017). For the Zebrafish embryo, there are 1,152 cells and 11,978 genes in the scRNA-seq dataset with all the 47 spatial genes included (Satija et al., 2015). For the murine frontal cortex, we used the scRNA-seq dataset provided by the Allen Institute (Tasic et al., 2016), generated with SMART-Seq2, which contains 14,249 cells and 34,617 genes, from which a set of 2,755 genes were found to be present in the top 3,000 highly variable genes in spatial atlas. These four datasets cover different situations. The follicle data contains a moderate number of locations, and the cells form well-defined layered structures such that there could be long and thin spatial regions that contain the same cells. The zebrafish embryo spatial data has a suboptimal resolution such that each spatial location consists of multiple cells. This data helps to evaluate the methods in treating coarse spatial atlases. The Drosophila embryo data contains rich spatial characteristics. There is a well-defined global ventral-dorsal/anterior-posterior coordinate system. Locally, there is also a striped pattern in the lateral side of the embryo. The frontal cortex data examines spatial gene expression at the transcriptomics level, and functions as a demonstration that DEEPsc is able to maintain a high performance on high-dimensional datasets.

For the baseline models, we linearly normalized each gene in the log-normalized scRNA-seq dataset onto the interval [0,1]. Continuous spatial atlases with expression values in the [0,1] range were used for the follicle, Drosophila embryo, and murine frontal cortex systems, the latter two having been linearly normalized to [0,1] in the same fashion as the scRNA-seq data. Since a continuous spatial atlas for Zebrafish embryo is lacking, we applied a spatial convolution to the binary atlas and added a small amount of Gaussian noise to simulate a continuous atlas. The 2-norm, Inf-norm, percent difference, and LMNN baseline models are then applied to the vectors of the commonly expressed genes in the spatial atlas and scRNA-seq data. For DEEPsc, we first applied a PCA reduction to the spatial atlas, and then applied the same linear transformation to the normalized expression values of the common genes in the scRNA-seq data. The feature vectors for the locations in the spatial atlas and the cells in the scRNA-seq data in the PCA space were then fed to the neural network. For the four existing methods, we followed the procedure as described in the associated original publications, scaling the resulting correspondence scores to [0,1] for direct comparison with baseline methods. For all the methods, we compute the predictive reproducibility by iterating over all common genes, attempting to reconstruct the expression of one gene using the k-fold cross validation scheme described in the previous section. We used k=4 for the follicle and Drosophila embryo dataset, and k=5 for the zebrafish embryo and cortex dataset.

DEEPsc has a comparable accuracy compared to other methods, and it also has a consistent performance across different systems (Table 2 and Figure 6). This consistent performance further demonstrates the system-adaptive advantage of DEEPsc and the benefit of using adaptive metrics over predefined ones. We also notice that similar to the simulated case, DEEPsc also achieves a balance between precision and robustness in the case of real scRNA-seq data. For example, while it exhibits high precision by mapping the example cell to a specific local spot in the Zebrafish embryo or a local strip in Drosophila embryo, it also robustly maps a cell to the entire outer bulge of the follicle instead of only part of it (Figure 7). The high precision ensures that we can resolve the heterogeneity in the spatial environment and further relate them to the heterogeneity in cell fates. The high robustness prevents the identification of false correlations. Overall, DEEPsc achieves a high predictive reproducibility across all cells in the scRNA-seq dataset.

Dimension reduction is a crucial initial step of DEEPsc. A dimension reduction method that can be trained on one dataset and deterministically applied to another is needed due to the separated training and predicting steps. Here, we explore two different representative dimension reduction methods in the linear and nonlinear categories, PCA and Uniform Manifold Approximation and Projection (UMAP; McInnes et al., 2018). To compare these two methods, we trained several networks with varying amounts of added noise on the reference atlases of the four studied biological systems (Figure 8). We compared PCA (8 principal components), UMAP30 (n_components = 8, n_neighbors = 30), and UMAP5 (n_components = 8, n_neighbors = 5). While on the follicle system all three reduction methods performed virtually identically, on all three other systems PCA outperformed the other reduction methods by achieving a higher robustness score while maintaining similar accuracy.

Given a matrix of scRNA-seq read counts where each row is a different gene and each column is a different cell, and a matrix representing a spatial reference atlas where each row is a different gene and each column is a different spatial position, we first select common genes by eliminating rows in each corresponding to genes not in the other matrix. Once we have eliminated genes not in common, we are left with a number of cells (C) × number of genes (G) matrix for the scRNA-seq data and a number of positions (P) × number of genes (G) matrix for the spatial reference atlas.

We then apply dimensionality reduction to the atlas in the form of a PCA projection, selecting a user-configurable number of principal components to serve as feature vectors. We find in our analysis that keeping the top eight principal components yields satisfactory results. The same PCA coefficients are used to project the scRNA-seq matrix into these principal components. After projection, both matrices are normalized by dividing by the largest element in each, so that the elements are all in [0,1].

For the comparisons in section “Comparison of Dimensionality Reduction Methods,” we use the UMAP implementation by Meehan et al. (2021), found on the MATLAB Central File Exchange at https://www.mathworks.com/matlabcentral/fileexchange/71902. Specifically, we ran the run_umap() function on the spatial reference atlas with n_dimensions = 8 and n_neighbors = 30 or n_neighbors = 5 for UMAP30 and UMAP5, respectively.

To train the DEEPsc network, we use the spatial position feature vectors themselves as simulated scRNA-seq data. The training data is a set of P² vectors of length 2N, where N is the reduced dimensionality of the reference atlas. The first N components correspond to a feature vector of one position in the reference atlas (functioning as a simulated cell) and the last N components correspond to some other position in the reference atlas. Each simulated cell is compared pairwise with every position in the spatial reference atlas; if the simulated cell is an exact match to a given position, the target output is chosen to be 1 (a high likelihood of origin), and if the simulated cell and chosen position are not an exact match, the target output is chosen to be 0 (a low likelihood of origin).

The DEEPsc architecture is an artificial neural network with 2N inputs, two fully connected hidden layers with N nodes each and a single node in the output layer. Sigmoid activation functions are attached to each node, including the output node, so that the resulting output is in [0,1] and can be interpreted as a likelihood that the input cell originated from the input spatial position. To preserve robustness and avoid overfitting the training data, a layer of Gaussian noise is added to the simulated cells so that the network is pushed to learn complex nonlinear relationships among the spatial positions in the reference atlas rather than simply activate when an exact match is encountered. This Gaussian noise layer allows the user to configure the standard deviation of the added noise, as well as to configure the probability that any noise will be added in a given training epoch. We find empirically that a noise level of about 0.10 and a probability of 0.5 yield reasonable robustness to noise, though this may vary from system to system.

Since the training data will naturally consist of many more non-matches than matches, and thus the target data will contain many more zeros than ones, we use a novel custom objective function,

where y_i is the network’s predicted output and t_i is the target output (t_i = 1 if exact match and t_i = 0 if not), to more heavily penalize the network when it gives a false negative (low likelihood when it should be high) than when it gives a false positive (high likelihood when it should be low). This acts to counteract the tendency of the network to “learn” to simply return 0 for every single input and “ignore” any comparably rare training data with t_i = 1.

To further account for the sparsity of exact matches in the training set, we split it into a test and validation set, the former consisting of a configurable fraction of the inputs corresponding to exact matches as well as a configurable multiple of the inputs corresponding to non-matches. If trainFrac = 0.9 and trainingMultiple = 99, for example, 90% of the exact matches will be added to the training set and 99x more non-matches will be added, so that the exact matches make up 1% of the training set. The rest of the inputs are reserved for the (generally much larger) validation set. This is beneficial in reducing training time because it allows us to train with a much smaller fraction of the P² input vectors, giving preference to the exact matches. Indeed, this reduces the size of the actual training set to scale linearly with the size of the atlas rather than quadratically.

Training is performed in MATLAB using the trainNetwork() function in the Deep Learning Toolbox (The Mathworks, Inc, 2019a), for which we implemented the above-described custom network layers. Since the input data is already normalized in preprocessing, we disable the default normalization of trainNetwork(). We use the default Glorot (Xavier and Yoshua, 2010) initialization of weights and biases in the fully connected layers. We then train each network for a maximum of 50,000 epochs of standard gradient descent with a learning rate of η=0.01, shuffle the order of the data each epoch, and use the ADAM optimization method (Kingma and Ba, 2014) with the default parameters β₁ = 0.9, β₂ = 0.999, and ε = 10⁻⁸. In addition to the custom objective function layer we describe above, trainNetwork() by default adds an L²-regularization term to the loss with a regularization factor of λ=0.0001. We monitor the RMSE of the validation set throughout training and manually stop training if it is no longer improving before the maximum number of epochs has been reached. The trainNetwork() function also allows for parallel computation via the Parallel Computing Toolbox (The Mathworks, Inc, 2019b), which is highly recommended but not strictly required for training.

To create a spatial reference atlas for the murine follicle system, we patterned the spatial coordinates of each position in the atlas off of a standard diagram of a mouse follicle found in Figure 1 of Joost et al. (2016). We constructed a Voronoi diagram around each of the cell centers and made manual adjustments to the vertices as we saw fit aesthetically. We then selected the eight genes in the atlas from the systematic staining catalog made available by Joost. We chose the genes based on a combination of high image quality and spatial diversity. Gene expression levels in [0,1] were chosen manually to best represent the images, though to eliminate any implicit bias we also added a small level of Gaussian noise to the atlas. For all methods requiring a binary atlas, we chose a cutoff of 0.2 to represent “on” expression in this atlas.

To implement a LMNN baseline for benchmarking comparison, we used code from the MATLAB Toolbox for Dimensionality Reduction found at https://lvdmaaten.github.io/drtoolbox/ and modified it for our uses. Specifically, we used the lmnn() function in the “techniques” subfolder, and modified the code to set mu = 1, i.e., to remove the “pull” term, as well as setting the number of targets to 1 (the point itself) and treating all other points as imposters. Further, we modified the slack variables to enforce a minimum separation of D, where D is the dimensionality of the space (D=G for our applications). For the numerical experiments of the LMNN method with the cortex dataset, a PCA dimension reduction (50 PCs) was performed before applying LMNN to accommodate the large number of genes.