Addressing the mean-variance relationship in spatially resolved transcriptomics data with spoon

Abstract

An important task in the analysis of spatially resolved transcriptomics data is to identify spatially variable genes (SVGs), or genes that vary in a 2D space. Current approaches rank SVGs based on either p-values or an effect size, such as the proportion of spatial variance. However, previous work in the analysis of RNA-sequencing identified a technical bias, referred to as the "mean-variance relationship", where highly expressed genes are more likely to have a higher variance. Here, we demonstrate the mean-variance relationship in spatial transcriptomics data. Furthermore, we propose spoon, a statistical framework using Empirical Bayes techniques to remove this bias, leading to more accurate prioritization of SVGs. We demonstrate the performance of spoon in both simulated and real spatial transcriptomics data. A software implementation of our method is available at https://bioconductor.org/packages/spoon.

Keywords: Gaussian process regression; empirical Bayes; mean-variance bias; spatial transcriptomics; spatially variable gene.

Figure 1:. Calculating precision weights for individual observations.

These data are from Invasive Ductal Carcinoma breast tissue analyzed with 10x Genomics Visium [39], hereafter referred to as “Ductal Breast”. (A-C) The square root of the residual standard deviations estimated using nearest neighbor Gaussian processes ($\sqrt{s_g}$ defined in Equation 3) are plotted against average logcount ($\tilde{\mathit{r}}$). (B) Same as A, except a spline curve (purple) is fitted to the data to estimate the gene-wise mean-variance relationship. (C) Using the fitted spline curve, each predicted count value $({\hat{\lambda }}_{gi})$ is mapped to its corresponding square root standard deviation value using $\mathrm{s}\mathrm{p}\mathrm{l}{({\hat{\lambda }}_{gi})}^{-4}$.

Figure 2:. Mean-variance relationship exists in spatially resolved transcriptomics.

Using data from different human tissues, in order from top to bottom: DLPFC [15], Ductal Breast cancer [39], HPC [16], LC [42], and Ovarian cancer [43], we quantified the mean-variance relationship. Each point is a gene colored by the likelihood ratio statistic for a test that compares the fitted model against a classical linear model for the spatial component of variance using a NNGP [9]. The likelihood ratio statistics (LR Stat) are scaled by the maximum likelihood ratio statistic for each dataset in order to have more uniform visualization. The x-axis is mean logcounts and the y-axes represent different components of variance, in order from left to right: total variance ${\sigma }^2+{\tau }^2$, spatial variance ${\sigma }^2$, nonspatial variance ${\tau }^2$, and proportion of spatial variance ${\sigma }^2/\left({\sigma }^2+{\tau }^2\right)$.

Figure 3:. Mean-rank relationship exists in spatial transcriptomics data.

Using three datasets, in order from top to bottom (DLPFC [15], Ovarian cancer [39], and Lobular Breast cancer [40]), we quantified the mean-rank relationship. The genes were binned into deciles based on mean logcounts. Decile 1 contains the lowest mean expression values. The x-axis represents the rank. Within each decile, the density of the top 10% ranks is plotted as the signal in blue, while the density of the remaining ranks is plotted as the background in orange. Each subfigure shows the mean-rank relationship that persists after applying each method, from left to right: Moran’s I [47], nnSVG [9], SPARK-X [11], SpaGFT [45], and SpatialDE2 [46].

Figure 4:. Spoon removes the mean-variance relationship when detecting spatially variable genes.

This dataset consists of 1,000 simulated genes across 968 spots using a lengthscale of 100. Separately for unweighted and weighted methods, the genes were binned into deciles based on mean logcounts. Decile 1 contains the lowest mean expression values. Ridge plots for the (A) unweighted ranks and (B) weighted ranks are shown. Within each decile (y-axis), the density of the top 10% of ranks is plotted as the signal, while the density of the remaining ranks is plotted as the background. (C) False discovery rate (FDR) as a function of Type I error ($\alpha$). As a function of FDR, we show the (D) true negative rate (TNR) and (E) true positive rate (TPR). The red represents weighted nnSVG and the blue represents unweighted nnSVG. These plots represent the average performance across five iterations of the same simulation, each with unique random seeds.

Figure 5:. Spoon helps to detect SVGs associated with cancer that are lowly expressed.

We used four datasets to evaluate the detection of cancer-related genes: Subtype Breast cancer [41], Ovarian cancer [43], Lobular Breast cancer [40], and Ductal Breast cancer [39]. Each bar contains the intersection of the set of genes of interest with genes within the set associated with cancer. For the first four rows, we defined low mean genes as those with means less than the 25th percentile in the dataset. Within the set of low mean genes, we found genes that were in the lowest 10% of ranks before weighting and then increased to the highest 10% of ranks after weighting. This is the set of genes of interest. The intersection in orange is the number of low mean and higher ranked genes that were found to be associated with the cancer of the dataset. For the last four rows, we defined small lengthscale genes as those with lengthscales between 40 to 90. Within the set of small lengthscale genes, we found genes that were ranked higher after weighting. This is the set of genes of interest. The intersection in orange shows the number of small lengthscale genes that were ranked higher and found to be associated with the cancer type of the dataset.