SMASH: Scalable Method for Analyzing Spatial Heterogeneity of genes in spatial transcriptomics data

In high-throughput spatial transcriptomics (ST) studies, it is of great interest to identify the genes whose level of expression in a tissue covaries with the spatial location of cells/spots. Such genes, also known as spatially variable genes (SVGs), can be crucial to the biological understanding of both structural and functional characteristics of complex tissues. Existing methods for detecting SVGs either suffer from huge computational demand or significantly lack statistical power. We propose a non-parametric method termed SMASH that achieves a balance between the above two problems. We compare SMASH with other existing methods in varying simulation scenarios demonstrating its superior statistical power and robustness. We apply the method to four ST datasets from different platforms revealing interesting biological insights.

1 Further details on the SMASH test statistic 1

.1 Choice of kernel covariance matrices
In the main text, we define the SMASH test statistic for a gene k to have the following form, tr(E k H(S)) N ; E k = y k (y T k y k ) −1 y k ; (1) where H(S) is any N × N kernel-based covariance matrix based on the locations S and y k is the mean-standardized expression vector of the gene k. As mentioned in the main text, when additional covariates are present, we would replace y k by y * k = [I − P X ]y k ; P X = X(X T X) −1 X T , where X is the design matrix. We consider H(S) to have following three forms: 1. The Gaussian kernel-based covariance matrix: for ten values of the lengthscale parameter l.
2. The cosine kernel-based covariance matrix of the form: for ten values of the period parameter p.
3. The linear kernel-based covariance matrix of the form: where three choices of the coordinate-wise transformation g are considered: • An identity transformation i.e., g(s i1 ) = s i1 , g(s i2 ) = s i2 for i = 1, . . . , N , or, g(S) = S.
The Gaussian and cosine kernel-based covariance matrices have earlier been used in SpatialDE [5] and SPARK [4]. We follow SpatialDE to choose a set of fixed grid points for both the lengthscale parameter l and the period p. In particular, we first obtain the minimum (d min ) and the maximum (d max ) values of the non-zero Euclidean distances across all pairs of spatial locations.
Then, we extract ten equally-spaced values between log 10 (d min /2) and log 10 (d max /2). These values are then converted to the original scale by taking the power of ten and subsequently used as the values of l and p. Let us denote the test statistics corresponding to these twenty covariance matrices (10 each for the Gaussian and cosine kernels) as, T SMASH kr , r = 1, . . . , 20. The linear kernelbased covariance matrices with transformed coordinates (g-transformation) have been earlier used in SPARK-X [6]. Following SPARK-X, we vary the transformation parameters, t 1 , t 2 , ϕ 1 and ϕ 2 to be the 20%, 40%, 60%, 80%, and 100% quantiles of the absolute values of the x and y coordinates in the data. Let us denote the test statistic corresponding to these eleven (1 for the identity transformation and five each for the Gaussian and cosine transformations) linear kernel-based covariance matrices as, T SMASH kr , r = 21, . . . , 31. We denote the p-values corresponding to each T SMASH kr as p kr .
Then, we can get SPARK-X's result by combining the p-values, p kr , r = 21, . . . , 31 using a Cauchy combination rule [2]. Note that the same Gaussian and cosine kernel-based covariance matrices (r = 1, . . . , 10) are also used in SpatialDE [5] but p k,f inal = 2 * min{p k,comb1 , p k,comb2 }. We do acknowledge that using two levels of minimum p-value combination can be conservative in some cases, which is partly observed in the null simulation studies performed using the real datasets in Section 3. However, our focus was on preventing the high degree of inflation in the p-values routinely observed in model-based methods like SpatialDE and SPARK.

Similarity with SpatialDE
From Equation (2) of the main text, recall that SpatialDE considers the following model, SpatialDE optimizes two likelihood functions: one corresponding to the 'full' model shown above and the other corresponding to the reduced 'null' model without the τ 2 k Σ component in the covariance.
Then, to test the null hypothesis, H 0 : τ 2 k = 0, it constructs a likelihood ratio test (LRT) comparing the optimized values of the above two likelihood functions. If y k is mean-standardized i.e., µ k = 0, the corresponding score test statistic [1] can easily be shown to have the form, (1).
Thus, we have derived a rough similarity between the SMASH test statistic with the framework of SpatialDE.

SPARK-X's equivalence with a multiple linear regression model
Let us consider the following linear regression model with y k as the dependent variable and the columns of the spatial location matrix S (x and y coordinates) as the predictors, The fixed effect coefficients vector β k is of length 2 and its OLS estimate has the form,β k = where D = S(S T S) −1 S T from the main text. To test the null hypothesis, H 0 : β k = 0, we would consider the following Wald test statistic, Thus, we have shown that the SPARK-X test statistic T SPARKX k is proportional to the simple Wald test statistic obtained from a multiple linear regression model. However, the asymptotic distributional assumptions under null are different in both cases. Note that SPARK-X also replaces S by g(S) where g is a coordinate-wise transformation function as discussed in Section 1.1.

QQ-plots under null simulations
As mentioned in the main text, we performed null simulation studies to construct an empirical null distribution of the p-values for every method. With each of the four datasets, we randomly permuted the cell/spot coordinates five times and then applied the three methods, SMASH, SPARK-X, and SpaGene to obtain the respective p-values. The p-values were then transformed as − log 10 (p-value) and displayed against the expected values as quantile-quantile plots (QQ-plots) in Figure S1. We noticed that SMASH showed no sign of p-value inflation and was rather slightly conservative. It is expected since the minimum p-value combination rule we use, is known to be conservative [3] (see Section (1.1)). SpaGene produced slightly inflated p-values in the SCCOHT dataset while SPARK-X did not show any sign of inflation.