Joint identification of spatially variable genes via a network-assisted Bayesian regularization approach

Consider an undirected network $G(V,E)$ that is constructed using biological information, where $V$ is the node set consisting of $p$ genes and $E=\left\{e(j,l),j,l\in\{1,\right.$ $\left.\cdots,p\}\right\}$ is the set of edges between nodes. For genes connected within the network, it is expected that they have similar biological functionalities, leading to potentially similar spatial variability, and thus are more likely to be SV or non-SV genes simultaneously. To induce such network-assisted identification, an adjacency matrix $\boldsymbol{A}=(a_{jl})_{p\times p}$ is first constructed based on $G(V,E)$, with $a_{jl}=1$ if there is an edge $e(j,l)$ between the $j$th and $l$th genes and $a_{jl}=0$ otherwise, and ${a}_{jj}=0$, for $j,l=1,2,\ldots,p$.

As for cellular compositions, denote $\boldsymbol{w}_{i}=\left(w_{i1},\dots,w_{iK}\right)^{\mathrm{T}}$ as the vector of cellular proportions for the $i$th spot, which satisfies the constraint that $0\leq w_{ik}\leq 1,k=1,\dots,K$ and $\sum_{k=1}^{K}w_{ik}=1$. Such information is typically available as ground truth or can be obtained using deconvolution methods such as RCTD (Cable et al., 2022), Redeconve (Zhou et al., 2023), and SONAR (Liu et al., 2023). Since the distributions of cellular proportions usually display spatial relatedness, which may confound SV gene detection, we correct for this potential confounding in SV gene identification.

We introduce latent binary variables to accommodate the zero-inflation in $\mathbf{Y}$ and consider the following zero-inflated negative binomial (ZINB) model:

where $r_{ij}=1$ indicates that $Y_{ij}$ is from a Dirac probability measure $\delta_{0}(\cdot)$ with a point mass at zero, and otherwise $Y_{ij}$ is from a NB distribution $\mathrm{NB}\left(Y_{ij}\mid\mu_{ij},\phi\right)$ with expectation $\mu_{ij}$ and dispersion $1/\phi$. The NB distribution has variance $\mu_{ij}+\mu_{ij}^{2}/\phi$, thus allowing modeling extra variation.

To accommodate spatial differential expression and cell-type-specific confounding, the logarithm of $\mu_{ij}$ is modeled as:

where $\mu_{0j}$ is the baseline expression parameter, $\alpha_{k}$ is the cell-type-specific coefficient to accommodate the potential influence of cellular distributions, and $\beta^{(1)}_{j}$ and $\beta^{(2)}_{j}$ reflect the degree of spatial differential expression, with $\mathcal{K}(\cdot)$ being the pre-specified spatialization function to measure the specified trend of spatial gene expression variation.

We adopt the ZINB model due to its superiority for simultaneously accommodating the count measure, over-dispersion, and excessive zeros caused by dropouts. For SV gene identification, significantly different from the previous studies that introduce spatial differential expression via a covariance matrix, which is usually space- and time-consuming and lacks efficiency and scalability for large-scale analysis, we adopt the mean-value-based formulation. The spatial variability is intuitively reflected by $\beta^{(1)}_{j}$ and $\beta^{(2)}_{j}$, where the genes with $\beta^{(1)}_{j}\neq 0$ or $\beta^{(2)}_{j}\neq 0$ are regarded as SV. This requires less storage space and also involves a simpler estimation procedure. Moreover, we innovatively utilize a set of cell-type-specific $\alpha_{k}$’s for eliminating the confounding impact of the latent cellular composition, where spots that are spatially closer are often observed to have similar cell-type proportions (Cable et al., 2022). Different from the mean-value-based strategy adopted by Yu and Luo (2022) for marginal cell-type-specific SV gene detection, we focus on the detection of global SV genes while accommodating the cell-type proportion confounding.

The proposed priors are defined as follows:

where $\boldsymbol{\beta}=\left((\boldsymbol{\beta}^{(1)})^{\mathrm{T}},(\boldsymbol{\beta}^{(2)})^{\mathrm{T}}\right)_{(2p)}^{\mathrm{T}}=\left(\beta^{(1)}_{1},\ldots,\beta^{(1)}_{p},\beta^{(2)}_{1},\ldots,\beta^{(2)}_{p}\right)^{\mathrm{T}}$, $\boldsymbol{\gamma}=\left((\boldsymbol{\gamma}^{(1)})^{\mathrm{T}},(\boldsymbol{\gamma}^{(2)})^{\mathrm{T}}\right)_{(2p)}^{\mathrm{T}}$ represents the effect size of the genes, and $\circ$ denotes the element-wise product. $\mathbf{t}_{\lambda,\boldsymbol{\rho}}(\boldsymbol{\gamma})=\left\{\operatorname{I}\left(\left|\gamma^{(1)}_{1}\right|>\lambda\cdot\rho^{(1)}_{1}\right),\ldots,\operatorname{I}\left(\left|\gamma^{(2)}_{p}\right|>\lambda\cdot\rho^{(2)}_{p}\right)\right\}^{\mathrm{T}}$ is a vector thresholding function with $\lambda$ being a parameter controlling model sparsity and $\boldsymbol{\rho}=\left((\boldsymbol{\rho}^{(1)})^{\mathrm{T}},(\boldsymbol{\rho}^{(2)})^{\mathrm{T}}\right)_{(2p)}^{\mathrm{T}}$ being the adaptive weights. Moreover, $\mathbf{L}=\left(\operatorname{sgn}(\tilde{\boldsymbol{\beta}})\operatorname{sgn}(\tilde{\boldsymbol{\beta}})^{\mathrm{T}}\right)\,\circ\,\widetilde{\mathbf{L}}$ is a block diagonal adaptive Laplacian matrix, where $\operatorname{sgn}(\tilde{\boldsymbol{\beta}})=\left({\operatorname{sgn}}\left(\tilde{\beta}^{(1)}_{1}\right),\ldots,{\operatorname{sgn}}\left(\tilde{\beta}^{(2)}_{p}\right)\right)^{\mathrm{T}}$ with $\tilde{\boldsymbol{\beta}}$ being a rough estimate of $\boldsymbol{\beta}$ and $\widetilde{\mathbf{L}}=\mathbf{I}-\widetilde{\boldsymbol{D}}^{-1/2}\widetilde{\boldsymbol{A}}\widetilde{\boldsymbol{D}}^{-1/2}$ with $\widetilde{\boldsymbol{A}}=\left(\begin{array}[]{cc}\boldsymbol{A}&\mathbf{0}\\ \mathbf{0}&\boldsymbol{A}\end{array}\right)=(\widetilde{a}_{jl})_{2p\times 2p}$ and $\widetilde{\boldsymbol{D}}=\operatorname{diag}\left(d_{1},\cdots,d_{p},d_{1},\cdots,d_{p}\right)$ with $d_{j}=\sum_{l=1}^{p}a_{jl}$. $\varepsilon$ is a small constant to make $\left({\mathbf{L}}+\varepsilon\mathbf{I}_{2p}\right)$ strictly positive-definite, which is set as $10^{-3}$ in our numerical work.

The proposed priors have been motivated by the following considerations. The identification of SV genes is achieved using the hard-thresholded Gaussian prior, where $\beta^{(d)}_{j}$ is shrunk to zero when $\mathrm{t}_{\lambda,\rho^{(d)}_{j}}\left(\gamma^{(d)}_{j}\right)=\operatorname{I}\left(\left|\gamma^{(d)}_{j}\right|>\lambda\cdot\rho^{(d)}_{j}\right)=0$ ($d=1,2$). The thresholded Gaussian prior has been shown as a useful alternative to shrinkage priors in Bayesian sparse analysis (Wu et al., 2024), and is favored here for its simplicity and flexibility as well as its appealing interpretability as a minimal detectable signal. Here, we further introduce a series of weights $\rho_{j}^{(d)}$’s in $\mathbf{t}_{\lambda,\boldsymbol{\rho}}(\boldsymbol{\gamma})$ to adjust the shrinkage of various $\gamma^{(d)}_{j}$’s to improve selection efficiency, where the genes with strong spatial variability are potentially assigned small weights and thus are more likely to have nonzero $\beta^{(d)}_{j}$’s.

Moreover, the network dependency is introduced via the graph Laplacian matrix in the covariance matrix of the hard-thresholded Gaussian prior. The proposed network-assisted strategy is motivated by the successes of Laplacian shrinkage in high-dimensional regression analysis (Chakraborty and Lozano, 2019; Cai et al., 2020). Different from these studies, we innovatively conduct SV gene detection with the network structures among genes incorporated. In particular, the Laplacian matrix $\widetilde{\mathbf{L}}$ is further modified with a pre-defined sign matrix $\operatorname{sgn}(\tilde{\boldsymbol{\beta}})\operatorname{sgn}(\tilde{\boldsymbol{\beta}})^{\mathrm{T}}$ to accommodate the scenario where two neighborhood genes are negatively correlated and have opposite directions of spatial variability. With the proposed priors, we have $\boldsymbol{\gamma}^{\mathrm{T}}\left(\mathbf{L}+\varepsilon\mathbf{I}_{2p}\right)\boldsymbol{\gamma}=\sum_{j\sim l}\left(\frac{\text{sgn}\left(\tilde{\beta}_{j}^{(1)}\right)\gamma_{j}^{(1)}}{\sqrt{d_{j}}}-\frac{\text{sgn}\left(\tilde{\beta}_{l}^{(1)}\right)\gamma_{l}^{(1)}}{\sqrt{d_{l}}}\right)^{2}+\sum_{j\sim l}\left(\frac{\text{sgn}\left(\tilde{\beta}_{j}^{(2)}\right)\gamma_{j}^{(2)}}{\sqrt{d_{j}}}-\frac{\text{sgn}\left(\tilde{\beta}_{l}^{(2)}\right)\gamma_{l}^{(2)}}{\sqrt{d_{l}}}\right)^{2}$, where for genes $j$ and $l$ with an edge $e(j,l)$, the absolute values of $\gamma^{(d)}_{j}$ and $\gamma^{(d)}_{l}$ are promoted to be similar, further inducing simultaneous SV or non-SV.

We assign a Bernoulli prior for the latent variable $r_{ij}$ with the hyperparameter $\pi_{j}\sim\text{Beta}(a_{\pi},b_{\pi})$, where $\pi_{j}$ is the probability that $Y_{ij}$ is a dropout zero. A Gamma distribution $\text{Ga}(a_{\phi},b_{\phi})$ is assumed for the dispersion parameter $\phi$. For the variance term $\sigma_{\boldsymbol{\gamma}}^{2}$, we use the conjugate prior by assigning the Inverse-Gamma distribution $\operatorname{IG}\left(a_{\gamma},b_{\gamma}\right)$. The non-negative uniform prior $\text{Unif}\left(\lambda_{l},\lambda_{u}\right)$ is assigned for the threshold parameter $\lambda$. For $\mu_{0j}$ and $\alpha_{k}$, the normal priors with mean 0 and variance $\sigma_{0j}^{2}$ and $\sigma_{\alpha_{k}}^{2}$ are assumed, respectively. These priors have been popular in the existing Bayesian studies.

The model parameter space consists of $\left(\boldsymbol{\gamma},\boldsymbol{R},\boldsymbol{\mu}_{0},\boldsymbol{\alpha},\phi,\lambda\right)$, where $\boldsymbol{R}=(r_{ij})_{(n\times p)}$, and $\boldsymbol{\mu}_{0}$ and $\boldsymbol{\alpha}$ are the vectors consisting of $\mu_{0j}$’s and $\alpha_{k}$’s, respectively. The posterior distribution is

where $\operatorname{Be-Bern}$ denotes the Beta-Bernoulli distribution with the probability mass function $\operatorname{Be-Bern}\left(r_{ij}\mid a_{\pi},b_{\pi}\right)=\frac{\operatorname{Beta}\left(a_{\pi}+r_{ij},b_{\pi}-r_{ij}+1\right)}{\operatorname{Beta}\left(a_{\pi},b_{\pi}\right)}$.

The posterior sampling is conducted based on the MCMC algorithm. We first introduce the sampling variances $\tau_{\mu_{0}}^{2},\tau_{\alpha}^{2},\tau_{\phi}^{2},\tau_{\gamma}^{2},$ and $\tau_{\lambda}^{2}$ for $\mu_{0j}$’s, $\alpha_{k}$’s, $\phi$, $\boldsymbol{\gamma}$, and $\lambda$, respectively, and then conduct the following steps at each MCMC iteration, where the symbol “$\boldsymbol{-}$” in the condition position denotes “the rest parameters”.

1. (a)

   Sequentially update $r_{ij}$ for $\left\{\left(i,j\right):Y_{ij}=0\right\}$ with the conditional distribution of $r_{ij}$ given by $f\left(r_{ij}\mid\boldsymbol{-}\right)\propto\left(\frac{\phi}{\mu_{ij}+\phi}\right)^{\phi(1-r_{ij})}\times\operatorname{Be-Bern}\left(r_{ij}\mid a_{\pi},b_{\pi}\right)$.

2. (b)

   Sequentially sample $\mu_{0j}^{*}\sim N(\mu_{0j},\tau_{\mu_{0}}^{2})$ for $j=1,\ldots,p$, and accept $\mu_{0j}^{*}$ with probability $\min\left\{1,\frac{N\left(\mu_{0j}\mid\mu_{0j}^{*},\tau_{\mu_{0}}^{2}\right)\times N\left(\mu_{0j}^{*}\mid 0,\sigma_{0j}^{2}\right)\times\prod_{\left\{i:r_{ij}=0\right\}}\text{NB}\left(Y_{ij}\mid\mu_{0j}^{*},\boldsymbol{-}\right)}{N\left(\mu^{*}_{0j}\mid\mu_{0j},\tau_{\mu_{0}}^{2}\right)\times N\left(\mu_{0j}\mid 0,\sigma_{0j}^{2}\right)\times\prod_{\left\{i:r_{ij}=0\right\}}\text{NB}\left(Y_{ij}\mid\mu_{0j},\boldsymbol{-}\right)}\right\}.$

3. (c)

   Sequentially sample $\alpha_{k}^{*}\sim N(\alpha_{k},\tau_{\alpha}^{2})$ for $k=1,\ldots,K$, and accept $\alpha_{k}^{*}$ with probability $\min\left\{1,\frac{N\left(\alpha_{k}\mid\alpha_{k}^{*},\tau_{\alpha}^{2}\right)\times N\left(\alpha_{k}^{*}\mid 0,\sigma_{\alpha_{k}}^{2}\right)\times\prod_{\left\{(i,j):r_{ij}=0\right\}}\text{NB}\left(Y_{ij}\mid\alpha_{k}^{*},\boldsymbol{-}\right)}{N\left(\alpha^{*}_{k}\mid\alpha_{k},\tau_{\alpha}^{2}\right)\times N\left(\alpha_{k}\mid 0,\sigma_{\alpha_{k}}^{2}\right)\times\prod_{\left\{(i,j):r_{ij}=0\right\}}\text{NB}\left(Y_{ij}\mid\alpha_{k},\boldsymbol{-}\right)}\right\}.$

4. (d)

   Sample $\phi^{*}$ from $N(\phi,\tau_{\phi}^{2})$ truncated at 0, and accept $\phi^{*}$ with probability $\min\left\{1,\right.$ $\left.\frac{N_{+}\left(\phi\mid\phi^{*},0,\infty,\tau_{\phi}^{2}\right)\times\mathrm{Ga}\left(\phi^{*};a_{\phi},b_{\phi}\right)\times\prod_{\left\{\left(i,j\right):r_{ij}=0\right\}}\text{NB}\left(Y_{ij}\mid\phi^{*},\boldsymbol{-}\right)}{N_{+}\left(\phi^{*}\mid\phi,0,\infty,\tau_{\phi}^{2}\right)\times\mathrm{Ga}\left(\phi;a_{\phi},b_{\phi}\right)\times\prod_{\left\{\left(i,j\right):r_{ij}=0\right\}}\text{NB}\left(Y_{ij}\mid\phi,\boldsymbol{-}\right)}\right\}$.

5. (e)

   Sample $\boldsymbol{\gamma}^{*}$ from $\mathrm{N}\left(\mu\left(\boldsymbol{\gamma}\right),\tau_{\boldsymbol{\gamma}}^{2}\cdot\sigma_{\gamma}^{2}\left(\mathbf{L}+\varepsilon\mathbf{I}_{2p}\right)^{-1}\right)$, where

   $$\mu\left(\boldsymbol{\gamma}\right)=\sqrt{1-\tau_{\boldsymbol{\gamma}}^{2}}\boldsymbol{\gamma}+\left(1-\sqrt{1-\tau_{\boldsymbol{\gamma}}^{2}}\right)\sigma_{\gamma}^{2}\left(\mathbf{L}+\varepsilon\mathbf{I}_{2p}\right)^{-1}\nabla_{\boldsymbol{\gamma}}\log f(\mathbf{Y}\mid\boldsymbol{\gamma},\boldsymbol{-}),$$

   (5)

   with $\nabla_{\boldsymbol{\gamma}}\log f(\mathbf{Y}\mid\boldsymbol{\gamma},\boldsymbol{-})$ being the first derivative of the log-likelihood function with respect to $\boldsymbol{\gamma}$. Then, accept $\boldsymbol{\gamma}^{*}$ with probability $\min\left\{1,\frac{\prod_{\left\{(i,j):r_{ij}=0\right\}}\text{NB}\left(Y_{ij}\mid\boldsymbol{\gamma}^{*},\boldsymbol{-}\right)}{\prod_{\left\{(i,j):r_{ij}=0\right\}}\text{NB}\left(Y_{ij}\mid\boldsymbol{\gamma},\boldsymbol{-}\right)}\right\}$.

6. (f)

   Update $\rho_{j}^{(d)}=\prod_{s=1}^{S}\left(\min_{\left\{l:l\in V_{s}\right\}}\frac{1}{\left|\boldsymbol{\gamma}^{(d)}_{l}\right|^{1/2}}\right)^{\operatorname{I}(j\in V_{s})}$ for $j=1,\ldots,p$ and $d=1,2$, and update $\mathbf{t}_{\lambda,\boldsymbol{\rho}}(\boldsymbol{\gamma})$ accordingly, where $V_{1},V_{2},\ldots,V_{S}$ are the index sets of the $S$ disconnected sub-networks in $G(V,E)$.

7. (g)

   Sample $\sigma_{\boldsymbol{\gamma}}^{2}$ from $\operatorname{IG}\left(\tilde{a}_{\boldsymbol{\gamma}},\tilde{b}_{\boldsymbol{\gamma}}\right)$ with shape parameter $\tilde{a}_{\boldsymbol{\gamma}}=a_{\boldsymbol{\gamma}}+p$ and scale parameter $\tilde{b}_{\boldsymbol{\gamma}}=b_{\boldsymbol{\gamma}}+\frac{\boldsymbol{\gamma}^{\mathrm{T}}\left(\mathbf{L}+\varepsilon\mathbf{I}_{2p}\right)\boldsymbol{\gamma}}{2}$.

8. (h)

   Sample $\lambda^{*}$ from $N\left(\lambda,\tau_{\lambda}^{2}\right)$ truncated at interval $[\lambda_{l},\lambda_{u}]$, and accept $\lambda^{*}$ with probability $\min\left\{1,\frac{N_{+}\left(\lambda\mid\lambda^{*},\lambda_{l},\lambda_{u},\tau_{\lambda}^{2}\right)\prod_{\left\{\left(i,j\right):r_{ij}=0\right\}}\text{NB}\left(Y_{ij}\mid\lambda^{*},\boldsymbol{-}\right)}{N_{+}\left(\lambda^{*}\mid\lambda,\lambda_{l},\lambda_{u},\tau_{\lambda}^{2}\right)\prod_{\left\{\left(i,j\right):r_{ij}=0\right\}}\text{NB}\left(Y_{ij}\mid\lambda,\boldsymbol{-}\right)}\right\}.$

Here, updating $r_{ij}$’s and $\sigma_{\boldsymbol{\gamma}}^{2}$ is achieved through the Gibbs sampler, while the Metropolis-Hasting (MH) algorithm is adopted for sampling $\mu_{0j}$’s, $\alpha_{k}$’s, $\phi$, and $\lambda$. For sampling $\boldsymbol{\gamma}$, we resort to the preconditioned Crank-Nicolson Langevin dynamics (pCNLD), which takes advantage of the gradient information of the target distribution to speed up convergence. Furthermore, with (5), pCNLD explicitly incorporates the network dependency into the sampling process. The details of the proposed pCNLD sampling are given in Supplementary Section S1. For $\rho_{j}^{(d)}$’s, since biological networks are often composed of multiple disconnected sub-networks, we propose adopting a set of group-wise weights for the $S$ sub-networks for more effectively utilizing network information. Specifically, a series of data-driven weights inversely proportional to the absolute effect sizes of the genes involved in the specific sub-networks are introduced, potentially facilitating the identification of SV genes with weak signals, which may be involved in the sub-networks consisting of SV genes with strong signals and small thresholds. The settings for the hyperparameters and sampling variances are provided in Supplementary Section S2.

For the spatial modeling function $\mathcal{K}(\cdot)$, we follow the published studies and consider three most popular functions. Specifically, for $x_{id}~{}(d=1,2)$ which has been transformed to have mean 0 and standard deviation 1, we consider one linear function $\mathcal{K}_{1}(x_{id})=x_{id}$, two exponential functions $\mathcal{K}_{2}(x_{id})=\exp\left(-\frac{x_{id}}{2\left(l_{d}^{(1)}\right)^{2}}\right)$ and $\mathcal{K}_{3}(x_{id})=\exp\left(-\frac{x_{id}}{2\left(l_{d}^{(2)}\right)^{2}}\right)$, and two periodic functions $\mathcal{K}_{4}(x_{id})=\cos\left(\frac{2\pi x_{id}}{l_{d}^{(1)}}\right)$ and $\mathcal{K}_{5}(x_{id})=\cos\left(\frac{2\pi x_{id}}{l_{d}^{(2)}}\right)$, where $l_{d}^{(1)}$ and $l_{d}^{(2)}$ are the 40% and 60% quantiles of $|x_{1d}|,\cdots,|x_{nd}|$, respectively. To accommodate the fact that results may be sensitive to the choice of scale parameter, two values are considered for the exponential and periodic patterns. Then, we conduct the MCMC sampling for each $\mathcal{K}_{s}(\cdot)(s=1,\cdots,5)$ and obtain the estimated posterior expectation $\hat{\boldsymbol{\beta}}^{{\mathcal{K}_{s}}}=\frac{\sum_{m=1}^{M}\boldsymbol{\gamma}^{(m)}\boldsymbol{\mathbf{t}}_{\lambda^{(m)},\boldsymbol{\rho}^{(m)}}\left(\boldsymbol{\gamma}^{(m)}\right)}{\sum_{m=1}^{M}\boldsymbol{\mathbf{t}}_{\lambda^{(m)},\boldsymbol{\rho}^{(m)}}\left(\boldsymbol{\gamma}^{(m)}\right)},\hat{\boldsymbol{\mu}}_{0}^{{\mathcal{K}_{s}}}=\frac{\sum_{m=1}^{M}\boldsymbol{\mu}_{0}^{(m)}}{M},\hat{\boldsymbol{\alpha}}^{{\mathcal{K}_{s}}}=\frac{\sum_{m=1}^{M}\boldsymbol{\alpha}^{(m)}}{M}$, $\hat{\phi}^{{\mathcal{K}_{s}}}=\frac{\sum_{m=1}^{M}\phi^{(m)}}{M}$, and $\hat{\boldsymbol{R}}^{{\mathcal{K}_{s}}}=\frac{\sum_{m=1}^{M}\boldsymbol{R}^{(m)}}{M}$, where $\left\{\boldsymbol{\gamma}^{(m)},\lambda^{(m)},\boldsymbol{\rho}^{(m)},\boldsymbol{\mu}_{0}^{(m)},\boldsymbol{\alpha}^{(m)},\phi^{(m)},\boldsymbol{R}^{(m)}\right\}_{m=1}^{M}$ denotes the $M$ samples obtained after burn-in and thinning (we omit the dependence on $\mathcal{K}_{s}$ to simplify notation).

To facilitate the combination of five models, instead of directly considering the values of $\hat{\boldsymbol{\beta}}^{{\mathcal{K}_{s}}}$ for SV gene identification, we further introduce a posterior inclusion probability vector estimated as ${\mathrm{PIP}}^{{\mathcal{K}_{s}}}=\frac{1}{M}\sum_{m=1}^{M}\boldsymbol{\mathbf{t}}_{\lambda^{(m)},\boldsymbol{\rho}^{(m)}}\left(\boldsymbol{\gamma}^{(m)}\right)$ and consider the posterior model probability $f\left(\mathcal{M}_{\mathcal{K}_{s}}\mid\mathbf{Y}\right)$ introduced in Quintana et al. (2011), where $\mathcal{M}_{\mathcal{K}_{s}}$ denotes the model with $\mathcal{K}_{s}(\cdot)$. Specifically, for each $\mathcal{M}_{\mathcal{K}_{s}}$, we calculate $f\left(\mathcal{M}_{\mathcal{K}_{s}}\mid\mathbf{Y}\right)=\frac{f\left(\mathbf{Y}\mid\mathcal{M}_{\mathcal{K}_{s}}\right)f\left(\mathcal{M}_{\mathcal{K}_{s}}\right)}{\sum_{s^{\prime}=1}^{5}f\left(\mathbf{Y}\mid\mathcal{M}_{\mathcal{K}_{s^{\prime}}}\right)f\left(\mathcal{M}_{\mathcal{K}_{s^{\prime}}}\right)}$. Here, $f\left(\mathcal{M}_{\mathcal{K}_{s}}\right)$ is the prior probability of model $\mathcal{M}_{\mathcal{K}_{s}}$, and we set a non-informative prior with $f\left(\mathcal{M}_{\mathcal{K}_{s}}\right)=1/5$. $f\left(\mathbf{Y}\mid\mathcal{M}_{\mathcal{K}_{s}}\right)$ is the likelihood function with the estimated parameters $\hat{\boldsymbol{\beta}}^{{\mathcal{K}_{s}}},\hat{\boldsymbol{\mu}}_{0}^{{\mathcal{K}_{s}}},\hat{\boldsymbol{\alpha}}^{{\mathcal{K}_{s}}},\hat{\phi}^{\mathcal{K}_{s}}$, and $\hat{\boldsymbol{R}}^{{\mathcal{K}_{s}}}$ under model $\mathcal{M}_{\mathcal{K}_{s}}$. We then utilize $f\left(\mathcal{M}_{\mathcal{K}_{s}}\mid\mathbf{Y}\right)$’s as the model-specific weights to obtain the combined $\operatorname{PIP}=\sum\limits_{s=1}^{5}f\left(\mathcal{M}_{\mathcal{K}_{s}}\mid\mathbf{Y}\right)\mathrm{PIP}^{{\mathcal{K}_{s}}}$. Moreover, as gene $j$ with at least one non-zero $\beta^{(d)}_{j}\left(d=1,2\right)$ is identified as SV, we introduce $\widetilde{\operatorname{PIP}}_{j}=\max\left(\operatorname{PIP}_{j},\mathrm{PIP}_{j+p}\right)$ for $j=1,\ldots,p$. Based on $\widetilde{\operatorname{PIP}}_{j}$’s, the Bayesian false discovery rate (BFDR) control strategy is adopted for controlling multiplicity. Specifically, $\operatorname{BFDR}(c)=\frac{\sum_{j=1}^{p}\left(1-\widetilde{\operatorname{PIP}}_{j}\right)\mathrm{I}\left(1-\widetilde{\operatorname{PIP}}_{j}<c\right)}{\sum_{j=1}^{p}\mathrm{I}\left(1-\widetilde{\operatorname{PIP}}_{j}<c\right)}$ with $\operatorname{BFDR}(c)$ being the desired significance level. We set $\operatorname{BFDR}(c)$ as 0.05 in our numerical studies. Then, the SV gene set is defined as $\left\{j:\widetilde{\operatorname{PIP}}_{j}\geq c\right\}$. To achieve improved stability and better false negative control, in numerical studies, we run five MCMC chains independently, and the genes identified in more than 80% of the chains are finally identified as SV.

To improve computational efficiency, parallelization is implemented with the R package RcppParallel. Specifically, the marginal sequential sampling is divided into parallel programming to reduce computer time. For $\boldsymbol{\gamma}$ with the dependency Laplacian matrix incorporated, benefiting from the sparse and block-wise nature of biological networks, the parallelism block-wise sampler is adopted to avoid sampling from a high-dimensional multivariate normal distribution, which further accelerates computation. Take a simulated dataset with $p=5,000$ and $n=1,024$ as an example. The average computer time is about 1.01 minutes for 100 iterations on a MacBook Pro with an Intel Core i5 CPU and 16GB RAM.

Simulation studies are conducted under the following settings. (a) $n=1,024$ spots located on a 32 by 32 square lattice, $p=5,000$ genes, and $K=6$ cell types. (b) The square lattice is partitioned into three regions as displayed in Figure S1, where the cellular compositions $\boldsymbol{w}_{i}$’s are independently sampled from Dirichlet distributions $\operatorname{Dirc}(1,1,1,1,1,1)$ (Region 1), $\operatorname{Dirc}(3,5,7,9,11,13)$ (Region 2), and $\operatorname{Dirc}(18,16,14,12,10,8)$ (Region 3). (c) Consider three types of spatial pattern $\mathcal{K}\left(x_{id}\right)$, including $x_{id}$ (Linear), $\exp\left(-\frac{x_{id}}{2}\right)$ (Exponential), and $\cos\left(2\pi x_{id}\right)$ (Periodic). (d) The networks are block-wise and composed of 100 disconnected sub-networks with 50 nodes each. For each sub-network, two types of network structure, namely Star and Scale-free, are considered to mimic the real-world transcription factor regulatory network and interaction network with scale-free properties, respectively. Illustrative examples of the sub-networks are presented in Figure S2(A) and S2(B). (e) All genes in the first ten sub-networks are SV, leading to a total of 500 SV genes. Both positive and negative signals are considered with various levels of magnitude. More detailed settings are provided in Supplementary Section S3.1. (f) The spatial transcriptomics count data is generated from model (1) with the dispersion parameter $\phi$ being 10. Two settings of dropout rates, 0.1 and 0.5, are considered, representing low and high sparsity. The baseline parameter ${\mu}_{0j}$ and cell-type-specific effect ${\alpha}_{k}$ are independently generated from $N(2,0.5^{2})$ and $N(0,3.5^{2})$, respectively. There are 12 scenarios, comprehensively covering a wide spectrum with different patterns of spatial expressions, different structures of networks and the corresponding spatially variable signals, and different degrees of sparsity.

In addition to the proposed approach, six alternatives are considered: (a) SpatialDE (Svensson et al., 2018), a likelihood ratio test method based on Gaussian process regression, (b) SPARK (Sun et al., 2020), a method based on a Poisson log linear model with a Gaussian process incorporated, (c) SPARK-X (Zhu et al., 2021), a scalable non-parametric test built on a robust covariance test framework, (d) MERINGUE (Miller et al., 2021), a test based on spatial cross-correlation analysis, (e) CTSV and CTSV-g (Yu and Luo, 2022), a test based on the ZINB model for cell-type-specific SV gene identification and global SV gene identification, respectively, and (f) HEARTSVG (Yuan et al., 2023), a distribution-free method that employs the exclusion of non SV genes to infer the presence of SV genes. Among them, SpatialDE, SPARK, CTSV, and CTSV-g are parametric, while SPARK-X, MERINGUE, and HEARTSVG are non-parametric. All the alternatives conduct tests marginally and implement multiple testing control. The implementation details for the competing methods are given in Supplementary Section S3.2.

For evaluating identification performance, we adopt Recall = $\frac{\text{TP+FN}}{\operatorname{TP}}$ , Precision =$\frac{\operatorname{TP+FP}}{\operatorname{TP}}$, and F1 score = $\frac{2\cdot\text{Precision}\cdot\operatorname{Recall}}{\text{Precision+Recall}},$ with TP, FP, and FN being the numbers of true positives, false positives, and false negatives, respectively. Under each scenario, we simulate 50 replicates. Summary results under the scenarios with a low dropout rate are given in Table 1. The rest of the results with a high dropout rate are provided in Supplementary Table S1. It is observed that the proposed approach achieves superior accuracy in identifying SV genes with higher F1 scores across all scenarios. Overall, CTSV-g achieves the second-best identification performance since it also accommodates dropouts and over-dispersion through the adoption of ZINB distribution. Improved performance over CTSV-g supports the validity of incorporating network information and accommodating cellular diversity. Furthermore, the majority of the alternatives exhibit good performance with the simple linear patterns while the proposed approach performs stably across different spatial expression patterns. Moreover, similar to those in some published studies, SPARK-X exhibits many false negatives, probably due to its covariance test framework and neglection of cellular composition confounding, and SPARK achieves excellent FP control across all scenarios, mainly attributable to the robust Cauchy combination rule. CTSV identifies the fewest global SV genes across all scenarios due to its concern with cell-type-specific SV genes and inability in accurately identifying genes with global spatial expression variations. Noticeably, under the scenarios with a scale-free network where the signals of the SV genes are larger, all approaches exhibit improved performance with higher F1 score values. However, the superiority of the proposed approach is again evidently observed. With a higher dropout rate, which is more common with practical ST data, all approaches have decayed performance, especially for SpatialDE which is not well-suited for sparse expression distribution. However, the superiority of the proposed approach becomes more prominent.

To better appreciate the operating characteristics of the proposed network-assisted strategy, in Figure S3, we take the scale-free network as an example and further examine the SV genes with large and small signals separately with two Recall indices (L:Recall and S:Recall). Here, besides the proposed approach, we also consider the corresponding approach without the assistance of the network (where $\mathbf{L}$ is a simple identity matrix). As shown in Figure S3, benefiting from network smoothness, the proposed approach has remarkably improved identification performance, especially for the SV genes with small signals. Such an improvement becomes more obvious for sparser data with a higher dropout rate. Moreover, selection stability is also improved.

We further evaluate performance of the proposed approach under scenarios where the data generation model is misspecified. In particular, we consider two types of model misspecification: (M1) the ZINB model considered in Yu and Luo (2022) where the logarithmic mean value is set as a mix of cell-type-specific expression levels; and (M2) the Poisson model considered in Sun et al. (2020) with the parameters inferred from a real mouse olfactory bulb study. More detailed settings are provided in Supplementary Section S3.4, and summary results are reported in Figure S6.

It is observed that the proposed approach can still maintain its superiority or at least perform comparably to the method that matches the data generation model. Specifically, under model (M1), which favors the CTSV method, the proposed approach achieves more accurate identification with a median F1 equal to 1, while CTSV exhibits the second-best performance. Moreover, even under model (M2) without cellular composition confounding, the proposed approach can still correctly identify SV genes with a median F1 score larger than 0.92. Here, spatialDE and SPARK also behave well since the expression pattern is inferred based on the analysis results of the two approaches.

We continue to examine scenarios with noisy network and cellular composition information. In particular, the noisy network is likely to involve edges that connect both SV and non-SV genes. Two specific generation models are considered, including the ZINB model considered in Section 3.1 with a star network, a linear spatial pattern, and two degrees of sparsity, and the misspecified model (M1) considered in Section 3.2, where in each informative sub-network, about 20% of the genes connected to the TF genes are uninformative. The summary results are reported in Figure S7. It can be seen that the proposed approach can still achieve an effective balance between the prior network information and observation likelihood, leading to certain robustness.

Moreover, to explore the impact of noisy cellular composition on identification performance, we consider the scale-free network with the linear pattern and high sparsity (dropout rate = 0.5), where the cell-type deconvolution estimates $\hat{\boldsymbol{w}}_{i}$ are sampled from $\operatorname{Dirc}(c\cdot\boldsymbol{w}_{i})$ with $\boldsymbol{w}_{i}$ being the true value and concentration parameter $c=100,80,$ and 50. The lower value of $c$ indicates less accuracy of the cellular composition estimates. The summary results are provided in Figure S8. As expected, the noisy cellular compositions tend to result in decreased identification performance with increased variance. However, even with $c=50$, the identification performance is still satisfactory, and the superiority over the alternatives maintains.

1. (a)

   SpatialDE: The source code for SpatialDE is publicly available on https://github.com/Teichlab/SpatialDE. Following the instructions, we filter the practically unobserved genes with total observations less than 3. The linear regression is conducted on the raw count matrix to account for the potential bias caused by library size or sequencing depth. The main function for SV gene detection is SpatialDE.run() with the arguments being the coordinate information $\left(x_{i1},x_{i2}\right)$ and the sample residual expressions corrected by the linear regression. To account for multiple testing, the genes with qval less than 0.05 are identified as SV.

2. (b)

   SPARK: The source code for SPARK is publicly available on https://github.com/xzhoulab/SPARK. Following the instructions, we filter genes that are expressed in less than 10% spots and spots with the total observations less than 10. The statistical model under the null hypothesis is first fit by employing the function spark.vc(), and the function spark.test() is subsequently employed for SV gene detection. Then, the genes with adjusted_pvalue less than 0.05 are identified as SV.

3. (c)

   SPARK-X: The source code for SPARK-X is publicly available on https://github.com/xzhoulab/SPARK. The SV gene detection is implemented with the main function sparkx() where the argument option is set as “mixture” for multiple kernels testing. To account for the impact of latent cellular compositions, as conducted in the original paper, each spot is assigned its major cell type, with the argument X_in being the matrix of binary indicators. Then, the genes with adjustedPval less than 0.05 are identified as SV.

4. (d)

   MERINGUE: The source code for MERINGUE is publicly available on https://github.com/JEFworks-Lab/MERINGUE. The non-expressed genes and spots are filtered following the tutorial. Normalization is conducted for the raw count matrix by employing the function normalizeCounts(). For SV gene detection, the neighborhood relationships are first constructed using the function getSpatialNeighbors() with the argument filterDist set as the default value 2.5. Then, the SV genes are identified with the function filterSpatialPatterns() with the arguments minPercentCells set as 0.05 to restrict that the SV genes are driven by more than 5% of spots. The adjusted significance threshold is set as 0.05 through setting the arguments adjustPv as TRUE and alpha as 0.05.

5. (e)

   CTSV and CTSV-g: The source code for CTSV is publicly available on https://github.com/jingeyu/CTSV. The SV gene detection is conducted through the main function ctsv(). Specifically, the cell-type-specific SV genes are identified with the argument $W$ being the $n\times K$ matrix composed of $w_{i}$’s ($i=1,\ldots,n$). The final SV gene set is the union of all cell-type-specific SV genes. We also conduct the global SV gene detection without accommodation for the cellular composition in the simulation studies (CTSV-g). Specifically, the argument $W$ is set as $(1,\ldots,1)^{\mathrm{T}}_{(n)}$. For both CTSV and CTSV-g, the SV genes are identified through function svGene() with the significance threshold thre.alpha set as 0.05.

6. (f)

   HEARTSVG: The source code for HEARTSVG is publicly available on https://github.com/cz0316/HEARTSVG. The detection of SV genes is conducted through the main function heartsvg() with the raw count matrix first scaled as recommended. The Holm method is further conducted for multiple testing control. The genes with adjusted P values p_adj less than 0.05 are identified as SV.