Cell2State: Learning Cell State Represen- tations From Barcoded Single-Cell Gene- Expression Transitions
Year: 2021

Abstract:
Genetic barcoding coupled with single-cell sequencing technology enables direct measurement of cell-to-cell transitions and gene-expression evolution over a long timespan. This new type of data reveals explicit state transitions of cell dynamics. Motivated by dimension reduction methods for dynamical systems, we develop a cell-to-state (cell2state) learning method that, through learning from such multi-modal data, maps single-cell gene expression profiles to low-dimensional state vectors that are predictive of cell dynamics. We evaluate the cell2state method using barcoded stem cell dataset (Biddy et al. (2018)) and simulation studies, compared with baseline approaches using features that are not dynamic-aware. We demonstrate the merits of cell2state in challenging downstream tasks including cell state prediction and finding dynamically stable clusters. Further, our method reveals potential latent meta-states of the underlying evolution process. For each of the meta-states, we identify a set of marker genes and development pathways that are biologically meaningful and potentially expand existing knowledge.

heading: 1 Introduction
text: With the explosive amount of data from single-cell genomics studies, one remaining major challenge is the lack of ability to understand cell transition on the individual level. Conventional methods for analyzing single-cell dynamics are mostly based on “ensemble” analysis (Kester & van Oudenaarden, 2018; Tanay & Regev, 2017; Bacher & Kendziorski, 2016; Stegle et al., 2015). Such analysis reveals population-level trends, but cannot reveal behaviors of individual cells.
Recent advances in genomic technology have enabled scientists to directly measure cell lineages and connect two cells that are far apart in the time course. This is opposed to inferring lineages from snapshots at nearby time points. The gene barcoding approach works by inserting into each cell a DNA sequence, i.e., a barcode, that randomizes across cells so that no two cells bear the same sequence (Woodworth et al., 2017). As the cell divides and differentiates, its descendants can be identified based on sequencing the label. This concept has been utilized recently in single-cell analysis of embryonic development (Yao et al., 2017; Wagner et al., 2018), stem cell reprogramming (Biddy et al., 2018), and fate determination in hematopoiesis (Weinreb et al., 2020). This genetic barcoding approach enables tracking of evolutionary trajectories across individual cell lineages.
In this paper, we focus on the new type of single-cell data enabled by the single-cell gene barcoding technology. The barcode can directly connect parent cells with their descendants over a long time span. Thus, gene barcoding coupled with RNA-seq generates pairs of gene-expression transitions {(X,X ′)}, where each (X,X ′) is the gene expression profiles for a parent cell at an early time point and one of its descendants at the later time point. One may view the gene-expression profile X as the raw state of a cell, which is a high-dimensional vector. Such state-transition data makes it possible to learn about a cell’s law of state transition. In other words, the new data type can let us decode the single-cell transition law in a way similar to system identification for dynamical systems.
We wish to learn mathematical abstractions of cell gene-expression states, i.e., a map Ψ from gene-expression profile X to a vector Ψ(X) of lower dimension. A good cell state abstraction should be low-dimensional and predictive, compressing predictive signals from gene-expression levels in a compact vector. In other words, we hope to maximize the mutual information I(Ψ(X), X ′) between the embedded parent cell state and its descendant. Ideally, we hope to achieve a nearly lossless encoding of states, i.e., I(Ψ(X), X ′) ≈ I(X,X ′). To learn such Ψ from noisy high-dimensional gene-expression data, we build on ideas from dynamical system theory, kernel machine learning, and low-rank optimization. In the area of molecular dynamics, Schütte et al. (2011) showed that the leading spectrum of transfer operator can be used to identify coresets for faster simulation. In reinforcement learning, computers need to figure out state abstraction of unknown transition systems in order to learn to control quickly (Sutton et al., 1998). Sun et al. (2019) developed a state embedding method to analyze game trajectories to significantly reduce the state dimension of one-player Atari games. See Appendix A for more discussion on related works.

heading: A summary of our work:
text: • Building on the spectral compression ideas from dynamic systems, we develop a cell-to-state (cell2state) representation learning method for analyzing gene-expression transition data made available by the single-cell barcoding technology. The cell2state algorithm is trained using gene-expression transition pairs; it finds a mapping to approximate and embed the transition distributions in a low-dimensional space. The embedding map is learned by first “lifting” the data’s dimension by random feature projection, then estimating a large matrix embedding of the transition distributions, and finally “compressing” the dimension by low-rank optimization.
• We provide information-theoretic analysis of the learnt cell2state embedding Ψ̂. In particular, we show that, upon appropriate quantization, the embedding map can be used to encode raw gene-expression data into a small number of bits. We show that this encoding can be nearly lossless, and we establish sample complexity bounds for preserving the mutual information between parent and descendants up to 1 bit.
• We apply cell2state to a published single-cell barcoded RNA-seq dataset for studying stem cells (Biddy et al., 2018). In this analysis, we used cell2state to map early-day cells to state vectors of dimension ≤ 8. Via the cell2state map, the cell populations demonstrated sharp polytope structures, where distinct vertices provide early signals that predict diverse dynamics and cell fates. To evaluate the learned cell state vectors, we test them in three downstream tasks: (i) finding dynamically stable cell cluster; (ii) early prediction of cell dynamics, such as cell descendants’ activities or fates, based on low-dimensional state vectors; (iii) subpopulation analysis to identify marker genes that signal distinct cell dynamics. Across these tasks, we observe substantially improved performance using the learned cell states, as compared with baselines that use either the raw data or features that are not dynamics-aware. In particular, our results show that cell2state achieves similar/better level of prediction accuracy using ≤7 dimensions, compared to neural networks that use raw gene expressions as input (up to 5000 dimensions).
• Further, we identify and examine subpopulations of cells that have the most representative low-dimensional states (in other words, subsets of cells that are close to likely meta-states, under the assumption of a latent-state model). These subpopulations are used to identify biologically relevant marker genes. These marker genes identified by cell2state are known to relate to stem cell reprogramming and epigenetic regulations.

heading: 2 Markov Branching Process Model for Single-Cell Dynamics
text: We model the time-course dynamics of gene expressions as a branching diffusion process (Edwards (1970); see Figure 1). Let Xt denote the gene-expression profile of a cell at a time point t, which is a high-dimensional vector. Each Xt has a random number of descendants Xit+1, i = 1, . . . , N with independent and identical distributions.
Definition 1. Define p(Xt+1|Xt) as the transition function for the gene expression profile Xt to evolve to a collection of descendants {X(i)t+1}Ni=1 in a fixed amount of time, i.e., for any measurable set S,
p(S|Xt) = E [ N∑ i=1 1{X(i)t+1 ∈ S} | Xt ] ,
where N is also a random variable.
Note that p is not necessarily a probability density function because a cell could have multiple descendants. If the growth rate is such that E[N |Xt] = ∫ p(y|x)dy > 1, we say the cell is actively growing.
Definition 2. Let P be the transition operator of the branching diffusion process with transition function p, given by
Pf(X) = E [ N∑ i=1 f{X(i)t+1 ∈ S} | Xt = X ] .
In single-cell analysis, the transition function, p, and operator, P, are infinite-dimensional, thus estimating them is largely intractable from finite noisy data. Like many other dynamic systems, cell dynamics is often driven by a small set of marker genes, and thus, it may admit an intrinsic low-dimension structure. We make the following assumption: Assumption 1. Let H be a space of functions. There exists a r-dimensional embedding map Φ∗ ⊂ H such that
Pf ∈ Span(Φ∗),∀f ∈ H.
Here, H will be specified later.
Such low-rank structure of Assumption 1 naturally exists in dynamical processes that admit latent states. Suppose that each x can be represented as a mixture over meta-states {z} such that
p(x′|x) = ∑ z p(x′|z)pZ(z|x).
This is a common latent-state model for stochastic processes; see Figure 2 for an illustration. In the single-cell context, a meta-state is often referred to as a “cell type”, which has a distinct “pathway” (i.e., future
dynamics). The “cell type” is defined as a function of the gene-expression profile, but the function is unknown and to be learnt. In this latent-state model , let Φ∗(x) = pZ(·|x). Then we can verify that Assumption 1 holds. In this case, finding the embedding map Φ∗ would make it possible to recover the set of meta-states {z} and aggregation distributions pZ(·|x).

heading: 3 Mapping Gene Expressions To Low-Dimensional Cell States
text: Recall that our goal is to find mathematical abstractions of high-dimensional expression profiles {Xt}. We will estimate an embedding map from gene-expression profiles to a low-dimensional vector space: Ψ : x ∈ Rd 7→ Ψ(x) ∈ Rr. Ideally, we hope Ψ(Xt) to be low-dimensional while still containing as much information about Xt+1 as possible.

heading: 3.1 Embedding the transition operator into a functional space
text: Consider a kernel mean embedding of the transition operator P: Q = ΠHPΠH, where the projections are with respect to appropriate norms. By assumption, we can verify that rank(Q) < r. We seek to estimate Q from cell transition data, perform singular value truncation, and then find the low-dimensional embedding map by transforming the left single functions of Q.
To guarantee the function space H is sufficiently expressive, we adopt a kernel composition approach for “lifting” the dimensions. First, we construct an initial kernel K0 that best fits the dataset’s topology and preserves neighborhood relations. To find such a K0, we can leverage existing dimension reduction methods for single-cell data analysis, such as PCA, manifold-based, and graph-based methods (see Kester & van Oudenaarden (2018); Tanay & Regev (2017); Bacher & Kendziorski (2016); Stegle et al. (2015) for reviews). Then, we construct the kernel function K = K0 ◦ K1 by taking the composition between K0 and another kernel function K1 (e.g., the Gaussian kernel) - this step would further lift the problem’s dimension and improve the function space’s expressibility. One can also take compositions of multiple kernels to mimic a multi-layer neural network (Cho & Saul, 2009).

heading: 3.2 Low-rank compression of cell states via random features
text: We propose a kernelized state embedding method based on random feature projection for computing an estimator Ψ̂ from transition data {(Xi, X ′i)}Ni=1, which can be obtained from cell trajectories. For analyzing single-cell sequencing data and embedding transition distributions, we will choose the function space H with the kernel function K tailored to the data’s geometry.
Suppose we have chosen a kernel function K (for example, a Gaussian kernel). We perform nonparametric estimation of Ψ∗ by generating a large number of random features to approximate the kernel space in large finite dimensions. Then, we downsize the estimator by using spectral decomposition. In the case where each parent cell has a single descendant, the cell2state method works by (informally):
(1) Generate random Fourier functions φ(·) = (φ1(·), . . . , φd(·))> by randomized decomposition of its kernel function K to approximately span H (Rahimi & Recht, 2008). (2) Estimate a finite matrix embedding of the scaled condition probability distribution
1√ p(x′)
p(x′|x) by P̂ = Σ−1/20 ( 1 N ∑N i=1 φ(Xi)φ(X ′ i) > ) Σ −1/2 1 , where Σ0,Σ1 are covariances
at the two time points.
(3) Let Ψ̂(·) = (ÛrΛ̂r)>Σ−10 Φ(·) where Ûr, Λ̂r are from top r truncation of the SVD of P̂ .
See Algorithm 1 in the Appendix for the full description of the cell2state algorithm, which also handles the case where cells have multiple descendants. Given a cell’s gene expression profile x, the vector Ψ̂(x) can be viewed as a low-dimensional mean embedding of the transition function p(·|x). Thus it should be predictive of this cell’s future dynamics. Runtime Complexity of Algorithm 1. The algorithm uses random features and singular value truncation, both designed for maximal computation efficiency. The overall runtime for training is at most O(n+ nD2 +D3), where n is number of cells, D is number of Fourier features (≤ n). This is the same complexity as computing covariances and PCA in the random feature space. After training, querying the embedding map Ψ̂ takes only O(rD) time. In our experiments, Algorithm 1 runs in seconds, while training an MLP (a deep neural network) for cell fate prediction takes 10-15 minutes.

heading: 3.3 Information-Theoretical Analysis
text: In this section, we analyze the information-theoretic property of the cell2state embedding map Ψ̂. Assume without loss of generality that the feature φ(x) is upper bounded by ‖φ(x)‖2 ≤ C,∀x. Let P = E[P̂ ]. Using an analysis similar to Sun et al. (2019), it can be shown that the distance distortion of state embedding map satisfies
Theorem 1. Let Assumption 1 hold. With probability 1− q, then for all x, y in the dataset,
|‖Ψ̂(x)− Ψ̂(y)‖ − ‖p(·|x)− p(·|y)‖H| ≤ 16Cκ
γ √ log 2dq N + 32C 3 2κ γ 3 2 log 2dq N
where d is the dimension of φ, and N is the number of sample transitions, γ−1 = max{‖Σ−10 ‖, ‖Σ −1 1 ‖}, κ = ‖P‖ σr(P ) .
Further, we will apply quantization to Ψ̂ and encode the parent cell data into finitely many bits. Suppose the state space admits a block structure, i.e., we have a mapping Ω∗ : X → [k], such that
p(·|x) = p(·|x′), ∀x, x′ s.t. Ω∗(x) = Ω∗(x′).
Suppose the SVD of P is P = UrΛrV >r and let Ψ(x) = (UrΛr)>Σ − 12 0 φ(x). Let A = {Aj} be an arbitrary quantization of Ψ̂(X ) such that ‖x − x′‖ ≤ δ8 ,∀x, x ′ ∈ Aj ,∀j, where δ = minx,x′:Ψ(x) 6=Ψ(x′) ‖Ψ(x)−Ψ(x′)‖. In addition, define the encoding map Ω̂Ψ̂ such that
Ω̂Ψ̂(x) = Aj , iff Ψ̂(x) ∈ Aj .
We can show that the learned encoding map largely preserves the mutual information between a parent cell’s raw gene expression and its descendants’ gene expression profiles. More precisely, we show that the loss of information can be bounded and estimated, as follows Theorem 2. Let Assumption 1 hold with the aforementioned block partition structure. The estimated encoding map Ω̂Ψ̂ is nearly loss-less in the following sense:
E[I(Ω̂Ψ̂(X), X ′)] ≥ I(X,X ′)− 2d exp
( − 3γ 2Nδ2
800C2κ2
) log k,
where the expectation is over the distribution that generates the transition data.
Therefore, in order to ensure EΨ̂[I(Ω̂Ψ̂(X), X ′)] ≥ I(X,X ′)−∆, we need the sample size N ≥ O ( C2κ2
γ2δ2 log d log k ∆
) . When ∆ = 1, we can quantify the sample size needed to preserve
mutual information up to 1 bit difference.

heading: 4 Experiment with Cell Reprogramming Data
text: We use the single-cell gene expression data with genetic barcoding tags from Biddy et al. (2018)to test our cell state embedding method. Biddy et al. (2018) tracked reprogramming of mouse embryonic fibroblasts (MEFs) to induced endoderm progenitors (iEPs) by using CellTagging technology. During this dynamic reprogramming process, DNA barcodes in the form of randomized nucleotides were delivered into pools of cells and remained in the cell via lentiviral genome integration. Through identifying cells with the same DNA barcode sequence, we can recover the lineage relationships between parental and descendent cells. Barcoded single-cell gene expressions were collected using scRNA-seq at 8 time points over a course of 28 days (Figure 3(a) for a tSNE visualization). A cell’s raw gene-expression profile is a vector of dimension 28,001.

heading: 4.1 Cell2state maps ancestral cells to a tetrahedron and implies the existence of at least 4 meta-states
text: For our experiment, we pick those cells profiled on day 12 and day 21, as parent cells and descendants respectively (Figure 3(b)), as these time points had a reasonably high number of sampled cells and pairs of genetic tags to establish lineage trajectories. After removing low-quality cells, we obtain a count matrix retaining 6233 cells (1997 on day 12 and 3509 on day 21) and spanning 17845 genes. We then process the single-cell data using their
cell tags and yield N =165716 cell-to-cell transition pairs. We pick the kernel K to be a composition between the principal component map (obtained by PCA) and a Gaussian kernel with tunable width γ. Then, we apply the cell2state method and compute the embedding map Ψ̂. We visualize in Figure 3(c) the top three features of the parent cells given by Ψ̂. These parent cells visibly form a tetrahedron with four vertices in the embedding space, implying at least four potential meta-states with distinct future pathways. Further, we validate that cells near the vertices and cells at the center have diverse cell fates based on gene expression profiles of their descendants on day 21 (Figure 3(d)).

heading: 4.2 Evaluation of low-dimensional cell states on downstream tasks
text: Next we investigate the predictive power and biological relevance of the low-dimensional cell state embeddings {Ψ̂(x)}. We consider three downstream tasks: clustering, prediction, and gene marker/pathway analysis.

heading: 4.2.1 Dynamics-stable cell clustering
text: A fundamental task in studying single-cell dynamics is to cluster cells into representative “cell types/states” that are biologically meaningful and stable across time. Many innovative tools have increasingly permitted the identification and classification of cell types/states (Kester & van Oudenaarden (2018); Wagner & Klein (2020)), but most of these tools do not account for cells’ temporal dynamics. With new barcoding data, we seek to integrate information from both gene expression and lineage trajectories to find dynamic-stable cell clusters.
Dynamic stability of embedding clusters. We apply k-means to the parent cell embeddings and identify 7 major clusters (Figure 3(f)). We evaluate the quality of the cluster assignment via examining their respective descendant distributions (Figure 3(g)). We clearly observe that the clusters are dynamically stable, i.e., descendants from the same parental cluster tend to stay near to one another (panel f-g). For comparison, the same cluster assignment in the raw data of parent cells (Figure 3(e)) is unstructured and mixed up. The contrast between Figure 3(e-f) suggests that these dynamically stable cluster structures were hidden in the raw data, but can be revealed by our low-dimensional embedding.
Comparison via computing the cluster assignment’s descendant inertia. We further compare the dynamic stability of cell2state clusters with the widely-used graph-based clustering method in Seurat that integrated Louvain algorithm (Blondel et al. (2008); Butler et al. (2018); Stuart et al. (2019)To quantitatively evaluate the dynamic stability, we compute the inertia of cluster assignments using descendant data. For the cluster assignment Ω = {Ω1, ...,Ωk}, we evaluate the inertia of Ω over all descendant cells, i.e., inertia(Ω) = ∑ i∈[k] minµi ∑ X1∈Ωi ||X1 − µi|| 2 2. The inertia measures the level of concentration of the clusters, thus if the clusters with smaller inertia across time are more dynamically stable. Figure 3(h) shows the descendant inertia values obtained from different cluster assignments. It shows that the cell2state k-means clusters from parent cells achieved a small inertia on descendant cells, similar to that from directly clustering descendants. Both are significantly smaller than assigning clusters on raw data without doing cell2state, which are not dynamics-stable. This validates that cell2state yields dynamics-stable cell clusters, i.e., inertia of the clusters remains small across time.

heading: 4.2.2 Early prediction of descendant cells’ proliferation activity
text: Next, we evaluate the predictive power of the learned cell state embeddings. One of the fundamental cell states that connects to the underlying biology of stem cells is a cell’s proliferation potential. Hence, we seek to use the cell state embedding learned for parent cells to predict proliferation activity of these cell’s descendants. Note that the parent and descendant in our experiments are 9 days apart. This is a relatively long time span, and to the best of our knowledge, there does not exist a prior attempt for such early prediction with mammalian cells.
Prediction targets. We measure the proliferation activities of cells on day 21 using two well-established metrics: (i) the G2M cell cycle score, defined by prior studies and widely used
in single-cell data analysis (Tirosh et al. (2016); Butler et al. (2018); Stuart et al. (2019)); (ii) the average expression level of 366 proliferation genes, designated by gene ontology (GO) annotations (Ashburner et al. (2000); Consortium (2019)).
Results and comparison. Figure 3(i,k) visualize the in-sample and out-of-sample prediction error as the model parameter (i.e., γ determines kernel width) varies. For comparison, we also trained predictors that use the principal components of raw data directly for the same prediction task, which we treat as the baselines. Observe that the cell2state-based predictors consistently outperforms the baseline across instances. Figure 3(j,l) illustrate the embedding dimension needed as the predication accuracy level varies. The results suggest that our cell embedding approach has substantially reduced dimensionality of the raw data, achieving similar or higher prediction accuracy using a fraction of dimensions.

heading: 4.2.3 Finding marker genes and cell development pathways
text: Finally, we seek to interpret the cell state embeddings and evaluate if these learned lowdimensional structures are biologically relevant. We examine the polytope structure of parent cells in the cell2state embedding space (Figure 3(c)). Based on the cluster assignment, we then select sample cells that are close to vertices (top 50% of cells in the cluster) and apply DESeq2 (Love et al. (2014)) to these representative cells from each vertex. DESeq2 is a tool that allows one to identify significantly enriched genes in a subpopulation as compared to the full data. This allows us to find marker genes that distinguish individual subpopulation from one another.
Finding top ranker marker genes Table in Figure 4(a) summarizes a list of top positive marker genes associated with each vertex as ranked by fold-change (logFC). We note that our results, from day12 cells only, already recover the majority of the genes/pathways implicated in Biddy et al. (2018), such as Apoa1, Col1a2, Peg3, as well as Wnt and Igf2 pathways (highlighted in bold).
Analysis of the development pathways of each set Additionally, our analysis demonstrates that each subpopulation has distinctive signature tissue programs via Genotype-Tissue Expression (GTEx)(Carithers & Moore (2015)), which reveals their putative cell fates and development pathways (Figure 4(b)).

heading: 5 Experiment with Synthetic Data
text: To further test cell2state, we artificially construct highly nonlinear, random cell dynamics for a simulated experiment. We use gene-expression data from Weinreb et al. (2020)but we added artificially nonlinear transitions between day2 and day6 cells. We generate simulated day2-day6 transitions with a Markov branching process as follows: We first partition the 2D UMAP embedding space for day2 cells into 3 regions/clusters using two concentric circles (Figure 5(a)). We also partition the 2D UMAP space for day6 cells into the 4 quadrant
regions, i.e, 4 cell fates (Figure 5(b)). For each day2 cell X, let c(X) ∈ {1, 2, 3} denote its cluster assignment. Then the number of descendants is given by N ∼ Poi(λc(X)), where λi is fixed for i = 1, 2, 3. We generate cell transitions using the cluster-to-cluster transition matrix P , by sampling descendants from the corresponding day6 cluster. This Markov branching process has highly nonlinear, discontinuous dynamics due to our artificial construction of the three initial regions and the four cell fates.
Application to cell fate prediction. We apply cell2state to the simulated cell transition data and learn low-dimensional state vectors of all parent cells. Figure 5(c) shows that cell2state maps the three parental clusters to a triangle, where each latent parental cluster is mapped to one vertex. We then use these learned cell states to predict the dominant cell fate of each parent cell, using a linear classifier. We consider both the cell2state based on random Fourier features and a variant of cell2state where exact kernel decomposition is used instead of random features. For comparison, we also tested the basic linear classifier and a multi-layer perceptron neural network, both taking raw gene-expression profiles as input for cell fate prediction (Appendix B for simulation experiment details).
Comparison with linear and neural network classifiers. Figure 5(d-e) illustrate the simulation results. We see that cell2state together with a basic linear classifier achieves the best out-of-sample accuracy across all experiments. The cell2state predictors using the low-dimensional states have comparable or slightly better performances than neural networks that take as input the raw data. This validates our theory that cell2state has largely compressed the high-dimensional information about the nonlinear cell dynamics into lower dimensions. Also note that, the basic linear classifier using raw data performs poorly in comparison to cell2state and neural networks, as it fails to capture the process’s nonlinearity. Overall, this simulation demonstrates cell2state’s ability to learn meaningful cell state representations even from highly irregular, nonsmooth dynamics.

heading: 6 Summary
text: We provide a random feature-based cell state embedding method for mapping single-cell gene expression profiles to low-dimensional representations based on barcoded gene-expression trajectories. Application to stem cell dataset and simulation studies suggests the learned cell state embedding carries predictive signals of cell dynamics.

heading: A Background and Related Works
text: In recent years, enabled by our ability to track and measure individual cells, single-cell analysis tools have demonstrated their significance in bringing next-level knowledge about cellular functions and regulations, thereby elucidating the mechanisms of human health and disease. As most biological or disease processes involve cells changing over time, a major need is to identify the lineage or state-transition of individual cells in single-cell analysis (Kester & van Oudenaarden, 2018). Using single-cell genomics (e.g. RNA-seq, proteomics), prior work demonstrated success based on mathematical inference and were broadly useful in understanding fundamental questions in cell evolution, differentiation, and development (Wolf et al., 2019; Kratochvil et al., 2019; Wang et al., 2018; Setty et al., 2018; Mukherjee et al., 2017; Welch et al., 2017; Qiu et al., 2017; Rizvi et al., 2017; Rashid et al., 2017; Ji & Ji, 2016; Shin et al., 2015; Marco et al., 2014). Nonetheless, as benchmarked in a recent meta-analysis on single-cell trajectory inference methods (Saelens et al., 2019), many existing tools have variabilities arising from the fact that, regardless of the underlying mathematics and statistical methodologies employed, inferred trajectories were constructed indirectly from genomics data such as gene-expression profiles. Thus, their performances are often sensitive to the qualities, dimensions, and underlying topologies of the original datasets (Saelens et al., 2019; Kester & van Oudenaarden, 2018). Additionally, most mathematical inference methods assume underlying “smoothness” of cell dynamics, yet discontinuous jumps and existence of “cycling” or “reversible” cell state transitions could confound these methods.
Hence, the latest single-cell studies began to employ new experimental methods to achieve genetic barcoding of individual cells for direct measurement of cell lineages (Woodworth et al., 2017; Kester & van Oudenaarden, 2018; Wagner & Klein, 2020). These new single-cell genomics methods and datasets, while more powerful, present significant challenges for computational and data analytical tools that we aim to take a substantial step to address in this proposal. Several recent single-cell studies leveraged DNA barcoding approaches, where each cell is given a static DNA sequence that randomizes across cells so that no two cells would bear the same sequence (Woodworth et al., 2017; Kester & van Oudenaarden, 2018). This sequence could then serve as a “tag” to identify cell lineage recoverable by single-cell sequencing. For example, this concept has been utilized in single-cell analysis of embryonic development (Yao et al., 2017; Wagner et al., 2018), stem cell reprogramming (Biddy et al., 2018), and fate determination in hematopoiesis (Weinreb et al., 2020). While static DNA barcoding is able to identify “clonal identity” across time, the emerging CRISPR barcoding technology can label cells using barcodes that evolve as the cell branches and reconstruct full lineage tree structures. An evolving barcoding system that gives each cell a unique tag will be useful for studying cell dynamics in a biological or disease process at finer resolution (Alemany et al., 2018; Kalhor et al., 2018; Raj et al., 2017; Spanjaard et al., 2018; Chan et al., 2019). These higher-resolution CRISPR barcoding tools let scientists track individual cell lineages more easily than ever. They present new open questions for single-cell data analysis. Hence,we are motivated to analyze individual cells’ gene-expression trajectories made available by these barcoding tools, and delve into the underlying single-cell dynamics. Our goal is to develop a streamlined analysis pipeline, connecting these genetic barcodes and single-cell profiles (e.g. gene-expression), to help uncover the hidden mechanisms of cell evolution and fate determination.
With the explosive amount of data from single-cell genomics studies, one remaining untapped problem is the lack of ability to understand cell transition on the individual level. Even if explicit barcodes are available and can define individual cell transitions, to decode single-cell dynamics from such data, we need to solve a highly complex system identification and state estimation problem whereas the true system dynamics is hidden under limited, noisy, and high-dimensional observations. The complexity of decoding single-cell dynamics is due to three multiplicative factors: (1) the stochastic nature of cellular dynamics, (2) nonlinearity or even discontinuity of single-cell gene expression across time, with possible existence of “cycling” and “reversible” transitions, and (3) high dimensionality of gene expression profiles. Hence, new types of dynamic single-cell datasets call for analytical methods that match such complexity and scale.
Conventional methods for analyzing single-cell data are mostly based on “ensemble” analysis, without barcodes to delineate individual transition paths (Kester & van Oudenaarden, 2018). Our work focuses on data made available by novel genetic barcoding approach that enables tracking evolutionary trajectories across individual cells while simultaneous performing genomics read-out of cell states (Woodworth et al., 2017; Kester & van Oudenaarden, 2018; Wagner & Klein, 2020). In these barcoding experiments, a pool of designed or randomized DNA sequences are synthesized and delivered into cells so that each cell bears a different sequence. These sequences, uniquely tagging individual cells, enable the recovery of cell clonality or lineage using genomics methods such as sequencing (Woodworth et al., 2017; Kester & van Oudenaarden, 2018). Recent new development in genetic barcoding tools hold the promise to enable tracking cells at higher capacity, with finer temporal resolution, and for more complex models (Wagner & Klein, 2020). Taking the latest CRISPR barcoding technology for example, a fixed-length DNA barcode is inserted into individual cells, and the barcode dynamically change as the cell branches by nature of the CRISPR editing process (McKenna et al., 2016; Kalhor et al., 2018; Chan et al., 2019). These CRISPR barcodes could allow the reconstruction of lineage history of cells conveniently with end-point data collection, revealing each cell’s transition paths at scale (Wagner & Klein, 2020). In this case, we have pairs of gene-expression transitions {(X,X ′)}, where each (X,X ′) are the gene expression profiles for a parent cell at one time point and its descent at the next time point. This new type of datasets, particularly datasets from recent CRISPR-based barcoding experiments, calls for continued innovation to invent suitable data methods that we are proposing here.
Our method is related to state representation learning and spectral dimension reduction in reinforcement learning, dynamical systems and scientific computing. In particular, our kernelized state embedding estimator generalizes the notation of diffusion map to reproducing kernel spaces. For molecular dynamics,(Schütte et al., 2011) showed that the leading spectrum of the transition operator contains information on slow dynamics of the system, and it can be used to identify coresets upon which a coarse-grained Markov state model could be built. Please see (Klus et al., 2018) for comprehensive surveys on various spectral dimension methods for dynamical systems. They did not study statistical properties of these methods which motivated our research. In control and reinforcement learning, computers need to figure out state abstraction (e.g, state features, state aggregation) of unknown transition systems in order to learn to control quickly (Singh et al., 1995; Venkatesh & Dahleh, 2001; Sutton et al., 1998; Sarkar et al., 2019). One recent work (Yang & Wang, 2019) showed that the statistical complexity of reinforcement learning scales linearly with respect to dimension of the state representation. (Sun et al., 2019) developed a statistically proven kernelized state learning method for Markov processes. These inspirations from dynamic system and reinforcement learning led us to develop single-cell analysis tools in order to map cell states to low-dimensional embeddings.
B Implementation Details
Algorithm 1 gives the full description of the Cell2state algorithm. It takes as input collections of cell transition pairs between two timepoints t0 < t1, organized into N clones. Each clone Ci = {Ci0, Ci1} consists of two sets, where Ci0 is a set of ancestor cells and Ci1 is a set of descendant cells. The ancestor/descendant cell sets, containing single-cell gene-expression profiles Y , are pairwise disjoint: ∀i 6= j ∈ [N ], Ci0 ∩ C j 0 = ∅, Ci1 ∩ C j 1 = ∅. For a geneexpression profile Y0 as an element in set Ci0, we also denote the vector representation X0 preprocessed from Y0 with X0 ∈ Ci0 for simplification; this notation is applied to all cell sets Ci0, Ci1, i ∈ [N ]. For every X0 ∈ Ci0, X1 ∈ Ci1, we consider them together as a cell-to-cell transition pair (X0, X1). In our experiments, we use the Gaussian kernel function, i.e., K(X,X ′) = exp(−γ||X −X ′||2), which is a positive semidefinite shift-invariant kernel. Steps 1-3 of Algorithm 1 corresponds to solving a regression formulation:
min P RSS = N∑ i=1 ∑ X0∈Ci0 ||Φ(X0)>Σ−1/20 P − 1 |Ci0| ∑ X1∈Ci1 Φ(X1) >Σ −1/2 1 ||2.
In Steps 4-5 of Algorithm 1, we approximate P̂ by Û Ŝ[1,..,r]V̂ and construct the Cell2state map Ψ̂(·).
Algorithm 1 Cell2state Algorithm (cell-to-cell) Input: Collection of cell-transition data {Ci = {Ci0, Ci1}}Ni=1, Rank r, Kernel function K Parameter: λ̃,γ,D,d,λ0,λ1 Notation: Let A0 = ⋃N i=1 C i 0, A1 = ⋃N i=1 C i 1.
1: Preprocess all gene-expression profiles Y ∈ A0 ∪ A1 by a selected function g(·) and normalize: X := ( Σ̃ + λ̃I )−1/2
[g(Y )− µ̃g] ∈ Rd̃, where Σ̃ =∑ Y ∈A0∪A1 [g(Y )− µ̃g] [g(Y )− µ̃g] > and µ̃g = 1|A0∪A1| ∑ Y ∈A0∪A1 g(Y ). // d̃ = 50,
λ̃ = 1 2: Generate random Fourier functions and keep top d PCs out of D random Fourier features:
Φ(·) = [Φ1(·),Φ2(·), ...,Φd(·)]> // D = 2000, d = 500 a: Draw D iid samples w1, ...wD ∈ Rd̃ from Fourier transform of chosen kernel K: p(w) = ∫ e−jw
′δK(δ)d∆. // p(w) gives Normal(0, 2γId̃) when K(δ) = exp(−γ||δ||2);γ = 0.4 b: Draw D iid samples b1, ...bD ∈ R from Uniform(0, 2π); c: Generate a randomized feature map Φ̃(X) = √ 2 D [cos(w > 1 X+b1), ..., cos(w > DX+bD)] >;
d: Run PCA on Φ̃(X) to generate map Φ(·) by top d leading eigenvectors ξ1, ..., ξd: Φ(X) = [ξ>1 Φ̃(X), ..., ξ > d Φ̃(X)]
>. Delete mean. 3: Estimate P̂ = Σ−1/20 [∑N i=1
1 |Ci0| ∑ X0∈Ci0 ∑ X1∈Ci1 Φ(X0)Φ(X1) > ] Σ −1/2 1 ,
where Σ0 = [∑ X0∈A0 Φ(X0)Φ(X0) >]+ λ0I and Σ1 = [∑X1∈A1 Φ(X1)Φ(X1)>]+ λ1I.
// λ0 = 0.1, λ1 = 0 4: SVD Û ŜV̂ = P̂ 5: Compute Cell2state features using top r singular values Ψ̂(x) = (Φ(x)>Σ−1/20 Û Ŝ[1,..,r])
Output: x 7→ Ψ̂(x)
The dataset we used (Biddy et al. (2018)) measured reprogramming activities of mouse embryonic fibroblasts to induced endoderm progenitors, containing 48,515 single-cell geneexpression profiles from day 0 to day 28, out of which 17,803 cells has transition information visualized on tSNE coordinates given by the original paper and colored by time points (Figure 3(a)). In particular, we input 1,997 day-12 ancestor cells and 3,509 day-21 descendent cells with 165,716 transition pairs in-between, colored on the same tSNE coordinates with downsampled transitions in Figure 3(b); these single-cell expression profiles measure over 28,001 genes. We choose the input K to be the composition of the principal component map and Gaussian kernel with rank r = 8. The parameter values producing cell2state features with low inertia, distinct descendants structure, and good prediction performance are preferred. The experiment results demonstrated in Figure 3(c-h) of Section 4 are associated with parameters λ̃ = 0, γ = 50, D = 2000, d = 900, λ0 = 0.1, λ1 = 0.
Data preprocessing. Our raw data is a collection of barcoded single-cell gene-expression profiles, with expression measurements of over 20000 genes. We use a standard single-cell data preprocessing pipeline, Seurat workflow(Butler et al. (2018); Stuart et al. (2019)). Firstly, we keep only genes expressed as non-zero in at least one cell and only cells expressed as non-zero in more than 200 genes. Then we filter out low-quality cells according to commonly used quality control metrics, such as the percentage of mitochondrial gene counts, and further log-transform the counts. Finally, we compare the barcode associated with each cell and arrange the cells into clones of transition pairs.
Predictions and Data splits. When testing the cell2state features of prediction of descendant activities, we use the ridge-regularized linear regressor/classifier from the ‘sklearn’ package. For computing out-of-sample errors, we repeated a random half-half data split 10 times and reported results with 95% confidence intervals.
Marker gene identification. Differential gene-expression analysis is used to find marker genes in each clusters. Marker genes are those genes that have significantly different levels of expressions in the cluster and could distinguish it from others. Specifically, we input the identity file of ancestor cells into R and run DESeq2(Love et al. (2014)) by the ‘FindMarkers’ function of Seurat(Butler et al. (2018); Stuart et al. (2019)) with parameter ‘min.pct’ = 0.25 and ‘logfc.threshold’ = 0.25.
Simulation experiment details. In Section 5, we highlighted a simulation setting where we take gene-expression data from Weinreb et al. (2020) and generate nonlinear transitions. In this section of the appendix, we provide further details as to how we go about our simulation experiments. First, we reduce the raw gene expression data to 50 dimensions via PCA as a preprocessing step. We then compare various embedding methods by providing these 50 PCs as input:
• Linear: Keep top 50 principal components of PCA.
• Cell2state: Run cell2state algorithm to get top 8 features.
• Cell2state KPCA: Compute the exact Gaussian kernel(not RFFs) in cell2state to get top 8 features.
• Neural network (NN): A multilayer perceptron (MLP) with Swish activations (Ramachandran et al., 2017) and cross entropy loss.
With the exception of NN, all other embedding methods were trained with logistic regression. We train NN with the ADAM optimizer (Kingma & Ba, 2015) for 400 epochs and batch size 128.
We also conduct a grid-search over hyperparameters/architecture settings: λ ∈ {1e-4, 1e-3, 1e-2}, layers ∈ {[128, 64], [128, 64, 32], [256, 128], [256, 128, 64]}, initial learning rate ∈ {1e-2, 5e-3, 1e-3, 5e-4, 1e-4}. We use 5-fold cross validation to train all models. Furthermore, in Table 1, we see that there is very little variation of validation accuracy across different architectural settings as long as we choose the best hyperparameters for each architecture.

heading: C Proof of Theorems
text: C.1 Assumptions
We first restate the assumptions needed for proofs of Theorems 1 and 2.
• ‖φ(x)‖2 ≤ C, ∀x; • The samples Xi are generated i.i.d. from some unknown distribution; • In the encoding map analysis (Theorem 2), we assume there exists a lossless encoding map Ω∗ : X 7→ [k], i.e., I(Ω∗(X);Y ) = I(X;Y ). It means that p(·|x) = p(·|x′), ∀x, x′ s.t. Ω∗(x) = Ω∗(x′). Assume w.l.o.g that r = rank(Φ∗) is no less than k. In this case, we may ensure ΠHp(·|x) 6= ΠHp(·|x′) when Ω∗(x) 6= Ω∗(x′). Otherwise we get |ΠHp(·|X )| < k, which implies a contradiction that rank(Φ∗) = rank(PH) < k. Then we are able to recover the block structure Ω∗ via estimating Φ∗ by Ψ̂.
In subsequent proofs, we focus solely on the case that each state generates exactly one descendant. Our proof can be easily extended to the case of a branching process where a state has a random number of independent descendants. When a state x has Nx descendants x′ := {x(1)′ , x(2)′ , . . . , x(Nx)′}, the only difference is to let Φ(x′) = ∑Nx j=1 Φ(x
(j)′) and p(·|x) = E[Nx]p′(·|x) where p′(·|x) is the distribution of a single descendant. In this case, we would still be able to establish concentration bounds for P̂ and the rest of the proof follows similarly.
C.2 Lemmas
First, we provide a matrix concentration bound for the estimated embedding matrix. Lemma 1. Let
P = Σ − 12 0 E[φ(X0)φ(X ′ 0) >]Σ − 12 1 .
For any ε > 0, we have P(‖P − P̂‖ > ε) ≤ 2d exp ( − 3γNε 2
2C2d(3 + ε)
) .
where γ−1 = max{‖Σ−10 ‖.
Proof. We have
P̂ − P = Σ− 1 2
0
[ 1
N N∑ i=1 (φ(Xi)φ(X ′ i) > − E[φ(Xi)φ(X ′i)>]) ] Σ − 12 1 .
Let Mi = Σ − 12 0 φ(Xi)φ(X ′ i) >Σ − 12 1 − E[Σ − 12 0 φ(Xi)φ(X ′ i) >Σ − 12 1 ], we have
‖Mi‖ = ‖Σ − 12 0 φ(Xi)φ(X ′ i) >Σ − 12 1 − E[Σ − 12 0 φ(Xi)φ(X ′ i) >Σ − 12 1 ]‖ ≤ 2C
γ ,
and ‖E[M>i Mi]‖ =‖Σ − 12 1 E[(φ(Xi)φ(X ′ i) > − E[φ(Xi)φ(X ′i)]>)>Σ−10 (φ(Xi)φ(X ′i)> − E[φ(Xi)φ(X ′i)>])Σ − 12 1 ]‖
≤‖E[Σ− 1 2 1 (φ(Xi)φ(X ′ i) >)>Σ−10 φ(Xi)φ(X ′ i) >Σ − 12 1 ]‖
=‖E[[φ(Xi)>Σ−10 φ(Xi)]Σ − 12 1 φ(X ′ i)φ(X ′ i) >Σ − 12 1 ]‖
≤C γ ,
‖E[MiM>i ]‖ =‖E[Σ − 12 0 (φ(Xi)φ(X ′ i) > − E[φ(Xi)φ(X ′i)>])Σ−11 (φ(Xi)φ(X ′i)> − E[φ(Xi)φ(X ′i)>])>Σ − 12 0 ]‖
≤‖E[Σ− 1 2 0 φ(Xi)φ(X ′ i) >Σ−11 (φ(Xi)φ(X ′ i) >)>Σ − 12 0 ]‖
=‖E[φ(X ′i)>Σ−11 φ(X ′i)Σ − 12 0 φ(X ′ i)φ(Xi) >Σ − 12 0 ]‖
≤C γ .
By matrix Bernstein’s inequality, we have P(‖P − P̂‖ > ε) ≤ 2d exp ( − 3γNε 2
2C(3 + ε)
) .
Next we show that the estimated embedding map preserves the correct feature space with high probability.
Lemma 2. Let Ψ̂ be the estimated mapping, we have
P( inf O:O>O=I sup x∈X ‖OΨ̂(x)−Ψ(x)‖ > ε) ≤ 2d exp
( − 3γ 2Nε2
8C 3 2κ(12 √ Cκ+ √ γε)
) ,
where κ = ‖P‖σr(P ) .
Proof. Suppose the SVD of P is P = UrΛrV >r . Let P̃ = ÛrΛ̂rV̂ >r . Let Ψ(x) = (UrΛr) >Σ − 12 0 φ(x). Let O be an arbitrary orthonormal matrix. For any x ∈ X , we have
‖Ψ(x)−OΨ̂(x)‖ = ‖(UrΛr)>Σ − 12 0 φ(x)−O(ÛrΛ̂r)>Σ − 12 0 φ(x)‖
≤ √ C √ γ ‖(UrΛr)> −O(ÛrΛ̂r)>‖
= √ C √ γ ‖(PVr)> −O(P̃ V̂r)>‖
= √ C √ γ ‖(PVr)> −O(P ′V̂r)> + O(PV̂r)> −O(P̃ V̂r)>‖ ≤ √ C √ γ ( ‖P‖‖Vr − V̂rO>‖+ ‖P − P̃‖ ) .
Taking infinum over O, we get
inf O ‖Ψ(x)−OΨ̂(x)‖ ≤ inf O √ C √ γ ( ‖P‖‖Vr − V̂rO>‖+ ‖P − P̃‖ ) ≤ √ C √ γ (√ 2‖P‖‖ sin Θ(Vr, V̂r)‖+ ‖P − P̃‖
) ≤ √ C √ γ ( √ 2‖P‖‖P − P̂‖ σr(P ) + ‖P − P̃‖ ) .
By Weyl’s inequality, we have
‖P − P̃‖ ≤ ‖P − P̂‖+ ‖P̂ − P̃‖ = ‖P − P̂‖+ σr+1(P ) ≤ 2‖P − P̂‖.
In conclusion, we get
inf O ‖Ψ(x)−OΨ̂(x)‖ ≤ √ C √ γ (√ 2‖P‖ σr(P ) + 2 ) ‖P − P̂‖,∀x,
which implies
inf O sup x∈X ‖Ψ(x)−OΨ̂(x)‖ ≤ 2 √ C √ γ ( ‖P‖ σr(P ) + 1 ) ‖P − P̂‖.
According to the result of Lemma 1, we get P (
inf O sup x∈X ‖Ψ(x)−OΨ̂(x)‖ > ε
) ≤ P ( ‖P − P̂‖ >
√ γε
2 √ C(1 + κ) ) ≤ 2d exp ( − 3γ 2Nε2
4C 3 2 (1 + κ)(6 √ C(1 + κ) + √ γε)
)
≤ 2d exp ( − 3γ 2Nε2
8C 3 2κ(12 √ Cκ+ √ γε)
) .
C.3 Proof of Theorem 1
Proof. For any orthonormal matrix O, we have
‖Ψ(x)−Ψ(y)‖ = ‖Ψ(x)−OΨ̂(x) + OΨ̂(x)−OΨ̂(y) + OΨ̂(y)− Ψ̂(y)‖ ≤ ‖Ψ(x)−OΨ̂(x)‖+ ‖OΨ̂(x)−OΨ̂(y)‖+ ‖OΨ̂(y)− Ψ̂(y)‖ = ‖Ψ(x)−OΨ̂(x)‖+ ‖Ψ̂(x)− Ψ̂(y)‖+ ‖OΨ̂(y)− Ψ̂(y)‖.
Similarly we have
‖Ψ̂(x)− Ψ̂(y)‖ ≤ ‖Ψ(x)−OΨ̂(x)‖+ ‖Ψ(x)−Ψ(y)‖+ ‖OΨ̂(y)− Ψ̂(y)‖.
Therefore, we have
|‖Ψ(x)−Ψ(y)‖ − ‖Ψ̂(x)− Ψ̂(y)‖| ≤ 2 sup x∈X ‖Ψ(x)−OΨ̂(x)‖.
Taking infimum over O, and notice that according to the result of Lemma 2, we have w.p. 1− q,
inf O sup x∈X ‖Ψ(x)−OΨ̂(x)‖ ≤ 8Cκ γ √ log 2dq N + 16C 3 2κ γ 3 2 log 2dq N .
We get w.p. 1− q,
|‖Ψ(x)−Ψ(y)‖ − ‖Ψ̂(x)− Ψ̂(y)‖| ≤ 16Cκ γ √ log 2dq N + 32C 3 2κ γ 3 2 log 2dq N .
Now, define inner products 〈f, g〉H = [∫
y
f(y)Σ − 12 1 φ(y)dy ]> [∫ y g(y)Σ − 12 1 φ(y)dy ] .
we have ∫ y Σ − 12 1 φ(y)[Σ 1 2 1 C −1φ(y)]>dy = ∫ y Σ − 12 1 φ(y)φ(y) >C−1Σ 1 2 1 dy = I.
Therefore, e(y) = Σ 1 2 1 C −1φ(y) is a set of orthonormal basis w.r.t. 〈·, ·〉H. And for any fixed x, we have
ΠH(p(·|x)) = [∫ p(y|x)(Σ− 1 2 1 φ(y)) >dy ] Σ 1 2 1 C −1φ(·).
According to assumption 1, we know for any function f ∈ H, we have
P f(x) = ∫ p(y|x)f(y)dy ∈ H.
Therefore, we know that each entry of ∫ p(y|x)(Σ− 1 2
1 φ(y)) >dy is in H. Denote∫
p(y|x)(Σ− 1 2
1 φ(y)) >dy = φ(x)>V , we have
V = Σ−10
∫ π(x)φ(x)p(y|x)(Σ− 1 2
1 φ(y)) >dxdy.
Therefore, we get ΠH(p(·|x)) = (Σ−10 Φ(x))> [∫ π(x)Φ(x)p(y|x)(Σ− 1 2 1 Φ(y)) >dxdy ] Σ 1 2 1 C −1Φ(·)
= (Σ−10 Φ(x)) > [ Σ − 12 0 P ] Σ 1 2 1 C −1Φ(·)
= Ψ(x)>Σ 1 2 1 C −1φ(·)
= Ψ(x)>e(·). Therefore,
‖p(·|x)− p(·|y)‖H = ‖Ψ(x)−Ψ(y)‖, which implies
|‖p(·|x)− p(·|y)‖H − ‖Ψ̂(x)− Ψ̂(y)‖| ≤ 16Cκ
γ √ log 2dq N + 32C 3 2κ γ 3 2 log 2dq N .
C.4 Proof of Theorem 2
Proof. According to the proof of Theorem 1, we know ‖p(·|x)− p(·|y)‖H = ‖Ψ(x)−Ψ(y)‖.
Therefore, when Ω∗(x) = Ω∗(x′), we must have 0 = ‖p(·|x)− p(·|x′)‖H = ‖Ψ(x)−Ψ(x′)‖,
i.e., Ψ(x) = Ψ(x′). Therefore, Ψ is the linear mapping such that |Ψ(X )| ≤ k. For any given ψ̂, define Oψ̂ as an orthonormal matrix such that
sup x∈X ‖Oψ̂ψ̂(x)−Ψ(x)‖ ≤ inf O:O>O=I sup x∈X ‖Oψ̂(x)−Ψ(x)‖+ δ 8 .
According to Lemma 2, we have
P(∃x, ‖OΨ̂Ψ̂(x)−Ψ(x)‖ > δ/4) ≤ P( inf O:O>O=I sup x∈X ‖OΨ̂(x)−Ψ(x)‖ > δ 8 )
≤ 2d exp ( − 3γ 2Nδ2
8× 64C 32κ(12 √ Cκ+ √ γδ/8)
) .
According to our setting, for any x1, x2 such that Ψ(x1) = Ψ(x2) = Ψi, we have p(·|x1) = p(·|x2) = pi(·).
Therefore, by definition, we have
I(X,X ′) = ∫ p(y|x)π(x) log p(y|x)
pY (y) dxdy
= ∑ i ∫ p(y|x)π(x)1Ψ(x)=Ψi log p(y|x) pY (y) dxdy
= ∑ i ∫ pi(y)P(Ψ(X) = Ψi) log pi(y) pY (y) dy
= ∑ i ∑ j ∫ pi(y)P(Ψ(X) = Ψi,Ωψ̂(X) = Aj) log pi(y) pY (y) dy.
On the other hand, we have
I(Ωψ̂(X), X ′) = ∑ j ∫ p(y|x)π(x)1Ωψ̂(x)=Aj log p(y|Ωψ̂(x) = Aj) pY (y) dxdy
= ∑ i ∑ j ∫ p(y|x)π(x)1Ψ(x)=Ψi,Ωψ̂(x)=Aj log p(y|Ωψ̂(x) = Aj) pY (y) dxdy
= ∑ i ∑ j ∫ p(y|x)π(x)1Ψ(x)=Ψi,Ωψ̂(x)=Aj log ∫ p(y|x)π(x)1Ωψ̂(x)=Ajdx pY (y)P(Ωψ̂(X) = Aj) dxdy
= ∑ i ∑ j ∫ p(y|x)π(x)1Ψ(x)=Ψi,Ωψ̂(x)=Aj log ∫ p(y|x)π(x)1Ωψ̂(x)=Ajdx∫ p(y|x)π(x)1Ψ(x)=Ψi,Ωψ̂(x)=Ajdx dxdy
+ ∑ i ∑ j ∫ p(y|x)π(x)1Ψ(x)=Ψi,Ωψ̂(x)=Aj log ∫ p(y|x)π(x)1Ψ(x)=Ψi,Ωψ̂(x)=Ajdx pY (y)P(Ωψ̂(X) = Aj) dxdy
= ∑ j ∫ p(y|x)π(x)1Ωψ̂(x)=Aj ∑ i P(Ψ(X) = Ψi|Y = y, Ω̂ψ̂(X) = Aj) log 1 P(Ψ(X) = Ψi|Y = y, Ω̂ψ̂(X) = Aj) dxdy
+ ∑ i ∑ j ∫ p(y|x)π(x)1Ψ(x)=Ψi,Ωψ̂(x)=Aj log ∫ p(y|x)π(x)1Ψ(x)=Ψi,Ωψ̂(x)=Ajdx pY (y)P(Ωψ̂(X) = Aj) dxdy
≥ ∑ i ∑ j ∫ p(y|x)π(x)1Ψ(x)=Ψi,Ωψ̂(x)=Aj log ∫ p(y|x)π(x)1Ψ(x)=Ψi,Ωψ̂(x)=Ajdx pY (y)P(Ωψ̂(X) = Aj) dxdy.
And we have∑ i ∑ j ∫ p(y|x)π(x)1Ψ(x)=Ψi,Ωψ̂(x)=Aj log ∫ p(y|x)π(x)1Ψ(x)=Ψi,Ωψ̂(x)=Ajdx pY (y)P(Ωψ̂(X) = Aj) dxdy
= ∑ i ∑ j ∫ pi(y)π(x)1Ψ(x)=Ψi,Ωψ̂(x)=Aj log ∫ pi(y)π(x)1Ψ(x)=Ψi,Ωψ̂(x)=Ajdx pY (y)P(Ωψ̂(X) = Aj) dxdy
= ∑ i ∑ j ∫ pi(y)P(Ψ(x) = Ψi,Ωψ̂(x) = Aj) log pi(y)P(Ψ(x) = Ψi,Ωψ̂(x) = Aj) pY (y)P(Ωψ̂(X) = Aj) dy.
which implies
I(Ωψ̂(X), X ′) ≥ ∑ i ∑ j ∫ pi(y)P(Ψ(x) = Ψi,Ωψ̂(x) = Aj) log pi(y)P(Ψ(x) = Ψi,Ωψ̂(x) = Aj) pY (y)P(Ωψ̂(X) = Aj) dy.
Combining the results so far, we get
I(X,X ′)− I(Ω̂ψ̂(X), X ′)
≤ ∑ i ∑ j ∫ pi(y)P(Ψ(X) = Ψi,Ωψ̂(X) = Aj) log pi(y) pY (y) dy
− ∑ i ∑ j ∫ pi(y)P(Ψ(x) = Ψi,Ωψ̂(x) = Aj) log pi(y)P(Ψ(x) = Ψi,Ωψ̂(x) = Aj) pY (y)P(Ωψ̂(X) = Aj) dy
= ∑ i ∑ j ∫ pi(y)P(Ψ(X) = Ψi,Ωψ̂(X) = Aj) log P(Ωψ̂(X) = Aj) P(Ψ(x) = Ψi,Ωψ̂(x) = Aj) dy = ∑ i ∑ j P(Ψ(X) = Ψi, Ω̂ψ̂(X) = Aj) log P(Ω̂ψ̂(X) = Aj) P(Ψ(X) = Ψi, Ω̂ψ̂(X) = Aj) .
Taking expectation on both sides, we get
EΨ̂[I(Y, Ω̂Ψ̂(X))] ≥I(Y,X) − ∫ ∑
i ∑ j P(Ψ(X) = Ψi, Ω̂ψ̂(X) = Aj) log P(Ω̂ψ̂(X) = Aj) P(Ψ(X) = Ψi, Ω̂ψ̂(X) = Aj)  dP(ψ̂). Note that∫ ∑
i ∑ j P(Ψ(X) = Ψi, Ω̂ψ̂(X) = Aj) log P(Ω̂ψ̂(X) = Aj) P(Ψ(X) = Ψi, Ω̂ψ̂(X) = Aj) dP(ψ̂)  =
∫ ψ̂:∀x,‖Oψ̂ψ̂(x)−Ψ(x)‖≤ δ 4 ∑ i ∑ j P(Ψ(X) = Ψi, Ω̂ψ̂(X) = Aj) log P(Ω̂ψ̂(X) = Aj) P(Ψ(X) = Ψi, Ω̂ψ̂(X) = Aj)  dP(ψ̂) +
∫ ψ̂:∃x,‖Oψ̂ψ̂(x)−Ψ(x)‖> δ 4 ∑ i ∑ j P(Ψ(X) = Ψi, Ω̂ψ̂(X) = Aj) log P(Ω̂ψ̂(X) = Aj) P(Ψ(X) = Ψi, Ω̂ψ̂(X) = Aj)  dP(ψ̂) ≤ ∫ ψ̂:∀x,‖Oψ̂ψ̂(x)−Ψ(x)‖≤ δ 4 ∑ i ∑ j P(Ψ(X) = Ψi, Ω̂ψ̂(X) = Aj) log P(Ω̂ψ̂(X) = Aj) P(Ψ(X) = Ψi, Ω̂ψ̂(X) = Aj)
 dP(ψ̂) +
∫ ψ̂:∃x,‖Oψ̂ψ̂(x)−Ψ(x)‖> δ 4 log k dP(ψ̂)
≤ ∫ ψ̂:∀x,‖Oψ̂ψ̂(x)−Ψ(x)‖≤ δ 4 ∑ i ∑ j P(Ψ(X) = Ψi, Ω̂ψ̂(X) = Aj) log P(Ω̂ψ̂(X) = Aj) P(Ψ(X) = Ψi, Ω̂ψ̂(X) = Aj)  dP(ψ̂) + 2d exp ( − 3γ 2Nδ2
8× 64C 32κ(12 √ Cκ+ √ γδ/8)
) log k.
It remains to consider the first term. Note that∑ i ∑ j P(Ψ(X) = Ψi, Ω̂ψ̂(X) = Aj) log P(Ω̂ψ̂(X) = Aj) P(Ψ(X) = Ψi, Ω̂ψ̂(X) = Aj) = ∑ i ∑ j:∃z∈Aj ,‖Oψ̂z−Ψi‖≤ δ 4 P(Ψ(X) = Ψi, Ω̂ψ̂(X) = Aj) log P(Ω̂ψ̂(X) = Aj) P(Ψ(X) = Ψi, Ω̂ψ̂(X) = Aj)
+ ∑ i ∑ j:∀z∈Aj ,‖Oψ̂z−Ψi‖> δ 4 P(Ψ(X) = Ψi, Ω̂ψ̂(X) = Aj) log P(Ω̂ψ̂(X) = Aj) P(Ψ(X) = Ψi, Ω̂ψ̂(X) = Aj) .
Therefore, when ψ̂ satisfies ∀x, ‖Oψ̂ψ̂(x) − Ψ(x)‖ ≤ δ 4 , then for any j such that ∀z ∈ Aj , ‖Oψ̂z −Ψi‖ > δ 4 , we have
P(Ψ(X) = Ψi, Ω̂ψ̂(X) = Aj) = P(Ψ(X) = Ψi, ψ̂(X) ∈ Aj) ≤ P(‖Ψ(X)−Oψ̂ψ̂(X)‖ > δ
4 ) = 0.
And for any j such that ∃z ∈ Aj , ‖Oψ̂z − Ψi‖ ≤ δ 4 and any l 6= i, since we have δ ≤ minl 6=i ‖Ψi −Ψl‖, we get ∀z′ ∈ Aj , ‖Ψl −Oψ̂z ′‖ = ‖Ψl −Ψi + Ψi −Oψ̂z ′‖ ≥ ‖Ψl −Ψi‖ − ‖Ψi −Oψ̂z
′‖ ≥ δ − ‖Ψi −Oψ̂z + Oψ̂z −Oψ̂z
′‖ ≥ δ − ‖Ψi −Oψ̂z‖ − ‖Oψ̂z −Oψ̂z ′‖
≥ δ − δ 4 − δ 4 = δ 2 .
Therefore,
P(Ψ(X) = Ψl, Ω̂ψ̂(X) = Aj) = P(Ψ(X) = Ψl, ψ̂(X) ∈ Aj) ≤ P(‖Ψ(X)−Oψ̂ψ̂(X)‖ > δ
2 ) = 0,∀l 6= i,
which implies
P(Ψ(X) = Ψi, Ω̂ψ̂(X) = Aj) = P(Ω̂ψ̂(X) = Aj)− ∑ l 6=i P(Ψ(X) = Ψl, Ω̂ψ̂(X) = Aj)
= P(Ω̂ψ̂(X) = Aj).
Combining the results so far, we get
∑ i ∑ j P(Ψ(X) = Ψi, Ω̂ψ̂(X) = Aj) log P(Ω̂ψ̂(X) = Aj) P(Ψ(X) = Ψi, Ω̂ψ̂(X) = Aj)
= ∑ i ∑ j:∃z∈Aj ,‖Oψ̂z−Ψi‖≤ δ 4 P(Ψ(X) = Ψi, Ω̂ψ̂(X) = Aj) log P(Ω̂ψ̂(X) = Aj) P(Ψ(X) = Ψi, Ω̂ψ̂(X) = Aj)
+ ∑ i ∑ j:∀z∈Aj ,‖Oψ̂z−Ψi‖> δ 4 P(Ψ(X) = Ψi, Ω̂ψ̂(X) = Aj) log P(Ω̂ψ̂(X) = Aj) P(Ψ(X) = Ψi, Ω̂ψ̂(X) = Aj)
= ∑ i ∑ j:∃z∈Aj ,‖Oψ̂z−Ψi‖≤ δ 4 P(Ω̂ψ̂(X) = Aj) log P(Ω̂ψ̂(X) = Aj) P(Ω̂ψ̂(X) = Aj) + ∑ i ∑ j:∀z∈Aj ,‖Oψ̂z−Ψi‖> δ 4 0
=0.
which implies
∫ ∑ i ∑ j P(Ψ(X) = Ψi, Ω̂ψ̂(X) = Aj) log P(Ω̂ψ̂(X) = Aj) P(Ψ(X) = Ψi, Ω̂ψ̂(X) = Aj) dP(ψ̂)  ≤2d exp ( − 3γ 2Nδ2
8× 64C 32κ(12 √ Cκ+ √ γδ/8)
) log k,
i.e.,
EΨ̂[I(Ω̂Ψ̂(X), X ′)] ≥ I(X,X ′)− 2d exp
( − 3γ 2Nδ2
8× 64C 32κ(12 √ Cκ+ √ γδ/8)
) log k.
Finally, noticing that
δ ≤ 2 sup x ‖Ψ(x)‖ ≤ 2
√ C‖P‖γ− 12 ≤ 2 √ Cκγ− 1 2 ,
we have finished the proof.

Review summaries:
The authors introduce cell2state, an algorithm that incorporates both genetic barcoding coupled with single-cell sequenced data to model explicit state transitions of cell dynamics over time. Single-cell gene expression profiles are mapped to low-dimensional state vectors that are predictive of cell dynamics. Cell2state is evaluated using barcoded stem cell dataset (Biddy et al. (2018)) and simulation studies are also presented. The model demonstrates better results for cell state prediction, finding dynamically stable clusters, and reveals potential latent meta-states of the underlying cellular evolution process.  

Authors proposed cell2state that could embed barcoded scRNA-seq trajectories into low-dimensional representation. Authors provided theoretic analysis of the embedding learnt by cell2state and demonstrated that the learnt embedding was almost lossless.  Authors applied this embedding framework on one barcoded scRNA-seq dataset (Biddy, et al., 2018) and demonstrated the learnt embeddings clearly distinguished different cell states. Furthermore, the learnt embeddings were able to substantially improve various downstream tasks other than identifying cell subpopulation.