The analysis of RNA data plays a crucial role in understanding cellular differentiation.
One widely-used methodology for analyzing RNA data is scVelo. However, in this paper, we show that, among other issues of scVelo, the current model formalization suffers from identifiability problems. We propose a Bayesian version of scVelo with modifications that address these issues

Bibbona,Enrico

Mastrantonio Gianluca

Sabbioni,Elena

Sanguinetti Guido

English

PORTO@iris (Publications Open Repository TOrino - Politecnico di Torino)

11 February 2024POLITECNICO DI TORINORepository ISTITUZIONALEAnalyzing RNA data with scVelo: identifiability issues and a Bayesian implementation / Sabbioni, Elena; Bibbona, Enrico;Mastrantonio, Gianluca; Sanguinetti, Guido. - ELETTRONICO. - (2023), pp. 538-543. (Intervento presentato alconvegno SIS 2023 - Statistical Learning, Sustainability and Impact Evaluation tenutosi a Ancona (ITA) nel 21/06/2023-23/06/2023).OriginalAnalyzing RNA data with scVelo: identifiability issues and a Bayesian implementationPublisher:PublishedDOI:Terms of use:Publisher copyright(Article begins on next page)This article is made available under terms and conditions as specified in the  corresponding bibliographic description inthe repositoryAvailability:This version is available at: 11583/2982276 since: 2023-10-04T09:25:25ZPearsonAnalyzing RNA data with scVelo: identifiabilityissues and a Bayesian implementationElena Sabbionia, Enrico Bibbonaa, Gianluca Mastrantonioa, and GuidoSanguinettibaPolitecnico di Torino; elena.sabbioni@polito.it, enrico.bibbona@polito.it,gianluca.mastrantonio@polito.it,bScuola Internazionale Superiore di Studi Avanzati (SISSA); gsanguin@sissa.itAbstractThe analysis of RNA data plays a crucial role in understanding cellular differentiation.One widely-used methodology for analyzing RNA data is scVelo. However, in this paper,we show that, among other issues of scVelo, the current model formalization suffers fromidentifiability problems. We propose a Bayesian version of scVelo with modifications thataddress these issues.Keywords: scVelo, RNA, Bayesian, identifiability1. IntroductionRNA velocity is a critical biological metric that facilitates the reconstruction of cellular differentia-tion at single-cell level. It provides insight into the future state of each cell, and it is closely associatedwith transcription from DNA to RNA, as well as the quantity of spliced mRNA in each cell. By ana-lyzing RNA velocity, researchers can gain a deeper understanding of the underlying mechanisms drivingcellular differentiation, which has important implications for fields such as developmental biology anddisease research. Single-cell RNA sequencing (scRNA-seq) techniques are commonly used to measurethe abundance of unspliced and spliced mRNA in each cell for each gene, which is essential for inferringRNA velocity. However, these techniques are destructive, as they permit only a single observation ofgene expression for each cell before it is destroyed. In this sector, one of the most influential worksis scVelo, presented in [1]. Despite its success in the scientific community, there are several criticismswhen it is analyzed from a mathematical and statistical point of view.As a primary contribution, we reframe the model under a Bayesian framework, which provides betterinsight into the parameters that can be estimated and identified. The use of Bayesian inference enablesus to compare the posterior estimates of the parameters with their corresponding priors and computecredible intervals. In contrast, such comparisons and interval estimates are not possible with the pointestimates provided by scVelo. Through simulated examples, we demonstrate that the “time” parameter isnot identifiable, which, to the best of our knowledge, has not been previously identified in the literature,representing another contribution of this paper. Furthermore, we propose modifications to the modelthat addresses some of the criticisms of the original scVelo. Collectively, these contributions have thepotential to enhance the accuracy and reliability of RNA velocity inference, as well as provide newinsights into cellular differentiation.11.5 2.0 2.5 3.0 3.5 4.0 4.50.51.01.52.02.5suSSgoffSSgon(s0goff , u0goff)InductionRepressionFigure 1: Solution of (2) in the space (s, u) for gene g. The solid line represents the gene’sbehavior if an early switch at time toff0g occurs, while the dashed line represents the potentialbehavior if the gene reaches the steady state SSong during the inductive phase. The blue uppercurve corresponds to the gene’s inductive phase, characterized by the rate αong , and the red lowercurve represents the repressive phase, associated with αoffg .2. Mathematical modelLet us consider a system with ng genes and nc cells. The model assumes a straightforward chemicalreaction network (CRN) to represent the processes of transcription, splicing and degradation, with gene-specific rates, according to mass-action kinetics. The CRN representation of the process, for a given g,isDNAαkg→ unspliced mRNA βg→ spliced mRNA γg→ ∅ (1)that is associated with the following ordinary differential equation (ODE) system:dug(t)dt = αkg − βgug(t),dsg(t)dt = βgug(t)− γgsg(t),ug(tk0g)= uk0g,sg(tk0g)= sk0g,(2)where ug(t) and sg(t) are the reads of unspliced and spliced mRNA at time t in a cell. The solution tothe ODEs isug(t) = uk0ge−βgτkg +αkgβg(1− e−βgτkg )sg(t) = sk0ge−γgτkg +αkgγg(1− e−γgτkg ) + αkg−βguk0gγg−βg (e−γgτkg − e−βgτkg ) with τkg = t− tk0g.(3)It should be noted that the variable t is not the real time, but a representation of the cell position in theODE dynamic, which is often called pseudotime in the literature, see, for example, [3].RNA velocity is defined asvg(t) :=dsg(t)dt= βgug(t)− γgsg(t).Accurately estimating the model parameters is crucial for obtaining a reliable estimator of this biologicalquantity.Parameters description For each gene, there exist two transcription rates, indicated as αong and αoffgwith αong > αoffg , represented in (2) as αkg , with k ∈ {on, off}. The rates regulate the conversion of DNAinto unspliced mRNA, as depicted in (1) (first and second components). This implies that a gene can existin two distinct states: an inductive phase, regulated by the transcription rate αong , and a repressive phase,where transcription either occurs at a lower rate or is absent altogether, dictated by αoffg . It is assumedthat each gene can be activated and then repressed only once, which is justified by the assumption that thetotal time length of the biological processes is sufficiently small. The rates βg and γg, illustrated in (1)(from the second to the fourth component), are responsible for the splicing and degradation mechanismsof mRNA. The gene time dynamic is depicted in Figure 1.The ODE system has two theoretical steady states, that are only gene-dependent and are identifiedby the coordinatesSSoffg =(αoffgβg,αoffgγg), SSong =(αongβg,αongγg),in the space (s, u), see Figure 1. After the cell is created, each gene remains at SSoffg for a time period ofton0g before being activated and entering the inductive phase, represented by the upper blue arc in Figure1. Time ton0g is not identifiable (see Section 3.) since there is no information in the data regarding the realtime point at which the cells are observed, and hence, without loss of generality, we assume ton0g = 0.This means that son0g = αoffg /βg and uon0g = αoffg /γg. However, before reaching the second steady stateSSong , the repressive phase is triggered at time toff0g + ton0g and the dynamic follows the evolution depictedby the solid line in Figure 1.In scVelo the pre-processed data (Yu,cg, Ys,cg)′ are assumed to be normally distributed and the un-spliced and spliced components to be independent, i.e.Yu,cg ∼ N (ug(tcg), σ2) Ys,cg ∼ N (sg(tcg), σ2) Yu,cg ⊥⊥ Ys,cg. (4)where tcg = τcg+ tkcg0g and ug(tcg) and sg(tcg) are evaluated in a time that is both cell- and gene-specific.The description of the data used to estimate the model is discussed in Section 3. The scVelo algorithmestimates the following parameters: (αoffg , αong , βg, γg, toff0g ) for each gene, and (τcg, kcg) for each cell andgene.There are several criticisms of this model, that, in our opinion, raise questions about the reliability ofthe results, which will be discussed in the next section.3. Critical issues of scVeloOne of the main concerns is related to the estimation of the cell- and gene-specific τcg. Since single-cell data only provides a single observation for each cell, inferring τcg is inherently difficult, if notimpossible. Despite this, the authors did not acknowledge this issue. As a first contribution, we demon-strate in Section 4.that τcg is at best weakly identifiable by showing that, the posterior distribution of τcgclosely resembles the one of τc′g, for c ̸= c′, and they are both very similar to the prior.Single-cell RNA-sequencing dataset contains discrete counts, describing the number of measuredRNA molecules in each cell. In scVelo, a series of pre-processing steps are applied to the raw data. Thisincludes filtering out genes that are not expressed, normalizing the counts to account for differences in se-quencing depth across cells, and smoothing the gene expression profiles among groups of cells with sim-ilar genetic expressions. Additionally, the logarithm of the pre-processed counts is taken. The variables(Yu,cg, Ys,cg)′ in equation (4), are the results of this pre-processing. While these pre-processing stepsare common in many biological pipelines, they significantly alter the nature of the data by transform-ing them from discrete counts to continuous values. This transformation can be problematic, especiallyin real data applications where the original counts are often very low (often in the range [0,10]). Thepipeline introduces dependence across genes which are not taken into account in the model, that assumesindependence, see (2) and (4). Additionally, the use of a logarithmic transformation is questionable sinceODE equations and solutions are not invariant under a non-linear transformation.Despite time-dynamic being dependent on four parameters (αoffg , αong , βg, γg)′, only three of themare identifiable due to the lack of information on tcg and its scale in the data. We can easily see that,if r ∈ R+, then the parameters (αoffg , αong , βg, γg)′ and (αoffg /r, αong /r, βg/r, γg/r)′ produce the samelikelihood if τcg is substituted with rτcg. This is because under both sets of parameters, the same value(ug(tcg), sg(tcg)) is obtained. While some of the issues discussed here have been previously addressedin the literature (e.g., [2; 5]), the non-identifiability of τcg and its impact on other parameter estimateshas not been adequately emphasized to the scientific community.In conclusion, scVelo has a further drawback in that it only provides point estimates of the parametersand does not compute any measure of their precision. This absence reduces the reliability of the resultsas it is not possible to assess the statistical differences among the parameters accurately.4. The Bayesian ImplementationIn our model formulation, we choose to use the original data without the non-linear pre-processingsteps applied in scVelo. The only step we keep is the filtration of the genes that are not sufficientlyexpressed. As a result, our data (Yu,cg, Ys,cg)′ is discrete, and (ug(tcg), sg(tcg))′, obtained as solution of(2), represents the mean of the original count data. A natural choice for modeling (Yu,cg, Ys,cg)′ is thePoisson, because this distribution arises from the chemical master equation (CME) associated with theCRN (1). Specifically, in the steady state, the Poisson distribution is the distribution of mRNA counts of asingle gene in a single cell, and in the transient part, CME distribution can be expressed as the convolutionof multinomial and product Poisson distributions [4]. On the other hand, to increase the model flexibilityand to take into account extra sources of variability we use Negative Binomial distribution, which is anoverdispersed version of the Poisson. Specifically, we assume thatYu,cg ∼ NB(ug(tcg), ηg) Ys,cg ∼ NB(sg(tcg), ηg) Yu,cg ⊥⊥ Ys,cgHere the Negative-Binomial is parameterized in terms of its mean µ and the overdispersion parameter η,such that if X ∼ NB(µ, η), then V(X) = µ(1+µη). As prior distributions we define αoffg , αong , γg, ηg ∼N[0,+∞)(0, 10000), where N[a,b] is a truncated Normal distribution with support in [a, b], and P (kcg =on) = 0.5. For τcg we define a mixed-type distribution with two masses of value 0.1 on τcg = 0 andτcg = ∞, respectively, and with probability 0.8 we have log τcg ∼ N(0, 100). The two masses are usedto locate the cell in the steady states. The prior on toff0g must depend on the set {τcg, kcg}ncc=1 sincetoff0g ≥ max{τcg|kcg = on, c = 1, . . . , nc} = τ ong,max (5)hence, we assume the following: log toff0g |{τcg, kcg}ncc=1 ∼ N[log τ ong,max,∞)(0, 100) To avoid identifiabilityissue, parameter βg is fixed to 1.Simulation setting We simulate data with ng = 5, nc = 3600, and γg randomly generated in[0.5, 0.8], αoffg in [0.05, 1], αong in [2, 5], and ηg in [0.01, 0.1]. These intervals have been chosen such thatthe empirical distribution of the raw data mimics the one of the real pancreatic dataset used in [1]. Thereare different issues when simulating parameters τcg, toff0g , and kcg. Indeed, we have to satisfy equation (5)and, to have realistic and diversified locations of (sg(tcg), ug(tcg)), as well as (sg(toff0g ), ug(toff0g )) close toSSong , SSoffg or in between. To achieve this, a sequence of if/else conditions were implemented, howeverfor the sake of brevity these details are omitted. We run the model for 100000 iterations, with thin 40and burning 10000, having then 2250 posterior samples.Discussion of the results The posterior distributions of different τcg, as shown in Figure 2, revealthat there are minimal differences between them. Additionally, the posterior distribution is similar to theprior distribution. Comparable results are obtained when changing the prior distribution, which are notpresented for the sake of brevity. It should be pointed out that, under this setting, for log τcg < −5 andlog τcg > 5 the coordinates (sg(tcg), sg(tcg)), for all c = 1, . . . , nc and g = 1, . . . , ng, are approximatelyequal to the steady states with a difference of order 10−3.−20 −10 0 10 200.000.010.020.030.040.05DensityP(τcg = 0) = 0.12P(τcg = ∞) = 0.07(a)−20 −10 0 10 200.000.010.020.030.040.05Cell 90DensityP(τcg = 0) = 0.14P(τcg = ∞) = 0.05(b)−20 −10 0 10 200.000.010.020.030.040.05DensityP(τcg = 0) = 0.14P(τcg = ∞) = 0.07(c)Figure 2: Gene- and cell-specific τcg model. Prior (red dashed line) and posterior (black solidline) of the logarithm of τcg for three cells. The vertical line represents the true value.0.5 1.0 1.5 2.0 2.5 3.0 3.50.51.01.52.02.5suRealEstimatedFigure 3: Phase plot in the space (s, u) for a given gene. The solid lines show the solutions(3) obtained with the real parameters (red) used to simulate the data and the posterior means(black). The dashed lines correspond to the potential dynamic if the steady state SSong is reached.For the estimated solution, the dashed and the solid line coincide.This illustrates the difficulty in estimating these parameters and emphasizes the importance of con-sidering the full distribution rather than just point estimates. This issue cannot be detected with theoriginal scVelo implementation, which only provides point estimates as output. As a consequence ofthis weak identifiability, all the other parameters are wrongly estimated, i.e., the associated 95% credibleintervals do not contain the true value, and the entire structure in the space (s, u) describing the timedynamic, is very different from the real one, as shown in Figure 3.Simulated examples demonstrate that estimating the parameter τcg and other unknowns in the modelis feasible when we have repeated measures for each (c, g). The results obtained with nc = 8, ng = 5and 450 repetitions are shown in Figure 4 as an example. However, in the case of single-cell data, it is notpossible to have true repetitions since the variable tcg is unobservable/unknown. Instead of repetitions,a mixture model can be used where data share a common coordinate in the space (s, u). These resultssuggest that this approach may be a viable direction.5. Conclusions and further developmentsThis study highlights the issues present in the current formalization of the scVelo model. Specifically,we focus on the weak identifiability of the variable τcg and its impact on the estimation of other parame-−20 −10 0 10 200.000.050.100.150.200.25DensityP(τcg = 0) = 0.11P(τcg = ∞) = 0(a)−20 −10 0 10 200.000.020.040.060.08DensityP(τcg = 0) = 0.14P(τcg = ∞) = 0.08(b)−20 −10 0 10 200.000.100.200.30DensityP(τcg = 0) = 0.07P(τcg = ∞) = 0(c)Figure 4: Repeated measurements model. Prior (red dashed line) and posterior (black solidline) of the logarithm of τcg for a three genes. The vertical lines represent the true values.ters. To address this, we introduce a new Bayesian version of scVelo and evaluate its performance usingsynthetic data. Upon inspection of the posterior distribution of τcg, we observe that, for a given gene, alldistributions are comparable and closely resemble the prior distribution, indicating weak identifiabilityof the parameters.In addition, we propose a potential solution to overcome the identifiability problem, which showspromising results in our initial investigations. We are pursuing this direction as a possible way forwardin improving the performance of scVelo.References[1] Bergen, V. and Lange, M. and Peidli, S. and Wolf, F. A. and Theis, F. J.: Generalizing RNAvelocity to transient cell states through dynamical modeling. Nature biotechnology 38.12, 1408–1414 (2020)[2] Gorin, G. and Fang, M. and Chari, T. and Pachter, L.: RNA velocity unraveled. PLOS Computa-tional Biology 18.9 (2022)[3] Gupta, R., Cerletti, D.Gut, G., Oxenius, A., Claassen, M.: Simulation-based inference of differen-tiation trajectories from RNA velocity fields. Cell Reports Methods, 2 1–15 (2022)[4] Jahnke, T. and Huisinga, W.: Solving the chemical master equation for monomolecular reactionsystems analytically. Journal of mathematical biology 54.1 1–26 (2007)[5] Marot-Lassauzaie, V. and Bouman, B. J. and Donaghy, F. D. and Demerdash, Y. and Essers, M. A.G. and Haghverdi, L.: Towards reliable quantification of cell state velocities. PLOS ComputationalBiology 18.9 (2022)

Analyzing RNA data with scVelo: identifiability issues and a Bayesian implementation

https://iris.polito.it/bitstream/11583/2982276/1/shortPaperSIS2023_Reviewed.pdf

Analyzing RNA data with scVelo: identifiability issues and a Bayesian implementation

Abstract

Similar works

Full text

Available Versions

PORTO@iris (Publications Open Repository TOrino - Politecnico di Torino)