Imputation of spatial gene expressions by tensor modeling

Let $$ \mathcal{T}\in {\mathbb{R}}_{+}^{n_p\times n_y\times n_x}$$ be the 3-way sparse tensor of the observed spatial gene expression data as show in Fig 1A, with the missing gene expressions represented as zeros, where n_p denote the total number of genes, n_x and n_y denote the dimensions of the x and y spatial coordinates of the spatial transcriptomics array. Our goal is to learn a complete spatial gene expression tensor $$ \widehat{\mathcal{T}}\in {\mathbb{R}}_{+}^{n_p\times n_y\times n_x}$$ from $$ \mathcal{T}$$ as illustrated in Fig 1C. However, it is computationally expensive and often infeasible to compute or store a dense tensor $$ \widehat{\mathcal{T}}$$, especially in high spatial resolutions with millions of spots. Therefore, we propose to compute an economy-size representation of $$ \widehat{\mathcal{T}}$$ via an equality constraint $$ \widehat{\mathcal{T}}=〚{\widehat{A}}_p,{\widehat{A}}_y,{\widehat{A}}_x〛$$, which is called Canonical Polyadic Decomposition (CPD) [28] of $$ \widehat{\mathcal{T}}$$ defined below

Graph regularized tensor completion model

The key ideas of modeling the task of spatial gene expression imputation are i) the inferred complete spatial gene expression tensor $$ \widehat{\mathcal{T}}$$ is regularized to integrate the spatial arrangements of the spots in the tissue array and the functional relations among the genes; ii) the observed part in $$ \mathcal{T}$$ is also required to be preserved in $$ \widehat{\mathcal{T}}$$ as the completion task requires; and iii) the inferred tensor $$ \mathcal{T}$$ is compressed as the economy-size representation $$ \widehat{\mathcal{T}}=〚{\widehat{A}}_p,{\widehat{A}}_y,{\widehat{A}}_x〛$$ for scalable space and time efficiencies. The novel optimization formulation is shown below in Proposition 1,

Proposition 1. The complete spatial gene expression tensor $$ \widehat{\mathcal{T}}\in {\mathbb{R}}^{n_p\times n_y\times n_x}$$ can be obtained by solving the following optimization problem:

There are two optimization terms in the model defined in Eq (2), consistency with the observations (the first term) and Cartesian product graph regularization (the second term), which are explained below,

Consistency with the observations

We introduce a binary mask tensor $$ \mathcal{M}$$ to indicate the indices of the observed entries in $$ \mathcal{T}$$. The (i, j, k)-th entry $$ {\mathcal{M}}_{i,j,k}$$ which is defined below, represents whether the (i, j, k)-th element in $$ \mathcal{T}$$ is observed or not.


By introducing the squared-error in $$ \mathcal{F}$$-norm $$ \left|\right|\mathcal{M}⊛(\mathcal{T}-\widehat{\mathcal{T}}){\left|\right|}_{\mathcal{F}}^2$$ in the model, we ensure the inferred spatial gene expression tensor $$ \widehat{\mathcal{T}}$$ is consistent with its observed counterparts in $$ \mathcal{T}$$.

Cartesian product graph regularization

Two useful assumptions to introduce prior knowledge for inferring the tensor are 1) the spatially adjacent spots should share similar gene expressions, and 2) the expressions of two genes are likely highly correlated if they share similar gene functions [29, 30]. We introduce a spatial graph and a protein-protein interaction (PPI) network into the model.

We first encode the spatial information in two undirected unweighted chain graphs G_x = (V_x, E_x) and G_y = (V_y, E_y), where V_x and V_y are vertex sets and E_x and E_y are edge sets. There are |V_x| = n_x vertices in G_x where n_x is the number of the spatial coordinates along the x-axis of the spatial array. Two vertices in G_x are connected by an edge if they are adjacent along the x-axis. The connections in G_y can be similarly defined to encode the y-coordinates of the tissue.

We also incorporate the topological information of a PPI network downloaded from BioGRID 3.5 [31] to use the functional modules in the PPI network. We denote the PPI network as G_p = (V_p, E_p) which contains |E_p| experimentally documented physical interactions among the |V_p| = n_p proteins. We then use the Cartesian product graph [32] $$ \mathfrak{G}(x,y,p)=(V,E)$$ of the three individual graphs G_x, G_y and G_p to regularize the elements in $$ \widehat{\mathcal{T}}$$, where |V| = n_x n_y n_p. The (v_x v_y v_p)-th vertex in V represents a triple of vertices {v_x ∈ V_x, v_y ∈ V_y, v_p ∈ V_p} from each of the three graphs. The (a_x, a_y, a_p)-th and (b_x, b_y, b_p)-th vertices in V are connected by an edge iff for any i, j ∈ {x, y, p}, (a_i, b_i) ∈ E_i and a_j = b_j ∈ V_j for all j ≠ i. For a graph G_i = (V_i, E_i) where i ∈ {x, y, p}, we denote its adjacency matrix as W_i, degree matrix as D_i, and graph Laplacian matrix as L_i = D_i − W_i. The adjacency and graph Laplacian matrices of $$ \mathfrak{G}(x,y,p)$$ are obtained as $$ \mathfrak{W}(x,y,p)=W_x\oplus W_y\oplus W_p$$ and $$ \mathfrak{L}(x,y,p)=L_x\oplus L_y\oplus L_p$$ respectively, where ⊕ denotes the Kronecker sum [33].

By introducing the term $$ \mathbf{\text{vec}}{\left(\widehat{\mathcal{T}}\right)}^T\mathfrak{L}(x,y,p)\mathbf{\text{vec}}\left(\widehat{\mathcal{T}}\right)$$ in Eq (2), the inferred gene expression values in $$ \widehat{\mathcal{T}}$$ are ensured to be smooth over the manifolds of the product graph $$ \mathfrak{G}(x,y,p)$$, such that a pair of tensor entries $$ {\widehat{\mathcal{T}}}_{a_p,a_y,a_x}$$ and $$ {\widehat{\mathcal{T}}}_{b_p,b_y,b_x}$$ share similar values if the (a_x, a_y, a_p)-th and (b_x, b_y, b_p)-th vertices are connected in $$ \mathfrak{G}(x,y,p)$$. A connection implies that the x-coordinate a_x and b_x is adjacent or y-coordinate a_y and b_y is adjacent or gene a_p and gene b_p are connected in the PPI, with the two other dimensions fixed. Using Cartesian product graph is a more conservative strategy to connect multi-relations in a high-order graph as we have shown in [34] since only replacing one of the dimensions by the immediate neighbors is allowed to create connections. Note that it also possible to use tensor product graph or strong product graph [34] but there could be too many connections to provide meaningful connectivity in the product graph for helpful regularization. It is known that genes’ connectivities in PPI network correlate with their co-expressions. We also justified this hypothesis on the spatial transcriptomics data by examining the relation between the connectivity in PPI network and the co-expression in spatial locations among the genes in the 10 different 10x Genomics Visium spatial genomics datasets used in this study. The results of this analysis are shown in S2 Fig. We observed higher co-expressions between the genes that are connected with less hops in the PPI, which clearly supports our assumptions.

FIST Algorithm

The model introduced in Eq (2) is non-convex on variables $$ \{{\widehat{A}}_p,{\widehat{A}}_y,{\widehat{A}}_x\}$$ jointly, thus finding its global minimum is difficult. In this section, we propose an efficient iterative algorithm Fast Imputation of Spatially-resolved transcriptomes by graph-regularized Tensor completion (FIST) to find its local optimal solution using the multiplicative updating rule [35], based on derivatives of $$ {\widehat{A}}_p$$, $$ {\widehat{A}}_y$$ and $$ {\widehat{A}}_x$$. Without loss of generality, we only show the derivations with respect to $$ {\widehat{A}}_p$$, and provide the FIST algorithm in Algorithm 1.

We first bring the equality constraint $$ \widehat{\mathcal{T}}=〚{\widehat{A}}_p,{\widehat{A}}_y,{\widehat{A}}_x〛$$ in Eq (2) into the objective function, and rewrite the objective function as

The partial derivative of $$ {\mathcal{J}}_1$$ with respect to $$ {\widehat{A}}_p$$ can be computed as

Following the derivations in [34], we obtain the partial derivatives of the second term $$ {\mathcal{J}}_2$$ as

Next, we combine $$ \frac{\partial {\mathcal{J}}_1}{\partial {\widehat{A}}_p}$$ and $$ \frac{\partial {\mathcal{J}}_2}{\partial {\widehat{A}}_p}$$ to obtain the overall derivative as

We then propose an efficient iterative algorithm FIST in Algorithm 1 to find the local optimum of the proposed graph regularized tensor completion problem with time complexity $$ O\left(r\right|\mathcal{M}|+{\sum }_{i\in \{x,y,p\}}(r^2{n}_i+r{n}_i^2))$$. FIST takes the incomplete spatial gene expression tensor $$ \mathcal{T}$$, PPI network and spatial chain graphs as input (line 1-2 in Algorithm 1), and outputs the inferred CPD representation of the complete spatial gene expression tensor $$ \widehat{\mathcal{T}}$$ (line 9 in Algorithm 1), via solving the optimization problem defined in Proposition 1 with the multiplicative updating rule (line 5-7 in Algorithm 1) based on the tensor calculus in Eqs (7)–(10). With the efficient tensor computation in Eqs (7)–(10), the algorithm can avoid computing the full Cartesian product graph and tensors, and break down the calculus into the computation on the individual graphs and the sparse tensors. Therefore, FIST is proven to be a scalable method, which only requires the space $$ O\left(\right|\mathcal{T}|+|\mathcal{M}\left|\right)$$ to store the sparse tensors, O(∑_i∈{x,y,p} |E_i|) to store the graphs, and O(∑_i∈{x,y,p} rn_i) to store the factor matrices. Thus, the overall space complexity is $$ O\left(\right|\mathcal{T}|+\left|\mathcal{M}\right|+{\sum }_{i\in \{x,y,p\}}\left(\right|E_i|+r{n}_i\left)\right)$$. The theoretical convergence analysis of FIST is also given in S1 File.

Algorithm 1 FIST: Fast Imputation of Spatially-resolved transcriptomes by graph-regularized Tensor completion

1: Input: 1) spatial gene expression tensor $$ \mathcal{T}\in {\mathbb{R}}_{+}^{n_p\times n_y\times n_x}$$, 2) binary mask tensor $$ \mathcal{M}\in {\mathbb{R}}_{[0,1]}^{n_p\times n_y\times n_x}$$ which indicates the observed part in $$ \mathcal{T}$$, 3) protein-protein interaction (PPI) network G_p and 4) hyper parameter λ.

2: Construct the spatial chain graphs G_x and G_y as described in the text.

3: Randomly initialize $$ {\widehat{A}}_p\in {\mathbb{R}}_{+}^{n_p\times r}$$, $$ {\widehat{A}}_y\in {\mathbb{R}}_{+}^{n_y\times r}$$ and $$ {\widehat{A}}_x\in {\mathbb{R}}_{+}^{n_x\times r}$$ as non-negative matrices.

4: while not converge do

5:  update $$ {\widehat{A}}_p$$ by Eq (11).

6:  update $$ {\widehat{A}}_y$$ by Eq (12).

7:  update $$ {\widehat{A}}_x$$ by Eq (13).

8: end while

9: Output: the low-rank matrices $$ {\widehat{A}}_p,{\widehat{A}}_y$$ and $$ {\widehat{A}}_x$$, which form the CPD representation of the inferred spatial gene expression tensor $$ \widehat{\mathcal{T}}=〚{\widehat{A}}_p,{\widehat{A}}_y,{\widehat{A}}_x〛\in {\mathbb{R}}_{+}^{n_p\times n_y\times n_x}$$.

Preparation of spatial gene expression datasets

We downloaded the spatial transcriptomic datasets from 10x Genomics (https://support.10xgenomics.com/spatial-gene-expression/datasets/), which is a collection of spatial gene expressions in 10 different tissue sections from mouse brain, mouse kidney, human breast cancer, human heart and human lymph node as listed in Table 2. All the sptRNA-seq datasets were collected with 10x Genomics Visium Spatial protocol (v1 chemistry) [14] to profiles each tissue section with a high density hexagonal array with 4,992 spots to achieve a resolution of 55 μm (1-10 cells per spot). To fit a tensor model on the spatial gene expression datasets, we organized each of the 10 datasets into a 3-way tensor $$ \mathcal{T}\in {\mathbb{R}}^{n_p\times n_y\times n_x}$$, where the (i, j, k)-th entry in $$ \mathcal{T}$$ is the UMI count of the i-th gene at the (k, j)-th coordinate in the array. Note that the spots are arranged in a perfect grid in earlier spatial transcriptomic arrays and the rows and columns in the grid can be used directly as the coordinates (n_x,n_y). In the Visium array slide, the spots are arranged in a hivegrid. To map the spatial coordinates (n_x,n_y), we shifted the odd-numbered rows by a half of a spot for a convenient arrangement of the spots in the tensor without loss of generality. We set the entries in $$ \mathcal{T}$$ to zeros if their UMI counts is lower than 3. We then removed the genes with no UMI counts across the spots, and removed the empty spots in the boundaries of the four sides in the H&E staining from $$ \mathcal{T}$$. The log-transformation is finally applied to every entry of $$ \mathcal{T}$$ as $$ {\mathcal{T}}_{i,j,k}\leftarrow \text{log}(1+{\mathcal{T}}_{i,j,k})$$. The sizes and densities of the 10 different spatial gene expression tensors after prepossessing are summarized in Table 2. Finally, we downloaded the full Homo sapiens and Mus musculus protein-protein interactions (PPI) networks from BioGRID 3.5 [31] as the gene network G_p to match with the genes in each dataset.

Table 2

10x Genomics spatial transcriptome data from 10 tissue sections.

Dataset	Tissue section	Tensor dimensions	Density
HBA1	Human Breast Cancer (Block A Section 1)	13, 426 × 60 × 77	0.093
HBA2	Human Breast Cancer (Block A Section 2)	13, 470 × 58 × 75	0.100
HH	Human Heart	7, 487 × 63 × 70	0.049
HLN	Human Lymph Node	12, 368 × 61 × 78	0.088
MKC	Mouse Kidney Section (Coronal)	12, 264 × 41 × 56	0.103
MBC	Mouse Brain Section (Coronal)	13, 570 × 49 × 74	0.110
MB1P	Mouse Brain Serial Section 1 (Sagittal-Posterior)	15, 404 × 62 × 67	0.115
MB2P	Mouse Brain Serial Section 2 (Sagittal-Posterior)	12, 497 × 63 × 65	0.077
MB1A	Mouse Brain Serial Section 1 (Sagittal-Anterior)	12, 658 × 59 × 66	0.105
MB2A	Mouse Brain Serial Section 2 (Sagittal-Anterior)	12, 295 × 63 × 66	0.082

Evaluations and performance measures

We applied 5-fold cross-validation to evaluate the performance of imputing spatial gene expressions by spatial spots or genes as follows:

Spot-wise evaluation: We chose 4-fold of the non-empty spatial spots for training and validation, and held out the rest 1-fold non-empty spatial spots as test data. When evaluating the expressions $$ {\mathcal{T}}_{:,j,k}\in {\mathbb{R}}^{n_p\times 1}$$ in the (k, j)-th spatial spot, we set the vectors $$ {\mathcal{T}}_{:,j,k}$$ and $$ {\mathcal{M}}_{:,j,k}$$ in the input tensor $$ \mathcal{T}$$ and mask tensor $$ \mathcal{M}$$ to zeros to indicate that the expressions in this spot are unobserved; next, we use the learned low-rank matrices $$ {\widehat{A}}_p,{\widehat{A}}_y$$ and $$ {\widehat{A}}_x$$ to construct the predicted gene expressions $$ {\widehat{\mathcal{T}}}_{:,j,k}$$ as described in Eq (1).

Gene-wise evaluation: For each gene, we chose 4-fold of its observed expressions (non-zeros in $$ \mathcal{T}$$) for training and validation, and held out the rest 1-fold of observed expressions as test data. When evaluating the 1-fold expressions in the i-th gene $$ {\mathcal{T}}_{i,:,:}\in {\mathbb{R}}^{n_y\times n_x}$$, we set the corresponding entries in $$ {\mathcal{T}}_{i,:,:}$$ and $$ {\mathcal{M}}_{i,:,:}$$ in the input tensor $$ \mathcal{T}$$ and mask tensor $$ \mathcal{M}$$ to zeros to indicate the expressions in this fold are unobserved; next, we use the learned low-rank matrices $$ {\widehat{A}}_p,{\widehat{A}}_y$$ and $$ {\widehat{A}}_x$$ to construct the predicted gene expressions $$ {\widehat{\mathcal{T}}}_{i,:,:}$$ as described in Eq (1).

The hyper-parameter λ is optimized by the validation set for FIST and the baseline methods. Denoting vectors $$ \mathbf{\text{t}}\in {\mathbb{R}}^{n\times 1}$$ and $$ \widehat{\mathbf{\text{t}}}\in {\mathbb{R}}^{n\times 1}$$ as the true and predicted expressions in the held-out spatial spot $$ {\mathcal{T}}_{:,j,k}$$ or the held-out entries in gene $$ {\mathcal{T}}_{i,:,:}$$, the imputation performance is evaluated by the following three metrics,

MAE (mean absolute error) $$ =\frac{1}{n}({\sum }_{i=1}^n|{\mathbf{\text{t}}}_i-{\widehat{\mathbf{\text{t}}}}_i\left|\right),$$

MAPE (mean absolute percentage error) $$ =\frac{1}{n}({\sum }_{i=1}^n|\frac{{\mathbf{\text{t}}}_i-{\widehat{\mathbf{\text{t}}}}_i}{{\mathbf{\text{t}}}_i}\left|\right),$$

R² (coefficient of determination) $$ =1-\left({\sum }_{i=1}^n{({\mathbf{\text{t}}}_i-{\widehat{\mathbf{\text{t}}}}_i)}^2\right){\left({\sum }_{i=1}^n{({\mathbf{\text{t}}}_i-\frac{{\sum }_{j=1}^n{\mathbf{\text{t}}}_j}{n})}^2\right)}^{-1}.$$

We expect a method to achieve smaller MAE and MAPE and larger R² for better performance.

FIST significantly improves the accuracy of imputing spatial gene expressions

The performances of FIST and the baseline methods except for ZIFA in the spot-wise evaluation are compared in Fig 2. ZIFA was excluded from this spot-wise evaluation as it does not allow empty rows (spots) in the implementation of its package, and thus is not applicable to the prediction of the held-out test spots. The average performances of all the spatial spots using each of the 10 sptRNA-seq datasets are shown as bar plots. FIST consistently outperforms all the baselines with lower MAE and MAPE, and larger R² in all the 10 datasets. We further applied right-tailed paired-sample t-tests on R² values to test against the alternative hypothesis that the R² produced by FIST has a larger mean than the mean of those produced by each of the baseline methods; and we also applied left-tailed paired-sample t-tests on MAE and MAPE values to test against the alternative hypothesis that the MAE and MAPE produced by FIST have a smaller mean than the mean of those produced by each of the baseline methods. The p-values in S1 Table show that compared with the baseline methods, FIST has significantly lower MAE and MAPE in each of the 10 datasets, and larger R² in all comparisons but one, in which FIST performed only slightly better than GWNMF on HB2P dataset by R².

Fig 2

Spot-wise cross-validation on 10x Genomics data.

The performances of the four compared methods on the 10 tissue sections are measured by 5-fold cross-validation. Each bar shows the mean of the imputation performance of one method on all the spatial spots. The result on each of the 10 datasets is shown in one vertical column separated by dashed lines. The means are also compared between FIST and each of the baseline methods in S1 Table by paired-sample t-tests.

The performances of FIST and the baseline methods in the gene-wise evaluation are compared in Fig 3. The average and standard deviation of the prediction performances across all the genes are shown as error bar plots in Fig 3. Similar to the spot-wise evaluation, FIST clearly outperforms all the baselines with more robust performances across all the genes, as the variances in all the three evaluation metrics are also lower than the other compared methods. To examine the prediction performance more closely, we also showed the distributions of MAE, MAPE and R² of individual genes in the 10 datasets in S3–S5 Figs, respectively. The result is consistent with the overall performance in Fig 3. The observations suggest that FIST indeed performs better than the other methods in the imputation accuracy informed by the spatial information in the tensor model. It is also noteworthy that GWNMF, the MF method regularized by the spatial graph applied to each individual gene slice in tensor $$ \mathcal{T}$$, outperforms the other baselines in almost all the datasets. This observation further confirms that the spatial patterns maintained in each gene slice is informative for the imputation task. It is clear that FIST outperformed GWNMF with better use of the spatial information coupled with the functional modules of the PPI network G_p and the joint imputation of all the genes in the tensor $$ \mathcal{T}$$.

Fig 3

Gene-wise cross-validation on 10x Genomics data.

The performances of the five compared methods on the 10 tissue sections are measured by 5-fold cross-validation. Each error bar shows the mean and variance of the imputation performance for one method on all the genes. The result on each of the 10 datasets is shown in one vertical column separated by dashed lines.

Cartesian product graph regularization plays a significant role

To demonstrate that the Cartesian product graph regularization in FIST significantly improves the imputation accuracy, we show in Fig 4 the performance of FIST in each of the 10 datasets by varying the graph hyper-parameter λ in the spot-wise evaluation. By increasing λ from 0 to 0.1 to put more belief on the graph information, we observe an appreciable reduction on the MAE and MAPE, and increase on R² across all the 10 datasets. The observation strongly suggests that the predictions by FIST are improved by leveraging the information carried in the CPG topology, and the belief on the graph information can be effectively optimized by using a validation set in the cross-validation strategy.

Fig 4

Analysis of Cartesian product graph regularization with varying network hyper-parameter in spot-wise evaluation.

The plots show the imputation performance of FIST on the ten 10x Genomics datasets with varying network hyper-parameters in {0, 0.1, 1} by MAE, MAPE and R².

To further understand the associations between the CPG regularization and characteristics of the expressions of the genes, we analyzed the genes that are benefiting most from the regularization by the CPG in the gene-wise evaluation. In particular, in the grid search of the optimal λ weight on the CPG regularization term by the validation set, we count the percentage of the genes with optimal λ = 0.01 rather than 0, which means completely ignoring the regularization. To correlate the improved imputations with the sparsity of the gene expressions, we divided all the genes into 4 equally partitioned groups (L1-4) ordered by their densities in the sptRNA-seq data, where L1 and L4 contain the sparsest and the densest gene slices, respectively. For each of the four density levels, we count the percentage of gene slices that benefit from the CPG regularization and plot the results in Fig 5A. In the plots, there is a clear trend that the sparser a gene slice, the more likely it benefits from the CPG regularization in all the 10 datasets. In the densest L4 group, as low as 20% of the genes can benefit from the CPG regularization versus more than 50% in the sparest L1 group. This is understandable that there is less training information available for sparsely expressed genes (with more dropouts) and the spatial and functional information in the CPG can play a more important role in the imputation by seeking information from the gene’s spatial neighbors or the functional neighbors in the PPI network. This observation is also consistent with the fact that the performance of tensor completion tends to degrade severely when only a very small fraction of entries are observed [43, 44], and therefore those sparser gene slices tend to benefit more from the side information carried in the CPG.

Fig 5

Analysis of Cartesian product graph regularization on gene-wise evaluation.

(A) The percentages of genes benefit from the CPG are plotted by their densities in four different ranges. Each colored line represents one of the 10 datasets. (B) The total reduction of MAE using the original and permuted graphs are compared across the 10 tissue sections.

We also compared the performance of FIST using the CPG of G_x, G_y and G_p with the one using a randomly permuted graph from the CPG. To generate the random CPG, we first generated three random graphs by permute G_x, G_y and G_p individually which also preserves the degree distributions of the original graphs, by randomly swapping the edges in each graph while keeping the degree of each node. Then we measured the performances of FIST using the original CPG and the CPG obtained from the permuted graphs by MAE reduction, which is the total reduction of MAE on all the genes by varying hyperparameter λ from 0 to 0.01 meaning not using the graph versus using the graph. The comparisons across all the 10 datasets are shown in Fig 5B. We observe that the FIST using the original graphs receives much higher MAE reduction than the FIST using the permuted graphs. This observation suggests that the topology in the original CPG carries rich information that is helpful for the imputation task beyond just the degree distributions preserved in the random graphs.

FIST imputations recover spatial patterns enriched by highly relevant functional terms

To demonstrate that imputations by FIST can reveal spatial gene expression patterns with highly relevant functional characteristics among the genes in the spatial region, we performed comparative GO enrichment analysis of gene clusters detected with the imputed gene expressions. We conducted a case study on the Mouse Kidney Section data to further analyze the associations between the spatial gene clusters and the relevance between their functional characteristics and three kidney tissue regions, cortex, outer stripe of the outer medulla (OSOM) and inner stripe of the outer medulla (ISOM).

To validate the hypothesis that the imputed sptRNA-seq tensor $$ \tilde{\mathcal{T}}$$ given below

Fig 6

Enrichment analysis on the sparse and imputed sptRNA-seq data.

The total number of significantly enriched clusters (with at least one enriched GO term with FDR adjusted p-value < 0.05) in the 10 tissue sections are shown.

Finally, we conducted a case study on the Mouse Kidney Section and present the highly relevant functional characteristics in different tissues in mouse kidney detected with the imputations by FIST. For each of the 100 gene clusters generated by K-means as described above, we collapsed the corresponding gene slices in $$ \tilde{\mathcal{T}}$$ into a n_x × n_y matrix by averaging the slices to visualize the center of the gene cluster. The enrichment results of all the 100 clusters are given in S2 Table. We focus on 3 kinds of representative clusters in Fig 7 which match well with three distinct mouse kidney tissue regions: cortex, ISOM (inner stripe of outer medulla and OSOM (outer stripe of outer medulla). By investigating the enriched GO terms by the clusters (p-values shown in Table 3), we found their functional relevance to cortex, ISOM and OSOM regions. We found that the spatial gene cluster 9 which is highly expressed in cortex specifically enriched biological processes for the regulation of blood pressure (GO:0008217, GO:0003073, GO:0008015 and GO:0045777) and transport/homeostasis of inorganic molecules (GO:0055067 and GO:0015672). The spatial gene cluster 23 and 28 which are also highly expressed in cortex enriched cellular pathways that are critical for the polarity of cellular membranes (GO:0086011, GO:0034763, GO:1901017, GO:0032413 and GO:1901380) and the transport of cellular metabolites (GO:1901605, GO:0006520, GO:0006790 and GO:0043648), respectively. These observations are consistent with previous studies reporting the regulation of kidney function by the above listed biological processes in cortex [46–49]. In contrast, the pattern analysis of spatial gene expressions in cluster 4, 8, 25 and 52 which are highly expressed in OSOM in kidney showed that catabolic processes of organic and inorganic molecules are specifically enriched such as GO:0015711, GO:0046942, GO:0015849, GO:0015718, GO:0010498, GO:0043161, GO:0044282, GO:0016054, GO:0046395, GO:0006631, GO:0072329, GO:0009062 and GO:0044242. These cellular processes are known to be active in renal proximal tubule which exists across cortex and OSOM [50–55]. Distinctively, the spatial gene clusters highly expressed in ISOM enriched pathways for nucleotide metabolisms (GO:0009150, GO:0009259 and GO:0006163) in cluster 3 and renal filtration (GO:0097205 and GO:0003094) in cluster 5. Collectively, these observations demonstrate that FIST could identify physiologically relevant distinctive spatial gene expression patterns in the mouse kidney dataset. Further, it suggests that FIST can provide a high-resolution anatomical analysis of organ functions in sptRNA-seq data.

Fig 7

FIST recovers spatial gene expression patterns on Mouse Kidney Section.

The H&E image of the mouse kidnesy section is shown in the middle with circles roughly separating the tissue area of Cortex, the outer stripe of the outer medulla (OSOM) and the inner stripe of the outer medulla (ISOM) from outer to inner regions. The gene expression patterns of the clusters in each of the three regions are grouped in the same box labeled by the region.

Table 3

Functional terms enriched by spatial gene clusters (most significantly relevant functions).

Region	Cluster	Significantly Enriched GO terms
Cortex	Cluster 9	GO:0003073—regulation of systemic arterial blood pressure (p = 9.1 × 10−6)
		GO:0008217—regulation of blood pressure (p = 1.0 × 10−4)
		GO:0055067—monovalent inorganic cation homeostasis (p = 4.3 × 10−4)
		GO:0008015—blood circulation (p = 5.3 × 10−3)
		GO:0045777—positive regulation of blood pressure (p = 5.8 × 10−3)
		GO:0015672—monovalent inorganic cation transport (p = 2.3 × 10−2)
	Cluster 23	GO:0086011—membrane repolarization during action potential (p = 2.2 × 10−3)
		GO:0034763—negative regulation of transmembrane transport (p = 2.2 × 10−3)
		GO:1901017—negative regulation of potassium ion transmembrane transporter activity (p = 2.4 × 10−3)
		GO:0032413—negative regulation of ion transmembrane transporter activity (p = 2.7 × 10−3)
		GO:1901380—negative regulation of potassium ion transmembrane transport (p = 3.4 × 10−3)
	Cluster 28	GO:1901605—alpha-amino acid metabolic process (p = 4.8 × 10−10)
		GO:0006520—cellular amino acid metabolic process (p = 6.4 × 10−9)
		GO:0006790—sulfur compound metabolic process (p = 3.1 × 10−6)
		GO:0043648—dicarboxylic acid metabolic process (p = 8.4 × 10−6)
OSOM	Cluster 4	GO:0015711—organic anion transport (p = 7.7 × 10−7)
		GO:0046942—carboxylic acid transport (p = 1.1 × 10−4)
		GO:0015849—organic acid transport (p = 1.1 × 10−4)
		GO:0015718—monocarboxylic acid transport (p = 5.0 × 10−3)
	Cluster 8	GO:0010498—proteasomal protein catabolic process (p = 1.3 × 10−3)
		GO:0006497—protein lipidation (p = 1.3 × 10−3)
		GO:0042158—lipoprotein biosynthetic process (p = 1.3 × 10−3)
		GO:0043161—proteasome-mediated ubiquitin-dependent protein catabolic process (p = 1.3 × 10−3)
	Cluster 25	GO:0044282—small molecule catabolic process (p = 5.5 × 10−19)
		GO:0016054—organic acid catabolic process (p = 1.0 × 10−18)
		GO:0046395—carboxylic acid catabolic process (p = 1.0 × 10−18)
		GO:0006631—fatty acid metabolic process (p = 2.9 × 10−16)
		GO:0072329—monocarboxylic acid catabolic process (p = 9.6 × 10−14)
		GO:0009062—fatty acid catabolic process (p = 1.0 × 10−13)
		GO:0044242—cellular lipid catabolic process (p = 4.7 × 10−11)
	Cluster 52	GO:0006732—coenzyme metabolic process (p = 1.2 × 10−10)
		GO:0006520—cellular amino acid metabolic process (p = 1.6 × 10−10)
		GO:1901605—alpha-amino acid metabolic process (p = 2.3 × 10−9)
		GO:0044282—small molecule catabolic process (p = 2.1 × 10−8)
		GO:0000096—sulfur amino acid metabolic process (p = 2.3 × 10−7)
ISOM	Cluster 3	GO:0009150—purine ribonucleotide metabolic process (p = 7.4 × 10−5)
		GO:0009259—ribonucleotide metabolic process (p = 7.4 × 10−5)
		GO:0006163—purine nucleotide metabolic process (p = 7.4 × 10−5)
		GO:0019693—ribose phosphate metabolic process (p = 7.4 × 10−5)
		GO:0072521—purine-containing compound metabolic process (p = 7.4 × 10−5)
	Cluster 5	GO:0048872—omeostasis of number of cells (p = 4.5 × 10−5)
		GO:0030218—erythrocyte differentiation (p = 3.2 × 10−3)
		GO:0034101—erythrocyte homeostasis (p = 3.2 × 10−3)
		GO:0003094—glomerular filtration (p = 3.2 × 10−3)
		GO:0097205—renal filtration (p = 3.2 × 10−3)