{'fraction_non_alphanumeric': 0.03778588807785888, 'fraction_numerical': 0.016374695863746958, 'mean_word_length': 4.488182667912939, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Bayesian Modified Ising Model for Identifying Spatially Variable Genes from Spatial Transcriptomics Data

A recent technology breakthrough in spatial molecular profiling has enabled the comprehensive molecular characterizations of single cells while preserving spatial information. It provides new opportunities to delineate how cells from different origins form tissues with distinctive structures and functions. One immediate question in spatial molecular profiling data analysis is to identify genes whose expressions exhibit spatially correlated patterns, called spatially variable genes. Most current methods to identify spatially variable genes are built upon the geostatistical model with Gaussian process to capture the spatial patterns, which rely on ad hoc kernels that could limit the models' ability to identify complex spatial patterns. In order to overcome this challenge and capture more types of spatial patterns, we introduce a Bayesian approach to identify spatially variable genes via a modified Ising model. The key idea is to use the energy interaction parameter of the Ising model to characterize spatial expression patterns. We use auxiliary variable Markov chain Monte Carlo algorithms to sample from the posterior distribution with an intractable normalizing constant in the model. Simulation studies using both simulated and synthetic data showed that the energy-based modeling approach led to higher accuracy in detecting spatially variable genes than those kernel-based methods. When applied to two real spatial transcriptomics datasets, the proposed method discovered novel spatial patterns that shed light on the biological mechanisms. In summary, the proposed method presents a new perspective for analyzing spatial transcriptomics data.


I. INTRODUCTION
Cellular and molecular spatial organizations play essential roles in biological functions.
In recent years, spatial molecular profiling (SMP) techniques have made significant breakthroughs, which enable transcriptome measurement in high spatial resolution (Zhang and others, 2020). Gene expression profiling approaches are no longer limited tissue-dissociation, which led to the loss of spatial context of the measured gene expression (Femino and others, 1998). Sequencing-based spatial molecular profiling platforms, such as spatial transcriptomics (ST) (Ståhl and others, 2016) and the improved Visium platform, use spatial barcodes to capture RNA molecules and then synthesize and sequence their complementary DNA molecules. Through this technology, the expression levels of thousands of genes can be measured across hundreds of spatial locations, namely spots. Those spots are usually arrayed on a two-dimensional grid. Particularly, ST and Visium spots are arranged on the square and triangular lattices, respectively. SMP techniques enable researchers to study the gene expressions together with their spatial and morphological contexts, which provide new opportunities to advance our understanding of both cellular and molecular spatial organizations (Crosetto and others, 2015), and their relationships with diseases (Shah and others, 2018). Many new questions can be explored with the emerging SMP techniques. One of the most immediate ones is to identify genes whose expressions exhibit spatially correlated patterns, referred to spatially variable (SV) genes. The study of spatial patterns in gene expression could reveal significant insights into many aspects, such as embryo development (Satija and others, 2015), tumor progression (de Bruin and others, 2014), and the clinical impact of intra-tumor heterogeneity (Bedard and others, 2013).
Several methods have been developed in recent years to address the above fundamental question in ST studies. Trendsceek (Edsgärd and others, 2018) is based on marked point processes, and it is computationally intensive (Dries and others, 2021) and has noticeably unsatisfying performance (Sun and others, 2020). BinSpect (Dries and others, 2021) is an easy and fast computational method based on statistical enrichment of spatial network neighbors after binarizing gene expression levels. Most model-based analyses, such as SPARK (Sun and others, 2020), SpatialDE (Svensson and others, 2018), and BOOST-GP (Li and others, 2021), are built upon the geostatistical model with Gaussian process (GP), where a kernel must be selected with caution. Among these methods, SpatialDE transforms the measured counts at different locations into normalized data before analysis, SPARK uses a Poisson distribution to model the count data directly, and BOOST-GP accounts for the zero-inflation and mean-variance relationship existed in raw counts under a Bayesian framework. However, the GP-based models rely on ad hoc kernels that limit the models' ability to identify complex spatial patterns. Furthermore, none of the existing approaches take advantage of the additional spatial structure of the ST data; that is, the gene expression levels are measured on a lattice grid.
To enable the model to identify complex spatial patterns, we developed a novel approach named Bayesian mOdeling Of Spatial Transcriptomics data via a Modified Ising model (BOOST-MI) to identify SV genes in ST studies. It makes computation efficient by taking advantage of the fact that the ST experiments measure gene expression on a lattice.
Ahead of fitting BOOST-MI, we need to normalize sequence count data to relative gene expression levels and dichotomize relative expression levels to a binary spatial pattern. Then, BOOST-MI characterizes the binary spatial pattern via inferring the Ising model interaction parameter under a Bayesian framework. The double Metropolis-Hastings (DMH) algorithm (Liang, 2010) is used to sample from the posterior distribution with an intractable normalizing constant in the Ising model. Compared with other existing approaches, BOOST-MI tests the interaction energy parameter in a modified Ising model, which is able to characterize a broader type of spatial patterns than kernel-based modeling approaches. In addition, the proposed method enriches the inference via incorporating priors if necessary and naturally quantifying uncertainties under a Bayesian framework. We demonstrated the advantages of BOOST-MI in a comprehensive simulation study using both simulated data with various spatial patterns and zero-inflation settings and synthetic data from real ST experiments.
The proposed model showed an outstanding performance compared to the alternatives. Finally, when applied to two real ST datasets, the proposed method discovered novel spatial patterns that shed light on the biological mechanisms.
The remainder of the paper is organized as follows. Section 2 introduces the two main components of the data preparation for BOOST-MI: sequence count data normalization and relative expression level dichotomization. In Sections 3 and 4, we introduce the proposed modeling framework and describe the Markov chain Monte Carlo (MCMC) algorithm and the resulting posterior inference. In Section 5, we compare BOOST-MI with existing approaches on simulated and two real ST datasets. Section 6 concludes the paper and discusses future research directions.

II. DATA PREPARATION
In this section, we introduce the two data preprocessing steps (see the dashed-line block in Figure 1) before fitting BOOST-MI to identify SV genes. We first summarize the ST data notations. Let an n × p matrix Y denote the gene expression count table (i.e. the molecular profile). Each entry y ij ∈ N, i = 1, . . . , n, j = 1, . . . , p is the read count for gene j collected at spot i. We use an n × 2 matrix T to represent the geospatial profile, where each row t i = (t i1 , t i2 ) indicates the spot location on the two-dimensional space. All nonboundary spots have the same neighborhood structure. Suppose the gene expression levels are measured on a L-by-W square lattice grid, the coordinate of each spot i can be written as (t i1 = l, t i2 = w), l, = 1, . . . , L, w = 1, . . . , W . If it is not a boundary spot, then its four neighboring spots are at locations (l ± 1, w) and (l, w ± 1).
Normalization is critical to sequence count data analysis. To counteract various artifacts and bias due to biological and technical reasons, we convert each read count to its relative gene expression level, denoted byỹ ij = y ij /s i , where s i is the size factor of sample i, capturing all nuisance effects. The most straightforward way is to set s i ∝ Y i = p j=1 y ij , i.e. the total number of counts across all genes in each sample (known as sequencing depth or library size), combined with some constraint, such as n i=1 s i = 1. Note that SPARK suggests this normalization, namely total sum scaling (TSS). We could also consider other estimations on s i 's, which have been introduced for mitigating the influence of extremely low and high counts when analyzing bulk RNA-seq data, such as upper-quartiles (Q75) (Bullard and others, 2010), relative log expression (RLE) (Anders and Huber, 2010), and weighted trimmed mean by M-values (TMM) (Robinson and Oshlack, 2010). In addition to the above normalization method based on size factor estimation, we could use the method based on variance-stabilizing transformation (VST) (detailed in Section S1 in the supplement). A following sensitivity analysis (see Figure S1 in the supplement) indicates that 1) the methods based on size factor estimation noticeably outperformed the VST-based methods; 2) the performance based on TSS was not significantly different from Q75, RLE, and TMM. Thus, we suggest using TSS as the default setting for the sake of simplicity.
We then denoise the relative expression levels by partitioning all spots into two groups.
This step outputs the suitable data type required in the subsequent analysis and makes our approach more robust in the face of over-dispersion and zero-inflation, which are the two important characteristics of ST data. For each gene, we introduce a binary vector p j = (p 1j , . . . , p nj ) to represent the dichotomization result based on its relative expression levels (ỹ 1j , . . . ,ỹ nj ), with p ij = 1 indicating gene j is highly expressed at spot i and p ij = 0 otherwise. Since there is no consensus on the dichotomization of spots based on either absolute or relative expression level, we propose to estimate p j via fitting a two-component Here µ j = (µ j0 , µ j1 ) and σ 2 j = (σ 2 j0 , σ 2 j1 ) are the group means and variances that need to be estimated. In addition to the model-based clustering, we also consider k-means, which is implicitly based on pairwise distances between relative expression levels. Section S2 in the supplement provides details about the GMM and k-means implementations (k = 2). A following sensitivity analysis (see Figure S2 in the supplement) indicates that BOOST-MI performed equally well between the two choices.

III. MODELS
In this section, we first review the Ising model and then introduce a modified Ising model with external fields to identify SV genes. The schematic diagram is shown in Figure 1.

A. A brief review of the Ising model
The Ising model, first introduced by (Lenz, 1920) and used in statistical mechanics, is a model of interacting binary states on a crystalline lattice (Cipra, 1987). Let G = (V, E) denote a graph with a finite set of vertices V and a set of edges E. In statistics, the Ising model is considered as an undirected graph such that each vertex is geometrically regular assigned on a lattice and each edge is of the same length (Tucker, 1994). There are three typical types of two-dimensional lattices: square, triangular, and honeycomb (see Figure 2).
The former two are a special case of Bravais lattices, which can be defined as an infinite array of discrete points generated by a set of discrete translation operations (Ashcroft and others, 1976). ST and Visium (an improved ST platform) spots are arranged on square and triangular lattices, where each non-boundary spot has four and six neighbors, respectively.
Every vertex will be assigned a binary state (also known as a spin where i ∼ i denotes the collection of all neighboring spot pairs, θ j denotes the interaction energy between highly and lowly expressed spots for gene j, and I(·) denotes the indicator function. According to the Hammersley-Clifford theorem (Clifford, 1990), we can write the probability of observing a particular configuration of p j as Pr(p j |θ j ) = exp(−H(p j |θ j ))/ p ∈P exp(−H(p |θ j )), where P denotes the set of all configurations of spins. An exact evaluation of the normalizing constant (i.e. the denominator of the above equation) requires us to sum over the entire space, consisting of 2 n configurations. Thus, Pr(p j |θ j ) is intractable even for a small lattice size.

B. A proposal of using a modified Ising model to identify SV genes
Although Ising models have a wide range of applications in many areas, they are often examined assuming that the underlying abundance of different spins is equivalent. In the context of statistical mechanics, we can view the model as a system of interacting particles in the absence of an external field. However, we found that low-expression spots are usually much higher than high-expression spots due to excessive zeros presented in the ST data.
For instance, the averaged proportions of low-expression spots in the two processed real ST datasets analyzed in the paper were 78.6% and 86.8%, respectively (see Figure S3 in the supplement). Ignoring this feature could decrease the identification accuracy due to the inclusion of more false positives (see Figure 4).
To remedy the oversimplified Ising model, we add an external force in the Hamiltonian, where ω j = (ω j0 , ω j1 ) and θ j represent the first and second-order intensities. The first term is proportional to the number of neighboring spot pairs with different spins, while the remaining part (in the square bracket) can be viewed as the weighted average of the numbers of spots with different spins. As the energy function in Equation (2) is still locally defined, we can write the joint probability on p j , up to its normalizing constant, as which serves as the full data likelihood of the proposed BOOST-MI model. To interpret ω j and θ j , we give the conditional probability of observing a high-expression level of gene j at spot i, where Nei(i) denotes the set of all neighbors to spot i. Equation (4) is essentially a logistic regression, and hence the parameters ω j and θ j can be interpreted in terms of conditional odds ratios in general. According to Equation (4), if θ j = 0, then Pr (p ij = 1|·) = exp(ω j1 )/ exp(ω j0 + ω j1 ), implying that each dichotomized expression level p ij is independently and identically sampled from a Bernoulli distribution. Thus, no spatial pattern should exhibit, and gene j is a non-SV gene. The underlying abundance of low and high-expression levels in p j are characterized by exp(ω j0 )/(exp(ω j0 ) + exp(ω j1 )) and exp(ω j1 )/(exp(ω j0 ) + exp(ω j1 )), respectively. Fixing ω j , Equation (4) reveals that the smaller the θ j (i.e. θ j → −∞), the more likely the dichotomized expression level at any spot is concordant with the majority of its neighboring spots' levels, resulting in a repulsion pattern (i.e. the clustering of spots with the same dichotomized expression level). In contrast, when θ j takes a large positive value, we expect an attraction pattern; that is, the exhibition of clustering among spots with different dichotomized expression levels. Thus, the spatial correlation between the low and high-expression levels of gene j can be quantified by θ j . Figure 1 shows the three typical patterns in terms of θ j conditional on ω j0 = ω j1 .
In conclusion, BOOST-MI uses θ j to characterize the binary spatial pattern defined by the dichotomized gene expression levels p j on the lattice grid T . It is noteworthy that the kernel-based methods, such as SPARK, SpatialDE, and BOOST-GP, are only able to identify SV genes with a small subset of repulsion patterns defined by the selected kernel, while BOOST-MI accounts for all attraction and repulsion patterns.
To complete the model specification, we impose θ j ∼ N(0, σ 2 θ ). As for the first-order intensity ω j , we notice that an identifiability problem arises from Equation(4). For example, adding a nonzero constant c into ω jk , k = 0, 1 does not change the conditional probability Pr (p ij = k|·). Thus, we force ω j1 = 1 and set a normal prior on ω j0 ∼ N(1, σ 2 ω ).

IV. MODEL FITTING
In this section, we describe the MCMC algorithm for model fitting and the posterior inference. Our inferential strategy allows for simultaneously estimating the first-order intensity ω j , which reveals the underlying abundance of the low and high-expression levels of gene j, and the second-order intensity θ j (also known as the interaction parameter), which captures the spatial correlation between the low and high-expression levels. We give the details of our MCMC algorithm and the resulting posterior inference as below. Note that each gene is tested independently by BOOST-MI.

A. MCMC algorithm
The full data likelihood function is given in Equation (3), which involves an intractable . Taking the two real ST datasets analyzed in the paper as examples, it needs to sum over 2 250 ≈ 1.8 × 10 75 and 2 260 ≈ 1.9 × 10 78 elements, respectively. This makes the Metropolis-Hastings (MH) algorithm infeasible in practice. To overcome this challenge, we employ the double MH (DMH) algorithm (Liang, 2010) to estimate both ω j and θ j for each gene. The DMH is an asymptotic algorithm, which has been shown to produce accurate results by various spatial models (Li and others, 2019a,b). Unlike other auxiliary variable MCMC algorithms (Møller and others, 2006;Murray and others, 2012) that also aim to have the normalizing constant ratio canceled, it is more efficient because it does not require drawing the auxiliary variables from a perfect sampler, which is usually computationally expensive.
To update θ j within each iteration, we first simulate a new sample θ j from π(θ j ) using the MH algorithm starting with θ j . Then, we generate an auxiliary variable p j through m MH updates starting with the current state p j based on the new value θ j and accept it with π(θ j ) . If the auxiliary variable p j is accepted, we update θ j to θ j ; otherwise, we keep the value of θ j . No improvement in the performance was noticed beyond m = 5 in both simulation and application studies of this paper.
Just as in updating θ j , we use the DMH algorithm to update ω j0 . Specifically, we first simulate a new sample ω j = (ω j0 , 1) from π(ω j0 ) using the MH algorithm starting with ω j0 . Then, we generate an auxiliary variable p j through m MH updates starting with the current state p j based on the new value ω j and accept it with probability min(1, R ω j ), where π(ω j0 ) . If the auxiliary variable p j is accepted, we update ω j0 to ω j0 ; otherwise, the value of ω j0 remains the same.

B. Posterior inference
Our primary interest lies in the identification of SV genes via making inferences on the interaction parameter θ j . We obtain the posterior inference by post-processing of the MCMC samples after burn-in. To validate if gene j exhibits a repulsion pattern, we set the null and alternative hypotheses as M 0 : θ j ≥ 0 and M 1 : θ j < 0; while testing an attraction pattern, we set M 0 : θ j ≤ 0 and M 1 : θ j > 0. We could select the model via calculating the Bayes factor (BF) in favor of M 1 over M 0 , which is defined as the ratio of posterior odds to prior odds, where the prior odds cancel out as we choose a normal prior on θ j centered at zero, and the posterior odds can be approximated using the MCMC samples {θ Here U denotes the total number of MCMC iterations after burn-in. The larger the BF j , the more likely gene j is an SV gene, integrating over the uncertainty in all model parameters. We suggest choosing the BF threshold based on the scale for interpretation (Kass and Raftery, 1995).

A. Simulation
We performed a series of simulation studies to evaluate the performance of BOOST-MI and compared it with that of four existing methods: SpatialDE, SPARK, BOOST-GP, and BinSpect. Because of the poor performance of Trendsceek reported in most literature (Dries and others, 2021; Sun and others, 2020), we did not include it here. In addition, we fitted the classic Ising model (with the Hamiltonian defined in Equation (1)) under the same Bayesian framework.
We generated simulated data from three artificial spatial patterns (see Figure 3(a)-(c)) and two real spatial patterns (see Figure 3(d) and (e)). The two real patterns were constructed from the mouse olfactory bulb (MOB) and human breast cancer (BC) datasets analyzed in this paper. The first two artificial patterns named spot and linear were on a 16 × 16 square lattice (n = 256 spots), while the remaining one named MOB I was on n = 260 spots. The MOB II and BC patterns were on n = 260 and 250 spots, respectively.
We set p = 100, among which 15 were SV genes. We followed the data generative schemes (Edsgärd and others, 2018;Li and others, 2021;Sun and others, 2020) to simulate the gene expression count table Y , which was substantially different from the model assumptions of BOOST-MI. For each gene j, the log relative expression level at spot i was generated via where β 0 denotes the baseline relative expression level and ij denotes the non-spatial errors following N(0, σ 2 ). We set β 0 = 2 and σ = 0.3. For a non-SV gene, the relative expression levels were from a log-normal (LN) distribution with mean and variance being 2 and 0.3 2 .
Consequently, no spatial correlation should be observed. For an SV gene with the spot pattern, the values of e i 's of the four center spots at (8, 8), (8,9), (9, 8), and (9, 9) were set to log 6, while all others were linearly decreased to zero within a radius of five spots. For an SV gene with the linear pattern, the value of e i of the most bottom-left spot at (1, 1) was set to log 6, while all others were linearly decreased to zero along the diagonal line. For an SV gene with the remaining patterns, each spot was dichotomized into low and high-expression levels with e i = 0 and log 3, respectively. To mimic the excess zeros and over-dispersion in the real ST datasets, we sampled each gene expression count y ij from a zero-inflated negative binomial (ZINB) model, y ij ∼ π i I(y ij = 0) + (1 − π i )NB(s iỹij , φ j ), where the size factor s i ∼ LN(0, 0.2 2 ) and the dispersion parameter φ j was from an exponential distribution with mean 10. For the choice of the false zero proportion π i , we randomly selected 10%, 30%, or 50% counts and forced their values to zero. Combined with the five patterns and three zero-inflation settings, there were 15 different scenarios. For each scenario, we repeated the above steps to generate 10 replicates.
We chose to normalize the raw counts Y using TSS and dichotomize the relative expression levels for each gene using GMM as the default setting. As for BOOST-MI, the prior specification are ω 0 ∼ N(1, σ 2 ω ) and θ ∼ N(0, σ 2 θ ). We set σ ω = 2.5 and σ θ = 1/3. The former indicated that the underlying proportion of the low-expression spots was expected between 1% and 99% with a probability of 95%, while the latter ensured about 99% of θ j 's value ranged from −1 to 1 a priori. A follow-up sensitivity analysis indicated that BOOST-MI was incredibly insensitive to the choice of these two hyperparameters (see Section S3 and Figure S4 in the supplement). As for the MCMC algorithm, we ran four independent MCMC chains for each gene with 10, 000 iterations, discarding the first half as burn-in. We started each chain from a model by randomly drawing all parameters from their prior distributions.
Results we report below were obtained by pooling together the MCMC outputs from the four chains. and BinSpect using the Benjamini-Hochberg method (Benjamini and Hochberg, 1995) and chose a significance level of 0.05 as the cutoff. We chose the Matthews correction coefficient (MCC) (Matthews, 1975) as the secondary performance metric, because SV genes is usually only a small subset of all genes, making other binary classification metrics not suitable.
AUC yields a value between 0 and 1, and MCC value ranges from −1 to 1. For both of them, the larger the value, the more accurate the identification.
According to Figure 4, which displays the boxplots of AUCs by different methods over ten replicated datasets under each scenario, we concluded as follows. inferences. Besides, we reported the summaries of all methods' performance in terms of MCCs in Figure S5 and we used the MOB replicate 11, which contains 16, 218 genes measured on 262 spots. We applied the quality control steps described below before applying BOOST-MI. First, we excluded spots with fewer than ten total counts across all genes. Then, those genes with more than 80% zero read counts across all spots were dropped. After these two steps, the MOB data had n = 260 spots and p = 9, 769 genes. For other methods, we applied the pre-processing procedures suggested in their papers. We used the same prior specification, algorithm setting, and significance criteria used in the simulation study. We ran four independent MCMC chains and used the potential scale reduction factor (PSRF) (Gelman and others, 1992) to diagnose MCMC convergence. PSRF is a statistic comparing the estimated between-chains and within-chain variances for a model parameter. Its value should be close to one if multiple chains have converged to the target posterior distribution. The PSRFs for all ω j0 's and θ j 's were below 1.1, clearly suggesting that the MCMC algorithms converged.
Then, for each dataset, we pooled together the outputs from the four chains and selected SV genes based on BFs. We only compared BOOST-MI with SPARK and BinSpect-rank due to the poor performance of SpatialDE and heavy computational burden of BOOST-GP.
BOOST-MI identified 734 SV genes, which was approximately the same number of SV genes detected by SPARK (772 SV genes) and around half of the number of SV genes detected by BinSpect (1, 415 SV genes). Figure 5 also included in the result of BOOST-MI or BinSpect, which indicates that SPARK is a relatively conservative method for SV gene identification. BinSpect was shown to be the most aggressive method that reported the most SV genes, more than half of which (746 out of 1, 415) were not founded by either alternative.
To further explore the spatial patterns we had found, we performed the agglomerative hierarchical clustering on the SV genes identified by each of the three methods. For better visualization, we followed the data preprocessing (Sun and others, 2020; Svensson and others, 2018) to normalize the raw read counts to relative expression levels using log-VST. Next, based on the distance matrix computed from relative expression levels of all pairs of SV genes, a hierarchical clustering dendrogram was constructed. We then determined the number of clusters by cutting the hierarchical clustering dendrogram at a height corresponding to a clear separation. Last, we summarized the expression patterns via the averaged relative expression levels within each cluster. As a result, the SV genes detected by BOOST-MI and BinSpect were clustered into five groups, while there were four groups for SPARK, shown in Figure 5(b). Consistent with results in a previous study (Sun and others, 2020), three major spatial patterns were shown in the first three columns in Figure 5 Last but not least, BOOST-MI identified 60 genes that had an attraction pattern with a positive interaction parameter in the Ising model. Table S4 in the supplement lists all of them. To analyze the potential biological functions of these SV genes, we performed functional enrichment analysis. There were 603 mouse GO terms and 43 Kyoto Encyclopedia of Genes and Genomes (KEGG) terms with at least one gene overlapping with those SV genes. We also found some statistically significant GO and KEGG terms, which had adjusted p-values less than a significance level of 0.05. For example, holo TFIIH complex (GO:0005675) and nucleotide excision repair KEGG term were significantly enriched (adjusted p-value 0.040 and 0.013, respectively). In all, these discoveries highlighted the advantage of BOOST-MI.

C. Application to human breast cancer ST dataset
The second dataset is an available ST dataset in a human breast cancer (BC) study, which is also accessible on the Spatial Research Lab (http://www.spatialresearch.org).
There are four layers available and we used the BC layer 2, which contains 14, 789 genes measured on 251 spots. We applied the same quality control steps described in Section V B. The BC data had n = 250 spots and p = 2, 280 genes. We applied the same preprocessing procedures for other methods and the same prior specification, algorithm setting, and significance criteria as mentioned. To check convergence, we run four chains and the PSRFs for all ω j0 's and θ j 's were below 1.1. We also compared our result with BinSpect and SPARK.
There were 302 SV genes identified by BOOST-MI, which was slightly larger than the number of SV genes detected by SPARK (292 SV genes) and around one-tenth of the number of SV genes detected by BinSpect (3, 278 SV genes). A Venn diagram, as shown in Figure   6(a), indicates that over 40% of the SV genes detected by BOOST-MI (124 out of 302) were in common with the result obtained by both SPARK and BinSpect. There were only 27 SV genes identified only by BOOST-MI. Once again, it demonstrated that SPARK is a conservative method and Binspect is an aggressive one, since nearly 90% of the SV genes detected by SPARK (260 out of 293) were also included in the results of BOOST-MI or BinSpect, while nearly 90% of the SV genes reported by BinSpect (2, 874 out of 3, 278) were found by neither BOOST-MI nor SPARK.
We reported more detailed analysis results, following the same procedure when studying the MOB dataset. Each row in Figure 6(b) shows the distinct expression patterns detected by BOOST-MI, SPARK, and BinSpect, respectively. There were five, four, and five groups obtained by performing the agglomerative hierarchical clustering. It is noteworthy that BOOST-MI detected all four patterns discovered by SPARK. However, BOOST-MI detected 102 SV genes with the first pattern, while SPARK reported only 52 SV genes with a similar pattern. BOOST-MI identified fewer SV genes than SPARK for the second and third patterns, while the number of SV genes was almost the same for the fourth pattern.
The last pattern of BOOST-MI was unique, and it was approximately the reversed pattern of the first one. Additionally, we repeated the same procedure on the SV genes only detected by BOOST-MI and SPARK, which is shown in Figure 6(c) and (d). In this BC study, BOOST-MI reported no SV genes with an attraction pattern.
Finally, we performed GO enrichment analysis of the 302 SV genes identified by BOOST-MI. A total of 2, 477 human GO terms had at least one gene overlap with those identified SV genes. At an FDR of 5%, 116 GO terms were found. Figure 6(f) shows the top 15 GO terms with the smallest adjusted p-values. SPARK discovered many enriched gene sets, which were related to extracellular matrix organization and immune responses (Sun and others, 2020). Although these GO terms were not shown in Figure 6(f), BOOST-MI did detect the same related terms with a significant result (e.g. the adjusted p-value for extracellular matrix organization (GO:0030198) was 1.53 × 10 −14 ). Furthermore, more virus-related GO terms were found to be significant in our analysis. There is strong evidence that many types of the virus may have a causal relationship with human breast cancers (Lawson and Heng, 2010). For example, virus life cycle GO term (GO:0019058) was significantly enriched in the reported SV genes by BOOST-MI, while not statistically significant for SV genes identified by SPARK (the adjusted p-value was only 0.829).

VI. CONCLUSION
In this paper, we develop a multi-stage method to identify SV genes from ST data.
Instead of characterizing an SV gene via a pre-specified kernel by most existing methods, we define a spatial pattern via the Hamiltonian energy in the Ising model with external fields. BOOST-MI offers the flexibility to choose different normalization and dichotomization methods, and it is considerably robust to various biological and technique noise and biases.
In the simulation study, BOOST-MI had a noticeable advantage over alternative methods, especially when there were a great number of zeros in the data. This is very encouraging since BOOST-MI does not directly model excess zeros. BOOST-MI led to more discoveries in real data analysis, such as novel spatial patterns that had never been reported and novel SV genes that kernel-based methods are unable to detect.
Several extensions of our model are worth investigating. First of all, the proposed model could be extended to model k discrete gene expression levels via a Potts model or its modified version to characterize finer spatial patterns. The number of components k can even be estimated (Green and Richardson, 2002). Second, based on the multi-stage BOOST-MI, a full hierarchical Bayesian framework could be developed to directly model the ST count data to sharpen inference. For example, with a joint Bayesian inference, BOOST-MI could further incorporate pathway information as prior knowledge to integrate the regulatory relationships between genes to perform a joint selection on SV genes. Moreover, since our approach requires the gene expression measured on a lattice grid, it is necessary to generalize our model to detect SV genes from SMP data produced by other platforms based on the single-molecule fluorescence in situ hybridization (FISH), such as sequential FISH

S1. Full details of data normalization
Normalization is critical to the analysis of sequence count data that suffer from various sequence artifacts and bias. We provide eight normalization methods in two categories.
The first type is based on size factor estimation. Let s i be the size factor of sample i, capturing all nuisance effects. Each relative gene expression level can be computed as If the main interest is in the absolute gene expression level, then all s i 's are set to the same value (e.g. s 1 = . . . = s n = 1); otherwise, we compute s i 's directly from the gene expression count data. The simplest way is to set s i ∝ Y i = p j=1 y ij , i.e., the total number of counts across all genes in each sample (known as sequencing depth or library size). Note that SPARK (Sun and others, 2020) suggests this normalization, namely total sum scaling (TSS). In practice, we could consider other estimations on s i 's, which have been introduced for mitigating the influence of extremely low and high counts when analyzing bulk RNA-seq data, such as upper-quartiles (Q75) (Bullard and others, 2010), relative log expression (RLE) (Anders and Huber, 2010), and weighted trimmed mean by M-values (TMM) (Robinson and Oshlack, 2010). Table S1 provides the definitions of the above size factor estimations. The size factor estimation is usually combined with some constraint, such as n i=1 s i = 1.
The other type of normalization method is based on variance-stabilizing transformation (VST), which aims to transform a random variable with a negative binomial distribution into one with an approximately normal distribution. There are three options: Naïve, Anscombe (Anscombe, 1948), and logarithm, namely N-VST, A-VST, and log-VST, all of which can be abstracted asỹ ij = g(y ij , φ), where g is the transformation-specific function (see Table   S2) and φ is the dispersion parameter estimated from the count data. Then, the relative gene expression levels are further adjusted for the log-scale total read counts, i.e. log Y i , via a linear regression model. Note that SpatialDE (Svensson and others, 2018) employs the log-VST normalization before fitting the geostatistical model.
We conducted a sensitivity analysis to investigate how different normalization methods affect the SV gene identification. We simulated ten replicated datasets following the data generating process described in Section 5.1 in the manuscript. We include only the scenario with the medium zero-inflation setting (30% false zeros) and the two real spatial patterns (MOB II and BC). We assessed the model performance in terms of the area under the curve (AUC). The result is summarized in Figure S1. We found that the methods based on size factor estimation significantly outperformed the VST-based methods. BOOST-MI was robust to the four size factor-based normalization methods. We conducted the analysis of variance (ANOVA) test on all pairs of the four size factor-based normalization methods. All p-values were greater than 0.05, confirming no significant difference among the four choices.

S2. Full details of data dichotomization
After correcting for sequence artifacts and bias, we denoise the relative expression levels by partitioning all spots into two groups. This step outputs the suitable data type required in the subsequent analysis and makes BOOST-MI more robust in the face of over-dispersion and zero-inflation.
There is no consensus on the dichotomization of spots based on either absolute or relative expression level. BinSpect (Dries and others, 2021) suggests allocating those spots with the top 30% relative expression levels to the high-expression group and the remaining to the low-expression group. Meanwhile, it also considers k-means (k = 2) as an alternative to avoid choosing a hard percentage rank cutoff.
We provide two choices of clustering methods. The first one is similar in spirit to BinSpect-k-means, searching the p j corresponding to the minimum within-cluster sum of squares, arg min where I(·) denotes the indicator function, m k = iỹ ij I(p ij = k)/n k and n k = i I(p ij = k) are the sample mean and size of each group. However, before applying k-means, we first exclude those spots whoseỹ ij 's are larger thanq 0.75 j + 3 q 0.75 j −q 0.25 j , whereq x j is defined as the x-th sample quantile of the relative expression levels of gene j, and then removes those spots with y ij = 0. Note that in the context of box-and-whisker plotting, a data point is defined as an extreme outlier if it stands outside this limit. We directly allocate those discarded spots in the first and second steps to the low and high-expression groups, respectively. Compared with BinSpect, this additional preprocessing leads to a more robust performance when excessive zeros and outliers are presented.
In addition to k-means that is implicitly based on pairwise distances between relative expression levels, we propose to estimate p j via fitting a two-component Gaussian mixture model (GMM) with unequal variances, subjecting to µ j0 < µ j1 . Here µ j = (µ j0 , µ j1 ) and σ 2 j = (σ 2 j0 , σ 2 j1 ) are the group means and variances that need to be estimated. To ensure the dichotomized expression levels are of the best quality to perform the subsequent modeling, we implement the above filtering steps as well.
To evaluate the performance between the two choices, we conducted a sensitivity analysis.
We simulated ten replicated datasets following the data generating process described in Section 5.1 in the manuscript. We include only the scenario with the medium zero-inflation setting (30% false zeros) and the two real spatial patterns (MOB II and BC). We assessed the model performance in terms of AUC. The result is summarized in Figure S2. BOOST-MI was robust to the two Dichotomization methods. We conducted the pairwise t-test. The resulting p-value= 0.585 and 0.850 for MOB II and BC pattern, respectively, confirming no significant difference among the two choices.

S3. Sensitivity analysis for the choices of BOOST-MI hyperparameters
We conducted a sensitivity analysis to investigate the sensitivity of BOOST-MI to the choice of σ ω and σ θ . We applied BOOST-MI to each of the ten replicated datasets under the scenario with the low zero-inflation setting (10% false zeros) and MOB II pattern. We varied values of σ θ from 1/4 to 2 and σ ω from 1 to 100. We chose five values for each hyperparameter, resulting in 25 combinations. We assessed the model performance in terms of both AUC and MCC, where the latter was based on a Bayes factor (BF) threshold of 150.
The result is summarized in Figure S4, clearly indicating that BOOST-MI was not sensitive to the choices of hyperparameters.  Note 1: q x i is defined as the x-th sample quantile of all the counts in sample i, i.e.
there are xp features in sample i whose y ij 's are less than q x i .
Note 2: The M -value M j (i, r) = log(y ij /Y i )/ log(y rj /Y r ) and A-value A j (i, r) = (log y ij /Y i + log y rj /Y r )/2 are the ratio and average of log-scaled counts between sample i and the reference sample r, respectively. G * denote a subset of genes whose