==== Front
Bioinformatics
Bioinformatics
bioinformatics
Bioinformatics
1367-4803
1367-4811
Oxford University Press

37387155
10.1093/bioinformatics/btad253
btad253
Regulatory and Functional Genomics
AcademicSubjects/SCI01060
Robust reconstruction of single-cell RNA-seq data with iterative gene weight updates
Sheng Yueqi School of Engineering and Applied Sciences, Harvard University, Boston, MA 02134, United States

Barak Boaz School of Engineering and Applied Sciences, Harvard University, Boston, MA 02134, United States

Nitzan Mor School of Computer Science and Engineering, Racah Institute of Physics, Faculty of Medicine, The Hebrew University of Jerusalem, Jerusalem 9190401, Israel

Corresponding authors. School of Engineering and Applied Sciences, Harvard University, Boston, MA 02134, USA. E-mail: ysheng@g.harvard.edu (Y.S.); School of Computer Science and Engineering, Racah Institute of Physics, Faculty of Medicine, The Hebrew University of Jerusalem, Jerusalem 9190401, Israel. E-mail: mor.nitzan@mail.huji.ac.il (M.N.)
6 2023
30 6 2023
30 6 2023
39 Suppl 1 ISMB/ECCB 2023 Proceedings i423i430
© The Author(s) 2023. Published by Oxford University Press.
2023
https://creativecommons.org/licenses/by/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Motivation

Single-cell RNA-sequencing technologies have greatly enhanced our understanding of heterogeneous cell populations and underlying regulatory processes. However, structural (spatial or temporal) relations between cells are lost during cell dissociation. These relations are crucial for identifying associated biological processes. Many existing tissue-reconstruction algorithms use prior information about subsets of genes that are informative with respect to the structure or process to be reconstructed. When such information is not available, and in the general case when the input genes code for multiple processes, including being susceptible to noise, biological reconstruction is often computationally challenging.

Results

We propose an algorithm that iteratively identifies manifold-informative genes using existing reconstruction algorithms for single-cell RNA-seq data as subroutine. We show that our algorithm improves the quality of tissue reconstruction for diverse synthetic and real scRNA-seq data, including data from the mammalian intestinal epithelium and liver lobules.

Availability and implementation

The code and data for benchmarking are available at github.com/syq2012/iterative_weight_update_for_reconstruction.

NSF 10.13039/100000001 DMS-2134157 DARPA 10.13039/100000185 W911NF2010021 DE-SC0022199 Azrieli Foundation Early Career Faculty Fellowship Israel Science Foundation 10.13039/501100003977 1079/21 European Union 10.13039/501100000780 101040660 European Union or the European Research Council European Union 10.13039/501100000780
==== Body
pmc1 Introduction

Single-cell RNA-sequencing (scRNA-seq) technologies have revolutionized our understanding of heterogeneous cellular populations and the underlying biological processes they encode (Wu et al. 2014; Potter 2018). However, the original spatial context of the cells (and the resulting spatial gene expression patterns) is lost due to the experimental procedure. While multiple computational methods have been developed to reconstruct spatial cellular context (Mages et al. 2023; Nitzan et al. 2019; Satija et al. 2015; Teves and Won 2020; Zeng et al. 2022), this task is challenging, partially because the cells simultaneously encode multiple spatiotemporal processes, and due to technical and biological noise. A common pre-processing step is considering only the highly variable genes as input for such reconstruction algorithms. However, these genes may not all be correlated with the same underlying process and thus may interfere with the reconstruction procedure.

The task of reconstructing cellular structure (or low-dimensional manifold in gene expression space) is tightly related to the task of identifying informative genes with respect to that manifold. On the one hand, a panel of manifold-informative marker genes can yield a metric similar to cell distances on the original structure, thus allowing for the clustering of cells that are close by in the original tissue (Droin et al. 2020). On the other hand, given an embedding of cells, one can recover informative genes with respect to that embedding by viewing individual gene expression as a function of cell location using, for example, regression-type algorithms (Droin et al. 2020). Here, we propose an algorithm that obtains enhanced reconstruction of single-cell data by alternating between these two tasks. By learning a representation of the subset of genes informative for a certain spatiotemporal process (or expressed smoothly along a certain manifold), we show that iterative gene weight update improves the reconstruction quality of existing algorithms, given the weighted gene expression as input, without relying on prior information.

 

We propose an algorithm that iteratively identifies manifold-informative genes, which can be integrated into existing reconstruction algorithms for single-cell data as a subroutine. Given single-cell RNA-seq data and a baseline learning algorithm that reconstructs gene expression along a spatial or temporal manifold, our algorithm iteratively updates the weight distribution over genes in the gene expression matrix taken as input to the reconstruction algorithm, to gradually increase the weight over genes related to the underlying reconstructed manifold. This leads to the identification of manifold-informative genes and a better reconstruction of the spatial or temporal manifold encoded by the cells. We demonstrate that such gene weight update improves the quality of reconstruction and identification of manifold-informative genes in several settings, including (i) diverse synthetic single-cell datasets when integrated with an Autoencoder as the baseline algorithm and (ii) over scRNA-seq datasets of mammalian intestine epithelium (Moor et al. 2018) and liver lobules (Droin et al. 2020), when integrated with either an Autoencoder or with novoSpaRc (Nitzan et al. 2019; Moriel et al. 2021) (an optimal-transport-based tissue reconstruction method) as the baseline algorithms.

2 Materials and methods

2.1 Model description and assumptions

Let {x1,…,xN}∈(Rn)N be gene expression profiles, where n:=number of genes and N:=number of cells. We assume that each cell, x, is represented by its locations on various low-dimensional manifolds M1,…,Mk, each representing a different biological structure or process. Here, we use location in a generalized way that includes time and additional cell labels belonging to other types of cell structures. However, the true cell locations are generally lost due to the scRNA-seq experimental procedure. For example, the spatial location of cells within the tissue of origin is lost due to cell dissociation. A gene is correlated with, or informative with respect to a manifold Mi if its expression at cell x can be written as a smooth function (e.g. a low-degree polynomial) of its location on Mi.Definition 1. Let M1,…,Mkbe manifolds where true cell locations are sampled from. For each cell, x, denote the true cell location as

c:=[c1⋮ck]

where ci∈Mi. Let Sibe the set of genes correlated to Mis.t. ∑i|Si|=n. Suppose there exists a set of functions f1,…,fks.t. fi:Mi→R|Si|and a permutation matrix π∈{0,1}n×n. The gene expression at cell x iswhere {yx}are i.i.d. samples from a distribution over Rn.Remark 1. In general, gene expressions are not deterministic functions and we do not assume a subset of genes correlated with the same manifold corresponds to a consecutive subset of coordinates of gene expression vector. In the above definition, we use πto include all gene orderings and {yx}to introduce randomness.

x:=π([f1(c1)⋮fk(ck)])+yx

Denote X=[x1,…,xN]∈Rn×N as the input matrix where each column is a cell. Given X, the goal is to recover Si and the corresponding true cell locations {ci}.

Note that some gene g may not depend on any cell manifold Mi. Those genes are not considered informative and we assume gene expression over cells is distributed as independent standard normal random variable.

Notation

Denote Πn={p∈ℝ+n:||p||1=1} to be the set of distributions over [n]={1,…,n}. Let un be the uniform distribution over [n]. For S⊂[n], let χS∈Πn be the uniform distribution over S and pS=|S|⟨p,χS⟩ be the weight over the set of genes S. Let Πnc⊂Πn be the set of distributions with ||p||∞≤c.

2.2 Description of algorithm

To enhance the reconstruction of single-cell data along a certain spatial structure or temporal process, we propose to “boost” the performance of a weak reconstruction algorithm by iteratively shifting the emphasis, or the weight distribution over genes, toward the group of manifold-informative genes using multiplicative weight update (Arora et al. 2012). A pseudo-code description is given in Algorithm 1. Suppose we have a baseline algorithm that can reconstruct cellular location and gene expression along an underlying manifold Mi when all the genes used as input to that algorithm correlate with that manifold and suppose we know the subset of genes Si. In that case, we can first project the input data X onto the subset of genes Si that correlate to Mi and use Xi as input to the baseline learning algorithm to reconstruct Mi. Generally, if we do not know Si and the input gene expression matrix contains genes that are not correlated with Mi, we expect the algorithm’s error to increase. Therefore, our goal is to generate a mask over the input gene expression matrix, which emphasizes, or places a higher probability density, over Si. This motivates the following definition.Definition 2 (Baseline algorithm). Given X and a distribution p over {1,…,n}, a weak algorithm Aklearns a k dimensional projection E:RN→Rkand D:Rk→RNs.t. for some ε>0, the reconstruction error of Ak, err(Ak(⋅)), is bounded by some monotone function fks.t. fk(χSS¯)=0

Xi=diag(χSi)X

err(Ak(diag(p)X))∈(1±ε)fk(pS¯)

where S is a subset of genes recoverable by Akand S¯=[n]/S.

In this work, we use round to refer to an iteration in step 3 of Algorithm 1, where we train the baseline algorithm on reweighted data based on the current weights. Maximum number of rounds can be chosen based on application, for example, either terminate the algorithm when the improvement in loss drops below a predetermined threshold (similarly to gradient descent), or when the relative entropy drops below a certain threshold (indicating gene weights concentration on a subset of genes) (Fig. 1c and g).

Figure 1. Gene weight updates enhance reconstruction quality over synthetic single-cell data. Experimental results for D1 (a–d) and D2 (e–h). (a, e) Mutual information between final encoding and gene expressions versus noise scale. Genes are partitioned into two sets S1,S2 based on the underlying manifold they depend on. We denote S2 as a subset of genes that are captured by the AE. For D1, S2 corresponds to the subset of genes dependent on Unif([0,1])2, and for D2, S2 corresponds to the subset of genes whose expression can be expressed as degree 2 polynomials. (b, f) MSE loss versus weight update rounds. (c, g) Relative entropy between algorithm output and ground truth distribution versus weight update round. (d, h) AUC scores for D1 and D2 versus weight update round over five experiments.

2.2.1 Baseline algorithms

A baseline reconstruction algorithm, as defined in Definition 2, is a subroutine to our proposed weight update algorithm. The goal of our algorithm is to enhance the performance of the baseline algorithm by focusing on manifold-informative genes for its respective reconstruction.

We test our weight update algorithm over two baseline reconstruction algorithms: feed-forward Autoencoder and novoSpaRc.

Autoencoder (AE). AE takes as input the reweighted data matrix diag(p)X, outputs an encoder E:RD→Rk, corresponding to a function that maps input gene expressions to a low-dimensional representation, and a decoder D:Rk→RD, corresponding to a function that computes the reconstructed gene expression from the low-dimensional representation. Autoencoders and their variants are widely used for single-cell reconstruction (Eraslan et al. 2019; Amodio et al. 2019; Tran et al. 2021). In this paper, we train a shallow AE with Tanh activation in each round to demonstrate the effectiveness of our weight update algorithm. The largest layer width is set to 64–128, and the code dimension is set to be the minimum cell manifold dimension. The number of epochs is set to 70–100. Our algorithm outputs the subset of genes that can be better reconstructed by AE, which generally corresponds to the subset of genes correlated with a lower dimensional cell manifold or is generated by a simpler class of functions.

novoSpaRc. novoSpaRc (Nitzan et al. 2019; Moriel et al. 2021) is a generalized optimal-transport-based method designed to spatially reconstruct single-cell data based on an interpolation between a structural correspondence assumption (corresponding to spatial smoothness of gene expression) and potentially a reference atlas (such as imaging data for a subset of marker genes). novoSpaRc was shown to successfully reconstruct gene expression patterns over diverse tissues, including the intestinal epithelium and liver lobules. We show that our weight update algorithm improves the reconstruction quality of novoSpaRc by making it more robust when input data contains noisy or less informative genes.

2.2.2 Alternative minimization and multiplicative update

Given a baseline algorithm as in Definition 2, we wish to find a distribution p that minimizes the reconstruction error of Ak while not putting too much weight on any single gene. We do this via the multiplicative weight update algorithm with restricted distributions (Kale 2007). Let l(p,X,D(E(X))) be a loss function that measures the quality of reconstruction. The problem can be rephrased as the following optimization problem: where |S| is the size of the subset of genes correlated with the current cell manifold. That is, we minimize the loss over distributions that put at most 1ε|S| weight on any particular gene.

minp∈Πn1ε|S|min  E,Dl(p,X,D(E(X)))

At each round t, we compute the optimal encoder and decoder for X reweighted by pt−1 and update pt by where D(p,q)=∑ipi log piqi is the relative entropy. The first term of Equation (1) increases pt[i] for genes with lower reconstruction loss. The second term prevents pt from concentrating on a small subset of genes. Solving for pt, where mt=∇l(pt,X,Dt(Et(X))). Finally, we project pt onto Πn1ϵ|S| as follows: Let B be the set of genes with pt[i]>1ϵ|S|, set pt˜[i]=1ϵ|S| and rescale pt[i] for i∈B¯ so that they sum to 1−|B|ϵ|S|.

(1) pt:=arg minp∈Πn{η∇l(pt−1,X,Dt-1(Et-1(X)))Tp+D(p,pt−1)}

pt[i]=pt−1[i] exp(−ηmt−1[i]))∑jpt−1[j] exp(−ηmt−1[j]))

2.2.3 Loss function

Generally, we measure the reconstruction error of Ak by a weighted MSE loss: where Y=D(E(X)) is the reconstructed data matrix. The reconstruction error for each gene is and we can rewrite the MSE loss as 〈mt,pt〉.

l(p,X,Y)=1N∑i=1np[i]∑j=1N(X[i,j]−Y[i,j])2

mtg=∇p lpt,X,DtEtX[g]=1N∑xx[g] - Dt(Et(x))[g]2

To check that pt is moving in the right direction, one observation is that the marginal of pt over S is non-decreasing (Supplementary Information).

For novoSpaRc, we measure reconstruction error by the variance of gene expression over the output cell distribution. Recall that a gene is strongly correlated with the cell distribution if cells that are mapped to similar true cell locations have similar expressions for this gene. Thus, a natural generalization of MSE loss would be the sum of the variance of the genes. In addition, we can also compute the variance of the derivative of the mean gene expressions to enforce the smoothness assumption of true gene expressions. Let μ1,…μh:[N]→[0,1] be the cell distribution computed by novoSpaRc where h is the number of cell locations.

l(p,X,μ1,…,μh)=∑i=1np[i][1-λ∑j=1hVarμj(X[i])+λVarj(Eμj[X[i]]−Eμj−1(X[i]))]

where Varj is the empirical variance over true cell locations and λ in [0, 1].

Algorithm 1. Alternate minimization

Require: X∈Rn×N , ϵ

1: p0⇐un t⇐0

2: fort=0,1,…do

3:    (Et,Dt)=Ak(diag(pt)X)

4:    mt=∇pl(pt,X,Dt(Et(X)))

5:    pt+1[g]∝ exp(−ηmt[g])pt[g]

6:    Project pt+1[g] onto Πn1ϵ|S|

2.2.4 Evaluation

On Synthetic datasets, we evaluate the performance of the algorithm in the following way:

The encoder captures the true cell location for the corresponding cell manifold: we compute the mean and standard mutual information between the encoding of cells and gene expressions for each subset of genes.

The algorithm achieves better reconstruction of gene expression after each update of p and D, E: reflected by a decrease in the MSE loss after each update.

The algorithm recovers the correct subset of informative genes: reflected by a decrease in the relative entropy between the ground truth distribution and pt, D(χS,pt).

Round ptto recover the correct subset of genes: AUC score increases. By giving all genes s.t. pT[i]≳p0[i] the same label, we show that the true positive rate remains ≥90% while the false-positive rate decreases after  10 iterations.

For the real scRNA-seq datasets, we compute the Pearson’s correlation coefficients between known and reconstructed mean gene expressions. We first run the baseline algorithm on pre-weighted (original) data (output of round 0) and show that the Pearson’s correlation coefficients improve after ∼10 rounds of weight update.

2.3 scRNA-seq simulations

Table 2 is a summary of the simple synthetic single-cell datasets when genes correlate with a single-cell manifold. We first generated 2×104 synthetic cell locations. For each cell manifold, we simulated single-cell RNA-seq data by generating random low-degree polynomials as functions for gene expression. Each gene is normalized to have a mean of 0 and a variance of 1. To obtain a mixture of datasets, we stack two datasets containing genes correlated with a single-cell manifold so that each cell contains genes correlated with two manifolds. Table 3 contains a summary of mixtures of datasets.

Table 2. Synthetic datasets.

Cell distribution	Function type	
Unif([0,1]) , Unif([0,1])2	Polynomials with degree ≤3	
N(0,1) , N(0,I2)	Polynomials with degree ≤3	

Table 3. Mixture datasets.

Type of mixture	Cell location	Gene expression	
Different cell locations distribution (D1)	Unif([0,1])2 , N(0,I2)	Degree 3 polynomials	
Different degree for gene expression (D2)	N(0,1)	Degree 2 and 3 polynomials	
Different cell manifolds dimension (D3)	Unif([0,1])2 and Unif([0,1])	Degree 3 polynomials	
Different cell manifold dimension and gene expression degree (D4)	N(0,1) and N(0,I2)	Degree 2 and 1 polynomials	

When adding Gaussian noise to the data, we say it has noise scale β if the input is of the form where X is the noise-less data as described above and N¯ is the normalized noisy data where each entry N[i,j]∼N(0,1). Each gene expression in Xβ has a mean of 0 and variance of 1.

Xβ=1−βX+βN¯

2.4 scRNA-seq datasets

We analyze two publicly available single-cell RNA-seq datasets. In both datasets, the reconstructed low-dimensional manifolds correspond to one-dimensional spatial axes within the tissues of origin (crypt-to-villus axis in the intestinal epithelium, portal-to-central vein in liver lobules). Generally, zonated genes are genes whose expression can be expressed as a function of the manifold corresponding to the spatial axes.

The first dataset contains cells sampled along the crypt-to-villus axis in the mammalian intestinal epithelium (Moor et al. 2018). In the original study, cell locations are partitioned into seven spatial layers. Here, we denote zonated genes as highly expressed zonated genes with Q value <0.05 as given by Moor et al. (2018).

The second scRNA-seq dataset contains hepatocytes sampled from mammalian liver lobules (Droin et al. 2020). In the original study, genes were classified as being correlated to either spatial zonation across the liver lobule axis (space), or to the circadian rhythm (time), or to their combination. The two cell manifolds, space and time, are denoted as c=(s,t)∈{0,…,7}×{0,…,3}, respectively. Specifically, genes in (Droin et al. 2020) were classified into the following disjoint groups:

Space-dependent genes (Z): zg(s)=as2+bs+c

Time-dependent genes (R): rg(t)=a sin(ωt+b)

Space + time genes (Z + R): fg(s,t)=azg(s)+brg(t)

Space × time genes (Z × R): fg(s,t)=azg(s)rg(t)

Unlabeled genes: noisy genes or low variance genes across space and time.

where a, b, c are constants that are different for each gene. Notice that there is no need to explicitly model the gene expression for the use of the algorithm we propose for iterative weight updates: for the liver data it is used only for evaluation; for the intestine data it is not used at all.

In both cases, gene expression was pre-processed and normalized, as commonly done for scRNA-seq reconstruction algorithms (Moor et al. 2018; Nitzan et al. 2019; Droin et al. 2020), by log-transform, and z-score normalization.

3 Results

3.1 Reconstruction of synthetic single-cell datasets

We first evaluate the performance of our algorithm on various mixtures of synthetic single-cell datasets to showcase the identification of manifold-informative genes and an encoding of the corresponding cell location when signal-to-noise ratio <1. The synthetic datasets (Table 3) are mixtures of basic single-cell datasets (Table 2) that differ in cell manifolds dimension, distribution, and degree of gene expression, with noise scale β∈[0,0.75]. Employing iterative gene weight updates over single-cell reconstruction using an Autoencoder output improves both relative entropy and accuracy of reconstruction (Fig. 1, Table 1 and Supplementary Fig. S1). Specifically, iterations of gene weight updates lead to the gradual decrease in relative entropy between the output distribution and ground truth, decrease in gene expression reconstruction error, increase in AUC score for gene expression reconstruction, and decrease in false-positive rate (and consistently high, >93%, true positive rate) for manifold-informative genes (where in each case of the synthetic mixtures, the correct subset of manifold-informative genes could be obtained by thresholding their weight by 1−δn for some small δ=0.2) (Fig. 1 andSupplementary Figs S1–S5). The quality of reconstruction and gene partition gradually decreases with the increasing noise level, as expected (Fig. 1b and f and Supplementary Figs S1b and f and S2b–S5b).

3.2 Reconstruction of scRNA-seq intestinal epithelium data

Here, we show that coupling reconstruction to iterative gene weight update as a subroutine can improve the quality of tissue reconstruction based on scRNA-seq data. We first focus on scRNA-seq data from the mammalian intestinal epithelium (Moor et al. 2018). To capture spatial gene expression along the crypt-to-villus axis, we generated a one-dimensional embedding of the cells, based on the top 1000 highly variable genes (approximately 50% of which are labeled as zonated by Moor et al. 2018), using both the baseline Autoencoder and novoSpaRc. In both cases, iterative weight updates gradually improve the quality of tissue reconstruction (Fig. 2). Specifically, when tracking the correlation between the reconstructed and measured gene expression values, we find that the quality of reconstruction of spatially zonated genes improves with each iteration of the gene weight update (Table 1). Following 10 rounds of gene weight update, the reconstruction quality matches that of training only on zonated, highly variable genes (Fig. 2).

Figure 2. Iterative gene weight updates enhance spatial reconstruction of the intestinal epithelium. (a, g) Quality of reconstruction measured by correlation coefficients following weight updates when using an Autoencoder (a) or novoSpaRc (g) as a baseline algorithm. z hvg corresponds to the output when trained on the subset of zonated, highly variable genes. (b–l) Mean gene expressions for sample genes from each of the five clusters of zonated genes in (Moor et al. 2018), reconstructed by novoSpaRc (h–l) and AE (b–f) after 0 and 10 rounds of update (novo or AE 0 or 10, respectively), compared to when only trained on zonated highly variable genes (novo or AE Z).

Furthermore, our algorithm gradually shifts the distribution of weights during training, as shown by the decrease in relative entropy between the current weight distribution and uniform distribution over zonated genes with weight update rounds (Supplementary Fig. S8). While novoSpaRc’s quality of reconstruction is higher than that of the baseline Autoencoder, the qualitative effects described in this section are shared by both algorithms (Fig. 2).

3.3 Reconstruction of scRNA-seq liver data

We next focus on a more challenging scenario, when genes change along two structures or processes, encoded by multiple cellular manifolds, simultaneously. This is the case for hepatocytes within liver lobules, whose expression changes across both physical space (along the lobule axis) and circadian time, as reflected in a recently collected scRNA-seq dataset (Droin et al. 2020).

We construct a one-dimensional embedding of the cells using both an Autoencoder and novoSpaRc as baseline algorithms. We use the top highly variable genes (200 for AE and 500 for novoSpaRc) out of 15,000 genes for training, approximately 50% of which are zonated (Z, Z + R, or Z × R) genes.  Applying our algorithm for iterative weight updates gradually enhances the quality of reconstruction of zonated genes (Z, Z + R, Z × R) in the de novo setting (without using a reference atlas) (Table 1), and achieves comparable results to those obtained by using a reference atlas as prior knowledge (Fig. 3). This occurs due to the increased weight density over spatially informative genes achieved by the iterative weight update.

Figure 3. Iterative gene weight updates enhance spatial reconstruction of liver lobules. (a–f) Example outputs of the recovered mean gene expression based on the iterative weight update algorithm, compared to the true mean gene expression. ZT are ground truth mean gene expressions and reconst ZT are reconstructed mean gene expressions after 10 rounds of weight update using AE as the baseline algorithm. (g) Percentage of zonated genes labeled by the weight update algorithm for Z, Z+R and Z×R genes as well as non-zonated genes. Error bars are computed over five runs of training. (i, h) The improvement in reconstruction quality was measured by correlation coefficients after 10 rounds of weight updates for AE (i) and novoSpaRc (h). z hvg corresponds to the output when trained on the subset of zonated highly variable genes.

Furthermore, we are able to recover the correct labels for most zonated genes (Z, Z + R, Z × R) using the final cell embedding output by the algorithm, as described below. Similarly to Droin et al. (2020), we recover gene labels from the cell embedding by using a linear regression model to approximate (reconstructed) mean gene expression by low-degree polynomials of s and t. In (Droin et al. 2020), a mixed regression model was used to approximate measured mean gene expression, however, since we focus here on spatially informative genes whose expression changes monotonically with space, the zonated part of gene expression could be approximated by a monotone function captured by linear regression. Thus, we fit each mean gene expression to a linear function and label genes that can be approximated well (R score >0.1) and are informative with respect to the spatial manifold (slope >0.001). Our algorithm achieves 93.7% true positive rate on purely zonated genes and >80% of true positive rate on mixed zonated genes (Fig. 3g).

4 Discussion

Our algorithm for iterative gene weight update improves the accuracy of single-cell RNA-seq reconstruction. Our algorithm can reconstruct gene expression, including the case when multiple cell manifolds are simultaneously encoded, without depending on prior information such as manifold-informative marker genes. We showcase these results over two baseline reconstruction algorithms, an Autoencoder and novoSpaRc, and over synthetic single-cell data as well as two scRNA-seq datasets of mammalian intestinal epithelium and liver lobules. Moreover, our algorithm can be incorporated into a general baseline reconstruction algorithm, treating it as a blackbox. Therefore, while in this manuscript we focus on enhancing spatial reconstruction algorithms that are based on scRNA-seq data alone, our algorithm for iterative weight update can also be applied to the reconstruction of underlying structures in the data corresponding to non-spatial processes such as temporal trajectories (see benchmarking of trajectory inference methods in Saelens et al. 2019), as well as algorithms that combine multiple types of experimental data, such as scRNA-seq and spatial transcriptomics data for spatial reconstruction. While our algorithm requires the selection of an appropriate loss function for the chosen baseline algorithm that effectively measures the reconstruction quality of each gene, it can generally be informed by the objective and characteristics of the baseline algorithm.

Supplementary Material

btad253_Supplementary_Data Click here for additional data file.

Supplementary data

Supplementary data is available at Bioinformatics online.

Conflict of interest

None declared.

Funding

This work was supported by Simons Investigator Fellowship, NSF grant DMS-2134157, DARPA grant W911NF2010021, and DOE grant DE-SC0022199 (B.B.), an Azrieli Foundation Early Career Faculty Fellowship, the Israel Science Foundation (Grant no. 1079/21), and the European Union (ERC, DecodeSC, 101040660) (M.N.). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them.

Data availability

The scRNA-seq datasets used for this study were acquired from the Gene Expression Omnibus (GEO) database with the following accession numbers: liver (GSE145197), intestine (GSE109413). Reconstructed mean gene expression for the liver data was downloaded from (https://github.com/naef-lab/Circadian-zonation). Zonated gene labels for the intestine were downloaded from (https://zenodo.org/record/1320734). The data and code for generating synthetic single-cell datasets are included in github.com/syq2012/iterative_weight_update_for_reconstruction.

Table 1. Summary of results.

Dataset	Baseline	With weight update subroutine	
Real(Cor Coef)	AE	novoSpaRc	AE	novoSpaRc	
Intestine	0.69±0.11	0.86±0.04	0.86±0.01	0.91±0.01	
liver	0.70±0.23	0.69±0.14	0.87±0.04	0.82±0.04	
	
Synthetic(MSE)	AE (β=0)	AE (β=0.4)	AE (β=0)	AE (β=0.4)	
	
Cell dist (D1)	0.76	0.86	0.10	0.54	
Gene degree (D2)	0.79	0.89	0.16	0.62	
Cell dim (D3)	0.64	0.77	0.14	0.48	
Cell dim + gene deg (D4)	0.7	0.82	0.15	0.49	
Summary of reconstruction results on synthetic and scRNA-seq datasets for the baseline algorithms (Autoencoder and novoSpaRc) without and with the integration of the weight update algorithm as a subroutine. For real scRNA-seq datasets, we record the mean and variance of correlation coefficients between reconstructed and ground truth zonated gene expressions. For Synthetic datasets, we measure the MSE loss on different noise levels. Higher correlation coefficients and lower MSE loss indicate better gene expression reconstruction. In both cases, the integration of the weight update algorithm enhances the quality of reconstruction.
==== Refs
References

Amodio M , van DijkD, SrinivasanK et al Exploring single-cell data with deep multitasking neural networks. Nat Methods 2019;16 :1139–45.31591579
Arora S , HazanE, KaleS. The multiplicative weights update method: a meta-algorithm and applications. Theory Comput 2012;8 :121–64.
Droin C , El KholteiJ, HalpernKB et al Space-time logic of liver gene expression at sublobular scale. Nat Metab 2020;3 :43–58.
Eraslan G , SimonLM, MirceaM et al Single-cell RNA-seq denoising using a deep count autoencoder. Nat Commun 2019;10 :390.30674886
Kale S. Boosting and hard-core set constructions: a simplified approach. Electron. Colloquium Comput. Complex 2007;14 :131.
Mages S , MorielN, Avraham-DavidiI et al TACCO unifies annotation transfer and decomposition of cell identities for single-cell and spatial omics. Nat Biotechnol 2023;16 :1–9.
Moor AE , HarnikY, Ben-MosheS et al Spatial reconstruction of single enterocytes uncovers broad zonation along the intestinal villus axis. Cell 2018;175 :1156–67.e15.30270040
Moriel N , SenelE, FriedmanN et al NovoSpaRc: flexible spatial reconstruction of single-cell gene expression with optimal transport. Nat Protoc 2021;16 :4177–200.34349282
Nitzan M , KaraiskosN, FriedmanN et al Gene expression cartography. Nature 2019;576 :132–7.31748748
Potter SS. Single-cell RNA sequencing for the study of development, physiology and disease. Nat Rev Nephrol 2018;14 (8 ):479–92.29789704
Saelens W , CannoodtR, TodorovH et al A comparison of single-cell trajectory inference methods. Nat Biotechnol 2019;37 :547–54.30936559
Satija R , FarrellJA, GennertD et al Spatial reconstruction of single-cell gene expression data. Nat Biotechnol 2015;33 (5 ):495–502.25867923
Teves JM , WonKJ. Mapping cellular coordinates through advances in spatial transcriptomics technology. Mol Cells 2020;43 :591–9.32507771
Tran D , NguyenH, TranB et al Fast and precise single-cell data analysis using a hierarchical autoencoder. Nat Commun 2021;12 :1029.33589635
Wu AR , NeffNF, KaliskyT et al Quantitative assessment of single-cell rna-sequencing methods. Nat Methods 2014;11 :41–6.24141493
Zeng Z , LiY, LiY et al Statistical and machine learning methods for spatially resolved transcriptomics data analysis. Genome Biol 2022;23 :83.35337374