RESCUE: recovery of unattributed expression patterns in spatial transcriptomics

Domain: neuroscience
Published: 2026-04-10
Editorial outcome: minor_revision

Abstract
Spatial transcriptomics (ST) enables gene expression profiling while preserving the spatial architecture of intact tissue. Analyzing ST data often proceeds by first extracting cell-level information, typically through cell segmentation or cell-type deconvolution. However, a critical oversight has been that a substantial portion of molecular expression is systematically lost or unannotated by these methods. This lost expression can arise from diverse and biologically important sources like fragile or underrepresented cell types, subcellular structures like neurites, and extracellular expression. These omissions can result in biased analyses and incorrect or incomplete biological interpretations. We describe a new computational method, RESCUE, that can recover the unattributed spatial expression patterns missed by existing ST analysis methods and enable robust inference even when reference is incomplete. We validate RESCUE using MERFISH data from the honey bee brain and apply it to multiple ST datasets to demonstrate how it can reveal novel insights into complex tissue biology.

Introduction
Profiling cellular expression within its spatial context has become possible with the recent development of spatial transcriptomics (ST) technologies, which measures gene expression in intact tissue sections1,2. Analyzing ST data typically proceeds by extracting cell-level information from the spatial expression data, as cells are natural units of analysis. With high-resolution platforms such as MERFISH3, 10x Visium HD4, or 10x Xenium5, this is usually done using cell segmentation6–8. With low-resolution platforms like 10x Visium9 and Slide-seq10,11, where each spatial location can contain multiple cells, it is common to deconvolve the observed expression into separate contributions from different cell types, based on reference single-cell/nuclei RNA-seq (sc/snRNA-seq) data12–16.

However, a major problem with these existing ST analysis methods is that they can fail to capture important expression patterns. First, they may omit entire cells or cell types. Cell segmentation algorithms for high-resolution data have difficulty with cells with complex morphologies, such as fibroblasts, or cells lacking nuclei, like erythrocytes17,18. Deconvolution methods for low-resolution data cannot account for cell types that are absent from the reference sc/snRNA-seq data because they are too fragile to survive tissue dissociation19–21. Second, existing ST approaches can overlook expression patterns unique to cellular substructures. Segmentation methods simply drop transcripts located in substructures that are difficult to segment, such as the axons and dendrites of neurons. These fragile structures are also typically destroyed during tissue dissociation and therefore are not represented in reference single-cell data used in deconvolution. Indeed, it has been estimated that over 40% of the total transcriptome of an adult mouse brain may be missed by sc/snRNA-seq22. Similarly, reference snRNA-seq data, by definition, ignores cytoplasmic expression. Finally, existing methods do not capture extracellular tissue expression patterns, which cannot be studied by segmentation or deconvolution.

These difficult-to-capture expression patterns can be critically important. First, they can bias analyses that do not take them into account. An example with dramatic implications is that deconvolution estimates of cell-type proportions can be invalid if the reference is missing any cell types23, as we illustrate later in this work. Second, there may be important biological processes that are unique to the overlooked expression patterns. For example, important missed cell types include adipocytes5, which play crucial roles in metabolic disease and cancer24, and neutrophils25, which are core components of the immune system. Important missed substructural expression includes axonal or dendritic transcripts that are central to synaptic development and plasticity26–30 and diseases such as amyotrophic lateral sclerosis31,32, as well as cytoplasmic transcripts that regulate key processes such as metabolism, signaling, and stress responses33–35 and are involved in Alzheimer's disease36,37. Important missed extracellular expression patterns include extracellular RNAs that mediate intercellular communication in cancer38–40.

We describe a new computational method, RESCUE, that can recover the spatial expression patterns missed by single-cell inference on ST data. Our approach is based on making a fundamental correction to the class of statistical models currently used for integrating ST data and reference data, such as sc/snRNA or cell segmentation results. Briefly, our method partitions observed gene expression into one component explainable by the reference data, which we refer to as “canonical” expression, and another component that is not, which we refer to as “unattributed” expression. Our purified canonical ST expression component can be used to correct existing cell-type deconvolution methods. Our recovered unattributed ST expression component gives us the new ability to specifically target biology that is distinct from our reference data, as discussed above. The unattributed expression recovery step in RESCUE is analogous to negative selection in experimental biology41–44, where we identify and remove canonical expression. We demonstrate that RESCUE accurately recovers unattributed expression, enables robust deconvolution of ST data despite incomplete reference data, and reveals biologically meaningful patterns overlooked by existing methods. While we focus here on ST data, we note that RESCUE can also be directly applied to non-spatial modalities, such as bulk RNA-seq data.

Results
Overview of RESCUE
RESCUE is predicated on a new model that we propose for integrating ST and reference single-cell data, and can be applied to diverse ST platforms, including MERFISH and Visium, across different organisms. Most existing models12–16 posit that ST expression at each location is composed of a linear combination of expression patterns present in the reference, where the coefficients of this linear combination are interpreted as the proportions of each reference cell type in that location. However, as described in the introduction, there are numerous sources of ST expression that are left out of this model; we will refer to these as unattributed expression patterns, because they are not attributed to canonical expression patterns reflected in the reference.

To address this issue, we propose a corrected model (Fig. 1) that allows the ST profile to have both canonical and unattributed expression components. While it appears minor, this change has important implications. Under this model, results from the statistics literatures45–47 imply that under certain conditions it is possible to simultaneously recover unattributed expression and estimate the linear combination of reference patterns that constitute canonical expression. The conditions are that only relatively few genes show non-zero, or sparse expression, in the unattributed component, which is reasonable given the nature of unattributed expression and the read depth/detection efficiency of ST technologies (Methods).

Leveraging this insight, we developed RESCUE, an algorithm that simultaneously deconvolves ST data while rescuing spatial unattributed expression. RESCUE solves a regularized robust multivariate regression problem and can be applied to both low- and high-resolution ST data. Unlike existing deconvolution methods12–16, it is robust to incomplete references, such as snRNA-seq data which lack cytoplasmic expression profiles. We emphasize that the interpretation of the unattributed component depends on both the biological context of the tissue and the choice of the reference dataset (Methods). As a result, the unattributed component may comprise multiple biological sources, including expression from cell types missing in the reference, substructural expression, and extracellular RNA. RESCUE is designed to recover these signals jointly rather than to assign them to predefined categories. The rescued unattributed expression can then be further decomposed using standard unsupervised approaches, such as clustering and non-negative matrix factorization (NMF), to identify distinct expression patterns within the unattributed component. In addition, differentially expressed gene (DEG) and spatially variable gene (SVG) analyses can be applied to characterize the biological and spatial features of the recovered unattributed signals. Careful consideration of biological structures and tissue organization is therefore essential when interpreting the unattributed expression. We illustrate this in a broad range of examples below.

Simulations
To evaluate the performance of RESCUE, we conducted simulation studies on synthetic data designed to mimic real ST measurements with varying degrees of unattributed expression. We simulated ST data using real scRNA-seq data from the honey bee brain48 by introducing synthetic unattributed expression of a small number of genes (Methods). Unlike other tissues where cell bodies and non-cell-body structures are spatially intermingled, honey bee brains provide a unique benefit for simulating unattributed expression due to the existence of the neuropil region, where axons and dendrites are spatially segregated from cell bodies. Molecular expression in these regions is largely absent in single-cell sequencing data due to the structural fragility of neurites. During the tissue dissociation step, a necessary step in single-cell sequencing, neurites are mostly sheared off and lost. In contrast, ST measurements performed on intact tissue sections preserve these structures and therefore capture molecular expression within neuropils.

To guide the choice of sparsity levels used in our simulations, we examined empirical expression patterns in the neuropil region of our in-house MERFISH data. We found that neuropil regions exhibit highly sparse gene expression, with an average sparsity level of 0.92 across spatial locations (Fig. 2a). Here, sparsity refers to the structure of the unattributed component (i.e., that only a subset of genes is present at a given spatial location), rather than to low expression levels (Methods). Based on this empirical distribution, we simulated unattributed expression across a range of sparsity levels by varying the proportion of zero entries (ranging from 0.8 to 0.95). In addition, we varied the proportion of unattributed expression (ranging from 0 to 0.8) to cover different levels of the overall magnitude of unattributed expression. Importantly, these parameter ranges fall within those observed in real biological samples and are not artificially chosen to favor RESCUE.

We first compared the deconvolution performance of RESCUE with existing methods. Fig. 2b summarizes the accuracy of cell-type proportion estimates across methods, measured using the average Lin’s concordance correlation coefficient49 (CCC), root mean squared error (RMSE), and Jensen-Shannon divergence (JSD) (Methods). As expected, the performance of standard deconvolution methods deteriorates sharply as the level of unattributed expression increases. In contrast, RESCUE maintains high accuracy across all conditions, demonstrating robustness to both high unattributed content and reduced sparsity. These results indicate that failing to explicitly model unattributed expression can lead to systematic biases in cell-type proportion estimates.

We also evaluated the accuracy of unattributed expression recovery (Fig. 2c) using the relative Frobenius error between the true and estimated unattributed expression matrices (Methods). Overall, RESCUE achieves the lowest recovery error even under high unattributed expression and low sparsity levels, whereas other methods deviate substantially from the ground truth, especially when unattributed signals become strong. These results indicate that explicitly modeling unattributed expression enables accurate recovery of hidden transcriptional signals across a wide range of sparsity and signal-strength regimes.

To assess robustness and generalizability, we expanded the simulation study by varying noise levels via library size variability across spots and varying the effective grid resolution via changes in the average total UMI count per spot (Supplementary Fig. 1). As expected, overall performance tends to decrease with increasing library size variability and increase with higher average total UMIs. Across these regimes, RESCUE remains among the top-performing methods, demonstrating robustness across a range of noise levels and grid sizes.

RESCUE accurately recovers the unattributed transcriptome in real tissues
Using a custom MERFISH platform50,51, we generated single molecule-resolution ST data from coronal sections of a honey bee brain, targeting 130 genes relevant to aggression (Supplementary Data 1)52. Transcripts in axons and dendrites, which are collectively termed as neurites, have been difficult to profile due to segmentation challenges. Cell segmentation methods, such as Cellpose6 and Baysor7, can only identify the boundaries of neuronal cell bodies, also called somata. Although individual axons and dendrites can be stained in 2D cell cultures53, in native tissues they extend across multiple imaging planes and therefore cannot be reliably segmented as individual cellular units. As a result, neurite-associated expression has traditionally been ignored or incorrectly merged into nearby cytosolic expression.

We first experimentally validated that RESCUE can accurately recover gene expression in the neuropil region of the honey bee brain (Fig. 3a). This validation is possible because the physical separation of neurites and somata in the brain makes it simple for existing algorithms to separate somatic from unattributed expression. Cell soma were segmented using Cellpose6, revealing that only 73.5% of transcripts were localized within cells, with the remaining 26.5% comprising the unattributed expression. Using honey bee brain scRNA-seq data48 as a reference, we then applied RESCUE and competing methods to our MERFISH data after dividing the brain into grids of size 20×20 𝜇𝑚2 (Methods). For each grid, we estimated the proportion of canonical expression relative to total expression and obtained a ground-truth estimate by identifying transcripts localized within cell bodies via cell segmentation. RESCUE’s estimates closely matched the ground truth (Fig. 3b) and were more accurate than those from other methods, as quantified by Lin’s CCC and RMSE (Methods). Notably, RESCUE relies only on reference scRNA-seq data and does not use cell segmentation information.

To further clarify the source of discrepancies among methods, we evaluated canonical and unattributed expression separately in anatomically well-defined regions. Specifically, we examined grids corresponding to pure cell bodies and compared the average canonical expression of each gene against the ground truth derived from segmentation (Fig. 3c). RESCUE closely matched the true somatic expression levels, whereas competing methods systematically underestimated canonical expression. This behavior reflects the fact that existing deconvolution approaches are designed to estimate cell-type proportions rather than expression levels, and under the subtraction-based framework (Methods), expression that is not well explained by the reference profiles is left as residual rather than attributed to the canonical component. Conversely, in grids containing no identifiable cell bodies, RESCUE accurately recovered unattributed expression in neuropil-only regions (Fig. 3d), while competing methods were significantly less accurate. Together, these results demonstrate that RESCUE effectively partitions the full transcriptome into canonical and unattributed components in real tissue, capturing valuable expression signals that would otherwise be lost using conventional approaches.

Importantly, transcripts outside segmented somata reflect genuine biological signals rather than technical artifacts: in our MERFISH sample (Fig. 3a), neuropil-associated transcripts form structured spatial patterns that are clearly segregated from DAPI-stained nuclei (Supplementary Fig. 2). Moreover, the MERFISH protocol includes gel embedding, expansion, and clearing steps that reduce autofluorescence and suppress nonspecific signal54,55. These signals, which would have been ignored by standard cell segmentation approaches, highlight the importance of RESCUE in capturing overlooked transcriptomic information.

Finally, to investigate whether the unattributed expression recovered by RESCUE reflects structured biological organization rather than unstructured noise, we performed clustering analysis of the unattributed expression profiles across spatial locations. Clustering revealed anatomically distinct transcriptional niches within neuropil regions (Fig. 3e), indicating that the recovered unattributed component exhibits coherent spatial patterns. Specifically, distinct clusters localize to known neuropil structures, including the mushroom body calyces, antennal lobes, optic lobes, and regions proximal to the central complex. In contrast, spatial locations corresponding to somatic regions collapsed into a single cluster, consistent with the fact that canonical expression is explicitly removed, and only minimal unattributed signal remains in these regions. Examination of representative genes within the neuropil-associated clusters showed spatially structured expression patterns that are largely absent from the canonical component (Fig. 3f). These findings further support that the unattributed expression recovered by RESCUE captures biologically meaningful and spatially organized signals, revealing hidden tissue structure that is not accessible through conventional analyses.

RESCUE enables accurate deconvolution of ST data with incomplete reference
Here, we illustrate how RESCUE’s estimated canonical expression component can be used to achieve accurate deconvolution of ST data when key cell types are missing from the reference, a challenge for all existing methods that has been noted in a recent study23. We analyzed paired single-cell FFPE-seq (scFFPE-seq) and 10x Visium ST data from human breast cancer tissue5. In this ST dataset, adipocytes are manually annotated in approximately 10% of spatial spots (Fig. 4a). However, adipocytes are absent from the scFFPE-seq reference, likely due to their high lipid content interfering with fixation and processing. These cells are prone to floating, rupturing, or sticking to plastic surfaces during tissue dissociation, which hinders their recovery in single-cell sequencing protocols5,17. As a result, reference-based methods that assume the observed expression is fully explained by the available reference cell types tend to misattribute adipocyte-derived transcripts to other cell types.

In adipocyte-rich regions, existing deconvolution methods fail to account for the missing adipocyte transcriptome and instead misattribute the signal to other reference cell types, such as stromal and endothelial cells (Fig. 4b). This systematic misassignment likely arises from transcriptional similarity, as shown by the strong correlation between adipocyte and stromal/endothelia profiles in pseudo-bulk measurements from Visium and scFFPE-seq (Fig. 4c). In contrast, RESCUE explicitly models unattributed expression and thus avoids forcing unassigned signal into existing reference cell types, which results in more accurate and biologically plausible estimates. As shown in Fig. 4d-e, RESCUE’s estimated proportions closely align with the scFFPE-seq composition, while other methods exhibit significant distortions. Together, these results demonstrate that RESCUE improves deconvolution accuracy in settings with incomplete references by preventing the systematic misassignment of unattributed transcriptomes to incorrect cell types.

To further interpret the biological content of the unattributed expression recovered by RESCUE, we applied NMF to the unattributed component. This analysis revealed multiple latent factors capturing structured expression patterns that are not explained by the canonical component (Fig. 4f). Notably, the NMF Factor 1 (Unattributed) captures the spatial distribution of adipocyte-rich regions, effectively recovering a missing cell type absent from the reference. Genes with the highest loadings in this factor include well-established adipocyte markers such as FABP4, CFD, and G0S25,56 (Fig. 4g). In addition, the NMF Factor 2 (Unattributed) is enriched for genes such as DCN, COL14A1, and C1R, which have been implicated in extracellular matrix (ECM) domains57–59. Together, these findings indicate that the unattributed component captured by RESCUE is biologically meaningful, thereby supporting that explicitly modeling it leads to improved deconvolution in the presence of incomplete references.

Finally, to assess whether these deconvolution improvements generalize beyond the human breast cancer example, we benchmarked RESCUE using additional datasets, including human dorsolateral prefrontal cortex (DLPFC)60 and cichlid fish telencephalon61 Visium samples. In both datasets, reference profiles were derived from snRNA-seq, which inherently underrepresent cytoplasmic and process-localized RNA, leading standard deconvolution methods to misattribute unexplained ST signals. Across these datasets, RESCUE consistently outperformed competing methods (Supplementary Fig. 3-4), indicating that the benefits of explicitly modeling unattributed expression extend across tissues, species, and reference modalities.

RESCUE uncovers layer-specific unattributed expression in the human DLPFC
In many neuroscience experiments, cell-level information is often characterized using snRNA-seq62,63, which can leave a substantial fraction of biologically relevant expression unattributed. Maynard et al.60 profiled the DLPFC of neurotypical adult donors using paired snRNA-seq and 10x Visium ST. We applied RESCUE to eight samples in which all six cortical layers and the white matter (WM) were annotated, to recover unattributed expression patterns not represented in the snRNA-seq reference (Fig 5a). As shown in Fig. 5b, RESCUE reveals that Layer 1 and the WM contain substantially higher levels of unattributed expression compared to other cortical layers. This is consistent with the fact that Layer 1 has sparser cell densities and contains neuropil-rich regions57, and that WM contains a substantial amount of axons and myelin64, whose transcriptomic content is mostly cytoplasmic and thus is absent from the reference snRNA-seq data.

RESCUE enables us to study biology specific to the unattributed transcriptome of Layer 1. We first identified DEGs by comparing Layer 1 to all other layers using limma65, based on either the raw observed ST expression or the unattributed expression recovered by RESCUE. These two analyses gave substantially different results (Fig. 5c). For example, the top five genes most enriched in Layer 1 based on unattributed expression had either fold-changes in the opposite direction (MT-ND1, MT-ND2, and MT-CO1) or negligible fold-changes (MT3 and CST3) in the differential expression analysis using raw data. We also identified genes such as CLU, GFAP, and MALAT1 that exhibit clear spatial enrichment within the unattributed component of Layer 1 (Fig. 5d). These genes are associated with local translation and neurite outgrowth, and their enrichment outside nucleus is consistent with known roles of astrocytic processes and non-somatic cellular structures in Layer 1 of the cortex66,67.

To explore the underlying biological processes, we performed gene set enrichment analysis separately for DEGs derived from the unattributed expression and from the raw ST data (Fig. 5e). Notably, gene ontology (GO) terms enriched among the unattributed DEGs include mitochondrial processes: oxidative phosphorylation, energy derivation by oxidation of organic compounds, and generation of precursor metabolites and energy. In contrast, no mitochondria-related pathways reached statistical significance in the analysis based on observed ST expression; only the broad term, response to stress, passed the FDR < 0.05 threshold. These results align with the elevated metabolic demands of a neurite-enriched Layer 1, which is densely populated with synapses and reliant on local energy production to support axonal and dendritic activity68,69. Collectively, these findings demonstrate that RESCUE can uncover previously obscured yet functionally important gene expression programs in Layer 1 of the human DLPFC.

RESCUE reveals spatially structured unattributed expression in the cichlid fish telencephalon
Johnson et al.61 collected paired snRNA-seq and Visium ST data from the dorsal telencephalon of individual cichlid fishes exhibiting distinct social behaviors. The dorsal telencephalon is characterized by prominent radial glial (RG) processes that originate in the ventricular zone and extend long projections toward superficial pallial regions70,71. These RG projections span nearly the entire depth of the pallium and may play a role in shaping spatial gene expression profiles across distinct layers and compartments, but they do not survive tissue dissociation and thus their transcriptomes are not captured by snRNA-seq.

We applied RESCUE to this dataset to recover unattributed transcriptomic signals using the snRNA-seq data as a reference. To assess spatial heterogeneity, we performed SVG analysis using SPARK-X72 on both the observed ST expression and the unattributed expression. As shown in Fig. 6a, a substantial number of SVGs are uniquely identified in the unattributed transcriptome. In Fig. 6b, we show the expression maps of those SVGs that were undetectable in the observed ST expression yet exhibited clear spatial structure once nuclear signals were removed. Among these, app, cntn1, dner, and eno1 mark distinct transcriptional programs that may operate along RG projections. These genes highlight transcriptional signals that are likely to reflect localized molecular activities along RG scaffolds: app encodes a transmembrane protein that promotes synapse formation and plasticity73,74, cntn1 (Contactin-1) is a GPI-anchored adhesion molecule essential for axon guidance and neurite outgrowth75, dner mediates neuron-glia interactions that are critical for astrocyte differentiation76, and eno1 (𝛼-enolase) promotes neuritogenesis and contributes to structural remodeling in the developing nervous system77. Indeed, Johnson et al.61 found that bower-building behavior in cichlid fishes may be related to neuronal rebalancing driven by a subpopulation of RG cells and their processes, highlighting the importance of studying transcriptional processes in RG projections, which RESCUE enables.

To further characterize the structure of the unattributed transcriptome, we applied NMF to the unattributed expression component (Fig. 6c). Several unattributed factors are enriched for genes associated with synaptic and process-localized functions, including slc1a2, calm2, and mitochondrial ATP synthase (Fig. 6d), consistent with transcriptional activity in extended RG scaffolds and neuronal processes. In addition, one factor shows strong enrichment for ribosomal protein genes, suggesting localized protein synthesis within these structures78. Notably, these structured transcriptional programs emerge only after the canonical expression is removed, indicating that the unattributed expression recovered by RESCUE captures biologically meaningful signals associated with radial glial and cytoplasmic processes that are not represented in snRNA-seq reference data.

Methods
RESCUE model
We assume that our ST data profiles 𝐺 genes in 𝐽 grids across the tissue section. In low-resolution ST technologies, these grids are predefined as multicellular spots, whereas in higher-resolution ST technologies, these can be constructed at any desired size, which allows inference at different biological scales. To introduce the statistical model behind RESCUE, we begin by reviewing existing linear deconvolution methods12–16. These methods decompose the ST expression matrix 𝐘 ∈ ℝ𝐺×𝐽 of 𝐺 genes in 𝐽 grids as
𝐘 = 𝐗𝛂+𝐄, (1)
where 𝐗 ∈ ℝ𝐺×𝐾 is a matrix of 𝐾 cell-type expression reference profiles, 𝛂 ∈ ℝ𝐾×𝐽 is a matrix of inferred cell-type proportions within each grid, and 𝐄 ∈ ℝ𝐺×𝐽 is assumed to be a mean-zero noise matrix. The issue is that (1) implies that expectation of ST expression is 𝔼(𝐘) = 𝐗𝛂, which ignores unattributed transcriptome.

Instead, we further decompose the error 𝐄, motivated by the robust regression framework45–47. Specifically, we consider the multivariate mean-shift model:
𝐘 = 𝐗𝛂+𝐑+𝓔, (2)
where 𝐗𝛂 represents the expected expression values of each gene in each grid due to the presence of cells whose types are represented in the reference sc/snRNA-seq data, 𝐑 ∈ ℝ𝐺×𝐽 is an unattributed expression that allows for outliers that cannot be captured by the reference data, and 𝓔 ∈ ℝ𝐺×𝐽 is a noise matrix that allows for residual variability between the ST sample and the reference data.

Following the mean-shift model in the robust regression framework, we assume that the unattributed component 𝐑 is sparse. Specifically, the sparsity assumption reflects the expectation that only a subset of genes exhibits unattributed expression within each grid. This assumption is required for identifiability of the decomposition and is biologically plausible, as we expect relatively few genes to exhibit unattributed expression within each grid. Consistent with this assumption, we observe empirically that unattributed expression in neuropil regions of the honey bee brain is highly sparse, with 92% of genes showing no unattributed expression on average. Moreover, even in extracellular environments, ECM proteins typically involve highly sparse, restricted gene sets, with only 1.3–1.5% of genes showing nonzero expression, rather than transcriptome-wide activation81.

Optimization problem
RESCUE aims to simultaneously recover the sparse non-negative matrix 𝐑 and non-negative cell-type proportions 𝛂 in our model (2). Specifically, we consider the following optimization problem:
min‖𝐖⊙𝐑‖1 subject to ‖𝐘−𝐑−𝐗𝛂‖𝐹2 ≤ 𝜆, (3)
𝛂,𝐑≥0
where 𝐖 is a weight matrix, ⊙ denotes the Hadamard product, ‖⋅‖𝐹 denotes the Frobenius norm, and 𝜆 > 0 is a tuning parameter. We set each entry of the weight matrix as 𝑊𝑔𝑗 = 1/√Var(𝑅𝑔𝑗) to account for the fact that different genes have different variances across different locations. These weights will be estimated through an iterative variance estimation procedure. The ℓ1 objective function enforces sparsity of our estimated 𝐑, while the constraint ensures that the residual variability 𝓔 is not too large, in other words, that the ST data do not deviate too much from the reference sc/snRNA-seq data.

In practice, we solve the Lagrangian form of this problem:
min 𝑓(𝛂,𝐑) = ‖𝐖⊙𝐑‖1 + 𝜇 ‖𝐘−𝐑−𝐗𝛂‖2𝐹, (4)
𝛂,𝐑≥0 2
where 𝜇 > 0 is a tuning parameter. We apply an efficient alternating minimization to solve two convex sub-problems in each iteration:
𝛂̂ = argmin‖𝐘−𝐑−𝐗𝛂‖2𝐹, (5)
𝛂≥0
𝐑̂ = argmin‖𝐖⊙𝐑‖1 + 𝜇 ‖𝐘−𝐑−𝐗𝛂‖2𝐹 = 𝑆+ (𝐘−𝐗𝛂, 𝐖/𝜇), (6)
𝐑≥0 2
where 𝑆+ (𝑥,𝜏) = max (𝑥 −𝜏,0) denotes the entry-wise non-negative soft-thresholding function.

Estimation of adaptive weights
Kong et al.47 proposed adaptive lasso-type weights to reduce masking and swamping effects in the presence of multiple outliers. Our weighting scheme plays a similar role but differs in that the weights are learned iteratively through variance estimation of the unattributed component. In ST deconvolution models, the observed count 𝑌̃𝑔𝑗 is often modeled under a Poisson assumption13,16:
𝑌̃𝑔𝑗 |𝜆𝑔𝑗 ∼ Poisson(𝑁𝑗𝜆𝑔𝑗) (7)
where 𝑁𝑗 = ∑𝑔 𝑌̃𝑔𝑗 represents the total transcript count, or number of UMIs for grid 𝑗, and 𝜆𝑔𝑗 denotes the expected expression rate of gene 𝑔 at grid 𝑗. Rather than directly modeling 𝑌̃𝑔𝑗, we assume that it can be decomposed as 𝑌̃𝑔𝑗 = 𝐿̃𝑔𝑗 +𝑅̃𝑔𝑗, where
𝐿̃𝑔𝑗 |𝜆𝐋𝑔𝑗 ∼ Poisson(𝑁𝐋𝑗𝜆𝐋𝑔𝑗) and 𝑅̃𝑔𝑗 |𝜆𝐑𝑔𝑗 ∼ Poisson(𝑁𝐑𝑗𝜆𝐑𝑔𝑗). (8)
Here, 𝑁𝐋𝑗 = ∑𝑔 𝐿̃𝑔𝑗 and 𝑁𝐑𝑗 = ∑𝑔 𝑅̃𝑔𝑗 represent the total transcript counts of the canonical and unattributed expression for grid 𝑗, respectively, and (𝜆𝐋𝑔𝑗,𝜆𝐑𝑔𝑗) denote the expected expression rates of these components. We can write 𝑁𝐋𝑗 = (1−𝜋𝑗)𝑁𝑗 and 𝑁𝐑𝑗 = 𝜋𝑗𝑁𝑗, where 𝜋𝑗 ∈ [0,1] denotes the proportion of unattributed expression for grid 𝑗.

When using counts per million (CPM) normalization, i.e., 𝑌𝑔𝑗 = 𝑌̃𝑔𝑗 /𝑁𝑗 ×106, we can estimate the variance of unattributed component as
Var(𝑅𝑔𝑗) = Var( 𝑅̃𝑔𝑗 /𝑁𝑗 ×106) = (106/𝑁𝑗)2 × Var(𝑅̃𝑔𝑗) ≈ 106/𝑁𝑗 ×𝑅𝑔𝑗, (9)
where the approximation follows from V̂ar(𝑅̃𝑔𝑗) = 𝑁𝐑𝑗𝜆𝐑𝑔𝑗 and 𝑅𝑔𝑗 = 𝑅̃𝑔𝑗 /106 ×𝑁𝑗. Here, the term 𝑁𝐑𝑗𝜆𝐑𝑔𝑗 can be iteratively updated using the plug-in estimator 𝑅̃𝑔𝑗 ≈ 𝑁𝐑𝑗𝜆𝐑𝑔𝑗. To avoid division by zero, we set 𝑊𝑔𝑗 = 1/max(√106/𝑁𝑗 ×𝑅𝑔𝑗,𝛿), where 𝛿 is a small positive constant. We summarize the RESCUE algorithm in Supplementary Algorithm 1.

Tuning parameter selection
RESCUE aims to recover unattributed expression by balancing the sparsity of 𝐑, which reflects the number of rejected outliers, against the level of residual variability 𝓔. Following standard practices in robust regression and robust PCA45,79, we set 𝜇 = 𝐺𝐽/𝑐‖𝐘‖1 and focus on selecting the optimal value of the tuning parameter 𝑐. However, choosing the optimal tuning parameter is challenging in the presence of outliers80,82.

To select the optimal tuning parameter, we adopt a data-driven strategy recently proposed for related deconvolution problems in bulk RNA-seq anaysis83. The key idea is to generate simulated data that closely resemble the observed ST data and use this for evaluation. Specifically, we use a surrogate ST dataset with ground truth proportions of canonical and unattributed expression. To construct this surrogate ST dataset, we first apply the non-negative least squares (NNLS) to the observed ST data to obtain initial estimates for 𝛂 and 𝐑:
𝐑̂NNLS = (𝐘−𝐗𝛂̂NNLS)+, where 𝛂̂NNLS = argmin𝛂≥0‖𝐘−𝐗𝛂‖2𝐹, (10)
where (𝑥)+ = max(𝑥,0) denotes the zero-thresholding function. We then simulate a surrogate ST data 𝐘̃sim from a Poisson distribution with rescaled expected expression rates, preserving the overall signal-to-noise characteristics observed in the real data. RESCUE is applied to 𝐘̃sim across a range of 𝜇 values, and the optimal 𝑐 is selected as the value that minimizes the recovery error of the canonical and unattributed components in the surrogate dataset. We summarize the tuning procedure in Supplementary Algorithm 2.

Platform effect normalization
To mitigate cross-platform effects between ST data and sc/snRNA-seq reference data, we apply platform effect normalization. This step is essential because differences in experimental platforms can introduce systematic shifts in gene expression levels that confound downstream analysis. Without proper alignment of scale, these differences may lead to biased estimation of cell-type contributions or unattributed signals. We adopt the normalization strategy from RCTD13, which rescales each gene in the reference data to match the overall expression level of the corresponding gene in the ST data. Specifically, we first compute pseudo-bulk measurements for both ST and reference datasets and then apply a gene-wise log-scale correction to harmonize them.

Data normalization and reference matrix
To account for technical variability and the effect of library size, we adopt CPM normalization but fit RESCUE in the linear space. We avoid log-transformation for two main reasons. First, log-transformed input tends to underestimate true signals in gene expression deconvolution84. Second, log-transformation requires adding pseudo-counts, which can introduce significant artificial distortion between zeros and non-zeros85,86. Instead, we work in the linear scale, employing an iteratively reweighted algorithm (Supplementary Algorithm 1) to stabilize variance adaptively. This approach allows the model to adjust for varying data scales and minimizes the influence of outliers without artificially altering data distributions. The normalized reference matrix is then used to construct average cell-type expression.

Reference dependence of unattributed expression
RESCUE is a computational negative selection method which relies on the choice of reference, and interpretation of the recovered unattributed expression therefore requires consideration of both biological context and reference completeness. To further assess the impact of reference selection within the same tissue, we extended the honey bee brain analysis using snRNA-seq87 and segmented MERFISH cells as alternative reference profiles (Supplementary Fig. 5). Using segmented MERFISH cells produced unattributed expression patterns highly consistent with those obtained using scRNA-seq, as both references capture cellular transcripts beyond the nucleus and therefore explain similar fractions of the ST signal. In contrast, the snRNA-seq reference yielded higher levels of unattributed expression on average, reflecting its enrichment for nuclear RNA and underrepresentation of cytoplasmic and neuritic transcripts. These results indicate that variation in the recovered unattributed component reflects biologically meaningful differences in reference coverage rather than technical artifacts, and that interpretation of unattributed expression should be guided by the biological properties of both the tissue and the reference dataset.

Choice of grid resolution for honey bee brain MERFISH data
For high-resolution ST platforms such as MERFISH, RESCUE operates on spatially aggregated expression profiles defined over user-specified grids. The choice of grid size reflects a trade-off between spatial resolution and sparsity: smaller grids preserve fine spatial structure but increase sparsity, whereas larger grids reduce sparsity at the cost of spatial detail. As a practical guideline, we recommend selecting a grid size such that the average total transcript count per grid is comparable to the typical total transcript count per cell. In the honey bee MERFISH analysis, we used a grid size of 20×20 𝜇𝑚2, which approximates cell-scale resolution. To assess sensitivity to grid resolution, we varied grid size across a range of spatial scales and found that RESCUE produces consistent results across a broad range of grid resolutions (Supplementary Fig. 6).

Subtraction-based estimation of unattributed expression for existing methods
Most existing ST deconvolution methods are designed to estimate cell-type proportions, rather than to recover expression levels directly. Because these methods do not explicitly model unattributed expression, we adopted a subtraction-based framework, following a strategy used in related deconvolution settings for bulk RNA-seq data23. Specifically, we first estimated cell-type proportions using each deconvolution method and then recovered the unattributed expression by subtracting the estimated cellular contributions from the observed expression, followed by zero-thresholding to avoid negative values. Under this framework, the unattributed expression is defined as the portion of the observed ST expression that cannot be explained by the reference-based cellular expression implied by each method’s estimated cell-type proportions.

Data generating process for simulation studies
For a total of 𝐽 = 2,500 grids, we first generated the total UMI count for grid 𝑗 from a negative binomial distribution 𝑁𝑗 ∼ NB(mean = 𝑀,size = 1/𝜎). We generated the raw count for each gene in each grid as 𝑌̃𝑔𝑗 = 𝐿̃𝑔𝑗 +𝑅̃𝑔𝑗, where 𝑅̃𝑔𝑗 denotes the unattributed expression and 𝐿̃𝑔𝑗 represents the canonical expression that could have been captured by sc/snRNA-seq reference data.

To generate 𝐿̃𝑔𝑗, we used the real scRNA-seq data from honey bee brains48 as the source of canonical expression profiles. We selected the top 𝐺 = 500 number of highly variable genes and constructed the reference matrix 𝐗 using 𝐾 = 20 number of cell types annotated via Louvain clustering in Seurat88. Latent mixing proportions were drawn as 𝛽𝑘𝑗 ∼ Beta(1,3) and normalized as 𝛼𝑘𝑗 = 𝛽𝑘𝑗 /∑𝑘 𝛽𝑘𝑗, with at most 10 cell types per spot. Defining 𝜆𝐋𝑔𝑗 = (𝐗𝛂)𝑔𝑗, we generated the cell-type-driven canonical expression as
𝐿̃𝑔𝑗 ∼ Multinomial((1−𝜋𝑗)𝑁𝑗,𝑃𝑔𝑗 = 𝜆𝐋𝑔𝑗 /∑𝑔 𝜆𝐋𝑔𝑗). (11)
To generate sparse unattributed expression, we sampled
𝑅̃𝑔𝑗 ∼ Multinomial(𝜋𝑗𝑁𝑗,𝑄𝑔𝑗 = 𝜆𝐑𝑔𝑗 /∑𝑔 𝜆𝐑𝑔𝑗), (12)
where 𝜆𝐑𝑔𝑗 ∼ Gamma(1,1/3) with probability (1−𝑠) and 𝜆𝐑𝑔𝑗 = 0 with probability 𝑠. Here, we introduced a sparsity parameter 𝑠, which makes (𝑠×100) percentage of intensity parameter 𝜆𝐑𝑔𝑗 zero-inflated. Under this setting, we can generate 𝐿̃𝑔𝑗 and 𝑅̃𝑔𝑗 that satisfy ∑𝑔 (𝐿̃𝑔𝑗 +𝑅̃𝑔𝑗) = 𝑁𝑗. In our simulations, we fixed 𝑀 = 10,000 and 𝜎 = 0.5, while varying 𝜋𝑗 ∈ {0,0.2,0.4,0.6,0.8} and 𝑠 ∈ {0.8,0.85,0.9,0.95} to control the proportion of canonical expression and the sparsity level of unattributed transcriptome, respectively (Fig. 2). We then examined the robustness of RESCUE under additional simulation settings by varying the average total UMI count per spot (𝑀), as well as varying noise levels by adjusting library size variability across spots (𝜎) (Supplementary Fig. 1).

Evaluation metrics for cell-type deconvolution in simulation studies
To evaluate the accuracy of cell-type deconvolution, we considered three complementary metrics: Lin’s CCC, RMSE, and JSD. We first measured Lin’s CCC between the true and estimated proportion of the 𝑘th cell type, defined as
CCC(𝜶𝑘,𝛂̂𝑘) = 2Cov(𝛂𝑘,𝛂̂𝑘) / (𝜎2𝛂𝑘 +𝜎2𝛂̂𝑘 +(𝜇𝛂𝑘 −𝜇𝛂̂𝑘)2), (13)
where 𝜇𝛂𝑘 and 𝜇𝛂̂𝑘 are the means of the two variables, and 𝜎2𝛂𝑘 and 𝜎2𝛂̂𝑘 denote the corresponding variances. We used Lin’s CCC because it captures shifts in both location and scale, making it robust to high leverage points83. The concordance increases as the CCC value gets closer to 1. We report the average CCC across all cell types.

To quantify estimation error at the spatial location level, we computed the RMSE between the true and estimated cell-type proportion vectors. For spatial location 𝑗, RMSE is defined as
RMSE𝑗 = √ 1/𝐾 ∑𝑘=1𝐾 (𝛼𝑘𝑗 −𝛼̂𝑘𝑗)2, (14)
where 𝐾 is the number of cell types, and 𝛼𝑘𝑗 and 𝛼̂𝑘𝑗 denote the true and estimated proportions of cell type 𝑘 at location 𝑗, respectively. We report the average RMSE across all spatial locations, with smaller values indicating better performance.

Finally, we computed the JSD between the true and estimated cell-type proportion distributions at each spatial location. For spatial location 𝑗, JSD is defined as
JSD𝑗 = 1/2 KL(𝛂𝑗 ‖𝐦𝑗)+ 1/2 KL(𝛂̂𝑗 ‖𝐦𝑗), 𝐦𝑗 = 1/2 (𝛂𝑗 +𝛂̂𝑗), (15)
where KL(⋅‖⋅) denotes the Kullback-Leibler divergence. We report the average JSD across all spatial locations, with smaller values indicating better performance.

Evaluation metrics for unattributed expression recovery in simulation studies
To evaluate the accuracy of unattributed expression recovery, we measured the relative Frobenius error, defined as
(Relative Frobenius Error) = ‖𝐑−𝐑̂‖𝐹 / ‖𝐑‖𝐹, (16)
where 𝐑 and 𝐑̂ denote the true and estimated unattributed expression matrices, respectively.

Evaluation metrics for honey bee brain MERFISH data
We measured Lin’s CCC and the RMSE between the segmentation-based (true) and estimated unattributed expression proportions, defined as follows:
CCC(𝛒,𝛒̂) = 2Cov(𝛒,𝛒̂) / (𝜎2𝛒 +𝜎2𝛒̂ +(𝜇𝛒 −𝜇𝛒̂)2) and RMSE(𝛒,𝛒̂) = √ 1/𝐽 ∑𝑗=1𝐽 (𝜌𝑗 −𝜌̂𝑗)2, (17)
where 𝜌𝑗 = ∑𝑔 𝑅̃𝑔𝑗 /∑𝑔 𝑌̃𝑔𝑗 denotes the true unattributed expression proportion in grid 𝑗, and 𝜌̂𝑗 denotes its estimated value. Here, 𝜇𝛒 and 𝜇𝛒̂ are the means of two variables, and 𝜎2𝛒 and 𝜎2𝛒̂ denote the corresponding variances.

Benchmarking computational runtime
We evaluated the computational runtime of each method across ST datasets of varying sizes (Supplementary Fig. 7). Across all benchmarks, RESCUE completed within practical time limits for contemporary ST datasets, demonstrating scalability to datasets with hundreds to thousands of spatial locations and genes. All analyses were performed using CPU-only computation on a standard laptop (Apple M3 Pro with 36GB of memory), without GPUs or high-performance computing clusters.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Identifying the sources of unassigned or missing expression in spatial transcriptomics (ST) remains a major analytical challenge. In this manuscript, Lee et al. introduce RESCUE, a computational method designed to decompose ST data into canonical (reference-aligned) and idiosyncratic (reference-free) components. The key innovation lies in its robust sparse recovery framework, which models unexplained ST expression as a sparse residual term while maintaining biological interpretability through non-negativity constraints. The authors demonstrate that RESCUE successfully captures expression patterns missed by traditional segmentation and deconvolution approaches, such as neurite-enriched or extracellular transcripts, and improves deconvolution accuracy when reference datasets are incomplete. Applications to honey bee brain MERFISH, human breast cancer, human DLPFC, and cichlid fish telencephalon data compellingly illustrate the method’s potential to uncover biologically meaningful “dark transcriptomes.” Overall, the concept of computational negative selection is original and timely. However, several important concerns must be addressed before the manuscript is suitable for publication.\nMajor comments:\n1. While the manuscript convincingly shows that RESCUE can partition canonical and idiosyncratic components, the validation of the latter is largely qualitative. For instance, the honey bee neuropil example estimates that 26.5% of total transcripts are idiosyncratic, but it remains unclear whether these transcripts represent genuine biological signals or technical artifacts (e.g., diffusion, autofluorescence, or mis-segmentation). The authors should benchmark RESCUE’s recovered signals against independent datasets or orthogonal measurements (such as smFISH or subcellular RNA localization assays) to verify that the idiosyncratic component corresponds to meaningful biology.\n2. Since the canonical component in RESCUE is reference-dependent, the choice of single-cell or single-nucleus reference can influence results. The authors should systematically evaluate how the inferred idiosyncratic expression changes when using different references from the same tissue, or from external atlases. This would strengthen the argument that the recovered idiosyncratic signals are not artifacts of a specific reference dataset.\n3. The manuscript primarily focuses on differential gene expression analysis within the idiosyncratic component, which effectively highlights region-specific or process-specific genes. However, clustering the idiosyncratic expression profiles across spatial locations or genes could provide deeper insights into hidden tissue organization or novel substructures. For instance, clustering may reveal distinct transcriptional niches associated with neuritic regions, extracellular domains, or unannotated cell states. Incorporating such an analysis—or at least discussing its potential—would strengthen the biological interpretation and demonstrate that RESCUE can uncover structured, not merely residual, expression heterogeneity.\nMinors:\n1. Given that the central innovation of RESCUE lies in decomposing spatial transcriptomic signals into canonical and idiosyncratic components, a dedicated visualization that directly contrasts these two layers would greatly improve conceptual clarity in each example.\n2. In the simulation experiments, the idiosyncratic component of other deconvolution methods is defined as the residual obtained by subtracting the estimated cellular contributions from the observed ST expression. However, this indirect definition may not be conceptually fair or methodologically consistent, since most benchmarked methods are not designed to explicitly separate residual signals. Such residuals might contain confounding effects from normalization, model mismatch, or scaling differences rather than true idiosyncratic expression. The authors should clarify the rationale behind this subtraction-based estimation and, if possible, test whether the comparative results hold under alternative definitions or normalization schemes.\n3. The manuscript mentions that grids are user-defined but does not provide guidelines for choosing the grid size. Since grid size affects both sparsity and spatial resolution, it would be useful to provide a heuristic or sensitivity analysis showing how RESCUE results vary with grid resolution.", "reviewer_2": "This manuscript introduces RESCUE, a computational framework to recover “idiosyncratic” expression in spatial transcriptomics (ST) data—signals missed by segmentation or deconvolution approaches. The authors apply RESCUE to simulations and several real datasets and demonstrate improved deconvolution performance and novel biological discovery. The idea has some novelty, however, some methodological and biological issues require clarification and further support before publication. Below I provide a list of specific comments.\n- The concept of “idiosyncratic expression” is ambiguous. According to the description, It includes missed cell types, missed substructural expression, missed extracellular expression. However, the method only uses one component to capture all these missed components. The manuscript should better delineate how these categories differ, and provide guidance on how users should interpret RESCUE output in practical analyses. Moreover, the term “idiosyncratic transcriptome” is unusual and might be off-putting. Consider alternative terminology such as “non-canonical” or “unattributed expression”.\n- The simulations are based on honey bee scRNA-seq data with synthetic sparsity. It would be important to test additional scenarios (e.g., varying grid sizes, different tissue types, higher noise levels) to show generalizability. The authors should also clarify whether sparsity assumptions used in simulations may bias in favor of RESCUE compared to other methods.\n- In the honey bee data results: Intuitively, one would expect competing methods to attribute more expression to canonical (somatic) transcripts. However, Figures 3A and 3C show the opposite pattern, other methods consistently estimate lower somatic expression. This counterintuitive result should be explained.\n- The comparison is limited to NNLS, WNNLS, RCTD, and CARD. More recent methods (e.g., SONAR, Redeconve) are not evaluated. Since the field is rapidly evolving, omitting them weakens the benchmarking.\n- The manuscript relies mainly on average Lin’s concordance correlation coefficient (CCC) to assess cell-type deconvolution performance. More widely used metrics such as root mean squared error (RMSE) and Jensen–Shannon divergence (JSD) should also be reported to provide a fuller evaluation of accuracy and robustness.\n- One important application of RESCUE is to improve the cell type deconvolution. The capability is demonstrated in one human breast cancer dataset. This result is too thin. The authors need to add more results to show that this improvement is robust across different scenarios.\n- The assumption that the idiosyncratic component is sparse may not hold in all tissues (e.g., in extracellular RNA-rich tumor microenvironments). The authors should discuss scenarios where RESCUE might misclassify biologically meaningful widespread signals as noise, and how users can detect such misbehavior.\n- The paper presents striking biological discoveries (e.g., mitochondrial gene enrichment in Layer 1 of DLPFC, RG-related programs in fish). However, these remain computational predictions. At minimum, the authors should (i) validate with orthogonal datasets (e.g., in situ hybridization or existing histology), or (ii) discuss potential pitfalls of relying solely on computational inference without experimental follow-up.\n- The Methods section is difficult to follow. Clearer explanation of the optimization problem, identifiability conditions, and tuning parameter selection is necessary. Readers less familiar with robust regression might struggle to implement or extend the model. Equation 5 has two plus signs (+), might be a typo there.\n- The manuscript claims broad applicability (MERFISH, Visium, Slide-seq, etc.), but computational runtime, scalability, and memory requirements are not systematically reported. Given the size of current ST datasets, performance benchmarks (runtime vs. data size) should be included.", "reviewer_3": "Summary\nThis manuscript presents RESCUE, a new method to recover expression from genes in missing or underrepresented cell types following cell type deconvolution or cell segmentation in spatial transcriptomics data. The authors refer to this type of expression as \"idiosyncratic\" expression. This is a useful method, which addresses a clear limitation in current methods.\nThe authors demonstrate performance on simulated and real datasets, and compare against several existing methods. The method is available as an R package from GitHub, including a tutorial, with data available from Figshare. The manuscript is clearly written and free of errors. The main comment below relates to providing additional information on computational requirements and scalability.\nMajor comments\n• No information is provided on computational requirements and scalability. This is important information for users, and should be added in additional text and/or figures - e.g. computational resources required (e.g. standard laptop or high-performance compute cluster / GPUs), runtime for the evaluations discussed in the manuscript, and scalability with the number of cells / grids and/or genes.\nMinor comments\n• The method is available as an R package from GitHub, which can be installed easily, and includes a tutorial. However, the tutorial requires manually saving and loading a data file in .RData format, which may be not be clear for some users. If this data file is included in the R package, it could be made easier to load by storing in 'inst/extdata/' and loading with 'system.file()'. Alternatively, the start of the tutorial could be updated with some additional text to explain how to load the data."}, "author_response": "RESPONSE TO REVIEWERS’ COMMENTS\nReviewer #1 (Remarks to the Author):\nIdentifying the sources of unassigned or missing expression in spatial transcriptomics (ST)\nremains a major analytical challenge. In this manuscript, Lee et al. introduce RESCUE, a\ncomputational method designed to decompose ST data into canonical (reference-aligned)\nand idiosyncratic (reference-free) components. The key innovation lies in its robust sparse\nrecovery framework, which models unexplained ST expression as a sparse residual term\nwhile maintaining biological interpretability through non-negativity constraints. The authors\ndemonstrate that RESCUE successfully captures expression patterns missed by traditional\nsegmentation and deconvolution approaches, such as neurite-enriched or extracellular\ntranscripts, and improves deconvolution accuracy when reference datasets are\nincomplete. Applications to honey bee brain MERFISH, human breast cancer, human\nDLPFC, and cichlid fish telencephalon data compellingly illustrate the method’s potential\nto uncover biologically meaningful “dark transcriptomes.” Overall, the concept of\ncomputational negative selection is original and timely. However, several important\nconcerns must be addressed before the manuscript is suitable for publication.\nResponse: Thank you very much for your positive comments. We have revised the\nmanuscript accordingly and believe that these changes have substantially improved the\npaper. Below, we outline the specific changes made in response to your comments.\nTerminology note. Throughout the revision, we have replaced the term “idiosyncratic”\nexpression (or component) with “unattributed” expression (or component) to better reflect\nthe biological interpretation of signals not explained by the reference.\nMajor comments:\n1. While the manuscript convincingly shows that RESCUE can partition canonical and\nidiosyncratic components, the validation of the latter is largely qualitative. For instance, the\nhoney bee neuropil example estimates that 26.5% of total transcripts are idiosyncratic, but\nit remains unclear whether these transcripts represent genuine biological signals or\ntechnical artifacts (e.g., diffusion, autofluorescence, or mis-segmentation). The authors\nshould benchmark RESCUE’s recovered signals against independent datasets or\northogonal measurements (such as smFISH or subcellular RNA localization assays) to\nverify that the idiosyncratic component corresponds to meaningful biology.\nResponse: We agree that validation of the unattributed (previously termed “idiosyncratic”)\ncomponent is crucial. We have updated the manuscript to clearly show that the honey bee\nexample provides the quantitative validation and the identified unattributed expression is\nnot likely to be technical artifacts. We also provide additional analyses to confirm the\nbiological relevance of unattributed expression.\n1. One advantage of the honey bee brain is that it allows for quantitative biological\nvalidation of the RESCUE results. In the honey bee brain, neurites and soma are\nphysically separated, and thus the unattributed expression identified by RESCUE\nwithin neuropil regions directly matches with the measured expression from\nMERFISH. To make this point clear, we added an example MERFISH sample (Figure\n3a) and a representative DAPI image (Supplementary Figure 2) with further details\n(lines 218-221).\n2. We further validated biological significance of transcripts in neuropil regions by\nclustering the unattributed expression profiles, which revealed coherent and\nbiologically meaningful patterns (lines 226-240). Please see our response to your\nMajor Comment 3 for details.\n3. We further clarified why the recovered signals are unlikely to arise from technical\nartifacts (lines 221-225). In the MERFISH experiment, sample preparation includes\ngel embedding, expansion, and clearing. Anchoring mRNAs to a gel matrix has been\nshown to reduce background without detectable signal loss1. The clearing step\nremoves proteins and lipids, resulting in minimal autofluorescence2. Finally,\nsegmentation-related errors are unlikely in the honey bee brain due to the clear\nphysical separation between neurites and soma, which minimizes ambiguity at cell\nboundaries (lines 190-192).\n2. Since the canonical component in RESCUE is reference-dependent, the choice of single-\ncell or single-nucleus reference can influence results. The authors should systematically\nevaluate how the inferred idiosyncratic expression changes when using different references\nfrom the same tissue, or from external atlases. This would strengthen the argument that the\nrecovered idiosyncratic signals are not artifacts of a specific reference dataset.\nResponse: As the reviewer noted, the reference-dependent interpretation of the\nunattributed expression is a very important point. Indeed, the four examples presented in\nthe manuscript are intended to illustrate how unattributed expression can be interpreted\ndifferently depending on biological context and reference choice. We further clarified in the\nrevised manuscript that RESCUE is reference-dependent, and that the interpretation of the\nunattributed expression does require careful consideration of the chosen reference dataset\n(lines 114-126).\nTo systematically evaluate how reference selection influences the recovered signals within\na given tissue, we extended the honey bee analysis to include two additional references:\nsnRNA-seq3 and segmented MERFISH cells, alongside the original scRNA-seq reference.\nNotably, using segmented MERFSH cells as a reference yields unattributed expression\npatterns highly consistent with those obtained using scRNA-seq, whereas the snRNA-seq\nreference results in higher unattributed expression levels on average (Supplementary\nFigure 5). Thus, RESCUE’s results correctly reflect the fact that scRNA-seq and MERFISH-\nsegmented references capture cytoplasmic transcripts, allowing them to explain similar\nfractions of the ST signal, while snRNA-seq primarily profiles nuclear RNA, leading to\ncytoplasmic transcripts being underrepresented in the reference and thus remaining in the\nunattributed component. We have now added the above analysis to the manuscript (lines\n546-561).\n3. The manuscript primarily focuses on differential gene expression analysis within the\nidiosyncratic component, which effectively highlights region-specific or process-specific\ngenes. However, clustering the idiosyncratic expression profiles across spatial locations or\ngenes could provide deeper insights into hidden tissue organization or novel substructures.\nFor instance, clustering may reveal distinct transcriptional niches associated with neuritic\nregions, extracellular domains, or unannotated cell states. Incorporating such an\nanalysis—or at least discussing its potential—would strengthen the biological\ninterpretation and demonstrate that RESCUE can uncover structured, not merely residual,\nexpression heterogeneity.\nResponse: Thanks for the great suggestion. We agree that clustering and related\nunsupervised analyses can provide deeper insight into whether the unattributed\nexpression recovered by RESCUE reflects structured biological organization rather than\nunstructured residual noise. To address this, we performed two additional analyses that\ncan explore the structure of expression heterogeneity identified by RESCUE:\n1. Clustering analysis: Clustering the unattributed expression profiles across spatial\nlocation in honey bee brain reveals anatomically distinct transcriptional niches\nwithin neuropil regions (Figures 3e-f). Specifically, distinct clusters localize to\nknown neuropil structures, including the mushroom body calyces, antennal/optic\nlobes, and regions proximal to the central complex, indicating that the recovered\nsignals are spatially organized and biologically meaningful (lines 226-240).\n2. Non-negative matrix factorization (NMF): Applying NMF to unattributed component\nfurther uncovers latent biological structure across multiple systems. In the human\nbreast cancer dataset, NMF identifies a missing cell type (adipocyte) as well as\nextracellular domains (Figures 4f-g; lines 267-278). In the cichlid fish dataset, NMF\nreveals novel substructures related to synaptic functions and neuronal processes\n(Figures 6d-e; lines 356-365).\nMinors:\n1. Given that the central innovation of RESCUE lies in decomposing spatial transcriptomic\nsignals into canonical and idiosyncratic components, a dedicated visualization that directly\ncontrasts these two layers would greatly improve conceptual clarity in each example.\nResponse: We agree with the reviewer and include dedicated visualizations in each\nexample that directly contrast the canonical and unattributed expression, improving\nconceptual clarity of the RESCUE’s decomposition (Figures 3f, 4g, 5d, and 6c).\n2. In the simulation experiments, the idiosyncratic component of other deconvolution\nmethods is defined as the residual obtained by subtracting the estimated cellular\ncontributions from the observed ST expression. However, this indirect definition may not be\nconceptually fair or methodologically consistent, since most benchmarked methods are\nnot designed to explicitly separate residual signals. Such residuals might contain\nconfounding effects from normalization, model mismatch, or scaling differences rather\nthan true idiosyncratic expression. The authors should clarify the rationale behind this\nsubtraction-based estimation and, if possible, test whether the comparative results hold\nunder alternative definitions or normalization schemes.\nResponse: Thank you for raising this important conceptual point. We agree that most\nexisting deconvolution methods are not designed to explicitly recover residual or\nunattribued expression, and this limitation is precisely the reason why RESCUE fills an\nimportant gap in the literature.\nImportantly, we note that the primary goal of most existing deconvolution methods is to\nestimate “cell-type proportions”, not “expression levels”. Because they do not provide an\nexplicit residual component, we compared RESCUE to subtraction-based extensions of\nthese methods, in which the estimated cellular contribution is subtracted from the\nobserved signal. While indirect, this approach provides the most straightforward and\nconceptually transparent baseline for comparison, and it reflects common practice in\nrelated settings. Notably, a similar subtraction-based strategy has been used to estimate\nresidual signals in bulk deconvolution4. We now clarify this rationale in the main text and\ninclude the corresponding citation (lines 575-585).\nWe agree that alternative extensions of existing methods may be possible and may lead to\nnew methods of residual recovery. Exploring such alternatives is beyond the scope of the\ncurrent study, but we view this as an important direction for future research.\n3. The manuscript mentions that grids are user-defined but does not provide guidelines for\nchoosing the grid size. Since grid size affects both sparsity and spatial resolution, it would\nbe useful to provide a heuristic or sensitivity analysis showing how RESCUE results vary\nwith grid resolution.\nResponse: Thank you for this important suggestion. We now include a sensitivity analysis\nin which we vary the grid size and demonstrate that RESCUE results are robust across a\nrange of grid resolutions (Supplementary Figure 6; lines 571-573). In addition, we provide\na practical guideline for selecting grid size based on the average spot counts per cell,\nbalancing the sparsity level and spatial resolution (lines 563-571).\nReviewer #2 (Remarks to the Author):\nThis manuscript introduces RESCUE, a computational framework to recover “idiosyncratic”\nexpression in spatial transcriptomics (ST) data—signals missed by segmentation or\ndeconvolution approaches. The authors apply RESCUE to simulations and several real\ndatasets and demonstrate improved deconvolution performance and novel biological\ndiscovery. The idea has some novelty, however, some methodological and biological issues\nrequire clarification and further support before publication. Below I provide a list of specific\ncomments.\nResponse: Thank you very much for your constructive comments, which have substantially\nimproved the manuscript. We have carefully revised the manuscript in response to your\ncomments, and the specific changes are detailed below.\nTerminology note. Throughout the revision, we have replaced the term “idiosyncratic”\nexpression (or component) with “unattributed” expression (or component) to better reflect\nthe biological interpretation of signals not explained by the reference.\n- The concept of “idiosyncratic expression” is ambiguous. According to the description, It\nincludes missed cell types, missed substructural expression, missed extracellular\nexpression. However, the method only uses one component to capture all these missed\ncomponents. The manuscript should better delineate how these categories differ, and\nprovide guidance on how users should interpret RESCUE output in practical analyses.\nMoreover, the term “idiosyncratic transcriptome” is unusual and might be off-putting.\nConsider alternative terminology such as “non-canonical” or “unattributed expression”.\nResponse: Thank you for pointing out this important point. First, we have replaced the term\n“idiosyncratic expression” with “unattributed expression” and further clarified in the\nrevised manuscript that interpretation of the unattributed component requires careful\nconsideration of biological structures (lines 114-126).\nWe completely agree that unattributed expression may comprise multiple sources,\nincluding missed cell types, substructural or process-level expression, and extracellular\nexpression. RESCUE is designed to recover these signals jointly. These can be further\ndecomposed using standard unsupervised approaches. To illustrate this, we have now\napplied non-negative matrix factorization (NMF) to RESCUE’s results, revealing a missing\ncell type (adipocytes) as well as extracellular domains in the human breast cancer dataset\n(Figures 4f-g; lines 267-278), and novel substructures specific to neuronal processes in\nthe cichlid fish dataset (Figures 6d-e; lines 356-365). Similarly, we have now clustered\nRESCUE’s results in the honey bee brain, uncovering distinct transcriptional niches within\nneuropil regions (Figures 3e-f), demonstrating that the recovered signals are spatially\norganized and biologically meaningful (lines 226-240). Together, these examples illustrate\nhow unattributed expression can be further partitioned in different biological contexts and\nshow that proper interpretation of the unattributed expression requires biological insights\ninto structural and expression components that are missed by the reference data in\ncorresponding tissues.\n- The simulations are based on honey bee scRNA-seq data with synthetic sparsity. It would\nbe important to test additional scenarios (e.g., varying grid sizes, different tissue types,\nhigher noise levels) to show generalizability. The authors should also clarify whether\nsparsity assumptions used in simulations may bias in favor of RESCUE compared to other\nmethods.\nResponse: We agree that assessing generalizability and clarifying the role of sparsity\nassumptions are essential. We have therefore expanded the simulation study to include\nadditional scenarios by (i) varying effective grid sizes through changes in the average total\nUMI count per spot, and (ii) varying noise levels by adjusting library size variability across\nspots. These analyses demonstrate that RESCUE remains robust across a range of grid\nresolutions and noise conditions (Supplementary Figure 1; lines 170-176).\nImportantly, we clarified in the revised manuscript that the sparsity levels used in our\nsimulations are anchored to empirical observations from neuropil-rich regions of the honey\nbee brain (Figure 2a; lines 141-152), while also extending beyond these observations to\nstress-test model behavior. We further note that the sparsity assumption is biologically\nplausible in other settings, such as extracellular environments (lines 453-455). Thus, the\nsimulated sparsity settings are not tuned to favor RESCUE, but are designed to evaluate\nrobustness across a wide and conservative range of conditions.\n- In the honey bee data results: Intuitively, one would expect competing methods to\nattribute more expression to canonical (somatic) transcripts. However, Figures 3A and 3C\nshow the opposite pattern, other methods consistently estimate lower somatic expression.\nThis counterintuitive result should be explained.\nResponse: Thank you for noticing and raising what turned out to be a very subtle but very\nimportant point! We do agree that intuitively one would expect competing methods to\nattribute more expression to somatic transcripts, compared to RESCUE. This intuition\narises because at first blush, competing methods seem to be trying to find the linear\ncombination of the reference somatic expression profiles that can best decompose the\nobserved expression, while RESCUE allows for an additional unattributed component.\nHowever, the reality is more nuanced, because existing methods enforce a constraint on\nthis linear combination, in that its coefficients must add up to one. This is because these\nmethods are fundamentally not designed to decompose observed expression, but rather to\nestimate cell-type proportions (lines 575-585). In contrast, RESCUE places no such\nconstraint on this linear combination, so it can more completely decompose observed\nexpression into its reference components. This results in a larger amount of attributed\nsomatic expression.\nTo make this distinction explicit, we now quantify performance separately in regions with\nclear anatomical interpretation. Specifically, Figure 3c shows the average canonical\nexpression of each gene across grids corresponding to anatomically pure cell bodies,\nwhere RESCUE closely matches the ground truth while competing methods systematically\nunderestimate expression (lines 203-212). Conversely, Figure 3d reports the average\nunattributed expression across grids containing no cell bodies, where RESCUE accurately\nrecovers the residual signal, whereas subtraction-based extensions of other methods\nshow increased error (lines 212-215).\n- The comparison is limited to NNLS, WNNLS, RCTD, and CARD. More recent methods (e.g.,\nSONAR, Redeconve) are not evaluated. Since the field is rapidly evolving, omitting them\nweakens the benchmarking.\nResponse: Thank you for this suggestion. We have updated the benchmarking to include\ntwo more methods, SONAR and Redeconve, as additional competing methods. The results\nare included in the revised figures and analyses.\n- The manuscript relies mainly on average Lin’s concordance correlation coefficient (CCC)\nto assess cell-type deconvolution performance. More widely used metrics such as root\nmean squared error (RMSE) and Jensen–Shannon divergence (JSD) should also be reported\nto provide a fuller evaluation of accuracy and robustness.\nResponse: Based on the reviewer’s comment, we now report RMSE and JSD alongside\nLin’s CCC for cell-type deconvolution performance (Figures 2b-c, 4d, and Supplementary\nFigures 3-4).\n- One important application of RESCUE is to improve the cell type deconvolution. The\ncapability is demonstrated in one human breast cancer dataset. This result is too thin.\nThe authors need to add more results to show that this improvement is robust across different\nscenarios.\nResponse: We agree that demonstrating robustness across different scenarios is\nimportant. To address this, we added additional results using two further datasets analyzed\nin the manuscript: the human DLPFC and the cichlid fish telencephalon (Supplementary\nFigures 3-4; lines 279-287). In both cases, RESCUE shows improved performance across\nmultiple evaluation metrics, including Lin’s CCC, RMSE, and JSD, compared with\ncompeting methods. These results indicate that the benefits of RESCUE for improving cell-\ntype deconvolution generalize across tissues and species.\n- The assumption that the idiosyncratic component is sparse may not hold in all tissues\n(e.g., in extracellular RNA-rich tumor microenvironments). The authors should discuss\nscenarios where RESCUE might misclassify biologically meaningful widespread signals as\nnoise, and how users can detect such misbehavior.\nResponse: We would like to clarify that the sparsity assumption in RESCUE refers to the\nstructure of the unattributed component (i.e., that only a subset of genes is present per\ngiven spot), rather than to low expression levels (lines 144-146; lines 447-451). Thus,\nsparsity does not imply that the unattributed signal is weak or negligible, but rather that it is\nnot uniformly distributed across the transcriptome.\nWe further clarify in the manuscript that this assumption remains biologically plausible\neven in extracellular environments (lines 453-455). It is known that ECM proteins typically\ninvolve restricted gene sets rather than transcriptome-wide activation, and thus remain\nhighly sparse in this structural sense5. In addition, we validate RESCUE across a range of\nsparsity levels empirically observed in biological samples from the honey bee brain,\ndemonstrating robustness within realistic regimes (Figure 2; lines 141-152).\nWhile this assumption is biologically plausible across many tissues, we acknowledge that\nthere may exist scenarios in which the unattributed component is not sparse but instead\nexhibits alternative structure, such as low-rankness. Extending RESCUE to accommodate\nsuch residual regimes is an important direction for future work (lines 411-416).\n- The paper presents striking biological discoveries (e.g., mitochondrial gene enrichment in\nLayer 1 of DLPFC, RG-related programs in fish). However, these remain computational\npredictions. At minimum, the authors should (i) validate with orthogonal datasets (e.g., in\nsitu hybridization or existing histology), or (ii) discuss potential pitfalls of relying solely on\ncomputational inference without experimental follow-up.\nResponse: Thank you for this comment. We have revised the manuscript to explicitly\ndiscuss the limitations and potential pitfalls of relying solely on computational inference\nwithout experimental follow-up (lines 395-404). In many tissues, orthogonal validation of\nunattributed expression is not practically feasible because cell bodies and non-cell-body\nstructures (e.g., processes, neuropil) are spatially intermingled. Directly measuring the\nunattributed expression is particularly difficult in lower-resolution ST data, and even at\nhigher resolution, accurately delineating cell-body boundaries (especially in three-\ndimensional tissues) remains challenging. These limitations motivate the need for\ndiscovery-oriented computational approaches such as RESCUE, which are intended to\ngenerate biologically plausible hypotheses rather than provide direct experimental\nconfirmation.\nTo provide direct validation where feasible, we benchmarked RESCUE using the honey bee\nbrain, which serves as a gold-standard system due to the physical separation of neurites\nand soma (lines 190–192). This strategy enables clear ground truth for evaluating\nunattributed expression, and we confirm that a substantial fraction of transcripts resides in\nnon-segmentable neuropil regions (lines 218-225) using a representative MERFISH sample\n(Figure 3a) and DAPI image (Supplementary Figure 2).\n- The Methods section is difficult to follow. Clearer explanation of the optimization problem,\nidentifiability conditions, and tuning parameter selection is necessary. Readers less familiar\nwith robust regression might struggle to implement or extend the model. Equation 5 has\ntwo plus signs (+), might be a typo there.\nResponse: Thank you for pointing this out. We have substantially revised the Methods\nsection to improve clarity, with a more explicit explanation of the optimization problem, the\nidentifiability conditions, and the tuning procedure. In addition, we corrected a typo in\nEquation 5 by removing the extra plus sign. We hope these revisions make the Methods\nsection more accessible to readers who may be less familiar with robust regression.\n- The manuscript claims broad applicability (MERFISH, Visium, Slide-seq, etc.), but\ncomputational runtime, scalability, and memory requirements are not systematically\nreported. Given the size of current ST datasets, performance benchmarks (runtime vs. data\nsize) should be included.\nResponse: Thank you for this great suggestion. We now include systematic benchmarks\nreporting runtime as a function of data size across multiple datasets (Supplementary\nFigure 7; lines 649-655).\nReviewer #2 (Remarks on code availability):\nThe code and software page appear to be adequate, although I didn't try to run the codes\nmyself.\nReviewer #3 (Remarks to the Author):\nSummary\nThis manuscript presents RESCUE, a new method to recover expression from genes in\nmissing or underrepresented cell types following cell type deconvolution or cell\nsegmentation in spatial transcriptomics data. The authors refer to this type of expression as\n\"idiosyncratic\" expression. This is a useful method, which addresses a clear limitation in\ncurrent methods. The authors demonstrate performance on simulated and real datasets,\nand compare against several existing methods. The method is available as an R package\nfrom GitHub, including a tutorial, with data available from Figshare. The manuscript is\nclearly written and free of errors. The main comment below relates to providing additional\ninformation on computational requirements and scalability.\nResponse: Thank you very much for your positive comments. We have revised the\nmanuscript accordingly. Below, we outline the specific changes made in response to your\ncomments.\nTerminology note. Throughout the revision, we have replaced the term “idiosyncratic”\nexpression (or component) with “unattributed” expression (or component) to better reflect\nthe biological interpretation of signals not explained by the reference.\nMajor comments\n• No information is provided on computational requirements and scalability. This is\nimportant information for users, and should be added in additional text and/or figures - e.g.\ncomputational resources required (e.g. standard laptop or high-performance compute\ncluster / GPUs), runtime for the evaluations discussed in the manuscript, and scalability\nwith the number of cells / grids and/or genes.\nResponse: Thank you for this important suggestion. We now report systematic runtime\nbenchmarks as a function of data size across multiple datasets (Supplementary Figure 7;\nlines 649-655).\nMinor comments\n• The method is available as an R package from GitHub, which can be installed easily, and\nincludes a tutorial. However, the tutorial requires manually saving and loading a data file\nin .RData format, which may be not be clear for some users. If this data file is included in\nthe R package, it could be made easier to load by storing in 'inst/extdata/' and loading with\n'system.file()'. Alternatively, the start of the tutorial could be updated with some additional\ntext to explain how to load the data.\nResponse: We updated the tutorial to include additional instructions explaining how to\nsave and load the example data, clarifying the use of the provided data file.\nReviewer #3 (Remarks on code availability):\nSee main comments section\nREFERENCES\n1. Wang, G., Moffitt, J. R. & Zhuang, X. Multiplexed imaging of high-density libraries of RNAs\nwith MERFISH and expansion microscopy. Sci. Rep.8, 4847 (2018).\n2. Moffitt, J. R. et al. High-performance multiplexed fluorescence in situ hybridization in\nculture and tissue with matrix imprinting and clearing. Proc. Natl. Acad. Sci.113, 14456–\n14461 (2016).\n3. Mu, X. et al. Single-nucleus and spatial transcriptomics identify brain landscape of gene\nregulatory networks associated with behavioral maturation in honeybees. Nat.\nCommun.16, 3343 (2025).\n4. Ivich, A. et al. Missing cell types in single-cell references impact deconvolution of bulk\ndata but are detectable. Genome Biol.26, 86 (2025).\n5. Manabe, R. et al. Transcriptome-based systematic identification of extracellular matrix\nproteins. Proc. Natl. Acad. Sci.105, 12849–12854 (2008)."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have properly addressed my concerns. The manuscript is suitable for publication.", "reviewer_2": "All my concerns have been satisfactorily addressed.", "reviewer_3": "The authors have addressed all my comments from the previous review. I have no further comments."}, "author_response": ""}]

Decision letter