Imputation of spatially-resolved transcriptomes by graph-regularized tensor completion

High-throughput spatial-transcriptomics RNA sequencing (sptRNA-seq) based on in-situ capturing technologies has recently been developed to spatially resolve transcriptome-wide mRNA expressions mapped to the captured locations in a tissue sample. Due to the low RNA capture efficiency by in-situ capturing and the complication of tissue section preparation, sptRNA-seq data often only provides an incomplete profiling of the gene expressions over the spatial regions of the tissue. In this paper, we introduce a graph-regularized tensor completion model for imputing the missing mRNA expressions in sptRNA-seq data, namely FIST, Fast Imputation of Spatially-resolved transcriptomes by graph-regularized Tensor completion. We first model sptRNA-seq data as a 3-way sparse tensor in genes (p-mode) and the (x, y) spatial coordinates (x-mode and y-mode) of the observed gene expressions, and then consider the imputation of the unobserved entries or fibers as a tensor completion problem in Canonical Polyadic Decomposition (CPD) form. To improve the imputation of highly sparse sptRNA-seq data, we also introduce a protein-protein interaction network to add prior knowledge of gene functions, and a spatial graph to capture the the spatial relations among the capture spots. The tensor completion model is then regularized by a Cartesian product graph of protein-protein interaction network and the spatial graph to capture the high-order relations in the tensor. In the experiments, FIST was tested on ten 10x Genomics Visium spatial transcriptomic datasets of different tissue sections with cross-validation among the known entries in the imputation. FIST significantly outperformed the state-of-the-art methods for single-cell RNAseq data imputation. We also demonstrate that both the spatial graph and PPI network play an important role in improving the imputation. In a case study, we further analyzed the gene clusters obtained from the imputed gene expressions to show that the imputations by FIST indeed capture the spatial characteristics in the gene expressions and reveal functions that are highly relevant to three different kinds of tissues in mouse kidney.

Introduction Dissection of complex genomic architectures of heterogeneous cells and how they are organized spatially in tissue are essential for understanding the molecular and cellular mechanisms underlying important phenotypes. For example, each tumor is a mixture of different types of proliferating cancerous cells with changing genetic materials [1]. The cancer cell sub-populations co-evolve in the micro-environment formed around their spatial locations. It is important to understand the cell-cell interactions and signaling as well as the functioning of each individual cell to develop effective cancer treatment and eradicate all cancer clones at their locations [2]. Conventional gene expression analyses have been limited to low-resolution bulk profiling that measures the average transcription levels in a population of cells. With singlecell RNA sequencing (scRNA-seq) [3][4][5], single cells are isolated with a capture method such as fluorescence-activated cell sorting (FACS), Fluidigm C1 or microdroplet microfluidics and then the RNAs are captured, reverse transcribed and amplified for sequencing the RNAs barcoded for the individual origin cells [6,7]. While scRNA-seq is useful for detecting the cell heterogeneity in a tissue sample, it does not provide the spatial information of the isolated cells. To map cell localization, earlier in-situ hybridization methods such as FISH [8], FISSEQ [9], smFISH [10] and MERFISH [11] were developed to profile up to a thousand targeted genes in pre-constructed references with single-molecule RNA imaging. Based on in-situ capturing technologies, more recent spatial transcriptomics RNA sequencing (sptRNA-seq) [12][13][14][15] combines positional barcoded arrays and RNA sequencing with single-cell imaging to spatially resolve RNA expressions in each measured spot in the spatial array [12,[16][17][18]. These new technologies have transformed the transcriptome analysis into a new paradigm for connecting single-cell molecular profiling to tissue micro-environment and the dynamics of a tissue region [19][20][21].
With in-situ capturing technology, RNAs are captured and sequenced in the spots on the spatial genomic array aligned to the locations on the tissue. For example, spatial transcriptomics technology based on 10x Genomics Visium kit reports the number of copies of RNAs by counting unique molecular identifiers (UMIs) in the read-pairs mapped to each gene [22]. There are still significant technical difficulties. First, in-situ capturing has a low RNA capture efficiency. The earlier spatial transcriptomics technology's detection efficiency is as low as 6.9% and 10x Genomics Visium has only a slightly improved efficiency [23]. In addition, the sample preparation requires highly specific handling of tissue sections. The spots in some tissue regions might entirely fail to fix and permeabilize RNAs due to various possible issues in preparing tissue sections. A few examples of such regions are shown in S1 Fig. Thus, sptRNA-seq data often only provides an incomplete profiling of the gene expressions over the spatial regions of the tissue. Similarly, in scRNA-seq data analysis, the missing gene expressions are called dropout events, which refer to the false quantification of a gene as unexpressed due to the failure in amplifying the transcripts during reverse-transcription [24]. It has been shown in previous studies on scRNA-seq data that normalizations will not address the dropout effects [22,25]. In the literature, many imputation methods such as Zero-inflated factor analysis (ZIFA) [26], Zero-Inflated Negative Binomial-based Wanted Variation Extraction (ZINB-WaVE) [27] and BISCUIT [25] have been developed to impute scRNA-seq. While these methods are also applicable to impute the spatial gene expressions, they ignore a unique characteristic of sptRNA-seq data, which is the spatial information among the gene expressions in the spatial array, and do not fully take advantage of the functional relations among genes for more reliable joint imputation.
To provide a more suitable method for imputation of spatially-resolved gene expressions, we introduce FIST, Fast Imputation of Spatially-resolved transcriptomes by graph-regularized Tensor completion. FIST is a tensor completion model regularized by a product graph as illustrated in Fig 1. FIST models sptRNA-seq data as a 3-way sparse tensor in genes (p-mode) and the (x, y) spatial coordinates (x-mode and y-mode) of the observed gene expressions ( Fig 1A). As shown in Fig 1B, a protein-protein interaction network models the interactions between pairs of genes in the gene mode, and the spatial graph is modeled by a product graph of two chain graphs for columns (x-mode) and rows (y-mode) in the grid to capture the spatial The input sptRNA-seq data is modeled by a 3-way sparse tensor in genes (p-mode) and the (x, y) spatial coordinates (x-mode and y-mode) of the observed gene expressions. H&E image is also shown to visualize the cell morphologies aligned to the spots. (B) A protein-protein interaction network and a spatial graph are integrated as a product graph for tensor completion. The spatial graph is also a product graph of two chain graphs for columns (x-mode) and rows (y-mode) in the grid. (C) After the imputation, the CPD form of the complete tensor can be used to impute any missing gene expressions, e.g. the entry (k, j, i) can be reconstructed as the sum of the element-wise multiplications of the three components ½Â p � k;: , ½Â y � j;: and ½Â x � i;: . https://doi.org/10.1371/journal.pcbi.1008218.g001

PLOS COMPUTATIONAL BIOLOGY
Imputation of Spatially-resolved Transcriptomes by Graph-regularized Tensor Completion relations among the (x, y) spots. The Cartesian product of these graphs with prior knowledge of gene functions and the spatial relations among the capture spots are then introduced as a regularization of tensor completion to obtain the Canonical Polyadic Decomposition (CPD) of the tensor. The imputation of the unobserved entries can then be derived by reconstructing the entries in the completed tensor shown in Fig 1C. In the experiments, we comprehensively evaluated FIST on ten 10x Genomics Visium spatial genomics datasets by comparison with widely used methods for single-cell RNA sequencing data imputation. We also analyzed a mouse kidney dataset with more functional interpretation of the gene clusters obtained by the imputed gene expressions to detect highly relevant functions in the clusters expressed in three kidney tissue regions, cortex, outer stripe of the outer medulla (OSOM) and inner stripe of the outer medulla (ISOM).

Materials and methods
In this section, we first describe the task of spatial gene expression imputation, and next introduce the mathematical model for graph-regularized tensor completion problem. We then present a fast iterative algorithm FIST to solve the optimization problem defined to optimize the model. We also provide the convergence analysis of proposed algorithm in S1 File. Finally, we provide a review of several state-of-the-art methods for scRNA-seq data imputation, which are also compared in the experiments later. The notations which will be used for the derivations in the forthcoming sections are summarized in Table 1.

Imputation of spatial gene expressions by tensor modeling
Let T 2 R n p �n y �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 T 2 R n p �n y �n x þ from T as illustrated in Fig 1C. However, it is computationally expensive and often infeasible to compute or store a dense tensorT , especially in high spatial resolutions with millions of spots. Therefore, we propose to compute an economy-size representation of

PLOS COMPUTATIONAL BIOLOGY
T via an equality constraintT ¼〚Â p ;Â y ;Â x 〛 , which is called Canonical Polyadic Decomposition (CPD) [28] ofT defined beloŵ where r is the rank ofT , and � denotes the vector outer product. Here, ½Â x � :;i is the i-th column of the low-rank matrixÂ x 2 R n x �r , which can be similarly defined for ½Â y � :;i and ½Â p � :;i . By utilizing the tensor CPD form, we replaced the optimization variables fromT toÂ p ,Â y andÂ x , reducing the number of parameters from n p n y n x to r(n p + n y + n x ). The advantage of the tensor representation is to incorporate the 2-D spatial x-mode and y-mode such that the grid structure is preserved within the columns and the rows of the spatial array in the tensor, which contains useful spatial information. Next, we introduce the tensor completion model overÂ p ,

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 tensorT 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 T is also required to be preserved inT as the completion task requires; and iii) the inferred tensor T is compressed as the economy-size representationT ¼〚Â p ;Â y ;Â x 〛for scalable space and time efficiencies. The novel optimization formulation is shown below in Proposition 1, Proposition 1. The complete spatial gene expression tensorT 2 R n p �n y �n x can be obtained by solving the following optimization problem: where λ 2 [0, 1] is a model hyperparameter; ⊛ denotes the Hadamard product; and jj:jj F denotes the Frobenius norm of a tensor.
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 M to indicate the indices of the observed entries in T . The (i, j, k)-th entry M i;j;k which is defined below, represents whether the (i, j, k)-th element in T is observed or not.
By introducing the squared-error in F -norm jjM ⊛ ðT ÀT Þjj 2 F in the model, we ensure the inferred spatial gene expression tensorT is consistent with its observed counterparts in 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 xaxis. 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 Bio-GRID 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] Gðx; y; pÞ ¼ ðV; EÞ of the three individual graphs G x , G y and G p to regularize the elements inT , where |V| = n x n y 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 Gðx; y; pÞ are obtained as Wðx; y; pÞ where � denotes the Kronecker sum [33].
By introducing the term vecðT Þ T Lðx; y; pÞvecðT Þ in Eq (2), the inferred gene expression values inT are ensured to be smooth over the manifolds of the product graph Gðx; y; pÞ, such that a pair of tensor entriesT a p ;a y ;a x andT 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 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 fÂ p ;Â y ;Â x g 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Â p ,Â y andÂ x . Without loss of generality, we only show the derivations with respect toÂ p , and provide the FIST algorithm in Algorithm 1.
We first bring the equality constraintT ¼〚Â p ;Â y ;Â x 〛in Eq (2) into the objective function, and rewrite the objective function as The partial derivative of J 1 with respect toÂ p can be computed as where T ð1Þ 2 R n p �n x n y denotes the matrix flattened from tensor T ; � denotes the Khatri-Rao product [28]. Note that the term M ð1Þ ⊛T ð1Þ in Eq (4) implies we only need to compute the entries inT which correspond to the non-zero entries (indices of the observed gene expression) in M. The rest of the computation in Eq (4) involves the well-known MTTKRP (matricized tensor times Khatri-Rao product) [36] operation, which is in the form of X ð1Þ ðÂ x �Â y Þ, and can be computed in OðrjXjÞ if X is a sparse tensor with jXj non-zeros, andÂ x andÂ y have r columns. Thus, the overall time complexity of computing Eq (4) is OðrjMjÞ.
Following the derivations in [34], we obtain the partial derivatives of the second term J 2 as where F i ¼Â T iÂi , and Y i ¼Â T i L iÂi , for all i 2 {x, y, p}. It is not hard to show that the complexity of computing Eq (5) is Oð P i2fx;y;pg ðr 2 n i þ rn 2 i ÞÞ. Next, we combine @J 1 @Â p and @J 2 @Â p to obtain the overall derivative as , which are defined below, where According to Eq (6), the objective function J objective will monotonically decrease under the following multiplicative updating rule, where ½Â p � a;b denotes the (a, b)-th element in matrixÂ p . Similarly, we can derive the update rule for ½Â x � a;b and ½Â p � a;b as follows, 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ðrjMj þ P i2fx;y;pg ðr 2 n i þ rn 2 i ÞÞ. FIST takes the incomplete spatial gene expression tensor 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 tensorT (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ðjT j þ jMjÞ to store the sparse tensors, O(∑ i2{x,y,p} |E i |) to store the graphs, and O(∑ i2{x,y,p} rn i ) to store the factor matrices. Thus, the overall space complexity is OðjT j þ jMj þ P i2fx;y;pg ðjE i j þ rn i ÞÞ. 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 T 2 R n p �n y �n x þ , 2) binary mask tensor M 2 R n p �n y �n x ½0;1� which indicates the observed part in T , 3) proteinprotein 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.

Related methods for comparison
To benchmark the performance of FIST, we compared it with three matrix factorization (MF)based methods (with graph regularizations) and a nearest neighbors (NN)-based method, which have been applied to impute various types of biological data including the imputation of dropouts in single-cell RNA sequencing (scRNA-seq) data. Note that NMF-based methods have been shown to be effective for learning latent features and clustering high-dimension sparse genomic data [37].
• ZIFA: Zero-inflated factor analysis (ZIFA) [26] factorizes the single cell expression data Y 2 R N�D where N and D denote the number of single cells and genes respectively, into a factor loading matrix A 2 R K�D and a matrix Z 2 R N�K which spans the latent low-dimensional space where dropouts can happen with a probability specified by an exponential decay associated with the expression levels. The imputed matrix can be computed asŶ ¼ ZA þ m, where m 2 R 1�D is the latent mean vector.
• REMAP: Since ZIFA is a probabilistic MF model which does not utilize the spatial and gene networks, we therefore also compare with REMAP [38], which was developed to impute the missing chemical-protein associations for the identification of the genome-wide off-targets of chemical compounds. REMAP factorizes the incomplete chemical-protein interactions matrix into the chemical and protein low-rank matrices, which are regularized by the compound similarity graph and protein sequence similarity (NCBI BLAST [39]) graph respectively.
• GWNMF: Both ZIFA and REMAP are only applicable to the spot-by-gene matrix which is a flatten of an input tensor T . Such flattening process assumes the spots are independent from each other, and thus loses the spatial information. To keep the spatial arrangements, we also apply MF to each n x × n y slice in tensor T . Specially, we adopt the graph regularized weighted NMF (GWNMF) [40] method to impute each n x × n y gene slice. We let GWNMF use the same x-axis and y-axis graphs G x and G y as described in the previous section to regularize the MF.
• Spatial-NN: It has been observed that in sparse high-dimensional scRNA-seq data, constructing a nearest neighbor (NN) graph among cells can produce more robust clusters in the presence of dropouts because of taking into account the surrounding neighbor cells [41]. Such rationale has be considered in the clustering methods such as Seurat [42] and shared nearest neighbors (SNN)-Cliq [41], and can also be adopted to impute the spatial gene expression data. We introduce a SNN-based baseline Spatial-NN using neighbor averaging to compare with FIST. Specifically, to impute the missing expression of a target spot, Spatial-NN first searches its spatially nearest spots with observed gene expressions, then assign their average gene expression to the target spot.
We used the provided Python package (https://github.com/epierson9/ZIFA) to experiment with ZIFA, and the provided MATLAB package (https://github.com/hansaimlim/REMAP) to experiment with REMAP. To apply both methods, we rearranged the data tensor T 2 R n p �n y �n x to a matrix T 2 R N�n p , where N = n x n y denotes the total number of spots. The spatial graph of REMAP is constructed by connecting two spots if they are spatially adjacent. REMAP adopts the same PPI network as the gene graph G p as used by FIST. We used MATLAB to implement GWNMF and Spatial-NN ourselves to impute each gene slice T i 2 R n x �n y in T . In the comparisons, the graph hyperparameter λ of FIST is selected from {0, 0.01, 0.1, 1}. The graph hyperparameters of REMAP and GWNMF are set by searching the grids from {0.1, 0.5, 0.9, 1} and {0, 0.1, 1, 10, 100} respectively as suggested in the original studies. Note that different methods use different scales of graph hyperparameters since the gradients of their variables with respective to the regularization terms are in different scales. The optimal hyper-parameters are selected by the validation set for each method. For FIST, REMAP and GWNMF, we applied PCA on matrix T 2 R N�n p to determine the rank r 2 [200, 300] of the low-rank factor matrices, such that at least 60% of the variance is accounted for by the top-r PCA components of T. The latent dimension K of ZIFA is set to 10 since it is time consuming to run ZIFA with a larger K. We also observed that increasing K from 10 to 50 does not show clear improvement on the imputation accuracy.

Results
In this section, we first describe data preparation and performance measure, and then show the results of spatial gene expression imputation. We also analyzed the results by the gene-wise density of the gene expressions and regularization by permuted graphs. Finally, we analyzed the imputed spatial gene expressions in the Mouse Kidney Section dataset to show several interesting gene clusters revealing functional characteristics of the three tissue regions, cortex, OSOM and ISOM.

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

PLOS COMPUTATIONAL BIOLOGY
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 T 2 R n p �n y �n x , where the (i, j, k)-th entry in 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 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 T . The log-transformation is finally applied to every entry of T as T i;j;k logð1 þ 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.

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 T :;j;k 2 R n p �1 in the (k, j)-th spatial spot, we set the vectors T :;j;k and M :;j;k in the input tensor T and mask tensor M to zeros to indicate that the expressions in this spot are unobserved; next, we use the learned low-rank matricesÂ p ;Â y andÂ x to construct the predicted gene expressionsT :;j;k as described in Eq (1).
• Gene-wise evaluation: For each gene, we chose 4-fold of its observed expressions (non-zeros in 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 T i;:;: 2 R n y �n x , we set the corresponding entries in T i;:;: and M i;:;: in the input tensor T and mask tensor M to zeros to indicate the expressions in this fold are unobserved; next, we use the learned low-rank matriceŝ A p ;Â y andÂ x to construct the predicted gene expressionsT i;:;: as described in Eq (1).
The hyper-parameter λ is optimized by the validation set for FIST and the baseline methods. Denoting vectors t 2 R n�1 andt 2 R n�1 as the true and predicted expressions in the heldout spatial spot T :;j;k or the held-out entries in gene T i;:;: , the imputation performance is evaluated by the following three metrics, • MAE (mean absolute error) ¼ 1 n ð P n i¼1 jt i Àt i jÞ; • MAPE (mean absolute percentage error) ¼ 1 n ð P n i¼1 j t i Àt i t i jÞ; We expect a method to achieve smaller MAE and MAPE and larger R 2 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 2 in all the 10 datasets. We further applied right-tailed paired-sample t-tests on R 2 values to test against the alternative

PLOS COMPUTATIONAL BIOLOGY
hypothesis that the R 2 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 2 in all comparisons but one, in which FIST performed only slightly better than GWNMF on HB2P dataset by R 2 .
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 2 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 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 T .

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 2 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.
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

PLOS COMPUTATIONAL BIOLOGY
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.
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 tensorT given below can better capture gene functional modules than the sparse sptRNA-seq tensor T does, we first rearranged both sptRNA-seq tensors into matricesT 2 R N�n p and T 2 R N�n p , where N = n x n y denotes the total number of spots. We then applied K-means on each matrix to partition the genes into 100 clusters. Next, we used the enrichGO function in the R package clusterProfiler [45] to perform the GO enrichment analysis of the gene clusters. The total number of significantly enriched gene clusters (FDR adjusted p-value < 0.05) in each of the 10 tissue sections are shown in Fig 6, which clearly tells that K-means on the imputed sptRNA-seq data produces much more significantly enriched clusters across all the 10 tissue sections than the sparse sptRNA-seq data without imputation.
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 inT 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 respectively. These observations are consistent with previous studies reporting the regulation of kidney function by the above listed biological processes in cortex [46][47][48][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][51][52][53][54][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.

Experiments on additional low-resolution spatial gene expression datasets
To demonstrate that FIST is broadly applicable to impute the spatial gene expression data generated with different platforms, we performed additional experiments on spatial transcriptomics datasets from 3 replicates of mouse tissue (olfactory bulb) provided from an earlier study [56]. Developed before 10x Genomics Visium Spatial protocol, the spatial transcriptomics technology [56] applies an aligned array to profile tissue with both lower spot density and larger spot size (1,007 spots in total, and 200 μm between spots). The design achieves a resolution of 100 μm (10-40 cells per spot). Similar to the experiments on the 10x Genomics data, we organized each of the 3 tissue replicates into a tensor T 2 R n p �n y �n x , where n x = 33 and n y = 35 in all the 3 replicates, and n p is 14,198, 13,818 and 138,40, respectively in replicate 1,2 and 3. The (i, j, k)-th entry in T is the RPKM value of the i-th gene at the (k, j)-th coordinate in the array.
We performed the spot-wise 5-fold cross-validation as we did in the 10x Genomics data to compare the performances of FIST and the same baseline methods. The distributions of MAE, MAPE and R 2 on all the spatial spots in each of the 3 tissue replicates are shown in   Consistent with the observations in the previous Figs 2 and 3, FIST clearly outperforms all the baselines with lower MAE and MAPE, and larger R 2 in all the 3 replicates. The results suggest that FIST has a potential to be applied to various spatial transcriptomics datasets of different resolution and sparcity to achieves better imputation performance by modeling the spatial data as tensors, and including the prior knowledge with the CPG regularization.

PLOS COMPUTATIONAL BIOLOGY
To confirm that the imputation accuracy of FIST is significantly improved by the CPG regularization, we showed in Fig 9 the performance of FIST in each of the 3 replicates by varying the graph hyper-parameter λ in the spot-wise evaluation. It is also consistent with the observation in the previous Fig 4, in which we can observe remarkable reduction on the MAE and MAPE and improvement on R 2 by increasing λ to 0.1. The observation also verifies that the CPG topology is informative for the imputation task.

Discussions
In this study, we propose to apply tensor modeling of multidimensional structure in spatiallyresolved gene expression data mapped by the 2D spatial array. To the best of our knowledge, this is the first work to model the imputation of spatially-resolved transcriptomes as a tensor completion problem. Our key observations in the experiments with the ten 10x Genomics Visium spatial transcriptomic datasets are that 1) the imputation accuracy is significantly improved by leveraging the tensor representation of the sptRNA-seq data, and 2) by incorporating the spatial graph and PPI network, the accuracy the imputation and the content of the functional information in the imputed spatial gene expressions can be further improved significantly.
We observed that the genes that are more sparsely expressed can benefit more from the adjacency information in the spatial graph and the functional information in the PPI network. These genes can be empirically detected with a validation set to tune the only hyper-parameter λ for deciding if the regularization by the product graph is needed for the imputation of a gene. Thus, we expect a low risk of overfitting in applying FIST to other datasets. In addition, the functional analysis of the spatial gene clusters detected on the Mouse Kidney Section data further confirms that FIST detects gene clusters with more spatial characteristics that are consistent with the physiological features of the tissue.
Overall, we concluded that FIST is an effective and easy-to-use approach for reliable imputation of spatially-resolved gene expressions by modeling the spatial relation among the spots in the spatial array and the functional relation among the genes. The imputation results by FIST are both more accurate and functionally interpretable. FIST is also highly generalizable to other spatial transcriotomics datasets with high scalability and only one hyper-parameter needed to tune.
Although our experiments mainly focused on medium density 10x Genomics Visium kit array (5000 spots or 1000 spots), we also plan to further develop variations of FIST for highresolution spatial transcriptomics datasets with millions of spots generated by high-definition spatial transcriptomics (HDST) [14]. The HDST datasets from the study includes 3 mouse tissue sections from olfactory bulb and 3 human tissue sections from breast cancer using hexagonal array to profile tissue with a high density (1,467,270 spots in total) to achieve a resolution of 2 μm. The imputation tasks on the HDST datasets are quite different for two reasons: First, each cell spans multiple spots. Simple imputation of each spot is not a

PLOS COMPUTATIONAL BIOLOGY
well-defined learning problem. Thus, the segmentation of the spots into each individual cell might be necessary as a pre-processing step. Second, the capture efficiency of HDST is as low as 1.3%, which leads to much sparser data for imputation. Our preliminary analysis indicate that the gene expression on the spots are too sparse to be meaningful unless they are aggregated among the spots in a larger region. We plan to develop variations of FIST to overcome these additional challenges.
Another interesting future direction is to develop a variation of FIST for imputing spatial gene expressions with additional information from matched single-cell RNA sequencing data. For example, probabilistic graphical models have been introduced for imputing spatial gene expressions by integration with scRNA-seq data [57]. With the integration of PPI network and tensor modeling, FIST has a great potential to achieve better scalability as well as better accuracy for imputation of transcriptome-wide spatial transcipromics data.  Table. p-values of paired-sample t-tests. The means of performance (measured by MAE, MAPE and R 2 as in Fig 2) for predicting all the spot fibers are compared between FIST and each of the baseline methods, using paired-sample t-tests. (XLSX) S2 Table. Enriched GO terms of spatial gene clusters. The GO terms significantly enriched by the genes in each spatial gene cluster (FDR adjusted p-value <0.05) are shown in the spreadsheet tables. (XLSX) S1 File. Convergence of FIST. (PDF)