Detecting Spatial Dependence in Transcriptomics Data using Vectorised Persistence Diagrams
Katharina Limbeck Bastian Rieck







































Detecting Spatial Dependence in Transcriptomics Data  
using Vectorised Persistence Diagrams
    
  Katharina Limbeck

  Helmholtz Munich, Technical University of Munich

  
  Bastian Rieck

  Helmholtz Munich, Technical University of Munich, University of Fribourg

  

    
=================================================================================================================================================================================



§ ABSTRACT 

Evaluating spatial patterns in data is an integral task across various domains, including geostatistics, astronomy, and spatial tissue biology.
The analysis of transcriptomics data in particular relies on methods for detecting spatially-dependent features that exhibit significant spatial patterns for both explanatory analysis and feature selection.
However, given the complex and high-dimensional nature of these data, there is a need for robust, stable, and reliable descriptors of spatial dependence.
We leverage the stability and multi-scale properties of persistent homology to address this task.
To this end, we introduce a novel framework using functional topological summaries, such as Betti curves and persistence landscapes, for identifying and describing non-random patterns in spatial data.
In particular, we propose a non-parametric one-sample permutation test for spatial dependence and 
investigate its utility across both simulated and real spatial omics data. 
Our vectorised approach outperforms 
baseline methods at accurately detecting spatial dependence. 
Further, we find that our method is more robust to outliers than alternative tests using `Moran’s I.' 




mkbothgobbletwo
oddheadempty
evenheadempty




§ INTRODUCTION


Omics data, that is measurements of molecules within an organism or a cell, provide crucial insights into biological processes and cellular function. 

Specifically, transcriptomics plays a vital role in studying RNA transcripts i.e. the expression of genes via the transcriptome and is thus capable of providing snapshots of cellular function and the mechanisms behind health and disease.

Recent advancements in spatial transcriptomics further enable the large-scale analysis of cells’ gene expression and their spatial locations across tissue samples. 
Leveraging these data, the detection of spatially variable genes (SVGs) is an important analysis step both for explanatory analysis and further downstream tasks, such as spatial domain identification and gene enrichment analysis, uncovering the spatial components of tissue biology <cit.>.  
Nevertheless, accurately identifying the spatial patterns in 
 omics data on both local and global scales remains challenging, particularly due to the high-dimensional, sparse and noisy nature of these data. The difficulty of modelling  
technical and systematic errors, the ongoing debate on even basic preprocessing stpdf, as well as an ever growing understanding of the fundamental processes in cellular biology, therefore necessitate the  development of robust and expressive analysis methods <cit.>. 
In particular, SVG detection methods should be robust to technical noise and sparsity, correctly identify spatially dependent genes, be independent of the sum of  expression values, and show low 
false positive error rates <cit.>. In practice, reliable SVG detection methods then help to identify the spatial patterns of disease e.g. by detecting spatial clusters of cancer cells or necrotizing tissue within an organ as well as identifying the genes driving these spatial changes <cit.>.


Meanwhile, topological data analysis (TDA) has given scholars a new toolbox for quantifying complex structures in biological networks at multiple scales motivated by the rigorous mathematical study of shapes <cit.>. 
Indeed, when adjusting for noise, molecular data in particular is equipped with meaningful coarse geometry 
<cit.>. 
Capturing both local and global characteristics of the shape of data, TDA is thus poised 
to quantify spatial patterns in molecular data. 
Computing persistent homology on spatial graphs gives stable and flexible summaries of the distribution of spatial features <cit.> summarising both local and global patterns in gene expression.  
In particular, our work demonstrates how toplogical descriptors can be used to distinguish spatial signals in transcriptomics data and identify genes whose expression varies across space.

Our contributions are as follows:
 

    
  * We propose a one-sample randomised permutation test for spatial dependence using functional topological summaries. 
    
  * We investigate the performance of this approach on both simulated and real data in comparison to alternative methods, Moran's I,  and .








§ RELATED WORK








 §.§ 
Spatially-Variable Gene Detection
 



The field of spatial omics data analysis has grown rapidly in the last 10 years fueled by advancements in sequencing technologies and data availability <cit.>. The term spatial omics denotes spatially-resolved molecular measurements in general, including features extracted from the genome, transcriptome, or proteome. These features are greatly relevant as they give insights into fundamental biological processes. As a key mechanism in biology, the set of DNA encoded in the genome is read out or transcribed by RNA molecules, which make up the transcriptome. Some of these transcripts later code for proteins, which are of great practical importance as they perform the actual biological functions in cells. Omics measurements then give key insights into these fundamental building blocks of biology. 
Depending on the sequencing method, observations from spatial datasets either correspond to individual cells, enabling the advent of single cell omics, aggregations of multiple cells or even sub-cellular measurements. Associated spatial locations are given by unstructured i.e. continuous coordinates on a tissue sample 
or spots on a pre-defined grid. 
Most frequently, features are transcriptomic measurements that  summarise the number of times a gene as measured by the number of RNA fragments detected at each capture location <cit.>, the focus of our investigation. 

Generally, SVG detection methods fall into numerous categories. Moran's I is the most established statistical measure for spatial auto-correlation across spatial statistics <cit.>. As well as its cousin Geary's C <cit.>, Moran's I is frequently used as a versatile descriptor of spatial patterns in gene expression values <cit.>. 
Despite its simplicity as a mean summary based on weighted variance and covariance estimates, and its well-established theoretical benefits, 
Moran's I is 
known to be sensitive to outliers and noise <cit.>, a disadvantage given the noisy  nature of spatial omics data <cit.>. 
Further, based on the rich tradition of spatial statistics, there also exists a variety of statistical methods that fit spatial regression models to asses the importance of spatial covariance on gene expression. The most  popular of these approaches include  and , which fit Gaussian regression processes, and  or  fitting mixed and non-parametric models respectively <cit.>. However, these models rely on fixed assumptions on the nature of the distribution of gene expression values, which remains a contested question in practice <cit.>. 
Lastly, there exist non-parametric methods, such as 
, which models gene expression as 
a marked point process and uses permutation testing to make inferences <cit.>. This approach benefits by allowing statistical reasoning under less stringent assumptions but is associated with increased computational costs. 
Further, 
model-free frameworks such as , which measures diffusion times, or , which models graph Laplacians, have also been proposed <cit.>. 
However, despite the appeal of such geometric and model-free approaches, multi-scale descriptors, such as persistent homology have not yet been employed in this context. 
Motivated by the use of geometric descriptors of spatial graphs, we thus aim to further develop persistent homology as a framework for spatial significance testing for spatial omics data analysis.




 §.§ Persistent Homology



Originating from algebraic topology, persistent homology (PH) enables practitioners to measure topological features from data <cit.> thus describing the shape of data <cit.>. 
This is of particular importance for molecular data, which after accounting for noise can be surmised to contain meaningful coarse geometry <cit.>. PH then gives one such a tool for encoding this geometry in a robust manner.

From a signal processing perspective, persistent homology has demonstrated its utility for detecting peaks in noisy signals guaranteeing stability properties and robustness to small perturbations <cit.>. The multi-scale nature of filtrations allows practitioners to not only summarise local but also global patterns in spatial data. 
Persistence-based methods can thus be applied to detect and encode both local and global extrema alike <cit.>.
Indeed, the use of persistent homology as spatial descriptor has been extended to investigating time-varying data and spatial clustering <cit.>. 
Thus, persistence has shown its utility for describing spatial patterns, such as local maxima and local minima, in various contexts. It is then only natural to consider the promising potential of using persistence for SVG detection in particular.

From a statistical perspective, 
the increasingly widespread use of PH for data analysis and budding applications to ML have also been driven by advancements in the statistical theory underpinning the theoretical foundation of PH. In particular, summaries of PH allow for the definition of parametric and non-parametric significance tests and uncertainty estimation. Amongst these approaches, non-parametric techniques, specifically permutation testing, are by far the most utilised and assumption-free methods used for statistical inferences <cit.>. 
Arguably, one-sample testing remains somewhat under-explored <cit.> with most common statistical approaches focusing on two-sample comparisons. 


Nevertheless, recent work by <cit.> demonstrate the use of persistent homology 
for the identification of spatially dependent features in geographical data. In particular, they propose a one-sample permutation testing approach using persistence summary statistics  
<cit.>. 
Their finding suggest that 
persistent homology performs well at testing for spatial dependence in the presence of outliers or on sparsely connected graphs resembling ladder graphs <cit.>. 
Further, they conclude that there is potential for future work to investigate the use of alternative functional summaries of persistent homology. 
Addressing this outlook and widening the application to spatial omics data, we now investigate the use of persistence curves for SVG detection.




§ BACKGROUND ON PERSISTENT HOMOLOGY





 §.§ Persistent Homology and Vertex-based Filtrations on Spatial Graphs
 

Throughout this paper, we set out to study a spatial dataset X = {x_1, x_2, ..., x_n} whose observations are assigned feature values f: X → and spatial locations s: X →^2. 
We then want to examine the spatial distribution of such a feature, which is modelled by constructing a spatial (neighbourhood) graph G. This approach is appealing because adjacency relationships play a central role in spatial data analysis and we can flexibly construct a graph that best describes the nature of the data as also commonly done for analysing spatial omics data. See <ref> for a more detailed description on our method.


 
Persistent homology can then be computed from such a spatial graph 
by first considering the clique-complex S=(G) defined as the simplicial complex made up of cliques of vertices in G. Here, a clique complex is a generalisation of a graph that consists of 0-simplices corresponding to vertices, 1-simplices corresponding to edges, 2-simplices corresponding to 3-cliques (i.e. triangles), and so on. In particular, this  simplicial complex is defined to be a set of simplices that is closed under taking intersections and face decomposition. 
In practice, we can consider S as a generalised triangulation of a spatial graph. 

As a multi-scale summary of geometry, persistent homology then studies a filtration on this complex. That is, we construct a sequence of topological spaces, and track the evolution of topological features given by homology groups. 
The added algebraic structure allows for the computation of boundaries and holes from the generalised graph. We can thus track the emergence and disappearance of connected components, that is 0th dimensional homology groups, H_0(S), and more generally k-dimensional homology groups, H_k(S), representing k-dimensional holes <cit.>. This process as well as the further details explained below are illustrated via an example in <ref>.

While there exist a multitude of methods for defining filtrations on graphs <cit.>, we are interested in one specific construction throughout this work. In particular, we track the evolution of superlevel-sets defined from the spatial vertex feature, f, to investigate the spatial distribution of this feature. 
That is, we filter through the simplicial complex as follows


    S_δ = Cl(G_δ)= {s ∈  Cl(G) | f(u_i) ≥δ for all  u_i ∈ s}


so that the filtration parameter δ decreases from max(f) to min(f). Here, s is a simplex in the clique complex S_δ of the subgraph G_δ∈ G defined as G_δ = {u_i ∈ G | f(u_i) ≥δ}. 

The collection of simplicial complexes S_max(f)⊂ ... ⊂ S_min(f) = Cl(G) then gives a Vietoris-Rips filtration. 
Intuitively, this vertex-based filtration can be constructed by 
filtering through 
the feature values starting with an empty simplicial set for δ > max(f). 
Then, as δ decreases, a vertex  x_i i.e a zero-dimensional simplex is included whenever its associated vertex feature value f(x_i) has been reached. Further, an edge i.e. a one-dimensional simplex, is added to the filtered simplicial complex S_δ whenever both its vertices have been included <cit.>, with the process generalising to higher dimensions. Note that to define a sublevelset-set filtration, this construction can be modified, so that δ increases from min(f) to max(f). Further, for data collected on a rectangular grid, cubical persistence, which filters through cubical complexes made up of cubical cells, offers a natural alternative to the triangle-based simplicial complex construction. See <cit.> for more details on how to generalise to this setting.

The output from computing these filtrations results in a set of persistence diagrams, {_0, _1}, one for each homology dimension, where each _· = { (p_1, q_1), (p_2, q_2), ..., (p_m, q_m)} is a multi-set of persistence pairs in ^2 encoding the birth and death, i.e. the emergence and merging scales of topological features in the filtration. 
In particular, peaks or local maxima in spatial feature values are represented by points in the 0-th dimensional persistence diagram above the diagonal <cit.>. 
For a deeper exploration of the mathematical framework and theoretical underpinnings of persistent homology, we refer the interested reader to  <cit.>. 
We however move on to consider the use of persistence homology for statistical analyses.




 §.§ Summarising Persistence Diagrams
 

Persistence diagrams consist of multi-sets of birth- and death pairs and are thus challenging to study from a statistical perspective. In fact, there exists no natural way to directly understand them in as a Hilbert space that allows for the application of standard statistical tools. Indeed, inherent distances between diagrams such as the Bottelneck or Wasserstein distances are expensive to compute and do not induce a norm on the space of persistence diagrams. 
Potential approaches for fitting PH into a statistical setting is the use of summary statistics, as for example explored by <cit.> in the context of spatial dependency detection. Moreover, to retain the multiscale properties of PH, functional summaries can be employed to understand PDs in a Hilbert space setting. 




  §.§.§ Summary Statistics

While various one-number descriptors, such as persistence entropy or lifetime-derived statistics are available, we choose to include total lifetime (TL) into our study based on the findings of <cit.>. 
In particular, the total lifetime is 
computed as 

    = ∑_i=1^m q_i - p_i
 i.e. the sum of the lifetimes of persistence pairs, defined as the difference between the death and birth scales of topological features. 



  §.§.§ Functional Summaries


Our analysis will focus on using functional summaries that map PDs to (piecewise) continuous functions from → V where V is a suitable vector space <cit.>. Note that in related literature these types of functional summaries are often called persistence curves. In particular, we compare Betti curves 
and persistence landscapes, which can be conveniently compared via suitable norms between these functions. We hypothesise that these functional summaries retain more complete information from the whole PD than single number summaries (as e.g. studied by <cit.>) while allowing for the application of statistical methods in a theoretically well-founded manner <cit.>. 



  
Betti curves
Betti curves are defined as 
: →ℕ where 

    (δ)=|{ (p_i, q_i) ∈|δ∈ [p_i, q_i]}|.
 That is, Betti curves sum up the number of topological features from the persistence diagram that persist at each scale of the filtration. 
Given one Betti curve, b, we can then take its L^p norm as a natural descriptor:

    L_p() = 
    ||||_p = (∫_min(f)^max(f) |(x)|^p  dx)^1/p

This L^1 norm then is closely related to computing

 total lifetime in terms of the information it summarises 
<cit.>.  
To compare multiple Betti curves, _1 and _2, we naturally extend this notion and use the L^p distance between Betti curves, that is 

    L^p(_1, _2) = 
    ||_2 - _1||_p = (∫_min(f)^max(f) |_1(x)-_2(x)|^p  dx)^1/p.




  
Persistence Landscapes

A frequently used tool for statistical analyses are persistence landscapes, which have been proven to follow a central limit theorem, even enabling the application of standard parametric tools <cit.>.
First, the rank functions are defined for each persistence pair given as 

    λ_(p_i, q_i)(δ) = 
    
    δ - p_i     if  p_i < δ≤p_i + q_i/2
    
     - δ + q_i     if p_i + q_i/2 < δ < q_i
    
     0      otherwise


and the k-th persistence landscapes function is then given by 
    Λ(k,δ) = kmax{λ_(p_i, q_i)(δ)}_i ∈ I
 where kmax is the function giving the k-th largest value of a set. The persistence landscape then is the collection of these functions. Further, we can define the L^p-norm for persistence landscapes as:


    L^p(Λ) = 
    ||Λ||_p = ∑_k=1^∞(∫_∞^∞ |Λ(k,x)|^p  dx)^1/p.











§ 
METHODS





 §.§ Computing Persistent Homology from Spatial Omics Data
 

Now, to develop our spatially variable gene detection method, we 

define a spatial graphs for omics data 
in the manner that best matches 
the underlying sequencing technology <cit.>. 
Amongst these technologies, spatial measurements are collected across varying resolutions. 
For example, for data with continuous spatial coordinates, such as high resolution measurements of single cells, we compute the Delaunay triangulation to find an appropriate triangular covering of the space. 
In comparison, Visium, one of the most popular sequencing methods, measures spots on a hexagonal grid with spots having diameters of 55 μ m and centers being a fixed distance of 100μ m apart <cit.>. Motivated by this, we choose a spatial graph connecting each grid point to its six directly adjacent neighbouring spots 
for further analysis. 

For both of these graph choices we then apply simplicial-complex based filtrations. 
Lastly, for data measured on rectangular grids, we connect each spot to its four adjacent neighbours and apply cubical persistence to capture the rectangular nature of these grid cells. 


Throughout this work we focus on computing zero-
dimensional persistent homology, which allows us to investigate the spatial patterns given by connected components. 
In particular, we are interested in encoding peaks in spatial gene expression and therefore construct super-levelset filtrations using gene expression as vertex feature as described in <ref>. Then, the resulting persistence diagrams encode key information on clustering behaviour. Here, we choose to use a super-levelset approach rather than a sub-levelset filtration because we reason that spatial maxima will be more suitable for revealing spatial trends than spatial minima for these data. This choice is informed by the fact that omics data is 
affected by technical noise and sparsity. 
Indeed, dropout is a well known measurement error describing the common phenomena of not observing a transcript even though it is present. In practice, this effect leads to a large proportion of excess zeros in gene expression values <cit.>, which could negatively affect the interpretation of local minima for our purpose. We therefore argue that spatial peaks will be more robust to such sparsity.

For each gene, whose spatial pattern we want to assess, we then compute persistence diagrams and persistence summaries as detailed in <ref>.
By convention, we bound death value of each persistence pair by minimum filtration value. That is, rather than allowing one persistence feature to achieve an infinite death value, we bound it by the lowest gene expression value observed in the data. 
Note that our construction is similar to <cit.>, but we further extend their methods to also consider cubical persistence for rectangular grids and investigate the use of functional summaries for spatial dependence testing as detailed below. 




 §.§ One-Sample Permutation Testing for Spatial Dependence


We use randomized permutation testing to test for the following hypotheses:

    
  * H_0: Expression values are randomly distributed in space
    
  * H_1: Expression values are not spatially random


We thus retain the original gene expression values as well as  the  spatial graph and randomly permute which expression value is assigned to each location to simulate spatial randomness, that is the lack of any spatial pattern. Then, significance is assessed as follows:



    p-value =  (||_̅0̅ - _i||_p ≤ ||_̅0̅ - _obs||_p) + 1/n_perm + 1
 where  is the indicator function summing up the number of times the L^p-norm, ||·||_p, between the observed Betti curve _obs and the empirical mean Betti curve under the null model _0 is more extreme than the difference between each individual Betti curve under the null model and their overall mean. 

Note that we compute _0, the mean Betti curve given the null hypothesis, as the mean across all n_perm permutations as well as the observed Betti curve while also adding 1 to both the nominator and denominator. 

This bias correction is included to avoid zero p-values and bias the test towards more conservative estimates by including this pseudo-count <cit.>. 
We thus conduct a two-sided, one-sample hypothesis test for each spatial feature. Whenever multiple genes are compared simultaneously, we further recommend to apply a suitable multiple testing adjustment such as  

Benjamini-Hochberg 
procedure to control the false discovery rate <cit.>, which we will apply throughout our experiments. Across our experiments, the number of permutations is set to 1000 giving a practical compromise with regards to computational efficiency considerations. 


Beyond statistical reasoning, there also exist 
more deterministic approaches to spatially variable gene discovery. 
Here, the goal is to either rank features by the strength of their spatial dependence, or to find a predefined number of the top most spatially variable features. 
In practice, these approaches are key for both explanatory analysis and for feature selection as a first step to uncover and understand spatial patterns in data <cit.>. 
 For further analysis, we thus choose to rank spatially expressed features by 
their empirical p-values following similar statistical approaches <cit.>. 



 §.§ Alternative Spatial Variability Detection Methods



In order to evaluate the performance of our spatial testing approach we compare our methods with three standard methods for detecting spatial dependence. 



  
Moran's I: Moran's I is one of the oldest and most widely-used measures for spatial auto-correlation. It is known to be a suitable baseline for SVG benchmarking <cit.> and can be computed as 


    I=n_obs/∑_i,j w_i,j∑_i=1^n_obs∑_j=1^n_obs w_i,j (f_i - f̅) (f_j - f̅)/∑_i=1^n_obs (f_i - f̅)^2
        .


We compute it using 
the spatial adjacency matrix to define the weights w_i,j. The value of Moran's I is directly interpretable, ranging from -1 indicating negative autocorrelation, to 0 corresponding to spatial randomness, and 1 indicating positive autocorrelation.
Further, for statistical analysis we use Moran's I in the same randomised permutation testing framework as the previously introduced topological summaries explained in <ref>.





  
sepal:   simulates diffusion processes, which models a stochastic process that tracks the spread of a feature on the spatial domain over time as it tends towards randomness. Diffusion is simulated on either a rectangular or hexagonal grid and  then measures diffusion times, which are used to assess spatial patterns and measure how long it takes the diffusion process to converge. The higher the diffusion time the stronger the spatial signal. Note that this method does not conduct any statistical testing and thus does not provide p-values. Instead, it provides a ranking of genes by their diffusion times. 
Moreover,  is only applicable to hexagonal or rectangular grids. Hence, whenever data is given with unstructured coordinates, we follow the implementation by <cit.> to convert spatial measurements  to a rectangular grid. 



  
SpatialDE: Finally,  applies spatial covariance testing via Gaussian process regression, an class of models originating from geostatistics. In a nutshell, gene expression is modelled via a multivariate normal model that includes both a non-spatial and a spatial covariance term.  
The model fitted with spatial covariance is then compared to a null model of spatial independence without this spatial component. 
Likelihood testing is used and p-values are computed analytically using the χ^2 distribution with one degree of freedom. As a multiple testing correction, we follow <cit.> and report q-values, which specifically control for the positive false discovery rate <cit.>.  assumes normally distributed residual noise, we follow the preprocessing stpdf by <cit.>. First, Anscombe’s transformation, a variance stabilising method for negative binomial data, is applied. Second, log total count values are regressed out  
to achieve independence from the total counts per spatial location.




 §.§ Datasets and Preprocessing




We use three real transcriptomics datasets for our investigation. First, we analyse spatial data of a mouse olfactory bulb (referred to as MOB, replicate 1) and a breast cancer tumour sample (breast, replicate 1) from <cit.>, standard datasets commonly used to test spatial variability detection methods. Both have been sequenced using 1k spatial transcriptomics arrays, which measure spots on rectangular grids. Each spot has a diameter of 100 μ m and encompasses 10-100 cells. 
Further, we study one sample from <cit.> (sample IZP10), which consists of spatial measurements of spots on a hexagonal grid. This sample has been taken from the heart of a patient that suffered from a heart attack from tissues with necrotic areas. 

Prior to analysis, peprocessing and quality control of omics data are crucial stpdf for accounting for technical errors as well as filtering out noise. Indeed, normalisation is a key step in spatial omics data analysis as many downstream tasks, including the detection of spatially variable features, rely on adequate preprocessing choices to account for technical effects <cit.>.

Following best practices <cit.> and alternative spatial variability studies <cit.>, we decide on the following pre-processing strategy and
select all genes that 
   have a minimum total expression count of 10, and
 
 are present in at least 1% of locations. 
Further, we remove all cells with fewer than 10 counts. These stpdf are done to avoid analysing cells or gene for which hardly any biological signals have been detected. Finally, 
we exclude all ribosomal and mitochondrial genes from analysis as these indicate technical artifacts and are not related to the mRNA, i.e. protein-coding genes that are the focus of transcriptomics data analysis <cit.>. 
For each real dataset, we then identify the 500 most variable genes using  
for further analysis. That is, we select a set of genes that vary most notably across cells. Our goal will then be to detect whether these genes that differ noticeably between dislocated cells also show distinct spatial patterns. Note that for the simulated data we skip these quality control stpdf as we want to compare all SVG detection methods across all simulated features. 


Further, as conventionally done in spatial omics data analysis, we decide to pre-process the expression values i.e. count data 
in the manner best fitting the SVG detection method <cit.>. For  we apply the normalisation strategy preferred for this method as detailed in <ref> <cit.>.
For , we apply the recommended shifted logarithm normalisation by applying the transformation log(f+2) to each feature f. This transformation is applied to 
stabilize the variance across the dataset to counteract variable sampling effects such as overdispersion and reduce the skewedness of the feature's distribution. The pseudo-count of 2 is added to enable the log-transformation also for zero-values counts. To ensure comparability across spatial variability tests, we similarly choose to compute our persistence-based test as well as Moran's I from the same log(f+2)-transformed features.  
Note however, that this pre-processing strategy is not the only viable choice and our method can be flexibly applied across varying pre-processing choices <cit.>.
In general, given that our persistence-based methods directly scales through each gene's expression values, we highly recommend to apply our framework to transformed and normalised features, but we leave the choice of the best data-dependent transformation to further investigation <cit.>. 





 §.§ Synthetic Data Generation





We create in-silico examples of spatial transcriptomics data using SRTsim <cit.>, an R package for the simulation of spatial patterns. In particular, we randomly sample spatial locations on a square domain and then simulate count data using (zero-inflated) Poisson or negative binomial distributions. 
The choice between these range of distributions of count data is motivated by the fact that the observed distribution and nature of transcriptomics data is highly debated and varies between sequencing technologies and datasets. 
Zero-inflation aims to account for the excess zeros often encountered in sparsely sampled gene expression data due to a high dropout probability. 
Across our simulations we thus consider zero-inflation 
to be a potential source of noise and assess the robustness of the proposed spatially variable gene detection methods against these technical errors. 
SRTsim allows us choose the following settings for simulations. Globally, we set the mean parameter (μ = 1), the dispersion parameter (s) and the zero-proportion parameter (z). 
We then define N spatial domains corresponding to distinct spatial shapes that are assigned an effect size x_i resulting in an x_i-fold increase in mean for the specified spatial domain i ∈0, 1, 2, .., N. The effect size of the background domain is set to x_0=1. Count data is then sampled separately for each domain with the thus specified parameters and assigned randomly to locations within the shape. 
We then simulate four distinct spatial patterns, a gradient, two lines, a set of clusters and a cellring illustrated in <ref> with the annotated effect sizes. Otherwise, if no spatial signal is simulated, gene expression is sampled randomly across the locations.





§ RESULTS


Across our experiments we consider two main goals: (i) using significance testing to determine spatially dependent features, (ii) 
rank features by their spatial dependence and identify a fixed number of the top most spatially variable genes. 
Throughout our empirical experiments we then demonstrate the following results:

    
  * PH-based approaches perform well at detecting SVGs in terms of the area under the precision recall curve and are more robust to zero-inflation i.e. high proportions of excess zeros in gene expression counts than Moran's I. 
    
  * PH-based significance testing leads to lower sensitivity but higher specificity identifying a smaller set of SVGs at a fixed 0.05 significance threshold. 
    
  * PH-based approaches are superior to alternatives at detecting a predefined number of the top most spatially variable genes. 
    
  * PH-based approaches capture orthogonal information to alternative methods.
    
  * PH-based spatial dependence test are less correlated to the total sum of feature values than Moran's I. 

We further explore the use of our PH-based approach to detect SVGs for three real spatial omics datasets, for which we plot and examine the top spatially variable features.  

We thus demonstrate that persistence-based methods, and Betti curves in particular, offer valid alternatives to existing SVG methods, such as Moran's I. 






 §.§ Simulation Study






  §.§.§ PH-based approaches are robust to zero-inflation


We simulate spatial data as described in <ref> and first vary the zero-inflation parameter z from 0.1 to 0.99 keeping the dispersion fixed at 0.3. For each shape and parameter choice, we then simulate 50 genes that possess the spatial signal 
and 50 genes that are randomly sampled without any spatial effect. 
Then, across each degree of zero-inflation we summarise the performance of each statistical test and scoring method in the following manner: 
<ref> 
reports the area under the precision-recall curve (AUPRC) as a function of the degree of zero-inflation for all methods. 
This comparison is chosen because it allows us to compare all methods, also the ranking-based  method, as a single summary statistic per degree of zero-inflation. 
An AUPRC of 1 indicates perfect performance, random guessing is expected to result in an AUPRC of 0.5, and 0 is the worst possible value. 





  §.§.§ Persistence-based testing performs well at detecting SVGs

We observe that all persistence-based methods achieve high AUPRC values very close to 1 until around 70% of zero-inflation across all pattern. 
This is followed by a decrease in AUPRC, which remains above 0.6, even for 99% of excess zeros. Overall, these trends show robustness and expressivity of persistent homology as a summary of spatial patterns even in the presence of high zero-inflation. 



  §.§.§ Persistence-based methods perform similarly to one another

We further observe that both Betti curves and total persistence 
follow very similar trends in <ref> with Betti curves performing better on the gradient pattern and total persistence showing higher scores on the streaks pattern.
This finding highlights that both approaches are alike in terms of the information they are summarising.
Persistent landscapes meanwhile achieve consistently higher AUPRC values of above 0.7 for higher degrees of zero-inflation across datasets and thus capture trends in spatial variability in a more robust manner for high degrees of dropout.
For this experiment, we thus observe that persistence landscape summarise different aspects of persistent homology than alternative summaries.



  §.§.§ Persistence-based methods outperform Moran's I

Our results thus indicate that persistent homology is more expressive at quantifying spatial signals in simulated expression values under the presence of zero-inflation than Moran's I. 
That is we observe that by tracking the evolution of connected components via superlevelset filtrations 
on spatial graphs, we get summary statistics with higher discriminative performance than Moran's I, which summarises spatial autocorrelation based on summarising local neighbourhoods and is the most commonly used baseline measure for spatial variability. Indeed, via this comparison in terms of AUPRC, we assess that persistence-based summaries offer a valid alternative to established tools from spatial statistics.



  §.§.§ Persistence detects spatial variability across all simulations
 
We note that our proposed persistence-based tests do well at detecting spatial variability, in particular for low levels of zero-inflation, as would be expected. However, in comparison, we see that some alternative methods do not follow the same clear trend across datasets.  starts decreasing in terms of AUPRC at comparably lower proportions of zero inflation. This difference is especially pronounced for the streaks pattern, where  performs notably worse than all alternative methods after 60% zero-inflation. Potentially this is due to this simulation only showing two relatively thin lines that are harder for  to detect as spatially variably under the presence of excess zeros. 
Meanwhile, 
 performs very well for higher values of z between 0.6 and 0.8 across datasets, showing higher AUPCRs than persistence landscapes. However,  fails at detecting the gradient pattern at low degrees of zero-inflation. For example,  receiving an AUPRC of only 0.6 at z=0.1 highlighting that diffusion times are not suitable for detecting this smooth trend in spatial variability. We thus demonstrate disadvantages of both  and  compared to our persistence-based approach. 
We conclude that PH-based methods achieve overwhelmingly higher or similar AUPRCs to the baseline method Moran's I highlighting their utility as robust SVG detection methods.




 




  §.§.§ PH-based methods are more specific than Moran's I

Extending our observations beyond AUPRC, we next decide to view spatially variable gene detection in the context of hypothesis testing. <ref> then compares the performance of each 
statistical test in terms of specificity and sensitivity at a fixed 0.05 significance  
threshold for corrected p-values.  
In agreement with results from <cit.> we observed that 
at this fixed threshold, total persistence is either similarly or slightly less powerful i.e. less sensitive than Moran's I. This is followed by
Betti curves, which detects a lower proportion of true positives, in particular on the streaks pattern, as well as persistence landscapes which give even more conservative estimate. However, we also observe that all persistence-based approaches show very high specificity across examples, whereas Moran's I demonstrates the highest false positive rate of around 5% to 25%, even after p-value correction, which is undesirable in practice, when the set of spatially variable genes should be correctly identified. 

Our findings thus agree with <cit.>, who find that persistence is less sensitive than Moran's I. 
However, <cit.> further claim that persistence cannot replace Moran's I, a findings we do not fully support. Arguably, specificity as well as sensitivity is important for interpreting spatial dependency tests in practice. Indeed, we observe that at a default 0.05 significance cutoff, corrected p-values from PH-based approaches achieve higher specificity than Moran's I, showing that PH-based results can more useful for tasks that require a low false positive rate and benefit from a more conservative assessment. Given the high AUPRC achieved by PH-based methods, we further assess that choosing a higher significance threshold could be a suitable adjustment to achieve higher sensitivity and identify a larger set of spatially variable features while retaining superior specificity.








  §.§.§ Agreement between methods








Further, we compute the Spearman correlation between each method's ranking 
and report the results in <ref>.
For persistence-based approaches we rank genes by their
p-values, else we use the values of the scores as defined in <ref> for  and Moran's I. 
Overall, we observe that total persistence and Betti curves perform most similarly across experiments with correlations between 0.64 and 0.82, which is likely driven by their common roots in accumulating persistence features <cit.>. 
Further, we find that there is varying but overall moderate levels of agreement between persistence and Moran's I (ρ≤ 0.53), likely explained by the fact that both are summaries of the same spatial graph and understand spatial variability via modelling clustering behaviour, as summarised by connected components in PH, or autocorrelation, as summarised by differences in neighbourhood scores. 
Interestingly, while  follows a fundamentally different modelling approach, correlation between  and the other approaches 
exist at comparable levels to other entries in the correlation matrix and the  score is actually most correlated with Moran's I on the clusters dataset (ρ = 0.58). Overall, we thus observe that persistence-based permutation testing captures complementary information not encoded by either Moran's I,  or . 




 §.§ PH-based approaches identify the top SVGs



Next, we keep the zero proportion at z=0.1 and the dispersion parameter at s=0.3, following estimations by  <cit.> from real data, but vary the effect size, by taking x̂_̂î = max(x_i / e, 1) for some e in { 6, 5, 4, 3, 2, 1} for each domain i effectively varying the strength of the spatial pattern compared to the background domain for which x_i=1. For each parameter choice and spatial shape we then simulate 50 features with spatial signals and 50 completely random patterns. 
We then investigate the identification of top spatially variable genes for different strengths of spatial effects.





  §.§.§ Correlation with the sum of gene expression counts

As reported in <ref>, we explore the correlation between each spatial variability score and the total count per gene / i.e. the sum of feature values. 
We want an ideal SVG detection method to not just correspond to these total counts, but score genes based solely on the existence of spatial dependencies in gene expression. 
Across datasets, we see that Moran's I and  show the highest correspondence to the sum of feature values. In comparison,  persistence-based ranking achieve consistently lower Spearman correlation of 0.6 to the number of total counts. Persistence landscapes reaching slightly higher correlation than Betti curves and total persistence. Overall, this rank correlation is similar to the one achieved by . Therefore, the comparison highlights that persistence offers a satisfactory alternative to established SVG tests.



  §.§.§ PH-based approaches correctly identify the top SVGs

Across methods, we then compare the ground truth 
given by the existence of spatial signal with each spatial variability score.  <ref> shows the proportion of truly positive features detected by each method. Ideally, a suitable spatial variability detection method would show a high proportion of true positives amongst the top 300 genes (given 300 out of 600 features show spatial patterns). 
Notably, Betti curves perform consistently well across examples correctly identifying more than 90% of top SVGs. This is followed by total persistence and persistence landscapes, which performs better on the gradient but worse on the cluster pattern. 
Interestingly, from these examples we see that while the proportion of genes correctly identified is high across  scenarios, there is no strong distinction between either persistence- or non persistence based methods based on these examples. However,  and Moran's I seem to perform slightly worse at identifying the most spatially variable features than persistence-based approaches, which together with their correlation to total counts, indicates that persistence-based testing is preferable in this scenario. 











 §.§ Analysing Transcriptomics Data






Having observed that persistence-based approaches give suitable spatial variability tests on simulated data, we next investigate their utility on real datasets as detailed in <ref>.




  §.§.§ PH-based methods identify spatially variable genes


First, we decide to focus on one of our proposed methods, namely Betti curves, which performed well across our simulated case study, and investigate spatially variable genes detected by this approach. Then, <ref> shows representative examples of genes that are both the most highly variable and the most spatially variable features reaching the lowest possible p-values in terms of the spatial dependence test. We see that Betti curves identify both genes with low and high expression counts to be spatially variable. 
For the mouse olfactory bulb (MOB), we see that there is a clear pattern of genes being more highly expressed and clustered in the center of this brain region, such as Mbp, Pcp4, Nrgn, or Gpsm1. What these genes have in common then is that they relate to neural activity in this region of the mouse's brain. For the beast cancer sample we see some more spares patterns that could indicate disease specific differences in tissue regions and immune response. Finally, for the heart sample the genes COL1A1 and COL1A2 are found to be spatially variable, which are coding for collagen type 1 mainly found in connective tissue and related to wound healing. Overall, studying the set of spatially variable genes in more depth can thus reveal important insights into tissue biology and disease response by specifically focusing on transcripts that vary across tissue regions.



  §.§.§ PH-based methods identify fewer SVGs than Moran's I

Across the 
real datasets considered for this study, we  report the number of spatially variable genes detected by each method. Results in <ref> then confirm that Moran's I reports the largest number of significantly spatially auto-correlated genes calling more than 50% of transcripts spatially variable. This is followed by Betti curves, which find more spatially variable patterns for the heart sample, but fewer for the mouse brain and breast cancer sample. Interestingly, across these real datasets, total persistence detects a low number of SVGs, which is comparable to Persistence landscapes, but  notably lower than Betti curves. This trend has not been observed across our simulations indicating that the results of using Betti curves or total persistence are not necessarily similar. 
Overall, the reported numbers of SVGs support our findings from the simulation study in <ref> on the sensitivity and specificity of each approach that lead us to conclude that persistence-based methods identify a smaller number of genes than Moran's I but do so with a lower type-1 error rate. 




  §.§.§ Persistence is robust to the sum of feature values

Next, we investigate the dependence of each score on the total counts measured for each gene as reported in <ref>. Ideally, spatial autocorrelation test evaluation also be somewhat robust of whether features have high or low values overall. That is, genes that are found to be spatially variable should not be baised towards genes with high expression levels, but also detect spatial patterns in lowly expressed genes.  
Using the Spearman correlation with total counts as a sanity check, we see that amongst all methods, Moran's I generally shows the highest correlation with total counts. This is followed by persistence-based methods which are more robust to the sum of feature values. Finally,  is designed to quantify spatial variability independently of total counts, but we find that in some examples it actually correlates more with total counts than persistence-based methods.



  §.§.§ Agreement between methods

Further, <ref> reports the rank correlations between the gene rankings reported by each spatially variable gene detection method. Overall, we once again see a high agreement between Betti curves and total persistence ρ∈ [0.7, 0.79]. This is followed by moderate correlation to persistence landscapes, Moran's I and  with ρ∈ [0.32, 0.69]. Persistence landscapes rank genes in a notably different order than other persistence-based approaches ρ∈ [0.42, 0.59]. 
Overall, the least agreement is reported between persistence-based methods and  with ρ∈ [0.24, 0.39]. Overall, our results point towards a relatively high disagreement between methods. Similar to the simulation study, we find that persistence summarises distinct information on spatial patterns not encoded by existing SVG detection methods. 




§ DISCUSSION


Throughout our study, we proposed randomised permutation testing procedures for spatial variability in omics data via Betti curves or persistence landscapes. We investigated these tests in comparison to baseline methods Moran'I, ,  and total persistence <cit.>. 


On simulated data, we showed that persistence is more robust to zero-inflation than Moran's I showing higher AUPCR and a better detection rate of truly variable features amongst the highest ranked genes. At a fixed significance cut-off we further find that persistence-based methods lead to slightly lower sensitivity but higher specificity making us conclude that persistence is potentially more useful than Moran's I in applications that require more conservative estimates. We thus demonstrate in-silico that persistence curves offer stable approaches to quantify spatial dependence in feature values.

On real data, we further observe that persistence is less dependent on the sum of feature values than Moran's I indicating that it summarises spatial patterns in a more robust manner. 
Finally, we visualise genes with representative spatial patterns that are both spatially and highly variable and describe their qualitative biological relevance. Thus, we describe the practical benefits of using persistence-based approaches over Moran's I and showcase examples where persistence curves are superior to 
 baseline methods.

We find that a non-parametric approach is less reliant on specific distributional assumptions about the spatial distribution of feature values, in contrast to parametric methods such as , making our approach preferable when such assumptions are not satisfied <cit.>, which arguably remains a contested question for omics data in particular. 

While our analysis focuses on illustrating examples on the use of persistent homology in the domain of spatial omics data analysis, our methodology can be generalised to other types of spatial data, such as spatial graphs with labelled vertices or edges. We thus propose a flexible method for detecting spatially variable features based on application-driven filtrations.




However, we want to also note the limitations of our study. The computational costs of non-parametric methods, especially the computation-heavy nature of permutation testing, remains a concern in practice <cit.>.
Computing a persistence curve for one spatial feature itself is extremely fast taking less than one second on a local CPU. 
However, the proposed spatial variablity tests pose some scalability constraints. Permutation testing requires thousands of repetitions resulting in a notable increase in computational costs <cit.>. While per feature this test might then take less than a minute, this scaling behaviour does not suffice when the goal is to study hundreds of thousands of genes simultaneously. Sequential permutation testing could be employed to reduce some of these costs <cit.>, but might not lead to sufficient computational improvements overall. There thus exists a great motivation for further work to develop alternative statistical tests using potentially parametric results, such as ones that could be derived from a universal null distribution for persistence diagrams <cit.> to speed up the proposed testing procedures. 


Further, we have observed a high degree of disagreement between SVG detection methods across our experiments, which is a known practical issue for benchmarking these methods <cit.>. Even across our simulation study, we find  no one spatial variability detection method that is preferable in all scenarios. Indeed, further investigation could be beneficial to gain a more complete understanding of the advantages and disadvantages of each approach. 
Nevertheless, we conclude that persistence-based approaches offer a valid alternative to conventional spatial variability detection methods, in particular Moran's I, as they achieve high AUPRC values even in the presence of zero-inflation, give more specific results than Moran's I, and show lower correlation to the total sum of feature values. 


In terms of the way persistence is used throughout this study, we investigated zeroth dimensional persistent homology 
because of its computational efficiency as well as the utility of encoding local peaks for detecting spatial patterns. However, we note that one could follow an approach similar to <cit.>, extend the filtration by its exterior and compute first dimensional persistent homology. This 
would allow to deduce information about the geographical location of local minima via tracking the generators of 1D homology classes, which could be combined with statistical methods to quantify exactly how spatial patterns vary from spatial randomness.




abbrvnat


