A Joint Spatiotemporal Differential Expression Modeling of Exposure Effects in Spatial Transcriptomics Data

A Joint Spatiotemporal Differential Expression Modeling of Exposure Effects in Spatial Transcriptomics Data Osafu Augustine Egbon1,2 and Benedict Anchang1,2 1National Institute of Environmental Health Sciences (NIEHS/NIH), Durham, NC 27709, United States 2National Cancer Institute, Bethesda, MD 20892, United States Abstract Differential expression analysis in spatial transcriptomics data is important for identifying localized biological processes and spatially structured gene regulation. While numerous methods have been developed for identifying differentially expressed genes (DEGs) in conventional single-cell RNA-seq data, approaches that explicitly account for spatial and temporal context remain scarce. A major challenge in spatial DEG analysis is the spatial misalignment between tissue sections under control and perturbed conditions, which renders coordinate-based comparisons unreliable. To address this gap, we developed SpatialDENet, a framework that identifies DEGs in spatially structured single-cell data by incorporating spatial information through a hierarchical network-based model. This approach accounts for spatial confounding and leverages tissue architecture to align biological signals across conditions. We validated SpatialDENet through simulation studies and comparative benchmarks with existing popular methods. Our findings demonstrate that SpatialDENet provides a powerful and interpretable alternative for spatial transcriptomics analysis and may inform the development of targeted therapies. Key Words: Bayesian inference, Network model, Differential expression analysis, Gaussian Markov Random Field, Zero-inflation. 1 Introduction Spatial transcriptomics technologies are increasingly popular for investigating the spatial architecture of tissues, uncovering cell-cell interactions, and analyzing gene expression patterns in situ with high resolution1,2. However, conducting Differential Expression (DE) analysis to identify differentially expressed genes (DEG) between two biological conditions on spatial transcriptomics data presents significant challenges due to its high dimensionality, noise, sparsity, and spatial structure. Methods for DE analysis in studying biological variations at the single-cell level are emerging quickly. For instance, the popularly used methods include DESeq23, EdgeR4, Lima5, MAST6, Zinbwave7, among others. However, most of these methods do not account for spatial variation, an important source of biological and technical variability that can confound differential expression analysis results8. Existing statistical and computational methods8–13 for analyzing spatially resolved 1 single-cell data are particularly designed to identify spatially variable genes (SVGs), and no adequate attention has been given to identifying DEGs between control and treatment scenarios in spatially resolved single-cell data. To address this concern, we proposed a statistical framework, SpatialDENet (Sptial Differential Expression Analysis through Network), to identify differential expression genes between control and perturbed conditions while correcting for spatial confounding (SC). SpatialDENet addresses three main analytical concerns in spatially resolved single-cell data. Firstly, SpatialDENet adjusts for the SC issue through spatial hierarchical modeling. This method allows the accurate estimation of the impact of a specific treatment on cells after controlling for the spatial effects. Secondly, spatially resolved single-cell data is associated with spatial misalignment due to the non-longitudinal and nonhomologous nature of the data across samples. That is, the spatial coordinates change from tissue samples to samples. To address this concern, SpatialDENet leverages network models to map cellular spatial neighborhoods between multiple tissues unto a latent space, creating a unified spatial field for adjusting for SC and identifying DEGs. Lastly, SpatialDENet accounts for the sparsity in the expression data. The zeros due to sparsity can arise from technical limitations in the sequencing process or reflect true biological absences of gene expression for specific cells. Distinguishing between these two sources of zeros is critical for accurate interpretation. SpatialDENet utilizes a two-part (count and activation) zero-inflated model to account for excess zeros. The remainder of this paper is organized as follows. We first briefly describe the statistical framework of SpatialDENet and its extension to accommodate zero-inflated data. We then benchmark the framework with popularly used differential expression analysis methods through a simulation study. Finally, we conclude with a summary of our contributions and potential directions for future research. 2 Statistical Framework 2.1 Background of differential analysis framework In single-cell DE analysis, the general framework adopted by existing methods can be described as follows. Let 𝑌𝑔𝑠 be the observed expression level (e.g., read count) for gene 𝑔in cell 𝑠, where 𝑔=1, . . . , 𝐺 and 𝑠=1, . . . , 𝑆. A generalized linear model (GLM) is typically used to model gene expression as: 𝑌𝑔𝑠 ∼NB(𝜇𝑔𝑠, 𝜙𝑔),(1) where NB denotes a Negative Binomial distribution, 𝜇𝑔𝑠 is the expected expression level of gene 𝑔 in cell 𝑠,𝜙𝑔is a gene-specific dispersion parameter. The parameter 𝜇𝑔𝑠 is then linked to covariates 2 through a log-linear model given as log(𝜇𝑔𝑠)=𝛽𝑔0+log(𝑙𝑠) +𝑐𝑠𝛽𝑔1+x⊤ 𝑠𝜷𝑔,(2) where 𝑙𝑠is a size factor or library size normalization for cell 𝑠,𝛽0is the model intercept, 𝑐𝑠∈ {0,1} is the treatment indicator, with effect size 𝛽1,x𝑠is a vector of covariates to adjust for observed confounding variables, and 𝜷𝑔is a corresponding vector of regression coefficients for gene 𝑔. To test for differentially expressed genes between control and treatment conditions, we define the hypothesis (Hypothesis 1) for a specific gene 𝑔as 𝐻0:𝛽𝑔1=0vs. 𝐻1:𝛽𝑔1≠0,(3) where 𝛽𝑔1is the coefficient corresponding to the condition of interest. In single-cell spatially resolved data with genes exhibiting spatial variation and zero inflation, Hypothesis 1 becomes problematic and inappropriate for finding differential expression genes in single-cell analysis. 2.2 SpatialDENet: The proposed framework To address the above concern, we proposed a unified framework that addresses spatial confounding, spatial misalignment, the temporal evolution of gene expression, and zero-inflation problems in single-cell spatially resolved data. 2.2.1 Response model Consider a scenario where we have pre-processed single-cell expression data collected over time. Let 𝑌𝑔𝑠𝑡, 𝑠 =1,2, ..., 𝑆;𝑔=1,2, ..., 𝐺,𝑡=1,2, ..., 𝑇 be a gene expression count at spatial location denoted by 𝑠at time 𝑡for gene 𝑔.𝑆is the total number of spots/cells across all tissues. Note here that the cell index 𝑠previously defined as a cell is now redefined as a spot index 𝑠in spatial transcriptomics data. We assume that 𝑌𝑔𝑠𝑡 |𝛽𝑔0, 𝛽𝑔1,𝜷𝑔∼NB(𝜇𝑔𝑠𝑡 , 𝜙𝑔), log(𝜇𝑔𝑠𝑡)=𝛽𝑔0+log(𝑙𝑠𝑡 ) +𝑐𝑠𝑡 𝛽𝑔1+x⊤ 𝑠𝑡 𝜷𝑔+𝜓𝑔𝑠𝑡, (4) where 𝜓𝑔𝑠𝑡 is a spatiotemporal latent variation used to adjust for spatial confounding in the gene expression data, 𝑙𝑠𝑡 is the size factor, 𝑐𝑠𝑡 ∈ {0,1}is the condition indicator of spatial point 𝑠at time 𝑡. The differential genes are then identified using Hypothesis 1. In standard spatial statistical modeling, where the spatial location domain remains fixed, 𝜓𝑔𝑠𝑡 in Equation (4) can be easily estimated by assuming a Gaussian Markov random field (GMRF) prior distribution. However, in single-cell analysis, the spatial domain (tissue sample) dissociates during data collection and thus, different tissues are used in the treatment and control conditions, 3 causing spatial misalignment due to non-homologous tissue samples. To address this challenge, we model 𝜓𝑔𝑠𝑡 using a Network-Based-GMRF prior distribution, which is discussed in the subsequent subsection. 2.2.2 Deriving network-based-GMRF To define a network-based Gaussian Markov Random Field (GMRF) prior for modeling spatiotemporal latent effects, we introduce a framework that constructs a dynamic cell-neighborhood network over a spatial domain, informed by gene expression summaries14. Let Ddenote a regular spatial domain for a specific tissue sample, such that 𝑠∈ D. For each time 𝑡or a specific tissue sample, the domain D𝑡⊂ D is non-linearly partitioned into non-overlapping spatial regions 𝐷𝑖𝑡 ,𝑖=1,2, ..., 𝐼, called a node, such that D𝑡= 𝐼𝑡 Ø 𝑖=1 𝐷𝑖𝑡 , 𝐷𝑖𝑡 ∩𝐷𝑖′𝑡=∅for 𝑖≠𝑖′.(5) For each partition 𝐷𝑖𝑡 , we calculate a representative centroid using a summary statistic, such as the mean expression of the principal components (PCs) derived from the gene expression data. Suppose the reduced gene expression space consists of 10 PCs; for each partition 𝐷𝑖𝑡 , we compute the mean value of each PC across all cells within that partition. Consequently, each 𝐷𝑖𝑡 is represented by a vector of 10 mean values, corresponding to the 10 PCs. The number of nodes is chosen to result in optimal precision. We then define a cell-neighborhood network Tby connecting each node 𝐷𝑖𝑡 to its 𝑘-nearest neighbors based on centroid proximity. Formally, define an undirected graph T=(G,E), where Gis the set of nodes (partitions across all 𝑡) and Eis the set of edges, such that (𝑖, 𝑡)∼(𝑗, 𝑡′) if and only if 𝐷𝑗𝑡′is among the 𝑘nearest neighbors of 𝐷𝑖𝑡 based on centroid distance. For example, if there exists a total of 1000 non-overlapping partitions of D, the network Twill have 1000 nodes (#G=1000), indicating that each cell/spot belongs to only one of the 1000 nodes. Given the nodes in G, the Network Tis used to incorporate the spot/cell interaction and dependencies informed by the relationship between the cell’s expression status across all the genes. To obtain the network-based GMRF, suppose network Thas 𝑛nodes, with 𝜓𝑔𝑖𝑡 ∈Rrepresenting the node latent effect for node 𝑖at time 𝑡. Let 𝝍𝑔𝑡 =(𝜓𝑔1𝑡, 𝜓𝑔2𝑡, ..., 𝜓𝑔𝑛𝑡 )𝑇and 𝝍𝑔= (𝝍𝑔1,𝝍𝑔2, ..., 𝝍𝑔𝑇 )𝑇. We assumed a Gaussian Markov Random Field (GMRF) on the latent variables 𝜓𝑔𝑖𝑡 , which captures the network structure. Let N𝑖be the collection of all the index sets of nodes having an edge with node 𝑖on network T, then the network-based GMRF is defined as 𝜓𝑔𝑖𝑡 |𝝍−𝑖 𝑔,N𝑖∼𝑁Í𝑗∈N𝑖𝜓𝑔 𝑗𝑡 #N𝑖 ,1 #N𝑖𝜏𝑔,(6) where #N𝑖is the total number of nodes connected to node 𝑖.𝝍−𝑖 𝑔is 𝝍𝑔but excluding component 𝜓𝑔𝑖𝑡 from the vector. Equation (6) shows that the expected value of the latent variable in node 𝑖is 4 the average of the values of the latent variables representing the nodes connected to node 𝑖, and the variance is scaled by the number of edges of node 𝑖. The equation captures the dependency between the nodes across the network. For example consider a network T=(G,E) such that the node set is G={1,2,3,4,5}and edge set E={1∼3,1∼2,2∼3,2∼4,3∼4,3∼5}, where 𝑖∼𝑗indicates that node 𝑖is connected to node 𝑗on the network. Therefore, N1={2,3},N2={1,3,4},N3= {1,2,4,5},N4={2,3},N5={3}. The joint distribution of the network-based-GMRF for the latent variable 𝝍𝑔|𝜏𝑔∼𝑁(0,𝚺(𝜏𝑔)), where 𝚺(𝝉𝑔)=(𝜏𝑔Q)−1and Qis a structured matrix that captures the spatial structure given in Equation (6), which is defined as 𝑄𝑖𝑖′=        |N𝑖|,if 𝑖=𝑖′, −1,if 𝑗∈ N𝑖, 0,otherwise. (7) 2.2.3 Zero-Inflation joint model A serious concern with single-cell data is its sparsity. The expression contains excessive zeros that can bias most modeling approaches, especially interpretable models. As a result, it is important to account for these zeroes in the modeling framework to differentiate between true zeroes and biological zeros. To address this concern, we adopted a “two-part" modeling technique. In the first part, we model the activation process of gene 𝑔through a logistic model. We model the probability that a specific gene is activated (expressed) or not. In the second step, the extent of activation is modeled through a different standard probability model. Mathematically, define a new variable 𝑍𝑔𝑠𝑡 such that 𝑍𝑔𝑠𝑡 =(1if 𝑌𝑔𝑠𝑡 >0, 0otherwise. (8) Thus, for a given gene 𝑔, we assume that 𝑍𝑔𝑠𝑡 ∼Bernoulli(𝜋𝑔𝑠𝑡), where 𝜋𝑔𝑠𝑡 denotes the probability of observing non-zero expression. Moreover, we define another variable 𝑉𝑔𝑠𝑡 =𝑌𝑔𝑠𝑡 only if 𝑍𝑔𝑠𝑡 =1, and 𝑉𝑔𝑠𝑡 is modeled using a zero-truncated negative binomial distribution. We define 𝜓𝑔𝑖𝑡 as the latent effect associated with node 𝑖at time 𝑡. Each spot 𝑠is assigned to a node 𝑖(𝑠), and thus the spot-level effect is given by 𝜓𝑔𝑠𝑡 =𝜓𝑔,𝑖(𝑠),𝑡 . The full model for the SpatialDENet is given as 5 𝑍𝑔𝑠𝑡 |.∼Bernoulli(𝜋𝑔𝑠𝑡), log 𝜋𝑔𝑠𝑡 1−𝜋𝑔𝑠𝑡 =𝛼𝑔0+𝑐𝑠𝑡 𝛼𝑔1+x𝑇 𝑠𝑡 𝜶𝑔+𝜓𝑔𝑖𝑡 , 𝑠 ∈node 𝑖 𝑉𝑔𝑠𝑡 |.∼NegativeBinomial(𝜇𝑔𝑠𝑡, 𝜙𝑔), log(𝜇𝑔𝑠𝑡 )=𝛽𝑔0+log(𝑙𝑠𝑡 ) +𝑐𝑠𝑡 𝛽𝑔1+x𝑇 𝑠𝑡 𝜷𝑔+𝜓𝑔𝑖𝑡 , 𝝍𝑔∼𝑁0,𝚺(𝜏𝑔), (𝛼𝑔0, 𝛼𝑔1,𝜶𝑔, 𝛽𝑔0, 𝛽𝑔1,𝜷𝑔)𝑇∼𝑀𝑉𝑁(0,104I), 𝚲∼𝑃(𝚲), (9) where 𝜋𝑔𝑠𝑡 is the probability that a specific gene 𝑔is biologically activated. Each spot 𝑠belongs to a node 𝑖(𝑠). Thus, in the linear predictor term, 𝜓𝑔𝑖𝑡 is included through 𝜓𝑔,𝑖(𝑠)𝑡, which is shared between the binary and the count components. As previously defined, 𝜓𝑔𝑖𝑡 adjusts for the spatial confounding in the activation component of the model using the network-based-GMRF. 𝚲is a vector of hyperparameters, 𝚲=(𝜏𝑔, 𝜙𝑔)𝑇and 𝑃(𝚲)is the prior distribution for 𝚲. We assigned loggamma prior for 𝜙𝑔,log(𝜙𝑔) ∼ 𝑙𝑜𝑔 −𝐺𝑎𝑚𝑚𝑎(1,0.00005). We adopted the PC prior distribution for 𝜏𝑔, which is among the most widely used prior distributions due to its robustness. The PC prior distribution for 𝜏=𝜏𝑔is given as 𝜋(𝜏)=𝜆 2𝜏−3/2exp −𝜆 𝜏−1/2,for 𝜏 > 0,(10) where 𝜆=−log 𝛼 𝑢, such that 𝑃(1 √𝜏> 𝑢)=𝛼. The prior is fully specified if 𝑢and 𝛼are fixed or known, which can be informatively derived from the empirical data. 3 Simulation study We conducted a simulation study to benchmark SpatialDENet with commonly used methods in the literature. To simulate the gene expression data, we adopted a zero inflated negative binomial distribution. To simulate the spatially random error, we adopted a Gaussian bump model. The data generation model is given as xcoord ∼uniform(0,1) ycoord ∼uniform(0,1) 𝑓𝑔(𝑠, 𝑡)=exp −2𝑡((xcoord −0.5)2+ (ycoord −0.5)2) log(𝜇𝑔𝑠𝑡)=𝛽𝑔0+𝑐𝑠𝑡 𝛽𝑔1+𝑓𝑔(𝑠, 𝑡) 𝑍𝑔𝑠𝑡 ∼Bernoulli(𝜋𝑔𝑠𝑡) 𝑌𝑔𝑠𝑡 =(0if 𝑍𝑔𝑠𝑡 =0, NB(𝜇𝑔𝑠𝑡 , 𝜙𝑔)if 𝑍𝑔𝑠𝑡 =1. (11) 6 where xcoord and ycoord represent the spatial coordinates of each location, and 𝑓𝑔(𝑠, 𝑡)denotes the spatially varying error term, which can potentially confound the estimation of 𝛽𝑔1. The error is higher towards the center of the domain and fades away towards the edges. We assumed an effect size of 𝛽𝑔1=[2,4], 𝑔 =1,2, ..., 𝐺 between the treatment and control conditions. We deliberately introduced this spatial error to evaluate how effectively the network-based model captures spatial structure, without directly simulating from the network-based GMRF. To simulate the differentially expressed genes, we set 𝑐𝑠𝑡 =1in the treatment condition and 𝑐𝑠𝑡 =0in the control condition. To simulate non-differentially expressed genes, we set 𝛽𝑔1=0for both treatment and control. In the simulation, 𝜋𝑔𝑠𝑡 =0.6,∀𝑔, 𝑠, 𝑡.𝛽0was fixed to 1and the dispersion parameter 𝜏𝑔=𝜙𝑔=1,∀𝑔. Given the simulated data, we compare SpatialDENet with commonly used differential analysis methods in the literature. Specifically, we considered DESeq23, EdgeR4, Limma5, MAST6, and monocle315. Figure 1a shows the randomly selected spatial location with the spatial domain (0,1)2. Using 𝑓(𝑠, 𝑡), the figure shows that the spatial pattern exhibits higher intensity moving closer to the center of the spatial domain. This indicates that genes that are more expressed in the center region tend to be relatively more affected by spatial confounding issues. Using the spatial effect 𝑓(𝑠, 𝑡)in the stochastic representation in Equation (11), differentially and non-differentially expressed genes are simulated. Figure 1b shows the SpatialDENet representation of the spatial domain. Each node is a cluster of neighboring points, and the size of the network determines the number of neighboring spatial locations present in the data. Larger nodes have a higher number of spatial locations. The graph is colored by the node-averaged simulated spatial effect. (a) (b) SpatialDENet EdgeR DESeq2 MAST Monocol3 False Discovery Rate Power (c) Figure 1: Simulation study: (a) The simulated spatial confounding variable distributed across (0,1) in both x and y axes. (b) The network representation of the simulated spatial region. (c) The power against false discovery rates of the competing models. We simulated 1000 genes in total. Of those, 100 (10%) differentially expressed genes between control and treatment conditions were simulated. The plot of the power against the false discovery rate between Monocle 3, MAST, DESeq2, EdgeR, and spatialDENet is shown in Figure 1c. The 7 power and the FDR were computed using Power =𝑇 𝑃 𝑇𝑃 +𝐹𝑁 ,FDR =𝐹𝑃 𝑇𝑃 +𝐹𝑃 ,(12) where 𝑇𝑃 is true positive, 𝐹𝑃 is false positive, and 𝐹𝑁 is false negative. The results showed that DESeq2 achieved a power of 0.80 at an FDR of 0.005, performing relatively well compared to other methods. However, spatialDENet outperformed DESeq2, attaining a power of over 0.82 at an FDR of 0.01 and exceeding 0.90 at an FDR of 0.02. In contrast, DESeq2 reached only 0.81 power at the same FDR level of 0.02. These findings highlight the importance of accounting for spatial confounding effects when identifying differentially expressed genes. We further conducted 30 independent simulation runs to compare spatialDENet with other competing models using FDR and power as primary evaluation metrics. For each simulation scenario, we calculated the FDR, power, and the corresponding composite metrics, F1 score, and Youden’s J statistic, to assess the overall performance of the models comprehensively. The F1 score and Youden’s J statistic are defined as F1 =2(1−FDR)Power (1−FDR) +Power,Youden’s J =Power −FDR.(13) Both the F1 score and Youden’s J statistic integrate power and FDR to provide a unified measure of a method’s true performance. Methods that achieve higher values for these metrics are considered to perform better overall, as they effectively balance sensitivity and false discovery control. 8 YoudenJ F1 FDR Power (a) (b) (d) (c) Figure 2: Benchmarking analysis in synthetic data: (a) False discovery rate, (b) Power, (c) unified F1 score, (d) Youden’s J statistic. Figure 2a & b show the box plot of the FDR and power of DESeq2, EdgeR, Limma, MAST, Monocle, and SpatialDENet. Findings showed that DESeq2 and SpatialDENet had similar levels of FDR, which is better than other competing models. However, judging by the power, SpatialDENet outperformed DESeq2 and other competing models (Figure 2b). Moreover, Monocle had the highest 9