Density estimation and function approximation is a fundamental problem in statistics and machine learning. Non-parametric methods such as spline regression ( Perperoglou et al., 2019), Fourier transformation ( Cochran et al., 1967) and Wavelet transformation ( Nason, 2008) have been used in such scenarios.

Consider a model with one predictor: y = f( x) + ϵ, where E ϵ 2 = σ 2 . We want to estimate f, the trend of the function. We assume the predictors are ordered x ₁ < x ₂ < … < x _n. We have a signal or frequency to estimate.

Consider the Haar mother wavelet: ψ x = 1 x ∈ 0 1 2 − 1 x ∈ 1 2 1 0 otherwise satisfying ∫ − ∞ ∞ ψ x dx = 0 . Unlike the Fourier basis, wavelets oscillate and decay fast, only contributing to a certain local area and zero elsewhere. We can generate wavelets from the Haar mother wavelet ψ j , k x = 2 j / 2 ψ 2 j x − k for integer j, k. These wavelets form orthonormal sets. We can decompose the trend as f x = ∑ j = − ∞ ∞ ∑ k = − ∞ ∞ d j , k ψ j , k x , where d j , k = ∫ − ∞ ∞ f x ψ j , k x dx = < f , ψ j , k > are the wavelet coefficients.

Wavelet based methods have several advantages among other non-parametric methods, especially dealing with sparse data. One property wavelet has is localization. If a sequence has a discontinuity, this will only influence the wavelet basis around it. In contrast, for a Fourier basis consisting of sine and cosine functions at different frequencies, every basis element will interact with this discontinuity, hence influencing every Fourier coefficient.

The simplest discrete wavelet transformation calculates the difference and sums between each adjacent pair. Suppose we have one vector with length n, where n is dyadic ( n = 2 ^J ). We computed d _{J−1, k} = y _{2 k} − y _{2 k−1} as finest-level detail and c _{J−1, k} = y _{2 k} + y _{2 k−1} as the finest-level averages. We have d k containing n/2 = 2 ^J−1 as our finest level coefficients. To obtain the next coarsest coeffients we set: d J − 2 , ℓ = c J − 1 , 2 ℓ − c J − 1 , 2 ℓ − 1 (1) c J − 2 , ℓ = c J − 1 , 2 ℓ + c J − 1 , 2 ℓ − 1 (2)

We continue this procedure until j reaches one. The set of details d k across levels are our wavelet coefficients.

We use J to denote the scale of wavelets, with larger J we have finer scale and approximation. However, finer scale sometimes introduces more parameters to capture minor details of the sequence (overfitting). The trade-off in choosing J is crucial in the wavelet transformation.

Regularization and smoothing can be used to prevent overfitting. Regularization is usually conducted by shrinking wavelet coefficients. The concept of wavelet shrinkage was first proposed by Donoho and Johnstone (1994). The motivation behind wavelet shrinkage is straightforward. Consider empirical wavelet coefficients, the large coefficients usually contain true signal and noise, whereas the small coefficients only contain noise. The shrinkage is often used when the wavelet coefficients are assumed to be a sparse vector.

Wavelet shrinkage is often conducted by setting a threshold (denoted by δ) and only keeping the coefficients above the threshold. To choose a threshold, a natural metric is the squared error between estimated function and the truth: M ̂ = ∑ i = 1 n f x i − f ̂ x i 2 and M = E M ̂ . Donoho et al. (1994) proposed a universal threshold δ = σ 2 log n , which induced M ≤ O log n σ 2 . Donoho and Johnstone (1995) also proposed Stein’s unbiased risk estimation (SURE) threshold based on Stein’s (1981) unbiased risk estimator. The optimal SURE threshold can be obtained in O n log n operations. Donoho and Johnstone (1995) also noted that SURE sometimes failed when the true signal coefficients are highly sparse, they proposed a hybrid scheme called SureShrink, combining the SURE and universal thresholds, using them depending on certain situations.

The extension of wavelet methods to 2D regularly spaced data (images) and such data in higher dimensions was proposed by Mallat (1989). We only consider 2-D wavelet transformation since we decompose 2-D spatial gene expression data. Suppose we have n × n matrix A where n = 2 ^J is dyadic. A simple discrete wavelet transformation on A first applies procedure (1) and (2) to the rows of the matrix. We then have two matrices of size n × n 2 , called H and G. Then we apply the same procedures to both the columns of H and G, resulting in four matrices HH, GH, HG, and GG each of size n 2 × n 2 . These are our finest level coefficients. HH is the local average of the original matrix used for the next level procedure.

Clustering and dimensionality reduction are widely used in genomics. In a gene-by-sample matrix, genes are often grouped into profiles, where genes from the same profile have a similar function. Statistically, we can treat them as correlated variables. We use the term factor genes as the principal component in our gene-by-sample data. Factor genes compose the linear combinations of genes. We consider each gene as a variable, and we aim to find variables ( f ₁, …, f _K ) such that each gene g _i has form g _i = a ₁ f ₁ + a ₂ f ₂ + ⋯ a _K f _K , where a _j are coefficients.

Suppose we have gene( P)-by-sample( N) matrix A ∈ ℝ N × P with sample covariance matrix S = 1 n à T à , where à is the column centered A. Consider the SVD on à = UΛ V ^T. Then we have SV = VΛ ². The columns of V are called eigenarrays ( Alter et al., 2000). The first K columns of Z = AV are the factor genes. The factor genes capture the mutual underlying information of genes.

Matrix Factorization is often used in capture factor genes. We have a formulation similar to Wang and Stephens (2021), consider the factorization model on observed samples by gene data Y ∈ ℝ N × P Y = FL T + E where F ∈ ℝ N × K denotes the factors, L ∈ ℝ P × K denotes the loadings for each factor, and E ∈ ℝ N × P denotes Gaussian noise with zero mean. Typical formulations would treat factors and loadings as fixed effects and use Maximum Likelihood Estimation (MLE). In our high-dimensional setting, penalty-based regularizers are often considered. Considering a prior for L and F under Bayesian setting has a similar effect and several advantages over adding a regularizer alone, such as simplifying hyperparameter search and selection of the number of factors K.

One feature of empirical Bayes approaches in matrix factorization proposed by Bishop (1999) is that the methods automatically select the number of factors K. In Wang and Stephens (2021), the author add factors with prior one at a time, estimating priors at each step until convergence. If the computed prior of the newly added factor is almost point mass on 0, the algorithm eliminates this factor and returns.

Following Wang and Stephens (2021), we consider Empirical Bayes in our setting, where we set a prior with unknown parameters of L and F. Empirical Bayes is not strictly Bayes, since the prior parameters are directly estimated by the MLE of the data. One example would be normal distribution D F of F and D L of L where coordinates are independent, which is conjugate prior. We estimate parameters in D F and D L by the maximize the marginal likelihood calculated by integral out L and F. Then we computed the posterior distribution of L and F.

We leverage a pre-processing step, which we call input generation, to combine spatial and expression data. We have both gene expression measurements and spatial coordinates for each sample. Each gene has an expression over each sample, and we can draw the sample over a 2D map by their coordinates.

As shown in Figure 2, the gene with the least sparsity in expression shows varying levels of expression across different spatial regions. Intuitively, we may consider the spatial expression pattern from this gene to be spatially related. However, most genes are sparse and do not show fluctuations in gene expression. Therefore, we first filter out those genes using a kOverA filter: we only keep genes that have an expression measure above A in at least k samples ( Gentleman et al., 2021). This has the effect of removing genes that are rarely active, though it is possible that strongly expressed genes still show no spatial relationships.

We want gene expression data to have a neat grid structure that can be expressed as a matrix form. However, at first, expression measurements are roughly staggered over the spatial locations, as shown in Figure 3. This is a consequence of the sampling strategy adopted by the 10x genomics Visium platform. To obtain an evenly resampled version of the expression pattern, we divide the two-dimensional space into several partitions and compute the local average of each partition. Detailes are shown in Algorithm 1. Each gene has an expression over the grid and we now can use a matrix Z ∈ ℝ D 2 × P to represent it. The matrix Z becomes the new input for analysis. We then obtain the same formulation as in Section Problem Setup, where we take n = D ² and each location corresponds to a new observation. We have f k ∈ ℝ D 2 , k = 1 , … , K as the factor genes, computed by vectorizing the gene expression matrix.

The choice of D depends on the wavelet scale used in the next section. We choose D to be dyadic to prevent handling edge cases in the wavelet transformation (recalling Subsection Wavelet Transformation). The dyadic size is also a natural choice in image analysis. In particular, we choose D = 2 ^J , where J is the level of scale in the wavelet transformation. Larger J allows finer recovery, but also requires more parameters and may result in overfitting.

We apply a wavelet transformation with shrinkage to denoise the gene expression matrix and smooth observed spatial expression patterns. In simulations, we find that this technique gives more accurate recovery as the signal-to-noise ratio decreases and yields less noisy visualizations for the spatially-related genes.

We implement wavelet transformation and shrinkage on the processed data matrix Z ∈ ℝ D 2 × P . Each column of Z is a vectorized expression matrix. For each column, we first reshape it into a D × D matrix and apply the wavelet transformation. We have a coefficient list with W coefficients associated with each gene. Then we conduct wavelet shrinkage on the coefficient list using the threshold strategies described in Subsection Wavelet Transformation. Finally, we vectorize the coefficient list of each gene into a length W vector for each gene. Stacking all vectors together, we have a coefficient matrix C ∈ ℝ W × P . Note that the wavelet scale is specified by the size D of the input matrix, i.e., we have D = 2 ^J , where J is the scale. The details are specified in Algorithm 2.

The matrix of the shrunk wavelet coefficients C gives a summary of the denoised matrix. This allows improved reconstruction of gene expression data. To obtain a low-dimensional approximation, we apply matrix factorization on the coefficient matrix C after transformation. We use the same notation in Section Problem Setup, with the subscript c to denote the decomposition on coefficient matrix C.

The resulting singular vectors can be used to estimate spatially structured factor genes (Subsection Factor Gene). We use SVD as a frequentist approach to estimating F c ̂ and L ̂ c . We also conduct EBMF (Section Empirical Bayes Matrix Factorization), using the posterior expectation of F _c and L _c to estimate F c ̂ and L c ̂ . EBMF can select the number of factors K by itself, whereas K must be manually specified in the SVD. The choice of K is informed by inspecting the scree plot, as in spectral clustering and PCA. We choose the number of factors when current factors explain sufficient information and the diminishing returns of additional factors are no longer worth the additional cost. Given F ̂ c and L ̂ C , we can compute the estimated coefficient matrix as C ̂ = F ̂ c L ̂ c T .

To transfer the coefficient matrix C ̂ back to the estimated location-by-gene matrix Z ̂ , we apply the inverse wavelet transformation on columns of C ̂ , each of which is a vectorized coefficient list for one gene. This results in a D ² expression matrix for each gene. Then we vectorize the matrix and stack all vectors together. This yields the reconstructed matrix Z ̂ . Details are given in Algorithm 3. For visualization, we also conduct a similar process for each column of F c ̂ . In this case, each column is associated with an factor gene. By applying the inverse wavelet transformation to each factor gene, we can build spatial gene expression matrices M ₁, …, M _K representing gene factors.

We evaluate our proposal both quantitatively and qualitatively. Our quantitative results measures reconstruction error between the estimated Z ̂ and Z using the Frobenius norm Z ̂ − Z F 2 .

This evaluation is also used for hyperparameter tuning, such as the scale J of the wavelet and the choice of wavelet thresholding method. Our qualitative result is given by visualization of the estimated factor gene expression matrices, with emphasis on capturing global spatial structure.

We use cross-validation in computing reconstruction error and calculating gene-wise errors. This evaluation is helpful in selecting spatially related genes. We also found a simple connection between the gradient of gene expression and spatial contribution, discussed in Section Real Data Experiment. This discussion requires a measure of spatial expression smoothness. To this end, a simple step computing the fluctuation of gene expression would give a spatial gene selection that coincides with the reconstruction error selection. For calculating successive differences, consider Z as the input matrix and let δ _jk = ( z _{j, k+1}− z _jk ) ², k = 1, …, D ². Then we define the gradient by ∑ jk δ jk .

We compute the reconstruction loss, 1 N ∑ i , j ∈ test Z ̂ ij − Z ij 2 where N is the number of test entries. We evaluate the loss only on test entries. We first tune parameters. We set up three settings: decompose the raw data with SVD or EBMF; wavelet transformation with hybrid thresholding; wavelet transformation with manual thresholding with threshold τ = 10, 20, …, 100. In each setting, we ran 100 replicates. The reconstruction error shown in Figure 8.

As shown in the upper panel of Figure 8, wavelet thresholding with manually set τ reduces error compared to decomposing raw data. The wavelet has the most positive effect when τ = 40, as shown in the bottom panel. We use τ = 40 in our following analysis.

We then evaluate how method performance varies for each gene by calculating the genewise reconstruction error. This reveals genes whose strong spatial expression structure leads to improved performance when using a wavelet basis. We still hold out 1 5 of the entries at random as a test set. We compute the reconstruction error of test entries on each column of Z ̂ − Z . We calculate entry-wise loss across 100 replicates and estimate the average loss. We compare genewise errors with or without wavelet transformation. The SVD and EBMF show the same result regarding the decomposition method. We only show result of EBMF here in Figure 9. We show the result of SVD in Figure 13b.

As we can see, most of the genes have lower reconstruction error when directly applying the SVD or EBMF. However, for some genes, wavelet smoothing reduces error. For example, this is seen in genes 712, 713, 716, 719, 715, 710, 711, 721, 717, 596, 373, and 289. We conjecture that these genes have acute fluctuations as well as clear spatial patterns. The expression matrix of these genes would have a larger gradient and distinct edges. To verify our conjecture, we calculated the sum of squares of the gradient of each gene expression matrix. We show the result in Figure 10.

The result verifies our conjecture: genes with larger gradients have better reconstruction under the wavelet-guided decomposition. This suggests a pre-processing step for selecting wavelet-suited genes by calculating the gradient of each gene. The spatially related genes would have a lower reconstruction error. Among these spatially related genes, one possibility is that we can divide them into two groups, one for decomposition directly and the other for the wavelet technique.

Now we visualize the top factor genes and genes with high loadings on these factor genes. In Figure 11a, we decompose Z using the SVD and plot the first factor gene and matrix slides of genes with the largest 5 loadings on that factor. The factor gene captured the same patterns as the original genes. To improve visualization, we find the analogous wavelet-based factors (we use manual thresholding with τ = 40), shown in Figure 11b.

The wavelet thresholding approach smooths over edges in the original visualization. The first factor gene captures the global spatial expression. We have a similar result for the second factor gene, as shown in Figure 12.

The second factor gene is orthogonal to the first one, capturing different spatial structures in gene expression. Different factor genes capture global and local structure, and using a wavelet decomposition provides denoised spatial expression visualization.

In conclusion, we can use this dimensionality reduction technique for spatial gene selection and extraction. We select wavelet-suited genes based on the calculated gradient. Then, we can select spatially related genes based on reconstruction error via cross-validation. Alternatively, we can extract factor gene to capture spatial information and use factor genes for visualization and further analysis.

We provide a small code block to show vignette of generating Figure 11 and Figure 12, shown in code block 1.

We have made a new package waveST containing the workflow we developed. The package is available at GitHub (see Software availability). The kOverA_ST function reduces data from an original ( Weber, 2021) class input using the kOverA technique. Then we use waveST function to construct a S4 class waveST, containing input generated by Algorithm 1. This object-oriented approach stores properties of the original spatial experiment and simplifies downstream calls, like decomposition and visualization. We use the decompose function to apply all decomposition methods. In line 6, we apply the SVD to our original data, setting the number of factors to 5. In line 11, we apply the wavelet-based reduction and apply a manual threshold with τ = 40. We use plot and k=1 to visualize the first factor gene and matrix slides of the genes with the largest 5 loadings on that factor.

This section we provide a comparison between the reconstruction results when using SVD and EBMF. In general, we find little difference between reduction using the two methods. We first provide an element-wise comparison between the reconstruction of SVD and of EBMF with wavelet technique in Figure 13a.

The results are similar when applied to real data. In particular, Figure 13b shows the reconstuction error per gene result for SVD. Our findings are similar to those for Figure 9. Similar findings are visible when using SVD and EBMF with and without an initial wavelet decomposition Figure 14. As we can see, SVD and Empirical Bayes Matrix Factorization (EBMF) have the very similar result. This is perhaps a consequence of EMBF being initialized usign the SVD.

Zenodo: OliverXUZY/waveST: waveST. https://doi.org/10.5281/zenodo.6983923 ( Zhuoyan, 2022)

This project contains the following underlying data:

raws.rda

simulation_data.rda

Data are available under the terms of the Creative Commons Zero “No rights reserved” data waiver (CC0 1.0 Public domain dedication).