# B Simulation Study

## B.1 Induced Correlation Study

In this section, we provide more details about the simulation study. First we consider $$3$$ different values of $$\rho$$ $$(0.1,0.15,0.2)$$ and make a induced correlation plot by using the squared exponential spatial kernel. The plots are generented for all cell types and cell type specific cases. The vertical line denotes the value of induced correlation at the distance $$0.01$$. For example the induced spatial correlations for all cell types (first Figure of B.1) w.r.t. $$0.01$$ distance are $$0.88$$, $$0.80$$, $$0.61$$ for $$\text S_{\text High}$$ ($$\rho = 0.2$$), $$\text S_{\text Med}$$ ($$\rho = 0.15$$), $$\text S_{\text Low}$$ ($$\rho = 0.1$$) methods respectively.

Estimation of parameter $$\rho$$ : We adapt the steps from the section 1.1 of supplementary information of S. Sun, Zhu, and Zhou (2020) to estimate the length scale parameter of the spatial kernel $$\rho_{c}$$ based on several grid points. The estimation procedure of $$\rho_{c}$$ uses the pair-wise distances among spatial locations in the data to ensure scale-invariance. First the maximum $$(m_{1})$$ and minimum $$(m_{2})$$ value of those pairwise spatial distances are obtained. Next, equidistant $$L$$ points are identified in the range from $$\text{log}(m_{1}/2)$$ to $$\text{log}(m_{2}*2)$$. In our case, we consider $$L=10$$ as a default value for all our simulations and spatial transcriptomics data applications. We consider the middle value of those $$10$$ grid points as the estimate of $$\rho_{c}$$. Table B.1 shows the true and estimated values of $$\rho_{c}$$ which shows the grid points based estimation procedure is effective.

 True value of $$\rho_{c}$$ 0.2 0.15 0.1 Estimated value of $$\rho_{c}$$ 0.18 0.16 0.11 Standard Errors 0.11 0.13 0.14

## B.2 Data generation with spatial correlation

Simulation Design: The details of the simulation design is provided in section 3 of the paper.

Comparative Models: We discussed the $$5$$ comparative models in section 3 of the paper. A summary of the comparative models is provided in Table B.2

SpaceX model Non-spatial Poisson model Gaussian model
Spatial information $$\checkmark$$ X X
Poisson likelihood $$\checkmark$$ $$\checkmark$$ X
Gaussian likelihood X X $$\checkmark$$
Table B.2: Overview of comparative models.

Metrics for comparison: To measure the co-expression estimation accuracy, we use the following metrics in Table B.3 such as RV coefficient and 4 Euclidean distance based norms (Frobenius, Log-Euclidean, Root Euclidean, Riemanian, Dryden et al. (2009)) as defined in Table B.3. These metrics are used to quantify the similarity between true and estimated covariance (co-expression) matrices. RV values close to $$1$$ ($$0$$) implying higher (lower) level of similarity. For rest of the norms on Table B.3, values closer to $$0$$ indicates higher level of similarity between true and estimated matrices.

 Name $\text{Notation}$ $\text{Form}$ RV $\text{RV}(S_{1},S_{2})$ $\frac{ tr( S_{1}^{T} S_{2} ) }{ \sqrt{ tr( S_{1}^{T} S_{1} ) tr( S_{2}^{T} S_{2} ) } }$ Euclidean (Frobenius) $d_{E}(S_{1},S_{2})$ $\mid \mid S_{1} - S_{2} \mid \mid$ Log-Euclidean $d_{L}(S_{1},S_{2})$ $\mid \mid \log(S_{1}) - \log(S_{2}) \mid \mid$ Root-Euclidean $d_{H}(S_{1},S_{2})$ $\mid \mid S_{1}^{1/2} - S_{2}^{1/2} \mid \mid$ Riemanian $d_{R}(S_{1},S_{2})$ $\mid \mid S_{1}^{-1/2} S_{2} S_{1}^{-1/2} \mid \mid$

### B.2.1 Comparative analysis with different norm measures

Figure B.2 represents the boxplot of distances between true $$(\Sigma_{True})$$ and estimated $$(\Sigma_{Est})$$ covariance matrices where the distances are measured in Euclidean, root Euclidean, log Euclidean and Riemanian norms respectively. In all the norms we observe that spatial settings are performing better in terms of estimation than the no-spatial settings. Among the spatial settings the estimation accuracy increase with an increment in induced spatial correlation.

### B.2.2 Estimation of latent factors

We follow same procedure from section 3.3 of De Vito et al. (2021) to estimate shared and cluster specific number of factors i.e. $$K$$ and $$K_{c}$$ $$(c=1,2, \dots C)$$. Figure B.3 shows shared and cluster specfic estimated factor loadings accross $$50$$ replicates for $$5$$ different methods. Figure B.4 shows the median estimate of shared and cluster specfic factor loadings for $$5$$ different methods. From both figures one can observe that spatial settings are estimating the loadings more precisely than the non-spatial settings.

## B.3 Data generation with no spatial correlation

Simulation design: We consider the same simulation design, comparative models and metrics for comparison. Now, we generate the from $$\text NS_{\text P}$$ model where we do not consider any spatial correlation. We model the spatial parameter with a multi-variate normal distribution with mean $$0$$ and identity as covariance matrix. We generate the data with no spatial correlation model $$\text NS_{\text P}$$ and fit all the comparative models as mentioned in Table B.2. We summarize our results based on $$50$$ replicated simulation study w.r.t co-expression estimation and network recovery.

Co-expression estimation: We display the boxplot of RV coefficients for shared ($$\bf G_{s}$$) and cluster-specific ($$\bf G_{c}$$, $$c=1, \dots, C$$) covariance matrices in Figure B.5A across 3 comparative models in Table B.2. Based on RV coefficients, the highest level of precision in estimation is obtained for the non-spatial setting in (IV) i.e. $$\text NS_{\text P}$$. For spatial settings, we do not observe a significant loss in estimation accuracy. For example, we only loose 1.3% accuracy in estimation for spatial settings with the SpaceX method in case of the shared network. A similar inference can be drawn based on different norm measures (Euclidean, log-Euclidean, root-Euclidean, Riemannian as defined in Table B.3 as provided in Figure B.6 (in clock-wise manner).

Network recovery: AUC based comparisons for shared and cluster specific networks are shown in Figure B.5B. The Figure leads us to infer that no significant reduction in network recovery while comparing the spatial and non-spatial settings while data is generated without spatial correlation.

In summary, we do not observe a significant loss in precision while applying the SpaceX model to the simulated data without any spatial correlation.

## B.4 Hub gene detection based simulation

Given the lack of ground truth for real data, we used a simulation study mimicking the real data structure, to evaluate the accuracy of the SpaceX method to recover true hub genes. To this end, we generate the data consisting of $$160$$ genes and $$5000$$ locations with $$7$$ clusters; additional design parameters and details about the simulation settings are provided in Section 3 of the paper. Next, we apply our SpaceX model to the simulated data to estimate the hub genes for each of the $$7$$ spatial clusters. We classify the hub gene into two categories: one with more than 40% connectivity (named “High”) and another with less than 40% connectivity (named “Low”). We choose this cut-off based on degree connectivity we found in the real data example (Section 4). Table B.4 shows the percentage of hub gene recovered for each cluster by setting $$3$$ different values of the spatial correlation parameter: $$\rho = 0.2, 0.15$$ and $$0.1$$ corresponding to low, medium and high levels of spatial correlation. The denominator of recovery rate is the total number of hub genes calculated from the true simulation settings whereas the numerator denotes the estimated number of hub genes obtained from the SpaceX method. From Table B.4, we observe that the percentage hub gene recovery increases with an increase as the level of spatial correlations increase $$\rho$$ and size of the cluster. For example, cluster 5 (size = 1700) has the following recovery rates 97.5% , 95.2%, 92.3% for highly connected hub genes corresponding to different level of spatial correlations i.e. $$\rho = 0.2,0.15$$ and $$0.1$$ respectively. The recovery rates are 97.5%, 93.2%, 90.3% for cluster 5, 4 and 1 respectively with different cluster sizes ($$N_{5}=1700$$, $$N_{4}=1000$$, $$N_{1}=700$$). As expected, the recovery rate of the hub genes with higher connectivity is more than in the ones with low connectivity.

### References

De Vito, Roberta, Ruggero Bellio, Lorenzo Trippa, and Giovanni Parmigiani. 2021. “Bayesian Multistudy Factor Analysis for High-Throughput Biological Data.” The Annals of Applied Statistics 15 (4): 1723–41.
Dryden, Ian L, Alexey Koloydenko, Diwei Zhou, et al. 2009. “Non-Euclidean Statistics for Covariance Matrices, with Applications to Diffusion Tensor Imaging.” The Annals of Applied Statistics 3 (3): 1102–23.
Robert, Paul, and Yves Escoufier. 1976. “A Unifying Tool for Linear Multivariate Statistical Methods: The RV-Coefficient.” Journal of the Royal Statistical Society: Series C (Applied Statistics) 25 (3): 257–65.
Sun, Shiquan, Jiaqiang Zhu, and Xiang Zhou. 2020. “Statistical Analysis of Spatial Expression Patterns for Spatially Resolved Transcriptomic Studies.” Nature Methods 17 (2): 193–200.