OSCR

Comprehensive RNA velocity by modeling the cascade of gene regulation, transcription, and splicing from single-cell RNA sequencing data with TSvelo.

Code ↔ Paper

8 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 8 matches
  1. [1] § Results › TSvelo can model 3D gene dynamics and predict cell fate on pancreas dataset ↔ FIG8_paper.r, lines 41–130 · score 0.68 · Ngn3 low EP, Ngn3 high EP, Ductal, Beta, pancreas, Alpha
  2. [2] § Methods › Lineages segmentation and pseudotime initialization ↔ TSvelo/TSvelo_branch.py, lines 27–47 · score 0.62 · PAGA graph, detect lineages, shortest, Edges, cluster
  3. [3] § Results › TSvelo can capture gene dynamics well and predict cell fate on mouse brain data ↔ TSvelo/TSvelo_pp.py, lines 16–50 · score 0.59 · Radial Glia, mouse brain, IPCs, RNA
  4. [4] § Results › Estimate RNA velocity with TSVelo ↔ functions.R, lines 753–786 · score 0.59 · analytical solutions, degradation rate, splicing rate, transcription rate, trajectory, ODE
  5. [5] § Methods › Preprocessing for scRNA-seq data ↔ dataSimulationScVelo.r, lines 43–107 · score 0.58 · highly variable genes, variance, scVelo, seq, clustering, cell
  6. [6] § Methods › Optimizing global time and Neural ODE in EM framework ↔ TSvelo_run.py, lines 37–77 · score 0.58 · ChEA, Neural ODE, ENCODE, databases, TF, trained
  7. [7] § Results › TSvelo can predict cell fate and model lineage-specific gene dynamics for multi-lineage tasks ↔ functions.R, lines 944–985 · score 0.54 · degradation rates, splicing rates, scVelo, transcriptional rates, trajectory, branches
  8. [8] § Methods › Optimizing global time and Neural ODE in EM framework ↔ TSvelo/TSvelo_model.py, lines 127–187 · score 0.53 · Adam, gradient, trained, optimizer, loss, zero

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

R · 1,298 lines · 58 KB · no license · 2 matches

  1. # ------------- ANALYTIC FUNCTIONS ----------------
  3. # -------------------------------------------------------------------------
  4. # Function: u
  5. # Purpose:
  6. # Computes the amount of *unspliced RNA* for multiple cells for a fixed gene,
  7. # based on RNA velocity ODE's model
  8. #
  9. # Parameters:
  10. # t : Current time (scalar or vector).
  11. # t0_off : Initial time when the gene switches to the "off" state.
  12. # t0_on : Initial time when the gene switches to the "on" state (default = 0).
  13. # u0_off : Initial unspliced RNA amount in the "off" state.
  14. # u0_on : Initial unspliced RNA amount in the "on" state (optional).
  15. # If NA, it is set to the OFF steady-state value alpha_off/beta.
  16. # k : scalar or vector of same length of t indicating the transcriptional state
  17. # for each cell (0 = off, 2 = on).
  18. # alpha : Two-element vector of transcription rates:
  19. # alpha[1] = rate when gene is off,
  20. # alpha[2] = rate when gene is on.
  21. # beta : Splicing rate constant (scalar). Equal to 1 by default.
  22. #
  23. # Returns:
  24. # A numeric vector of the same length of t containing the unspliced RNA values
  25. # for each cell at time t.
  26. # -------------------------------------------------------------------------
  28. u <- function(t, t0_off, t0_on = 0, u0_off, u0_on = NA, k, alpha, beta = 1) {
  30. # If not provided, initialize u0_on as the steady-state level for the "off" state
  31. if (any(is.na(u0_on))) {
  32. u0_on <- alpha[1] / beta
  33. }
  35. n_cells <- length(k) # number of cells
  37. # Compute elapsed time (tau) depending on cell state
  38. tau <- ifelse(k == 2, t - t0_on, t - t0_off)
  40. # Precompute exponential decay terms for both states
  41. expBetaOFF <- exp(-beta * tau[k == 0])
  42. expBetaON <- exp(-beta * tau[k == 2])
  44. # Initialize result vector
  45. res <- rep(NA, n_cells)
  47. # Assign appropriate initial unspliced RNA based on cell state
  48. u0 <- ifelse(k == 2, u0_on, u0_off)
  50. # Compute unspliced RNA
  51. res[k == 0] <- u0[k == 0] * expBetaOFF + alpha[1]/beta * (1 - expBetaOFF)
  52. res[k == 2] <- u0[k == 2] * expBetaON + alpha[2]/beta * (1 - expBetaON)
  54. return(res)
  55. }
  57. # -------------------------------------------------------------------------
  58. # Function: s
  59. # Purpose:
  60. # Computes the amount of *spliced RNA* for multiple cells for a fixed gene,
  61. # based on RNA velocity ODE's model
  62. #
  63. # Parameters:
  64. # t : Current time (scalar or vector).
  65. # t0_off : Initial time when the gene switches to the "off" state.
  66. # t0_on : Initial time when the gene switches to the "on" state (default = 0).
  67. # u0_off : Initial unspliced RNA amount in the "off" state.
  68. # u0_on : Initial unspliced RNA amount in the "on" state (optional).
  69. # If NA, it is set to the OFF steady-state value alpha_off/beta.
  70. # s0_off : Initial spliced RNA amount in the "off" state.
  71. # s0_on : Initial spliced RNA amount in the "on" state (optional).
  72. # If NA, it is set to the OFF steady-state value alpha_off/gamma.
  73. # k : scalar or vector of same length of t indicating the transcriptional state
  74. # for each cell (0 = off, 2 = on).
  75. # alpha : Two-element vector of transcription rates:
  76. # alpha[1] = rate when gene is off,
  77. # alpha[2] = rate when gene is on.
  78. # beta : Splicing rate constant (scalar). Equal to 1 by default.
  79. # gamma : Degradation rate constant (scalar).
  80. #
  81. # Returns:
  82. # A numeric vector of the same length of t containing the spliced RNA values
  83. # for each cell at time t.
  84. # -------------------------------------------------------------------------
  85. s <- function(t, t0_off, t0_on = 0, u0_off, u0_on = NA, s0_off, s0_on = NA, k, alpha, beta = 1, gamma){
  87. # If not provided, initialize u0_on and s0_on as the steady-state level for the "off" state
  88. if(is.na(u0_on)){
  89. u0_on = alpha[1]/beta
  90. }
  91. if(is.na(s0_on)){
  92. s0_on = alpha[1]/gamma
  93. }
  95. n_cells <- length(k) # number of cells
  97. # Compute elapsed time (tau) depending on cell state
  98. tau <- ifelse(k == 2, t - t0_on, t - t0_off)
  99. # Precompute exponential decay terms for both states
  100. exp_gammaOFF <- exp(-gamma*tau[k == 0])
  101. exp_gammaON <- exp(-gamma*tau[k == 2])
  102. exp_betaOFF <- exp(-beta*tau[k == 0])
  103. exp_betaON <- exp(-beta*tau[k == 2])
  104. # Initialize result vector
  105. res <- rep(NA, n_cells)
  106. # Assign appropriate initial unspliced and spliced RNA based on cell state
  107. u0 <- ifelse(k == 2, u0_on, u0_off)
  108. s0 <- ifelse(k == 2, s0_on, s0_off)
  110. # Compute unspliced RNA
  111. res[k == 0] <- s0[k == 0]*exp_gammaOFF + alpha[1]/gamma*(1-exp_gammaOFF)+ (alpha[1]-beta*u0[k == 0])/(gamma-beta)*(exp_gammaOFF - exp_betaOFF)
  112. res[k == 2] <- s0[k == 2]*exp_gammaON + alpha[2]/gamma*(1-exp_gammaON) + (alpha[2]-beta*u0[k == 2])/(gamma-beta)*(exp_gammaON - exp_betaON)
  113. return(res)
  114. }
  116. # -------------------------------------------------------------------------
  117. # Function: s_di_u
  118. # Purpose:
  119. # Computes the amount of *spliced RNA* for multiple cells for a fixed gene,
  120. # given the amount of the *unspliced RNA*.
  121. # Note that this formula holds when beta is different from gamma.
  122. #
  123. # Parameters:
  124. # u : Current value of unspliced RNA for he different cells (scalar or vector).
  125. # u0_off : Initial unspliced RNA amount in the "off" state.
  126. # u0_on : Initial unspliced RNA amount in the "on" state (optional).
  127. # If NA, it is set to the OFF steady-state value alpha_off/beta.
  128. # k : scalar or vector of same length of t indicating the transcriptional state
  129. # for each cell (0 = off, 2 = on).
  130. # alpha : Two-element vector of transcription rates:
  131. # alpha[1] = rate when gene is off,
  132. # alpha[2] = rate when gene is on.
  133. # beta : Splicing rate constant (scalar). Equal to 1 by default.
  134. # gamma : Degradation rate constant (scalar).
  135. #
  136. # Returns:
  137. # A numeric vector of the same length of t containing the spliced RNA values
  138. # for each cell with unspliced RNA u.
  139. # -------------------------------------------------------------------------
  140. s_di_u <- function(u, u0_off, s0_off, u0_on = NA, s0_on = NA, k, alpha, beta = 1, gamma){
  141. # check that beta not equal to gamma
  142. if(beta == gamma){
  143. warning("WARNING: this formula holds for beta not equal gamma, so in this case you can not use it!")
  144. }
  146. # If not provided, initialize u0_on and s0_on as the steady-state level for the "off" state
  147. if(is.na(u0_on)){
  148. u0_on = alpha[1]/beta
  149. }
  150. if(is.na(s0_on)){
  151. s0_on = alpha[1]/gamma
  152. }
  154. # set the parameters depending of the state
  155. if(k == 0){
  156. a <- alpha[1]
  157. u0 <- u0_off
  158. s0 <- s0_off
  159. }else{
  160. a <- alpha[2]
  161. u0 <- u0_on
  162. s0 <- s0_on
  163. }
  165. # compute the spliced RNA
  166. res <- (s0 - a/gamma + (a-beta*u0)/(gamma-beta))*((beta*u-a)/(beta*u0 - a))^(gamma/beta) - (a*beta)/(gamma*(gamma-beta)) + beta/(gamma-beta)*u
  168. return(res)
  169. }
  171. # -------------------------------------------------------------------------
  172. # Function: u0
  173. # Purpose:
  174. # Computes the amount of u-coordinate of the initial point (for on or off phase) for a fixed gene.
  175. #
  176. # Parameters:
  177. # t0_off : Initial time when the gene switches to the "off" state (optional).
  178. # If NA, it is set to Infinity, meaning that we reach the upper steady state and do not enter in the repressive phase.
  179. # t0_on : Initial time when the gene switches to the "on" state (default = 0).
  180. # u0_off : Initial unspliced RNA amount in the "off" state (optional).
  181. # If NA and k == 2, it is set to the ON steady-state value alpha_on/beta.
  182. # u0_on : Initial unspliced RNA amount in the "on" state (optional).
  183. # If NA and k == 0, it is set to the OFF steady-state value alpha_off/beta.
  184. # k : scalar indicating if we want to compute the initial point for the on (k = 2) or
  185. # off (k = 0) phase.
  186. # alpha : Two-element vector of transcription rates:
  187. # alpha[1] = rate when gene is off,
  188. # alpha[2] = rate when gene is on.
  189. # beta : Splicing rate constant (scalar). Equal to 1 by default.
  190. #
  191. # Returns:
  192. # A scalar with the computed u0.
  193. # -------------------------------------------------------------------------
  194. u0 <- function(t0_off = NA, t0_on = 0, u0_off = NA, u0_on = NA, k, alpha, beta){
  196. # we want to compute the initial u-coordinate for the repressive phase
  197. if(k == 0){
  198. if(any(is.na(u0_on))){
  199. u0_on <- alpha[1]/beta # the initial u-coordinate for the inductive phase coincides with the lower steady state
  200. }
  202. if(any(is.na(t0_off))){
  203. t <- rep(NA, length(k))
  204. t[which(is.na(t0_off))] <- rep(Inf, length(which(is.na(t0_off)))) # the cells do not switch and remain in the inductive state, such that the upper steady state is reached
  205. }else{
  206. t <- t0_off # the upper steady state is not reached and the dynamic switch to the repressive phase
  207. }
  208. res <- u(t = t, t0_off = t0_off, t0_on = t0_on, u0_off = NA, u0_on = u0_on, k = rep(2, length(t)), alpha = alpha, beta = beta)
  209. }
  211. # we want to compute the initial u-coordinate for the inductive phase
  212. if(k == 2){
  213. if(any(is.na(u0_off))){
  214. u0_off <- alpha[2]/beta # the initial u-coordinate for the inductive phase coincides with the lower steady state
  215. }
  216. res <- u(t = Inf, t0_off = 0, t0_on = NA, u0_off = u0_off, u0_on = NA, k = rep(0, length(t)), alpha = alpha, beta = beta)
  217. }
  218. return(res)
  219. }
  221. # -------------------------------------------------------------------------
  222. # Function: s0
  223. # Purpose:
  224. # Computes the amount of s-coordinate of the initial point (for on or off phase) for a fixed gene.
  225. #
  226. # Parameters:
  227. # t0_off : Initial time when the gene switches to the "off" state (optional).
  228. # If NA, it is set to Infinity, meaning that we reach the upper steady state and do not enter in the repressive phase.
  229. # t0_on : Initial time when the gene switches to the "on" state (default = 0).
  230. # u0_off : Initial unspliced RNA amount in the "off" state (optional).
  231. # If NA and k == 2, it is set to the ON steady-state value alpha_on/beta.
  232. # If NA and k == 0, the initial point is computed with the function u0()
  233. # u0_on : Initial unspliced RNA amount in the "on" state (optional).
  234. # If NA and k == 0, it is set to the OFF steady-state value alpha_off/beta.
  235. # s0_off : Initial spliced RNA amount in the "off" state (optional).
  236. # If NA and k == 2, it is set to the ON steady-state value alpha_on/gamma
  237. # s0_on : Initial spliced RNA amount in the "on" state (optional).
  238. # If NA and k == 0, it is set to the OFF steady-state value alpha_off/gamma
  239. # k : scalar indicating if we want to compute the initial point for the on (k = 2) or
  240. # off (k = 0) phase.
  241. # alpha : Two-element vector of transcription rates:
  242. # alpha[1] = rate when gene is off,
  243. # alpha[2] = rate when gene is on.
  244. # beta : Splicing rate constant (scalar). Equal to 1 by default.
  245. # gamma : Degradation rate constant (scalar).
  247. # Returns:
  248. # A scalar with the computed s0.
  249. # -------------------------------------------------------------------------
  250. s0 <- function(t0_off = NA, t0_on = 0, u0_off = NA, u0_on = NA, s0_off = NA, s0_on = NA, k, alpha, beta, gamma){
  251. if(k == 0){
  252. if(any(is.na(s0_on))){
  253. s0_on <- alpha[1]/gamma
  254. }
  255. if(any(is.na(u0_on))){
  256. u0_on <- alpha[1]/beta
  257. }
  259. if(any(is.na(t0_off))){
  260. t <- rep(NA, length(k))
  261. t[which(is.na(t0_off))] <- rep(Inf, length(which(is.na(t0_off)))) # the cells do not switch and remain in the inductive state, such that the upper steady state is reached
  262. }else{
  263. t <- t0_off # the upper steady state is not reached and the dynamic switch to the repressive phase
  264. }
  265. if(any(is.na(u0_off))){
  266. # compute the initial u-coordinate
  267. u0_off <- u0(t0_off = t0_off, t0_on = t0_on, u0_off = NA, u0_on = u0_on, k = 0, alpha = alpha, beta = beta)
  268. }
  270. # compute s0
  271. res <- s(t = t, t0_off = t0_off, t0_on = t0_on, u0_off = u0_off, u0_on = u0_on, s0_off = s0_off, s0_on = s0_on, k = rep(2, length(t)), alpha, beta, gamma)
  272. }
  274. if(k == 2){
  275. if(any(is.na(s0_off))){
  276. s0_off <- alpha[2]/gamma
  277. }
  278. if(any(is.na(u0_off))){
  279. u0_off <- alpha[2]/beta
  280. }
  282. res <- s(rep(Inf, length(t0_on)), t0_off = 0, t0_on = NA, u0_off = u0_off, u0_on = u0_on, s0_off = s0_off, s0_on = s0_on, k = rep(0, length(t0_on)), alpha = alpha, beta = beta, gamma = gamma)
  284. }
  286. return(res)
  287. }
  289. # -------------------------------------------------------------------------
  290. # Function: u_and_s_withTauMCMC
  291. # Purpose:
  292. # Compute unspliced (u) and spliced (s) RNA levels for multiple genes,
  293. # cells and MCMC samples (or other 3D structure) given per-cell elapsed
  294. # times (tau) and kinetic parameters. This implements the closed-form
  295. # solutions of linear ODEs for the 2-state transcription model
  296. # (states: off (k==0) and on (k==2)).
  297. #
  298. # Inputs (expected shapes):
  299. # tau : array (dim: [subtype, n_genes, n_samples] or similar)
  300. # elapsed time for each subgroup, gene and MCMC sample.
  301. # u0_off : array with initial unspliced RNA in the "off" state
  302. # u0_on : array with initial unspliced RNA in the "on" state.
  303. # s0_off : array with initial spliced RNA in the "off" state.
  304. # s0_on : array with initial spliced RNA in the "o" state.
  305. # k : integer array with same dim as tau indicating
  306. # transcriptional state per element (0 = off, 2 = on).
  307. # alpha : numeric array of transcription rates. Expected to
  308. # have dimensions [genes, state_index(1/2), samples]
  309. # (uses alpha[,1,] for off and alpha[,2,] for on).
  310. # beta : array for splicing rate ([n_genes, n_samples]).
  311. # gamma : array for degradation rate ([n_genes, n_samples]).
  312. # subtypeCell : vector indicating subtype index for each cell
  313. # (used to assign off-state initial values by subtype).
  314. # typeCellT0_off : vector indicating switching time labels
  315. #
  316. # Returns:
  317. # A list with:
  318. # $u : array of same shape as tau with computed unspliced RNA values.
  319. # $s : array of same shape as tau with computed spliced RNA values.
  320. # -------------------------------------------------------------------------
  321. u_and_s_withTauMCMC <- function(tau, u0_off, u0_on = NA, s0_off, s0_on = NA, k, alpha, beta, gamma, subtypeCell, typeCellT0_off){
  323. # Set inductive initial values equal to off steady state coordinates, if not assigned.
  324. if(any(is.na(u0_on))){
  325. u0_on <- alpha[, 1, ]/beta
  326. }
  327. if(any(is.na(s0_on))){
  328. s0_on <- alpha[, 1, ]/gamma
  329. }
  330. # expand the elements, such that they will have the same dimension of tau
  331. u0_on = sweep(k == 2, MARGIN = c(2, 3), STATS = u0_on, FUN = "*")
  332. s0_on = sweep(k == 2, MARGIN = c(2, 3), STATS = s0_on, FUN = "*")
  334. # expand also the off initial coordinates, such that they will have the same dimension of tau
  335. u0_off_new <- array(NA, dim = dim(tau))
  336. s0_off_new <- array(NA, dim = dim(tau))
  337. for(sty in unique(subtypeCell)){
  338. # assign the corresponding off-state initial arrays to all elements
  339. # in each subtypes
  340. tyT0_off <- typeCellT0_off[which(subtypeCell == sty)[1]]
  341. u0_off_new[sty, , ] = u0_off[tyT0_off, , ]
  342. s0_off_new[sty, , ] = s0_off[tyT0_off, , ]
  343. }
  344. # Compose the full initial conditions: choose on vs off initial values
  345. # depending on k (state indicator).
  346. u0 = u0_on * (k == 2) + u0_off_new * (k == 0)
  347. s0 = s0_on * (k == 2) + s0_off_new * (k == 0)
  349. # upper-steady state u-coordinate
  350. uSS_on <- alpha[,2,]/beta
  352. # Exponential decay factors for splicing (beta) and degradation (gamma) rates.
  353. expBeta <- exp(-sweep(tau, MARGIN = c(2, 3), STATS = beta, FUN = "*"))
  354. expGamma <- exp(-sweep(tau, MARGIN = c(2, 3), STATS = gamma, FUN = "*"))
  355. p2 <- sweep(k == 0, MARGIN = c(2, 3), STATS = alpha[, 1, ]/beta, FUN = "*")*(1-expBeta)
  357. # Transcription rate depending on the state
  358. alpha <- sweep(k == 2, MARGIN = c(2, 3), STATS = alpha[, 2, ], FUN = "*") + sweep(k == 0, MARGIN = c(2, 3), STATS = alpha[, 1, ], FUN = "*")
  360. # unspliced and spliced values
  361. resS <- s0*expGamma + sweep(alpha, MARGIN = c(2, 3), STATS = gamma, FUN = "/")*(1-expGamma) + sweep((alpha-sweep(u0, MARGIN = c(2, 3), STATS = beta, FUN = "*")), MARGIN = c(2, 3), STATS = (gamma-beta), FUN = "/")*(expGamma - expBeta)
  362. resU <- u0*expBeta + p2 + sweep(k==2, MARGIN = c(2, 3), STATS = uSS_on, FUN = "*")*(1-expBeta)
  364. return(list(u = resU, s = resS))
  365. }
  367. # -------------------------------------------------------------------------
  368. # Function: u0_MCMC
  369. # Purpose:
  370. # Compute initial unspliced RNA (u0) for MCMC-shaped arrays.
  371. #
  372. # Arguments:
  373. # t0_off : 3D array [n_switching_clusters, n_genes, n_samples] with time of OFF-switch.
  374. # t0_on : scalar or array; default 0, for the starting time of the on phase.
  375. # k : scalar or vector of states (0 or 2).
  376. # alpha : array [n_genes, 2, n_samples] for transcription rate (alpha[,1,] = off, alpha[,2,] = on).
  377. # beta : scalar or array [n_genes, n_samples] for splicing rate. By default equal to 1.
  378. #
  379. # Returns:
  380. # Array of same shape as t0_off with u0 values for each element.
  381. # Note that we can compute just one between u0_off and u0_on (for different genes) for each call of this function.
  382. # -------------------------------------------------------------------------
  383. u0_MCMC <- function(t0_off = NA, t0_on = 0, k, alpha, beta = 1){
  384. # If multiple k values provided, keep only unique values.
  385. if(length(k) > 1){
  386. k <- unique(k)
  387. }
  389. # initial point for repressive phase
  390. if(k == 0){
  391. # transform u0_on according to the dimension of t0_off
  392. u0_on <- alpha[, 1, ]/beta
  393. u0_on <- array(rep(as.vector(u0_on),dim(t0_off)[1]), dim = c(dim(t0_off)[2], dim(t0_off)[3], dim(t0_off)[1]))
  394. u0_on <- aperm(u0_on, c(3, 1, 2))
  396. t <- t0_off # we assume the switch from on to off phase occurs
  397. t0_on <- array(t0_on, dim = dim(u0_on))
  398. tau <- t - t0_on # elapsed time in the dynamic
  400. # coordinate of the on steady state
  401. uSS_ON <- alpha[, 2,]/beta
  402. uSS_ON <- array(rep(as.vector(uSS_ON), dim(t0_off)[1]), dim = c(dim(t0_off)[2], dim(t0_off)[3], dim(t0_off)[1]))
  403. uSS_ON <- aperm(uSS_ON, c(3, 1, 2))
  405. # modify the splicing rate accordingly to the dimension of t0:off
  406. beta <- array(rep(as.vector(beta), dim(t0_off)[1]), dim = c(dim(t0_off)[2], dim(t0_off)[3], dim(t0_off)[1]))
  407. beta <- aperm(beta, c(3, 1, 2))
  409. expBetaOFF <- exp(-beta*tau)
  410. # u0_off
  411. res <- u0_on*expBetaOFF + uSS_ON*(1-expBetaOFF)
  412. }
  414. # initial point for inductive phase
  415. if(any(k == 2)){
  416. res <- alpha[,1,]/beta
  417. }
  418. return(res)
  419. }
  421. # -------------------------------------------------------------------------
  422. # Function: s0_MCMC
  423. # Purpose:
  424. # Compute initial spliced RNA (s0) for MCMC-shaped arrays.
  425. #
  426. # Arguments:
  427. # t0_off : 3D array [n_switching_clusters, n_genes, n_samples] with time of OFF-switch.
  428. # t0_on : scalar or array; default 0, for the starting time of the on phase.
  429. # k : scalar or vector of states (0 or 2).
  430. # alpha : array [n_genes, 2, n_samples] for transcription rate (alpha[,1,] = off, alpha[,2,] = on).
  431. # beta : scalar or array [n_genes, n_samples] for splicing rate. By default equal to 1.
  432. # gamma : scalar or array [n_genes, n_samples] for degradation rate.
  433. #
  434. # Returns:
  435. # Array of same shape as t0_off with s0 values for each element.
  436. # Note that we can compute just one between s0_off and s0_on (for different genes) for each call of this function.
  437. # -------------------------------------------------------------------------
  438. s0_MCMC <- function(t0_off = NA, t0_on = 0, k, alpha, beta, gamma){
  439. # If multiple k values provided, keep only unique values.
  440. if(length(k) > 1){
  441. k <- unique(k)
  442. }
  444. # initial point for repressive phase
  445. if(k == 0){
  446. # transform s0_on according to the dimension of t0_off
  447. s0_on = alpha[, 1, ]/gamma
  448. s0_on <- array(rep(as.vector(s0_on),dim(t0_off)[1]), dim = c(dim(t0_off)[2], dim(t0_off)[3], dim(t0_off)[1]))
  449. s0_on <- aperm(s0_on, c(3, 1, 2))
  451. # transform s0_on according to the dimension of t0_off
  452. u0_on = alpha[,1,]/beta
  453. u0_on <- array(rep(as.vector(u0_on),dim(t0_off)[1]), dim = c(dim(t0_off)[2], dim(t0_off)[3], dim(t0_off)[1]))
  454. u0_on <- aperm(u0_on, c(3, 1, 2))
  456. t <- t0_off # we assume the switch from on to off phase occurs
  457. t0_on <- array(t0_on, dim = dim(u0_on))
  458. tau <- t - t0_on # elapsed time in the dynamic
  460. # coordinate of the on steady state
  461. sSS_ON <- alpha[,2,]/gamma
  462. sSS_ON <- array(rep(as.vector(sSS_ON), dim(t0_off)[1]), dim = c(dim(t0_off)[2], dim(t0_off)[3], dim(t0_off)[1]))
  463. sSS_ON <- aperm(sSS_ON, c(3, 1, 2))
  465. # transform the rates according to the dimension of t0:off
  466. beta <- array(rep(as.vector(beta), dim(t0_off)[1]), dim = c(dim(t0_off)[2], dim(t0_off)[3], dim(t0_off)[1]))
  467. beta <- aperm(beta, c(3, 1, 2))
  468. gamma <- array(rep(as.vector(gamma), dim(t0_off)[1]), dim = c(dim(t0_off)[2], dim(t0_off)[3], dim(t0_off)[1]))
  469. gamma <- aperm(gamma, c(3, 1, 2))
  470. alpha2 <- array(rep(as.vector(alpha[,2,]), dim(t0_off)[1]), dim = c(dim(t0_off)[2], dim(t0_off)[3], dim(t0_off)[1]))
  471. alpha2 <- aperm(alpha2, c(3, 1, 2))
  473. exp_gammaOFF <- exp(-gamma*tau)
  474. exp_betaOFF <- exp(-beta*tau)
  476. res <- s0_on*exp_gammaOFF + sSS_ON*(1-exp_gammaOFF)+ (alpha2-beta*u0_on)/(gamma-beta)*(exp_gammaOFF - exp_betaOFF)
  477. }
  479. # initial point for inductive phase
  480. if(any(k == 2)){
  481. res <- alpha[,1,]/gamma
  482. }
  483. return(res)
  484. }
  487. # ------------- LOAD THE DATA ----------------
  489. # -------------------------------------------------------------------------
  490. # Function: nameGenes
  491. # Purpose:
  492. # Load the name of the considered genes for the real data.
  493. #
  494. # Parameters:
  495. # pathData: path where the file (named as "var.csv") with the name of the genes is stored.
  497. # Returns:
  498. # A vector with the name of the genes.
  499. # -------------------------------------------------------------------------
  500. nameGenes <- function(pathData){
  501. var <- read.csv(paste(pathData, "/var.csv", sep = ""))
  502. nameG <- var$index
  503. return(nameG)
  504. }
  506. # -------------------------------------------------------------------------
  507. # Function: nameCells
  508. # Purpose:
  509. # Load the name of the considered cells for the real data.
  510. #
  511. # Parameters:
  512. # pathData: path where the file (named as "obs.csv") with the name of the cells is stored.
  514. # Returns:
  515. # A vector with the name of the cells.
  516. # -------------------------------------------------------------------------
  517. nameCells <- function(pathData){
  518. obs <- read.csv(paste(pathData, "/obs.csv", sep = ""))
  519. nameC<- obs$index
  520. return(nameC)
  521. }
  523. # -------------------------------------------------------------------------
  524. # Function: loadRealData
  525. # Purpose:
  526. # Load the values of spliced and unspliced RNA for the real data.
  527. #
  528. # Parameters:
  529. # typeProcessing: type of pre-processing that have been previously applied to the data.
  530. # It can be "filter", "filter_and_normalize", "filter_and_normalize_noLog", filter_and_normalize moments"
  531. # pathData: path where the data are stored.
  532. # names: boolean if we want to keep the name of cells and genes or not. By default it is FALSE.
  533. # sparse: boolean if we want to transform the loaded matrices in sparse objects or not. By default it is FALSE.
  535. # Returns:
  536. # A list with
  537. # -"unspliced": matrix unspliced RNA values.
  538. # -"spliced" : matrix spliced RNA values.
  539. # -"typeCell" : vector with the group label associated to each cell.
  540. # -------------------------------------------------------------------------
  541. loadRealData <- function(typeProcessing, pathData, names = FALSE, sparse = FALSE){
  542. if(is.null(pathData)){
  543. error("Provide the directory where the real data are stored")
  544. }
  546. pathData <- paste0(pathData, "/", typeProcessing)
  548. if(typeProcessing == "moments"){ # import pre-processed continuous data
  549. unspliced <- loadMu(pathData, names, sparse)
  550. spliced <- loadMs(pathData, names, sparse)
  551. }else{ # import raw counts
  552. unspliced <- loadUnspliced(pathData, names, sparse)
  553. spliced <- loadSpliced(pathData, names, sparse)
  554. }
  556. # import the type of cells
  557. typeCell <- loadObs(pathData)$clusters
  559. res <- list("unspliced" = unspliced, "spliced" = spliced, "typeCell" = typeCell)
  560. return(res)
  561. }
  563. # -------------------------------------------------------------------------
  564. # Function: loadVar
  565. # Purpose:
  566. # Load the var.csv file, containing the name of the genes for the real data.
  567. #
  568. # Parameters:
  569. # pathData: path where the file (named as "var.csv") is stored.
  570. #
  571. # Returns: a matrix with the content of the var file.
  572. # -------------------------------------------------------------------------
  573. loadVar <- function(pathData){
  574. var <- read.csv(file = paste(pathData, '/var.csv', sep =""))
  575. return(var)
  576. }
  578. # -------------------------------------------------------------------------
  579. # Function: loadObs
  580. # Purpose:
  581. # Load the obs.csv file, containing the name of the cells and the type of cells'labels for the real data.
  582. #
  583. # Parameters:
  584. # pathData: path where the file (named as "obs.csv") is stored.
  585. #
  586. # Returns: a matrix with the content of the obs file.
  587. # -------------------------------------------------------------------------
  588. loadObs <- function(pathData){
  589. obs <- read.csv(file = paste(pathData, '/obs.csv', sep = ""))
  590. return(obs)
  591. }
  593. # -------------------------------------------------------------------------
  594. # Function: loadUnspliced
  595. # Purpose:
  596. # Load the matrix with unspliced RNA counts.
  597. #
  598. # Parameters:
  599. # pathData: path where the file with unspliced counts is stored.
  600. # names: boolean if wee want to store the name of cells and genes. By default it is FALSE.
  601. # sparse: boolean if we want to transform the loaded matrix in sparse object or not. By default it is FALSE.
  603. # Returns: a matrix with the content of the unspliced data.
  604. # -------------------------------------------------------------------------
  605. loadUnspliced <- function(pathData, names = FALSE, sparse = FALSE){
  606. # import unspliced raw data and transform them into a sparse matrix, if required
  607. unspliced <- read.csv(file = paste(pathData, '/unspliced.csv', sep = ""))
  609. if(names){
  610. nameG <- nameGenes(pathData)
  611. nameC <- nameCells(pathData)
  612. colnames(unspliced) <- nameG
  613. rownames(unspliced) <- nameC
  614. }
  616. unspliced <- as(unspliced, "matrix")
  618. if(sparse){
  619. unspliced <- as(unspliced, "dgCMatrix")
  620. }
  622. return(unspliced)
  623. }
  625. # -------------------------------------------------------------------------
  626. # Function: loadSpliced
  627. # Purpose:
  628. # Load the matrix with spliced RNA counts.
  629. #
  630. # Parameters:
  631. # pathData: path where the file with spliced counts is stored.
  632. # names: boolean if wee want to store the name of cells and genes. By default it is FALSE.
  633. # sparse: boolean if we want to transform the loaded matrix in sparse object or not. By default it is FALSE.
  635. # Returns: a matrix with the content of the spliced data.
  636. # -------------------------------------------------------------------------
  637. loadSpliced <- function(pathData, names = FALSE, sparse = FALSE){
  638. # import spliced raw data and transform them into a sparse matrix, if required
  639. spliced <- read.csv(file = paste(pathData, '/spliced.csv', sep = ""))
  641. if(names){
  642. nameG <- nameGenes(pathData)
  643. nameC <- nameCells(pathData)
  644. colnames(unspliced) <- nameG
  645. rownames(unspliced) <- nameC
  646. }
  648. spliced <- as(spliced, "matrix")
  650. if(sparse){
  651. spliced <- as(spliced, "dgCMatrix")
  652. }
  654. return(spliced)
  655. }
  657. # -------------------------------------------------------------------------
  658. # Function: loadMu
  659. # Purpose:
  660. # Load the matrix with unspliced RNA moments (previously computed by scVelo)
  661. #
  662. # Parameters:
  663. # pathData: path where the file with unspliced moments is stored.
  664. # names: boolean if wee want to store the name of cells and genes. By default it is FALSE.
  665. # sparse: boolean if we want to transform the loaded matrix in sparse object or not. By default it is FALSE.
  667. # Returns: a matrix with the the unspliced moments.
  668. # -------------------------------------------------------------------------
  669. loadMu <- function(pathData, names = FALSE, sparse = FALSE){
  670. # import unspliced pre-processed moments and transform them into a sparse matrix, if required
  671. Mu <- read.csv(file = paste(pathData, '/Mu.csv', sep = ""))
  673. if(names){
  674. nameG <- nameGenes(pathData)
  675. nameC <- nameCells(pathData)
  676. colnames(Mu) <- nameG
  677. rownames(Mu) <- nameC
  678. }
  680. Mu <- as(Mu, "matrix")
  682. if(sparse){
  683. Mu <- as(Mu, "dgCMatrix")
  684. }
  686. return(Mu)
  687. }
  689. # -------------------------------------------------------------------------
  690. # Function: loadMs
  691. # Purpose:
  692. # Load the matrix with spliced RNA moments (previously computed by scVelo)
  693. #
  694. # Parameters:
  695. # pathData: path where the file with spliced moments is stored.
  696. # names: boolean if wee want to store the name of cells and genes. By default it is FALSE.
  697. # sparse: boolean if we want to transform the loaded matrix in sparse object or not. By default it is FALSE.
  699. # Returns: a matrix with the the spliced moments.
  700. # -------------------------------------------------------------------------
  701. loadMs <- function(pathData, names = FALSE, sparse = FALSE){
  702. # import spliced pre-processed moments and transform them into a spare matrix, if required
  703. Ms <- read.csv(file = paste(pathData, '/Ms.csv', sep = ""))
  705. if(names){
  706. nameG <- nameGenes(pathData)
  707. nameC <- nameCells(pathData)
  708. colnames(Ms) <- nameG
  709. rownames(Ms) <- nameC
  710. }
  712. Ms <- as(Ms, "matrix")
  714. if(sparse){
  715. Ms <- as(Ms, "dgCMatrix")
  716. }
  718. return(Ms)
  719. }
  721. # -------------------------------------------------------------------------
  722. # Function: nameReal_Pancreas
  723. # Purpose:
  724. # Convert the name of thee simulations into the names used in Table 5 of bayVel's paper
  725. #
  726. # Parameters:
  727. # name: name of the considered simulation
  729. # Returns: string with the corresponding name used in Table 5.
  730. # -------------------------------------------------------------------------
  731. nameReal_Pancreas <- function(name){
  732. label <- ""
  733. if(grepl("SW1", name)){
  734. label <- paste0(label,"K = 1,$")
  735. }else{
  736. label <- paste0(label,"K = 8,$")
  737. }
  739. if(grepl("T1", name)){
  740. label <- paste(label,"$R = 1$", sep = "-")
  741. }else if(grepl("T2", name)){
  742. label <- paste(label,"$R = 9$", sep = "-")
  743. }else if(grepl("T3", name)){
  744. label <- paste(label,"$R = 38$", sep = "-")
  745. }
  747. return(label)
  748. }
  751. # ------------- PLOT THE DATA ----------------
  753. # -------------------------------------------------------------------------
  754. # Function: plot_sVSu
  755. # Purpose:
  756. # Plot the phase trajectory (s vs u) for a single gene using the
  757. # analytical solutions of the ODE model (functions `u()` and `s()`).
  758. # The function draws both the dynamic when the switch occurs and the
  759. # potential dynamic when the switch does not occur. The function admit to
  760. # overlay observed points for cell types.
  761. #
  762. # Arguments:
  763. # t0_off : numeric scalar or Inf, time of the OFF switch (if Inf,
  764. # the system stays in the ON steady-state before switching).
  765. # t0_on : numeric scalar (default 0), time when ON state starts.
  766. # alpha : transcription-rate vector.
  767. # beta : splicing rate, equal to 1 by default.
  768. # gamma : degradation rate.
  769. # pos_u : vector of u positions (one per cell/subgroup). Omit if you want just the model dynamic.
  770. # pos_s : vector of s positions (one per cell/subgroup). Omit if you want just the model dynamic.
  771. # subGrLabels : vector of subgroup-labels for each cell.
  772. # g : gene identifier (used in the plot title).
  773. # add : boolean; if FALSE create a new ggplot, if TRUE add to `gg`. By default it is equal to FALSE.
  774. # gg : an existing ggplot object to add layers to when add = TRUE.
  775. # colCell : colors for subgroup points (vector or NA to auto-pick).
  776. # colDyn : color for the model dynamic (optional, default "red").
  777. # xlim, ylim : numeric scalars for axis limits (optional, if NA they are not set).
  778. # axisTitle.size, axisText.size, title.size : sizes for theme elements.
  779. # lineSize : line width for model dynamic.
  780. # shapePoint : shape of subgroup points (optional).
  781. # sizePoint : size of subgroup points (optional).
  782. #
  783. # Returns:
  784. # A ggplot object with the s vs u phase plot for the given gene.
  785. # -------------------------------------------------------------------------
  786. plot_sVSu <- function(
  787. t0_off,
  788. t0_on = 0,
  789. alpha,
  790. beta = 1,
  791. gamma,
  792. pos_u,
  793. pos_s,
  794. subGrLabels,
  795. g,
  796. add = FALSE,
  797. gg = NA,
  798. colCell = NA,
  799. colDyn = "red",
  800. xlim = NA,
  801. ylim = NA,
  802. axisTitle.size = NA,
  803. axisText.size = NA,
  804. title.size = NA,
  805. lineSize = 1,
  806. shapePoint = 21,
  807. sizePoint = 2,
  808. ...
  809. ){
  812. # Compute the initial points of the on and of the off phase
  813. u0_off <- u0(k = 0, alpha = alpha, beta = beta)
  814. u0_on <- u0(k = 2, alpha = alpha, beta = beta)
  815. s0_off <- s0(k = 0, alpha = alpha, beta = beta, gamma = gamma)
  816. s0_on <- s0(k = 2, alpha = alpha, beta = beta, gamma = gamma)
  818. # Compute the points lying on the on branch of the dynamic. We assume here the switching point is infinite, in order to plot the potential behavior up to the upper steady state.
  819. t_seq <- seq(t0_on, 500, 0.1)
  820. u_plot <- u(t_seq, t0_off = Inf, t0_on = t0_on, u0_off = u0_off, u0_on = u0_on, k = rep(2, length(t_seq)), alpha, beta)
  821. s_plot <- s(t_seq, t0_off = Inf, t0_on = t0_on, u0_off = u0_off, u0_on = u0_on, s0_off = s0_off, s0_on = s0_on, k = rep(2, length(t_seq)), alpha, beta, gamma)
  823. df_on <- data.frame(s_plot = s_plot, u_plot = u_plot)
  825. # Title for the plot
  826. main <- paste("Gene", g)
  827. if(length(unique(subGrLabels)) == 1){
  828. main <- paste(main, ", typeCell", unique(subGrLabels)) # we are plotting just the position of one subgroup
  829. }
  831. # Initialize or extend ggplot with the ON branch (dashed line)
  832. if(!add){
  833. gg <- ggplot(data = df_on, aes(x = s_plot, y = u_plot)) +
  834. geom_path(lty = 2, col = colDyn, linewidth = lineSize) +
  835. xlab("s") +
  836. ylab("u") +
  837. labs(title = main)
  838. }else{
  839. gg <- gg +
  840. geom_path(data = df_on, aes(x = s_plot, y = u_plot), lty = 2, col = colDyn, linewidth = lineSize)
  841. }
  842. # adjust graphical parameters
  843. if(is.numeric(axisText.size)){
  844. gg <- gg + theme(axis.text = element_text(size = axisText.size))
  845. }
  846. if(is.numeric(axisTitle.size)){
  847. gg <- gg + theme(axis.title = element_text(size = axisTitle.size))
  848. }
  849. if(is.numeric(title.size)){
  850. gg <- gg + theme(plot.title = element_text(hjust = 0.5, size = title.size))
  851. }
  852. if(is.numeric(xlim)){
  853. gg <- gg + xlim(0, xlim)
  854. }
  855. if(is.numeric(ylim)){
  856. gg <- gg + ylim(0, ylim)
  857. }
  860. # Compute the points lying on the off branch of the dynamic, describing the potential behavior from to the upper steady state.
  861. u_plot <- u(t_seq, t0_off = 0, t0_on = t0_on, u0_off = u0_off, u0_on = u0_on, k = rep(0, length(t_seq)), alpha, beta)
  862. s_plot <- s(t_seq, t0_off = 0, t0_on = t0_on, u0_off = u0_off, u0_on = u0_on, s0_off = s0_off, s0_on = s0_on, k = rep(0, length(t_seq)), alpha, beta, gamma)
  864. df_off <- data.frame(s_plot = s_plot, u_plot = u_plot)
  865. gg <- gg + geom_path(data = df_off, aes(x = s_plot, y = u_plot), lty = 2, col = colDyn, linewidth = lineSize)
  867. # If a finite t0_off is provided, compute the transient branches:
  868. # - ON branch from t0_on to t0_off (solid line)
  869. # - OFF branch from t0_off to t0_off + 500 (solid line)
  870. # Else (t0_off == Inf) we reuse the previously computed branch(s).
  872. if(t0_off != Inf){
  873. # compute initial conditions at the switching times
  874. u0_off <- u0(t0_off = t0_off, k = 0, alpha = alpha, beta = beta)
  875. u0_on <- u0(t0_on = t0_on, k = 2, alpha = alpha, beta = beta)
  876. s0_off <- s0(t0_off = t0_off, k = 0, alpha = alpha, beta = beta, gamma = gamma)
  877. s0_on <- s0(t0_on = t0_on, k = 2, alpha = alpha, beta = beta, gamma = gamma)
  880. # ON transient: from t0_on up to t0_off (solid line)
  881. t_on_trans <- seq(t0_on, t0_off, 0.1)
  882. u_plot <- u(t_on_trans, t0_off = t0_off, t0_on = t0_on, u0_off = u0_off, u0_on = u0_on, k = rep(2, length(t_on_trans)), alpha, beta)
  883. s_plot <- s(t_on_trans, t0_off = t0_off, t0_on = t0_on, u0_off = u0_off, u0_on = u0_on, s0_off = s0_off, s0_on = s0_on, k = rep(2, length(t_on_trans)), alpha, beta, gamma)
  885. df_on_trans <- data.frame(s_plot = s_plot, u_plot = u_plot)
  886. gg <- gg +
  887. geom_path(data = df_on_trans, aes(x = s_plot, y = u_plot), lty = 1, col = colDyn, linewidth = lineSize)
  889. # OFF transient: from t0_off onward (solid line)
  890. t_off_trans <- seq(t0_off, t0_off + 500, 0.1)
  891. u_plot <- u(t_off_trans, t0_off = t0_off, t0_on = t0_on, u0_off = u0_off, u0_on = u0_on, k = rep(0, length(t_off_trans)), alpha, beta)
  892. s_plot <- s(t_off_trans, t0_off = t0_off, t0_on = t0_on, u0_off = u0_off, u0_on = u0_on, s0_off = s0_off, s0_on = s0_on, k = rep(0, length(t_off_trans)), alpha, beta, gamma)
  894. df_off_trans <- data.frame(s_plot = s_plot, u_plot = u_plot)
  895. gg <- gg +
  896. geom_path(data = df_off_trans, aes(x = s_plot, y = u_plot), lty = 1, col = colDyn, linewidth = lineSize)
  897. }else{
  898. # If infinite off-time: overlay the earlier computed branches.
  899. gg <- gg +
  900. geom_path(data = df_on, aes(x = s_plot, y = u_plot), lty = 1, col = colDyn, linewidth = lineSize) +
  901. geom_path(data = df_off, aes(x = s_plot, y = u_plot), lty = 1, col = colDyn, linewidth = lineSize)
  902. }
  904. # If pos_* vectors are provided (no NAs) add points (pos_s, pos_u) for each subgroup
  905. # Colors are auto-chosen if colCell is NA.
  906. subGrC <- unique(subGrLabels)
  908. if(sum(is.na(pos_s)) == 0){
  909. # If lengths already match unique types, use directly; otherwise pick
  910. # the first observed position for each subtype.
  911. if(length(pos_s) == length(subGrC)){
  912. df <- data.frame(s = pos_s, u = pos_u)
  913. }else{
  914. tmpS <- c()
  915. tmpU <- c()
  916. for(i in 1:length(subGrC)){
  917. x <- subGrC[i]
  918. tmpS <- c(tmpS, pos_s[which(subGrLabels == x)[1]])
  919. tmpU <- c(tmpU, pos_u[which(subGrLabels == x)[1]])
  920. }
  921. df <- data.frame(s = tmpS, u = tmpU)
  922. }
  925. # Choose palette if colCell not provided
  926. if(sum(is.na(colCell)) == 1){
  927. if(length(subGrC) >= 3){
  928. colCell <- brewer.pal(length(subGrC), "Spectral")
  929. }else if(length(subGrC) == 2){
  930. colCell <- c("darkgreen", "blue")
  931. }else if(length(subGrC) == 1){
  932. colCell <- c("darkgreen")
  933. }
  934. }
  936. df$colCell <- colCell
  937. # Add points to the plot
  938. gg <- gg + geom_point(data = df, aes(x = s, y = u, fill = colCell), fill = colCell, col = colDyn, size = sizePoint, shape = shapePoint)
  939. }
  941. return(gg)
  942. }
  944. # -------------------------------------------------------------------------
  945. # Function: plot_sVSu_scVelo
  946. # Purpose:
  947. # Plot the scVelo-style s vs u phase trajectory for a single gene using
  948. # analytical ODE solutions (functions u() and s()). It takes as input
  949. # the output of scVelo's inference. The function draws:
  950. # - ON and OFF branches (dashed)
  951. # - transient branches (solid) if fit_t0_off is finite
  952. # - cell positions computed from fit_t (scaled/translated to model fit)
  953. #
  954. # Notes:
  955. # - scVelo uses beta_scaled = fit_beta * fit_scaling; alpha_off is assumed 0 and
  956. # alpha_on = fit_alpha
  957. # - Trajectories are rescaled and translated with fit_scaling, fit_u0_offset,
  958. # fit_s0_offset to match scVelo plotting conventions.
  959. #
  960. # Arguments:
  961. # fit_t0_off : numeric scalar or Inf, fitted switch time to OFF state.
  962. # fit_t0_on : numeric scalar (default 0), fitted ON start time.
  963. # fit_alpha : on transcription rate.
  964. # fit_beta : splicing rate.
  965. # fit_gamma : degradation rate-
  966. # fit_scaling : scalar used to scale u (beta is multiplied by this).
  967. # fit_u0_offset : scalar added to u after scaling (translation).
  968. # fit_s0_offset : scalar added to s (translation).
  969. # fit_t : numeric vector of per-cell times used to place cells on the curve.
  970. # subGrLabels : vector of group labels for each element of fit_t (used for colors).
  971. # g : gene identifier used in the plot title.
  972. # add : boolean; if FALSE create new ggplot, if TRUE add to gg.
  973. # gg : existing ggplot object to add layers to (when add = TRUE).
  974. # colCell : color for points (vector or NA to auto-pick).
  975. # colDyn : color for the model dynamic and points border (optional, default "red").
  976. # xlim, ylim : numeric scalars for axis limits (optional, if NA they are not set).
  977. # axisTitle.size, axisText.size, title.size : sizes for theme elements.
  978. # lineSize : line width for model dynamic.
  979. # shapePoint : shape of points (optional).
  980. # sizePoint : size of points (optional).
  981. #
  982. # Returns:
  983. # A ggplot object with the s vs u phase plot for the given gene, accordingg to scVelo plots.
  984. # -------------------------------------------------------------------------
  985. plot_sVSu_scVelo <- function(
  986. fit_t0_off,
  987. fit_t0_on = 0,
  988. fit_alpha,
  989. fit_beta,
  990. fit_gamma,
  991. fit_scaling,
  992. fit_u0_offset,
  993. fit_s0_offset,
  994. fit_t,
  995. subGrLabels,
  996. g,
  997. add = FALSE,
  998. gg = NA,
  999. colCell = NA,
  1000. colDyn = "red",
  1001. xlim = NA,
  1002. ylim = NA,
  1003. axisTitle.size = NA,
  1004. axisText.size = NA,
  1005. title.size = NA,
  1006. lineSize = 1,
  1007. shapePoint = 21,
  1008. sizePoint = 1,
  1009. ...
  1010. ){
  1012. # Prepare kinetic parameters according to scVelo convention
  1013. # scVelo uses beta_scaled = fit_beta * fit_scaling
  1014. # alpha is c(0, fit_alpha) so that alpha[1] = 0 (off), alpha[2] = fit_alpha (on)
  1015. beta <- fit_beta * fit_scaling
  1016. alpha <- c(0, fit_alpha)
  1018. # Compute the initial points of the on and of the off phase
  1019. u0_off <- u0(k = 0, alpha = alpha, beta = beta)
  1020. u0_on <- u0(k = 2, alpha = alpha, beta = beta)
  1021. s0_off <- s0(k = 0, alpha = alpha, beta = beta, gamma = fit_gamma)
  1022. s0_on <- s0(k = 2, alpha = alpha, beta = beta, gamma = fit_gamma)
  1024. # Compute the points lying on the on branch of the dynamic. We assume here the switching point is infinite, in order to plot the potential behavior up to the upper steady state.
  1025. t_on_seq <- seq(fit_t0_on, 500, 0.1)
  1026. u_plot <- u(t_on_seq, t0_off = Inf, t0_on = 0, u0_off = u0_off, u0_on = 0, k = rep(2, length(t_on_seq)), alpha, beta)
  1027. s_plot <- s(t_on_seq, t0_off = Inf, t0_on = 0, u0_off = u0_off, u0_on = 0, s0_off = s0_off, s0_on = 0, k = rep(2, length(t_on_seq)), alpha, beta, fit_gamma)
  1029. # Apply scaling and translation to match scVelo plotting
  1030. u_plot <- u_plot * fit_scaling + fit_u0_offset
  1031. s_plot <- s_plot + fit_s0_offset
  1033. df_on <- data.frame(s_plot = s_plot, u_plot = u_plot)
  1035. # Title for the plot
  1036. main <- paste("Gene", g)
  1037. if(length(unique(subGrLabels))== 1){
  1038. main <- paste(main, ", typeCell", unique(subGrLabels))
  1039. }
  1041. # Initialize or extend ggplot with the ON branch (dashed line)
  1042. if(!add){
  1043. gg <- ggplot(data = df_on, aes(x = s_plot, y = u_plot)) +
  1044. geom_path(lty = 2, col = colDyn, linewidth = lineSize) +
  1045. xlab("s") +
  1046. ylab("u") +
  1047. labs(title = main)
  1049. # adjust graphical parameters
  1050. if(is.numeric(axisText.size)){
  1051. gg <- gg + theme(axis.text = element_text(size = axisText.size))
  1052. }
  1053. if(is.numeric(axisTitle.size)){
  1054. gg <- gg + theme(axis.title = element_text(size = axisTitle.size))
  1055. }
  1056. if(is.numeric(title.size)){
  1057. gg <- gg + theme(plot.title = element_text(hjust = 0.5, size = title.size))
  1058. }
  1059. if(is.numeric(xlim)){
  1060. gg <- gg + xlim(0, xlim)
  1061. }
  1062. if(is.numeric(ylim)){
  1063. gg <- gg + ylim(0, ylim)
  1064. }
  1065. }else{
  1066. gg <- gg +
  1067. geom_path(data = df_on, aes(x = s_plot, y = u_plot), lty = 2, col = colDyn, linewidth = lineSize)
  1068. }
  1070. # Compute the points lying on the off branch of the dynamic, describing the potential behavior from to the upper steady state.
  1071. t_off_seq <- seq(0, 5000, 0.05)
  1072. u_plot <- u(t_off_seq, t0_off = 0, t0_on = 0, u0_off = u0_off, u0_on = u0_on, k = rep(0, length(t_off_seq)), alpha, beta)
  1073. s_plot <- s(t_off_seq, t0_off = 0, t0_on = 0, u0_off = u0_off, u0_on = u0_on, s0_off = s0_off, s0_on = s0_on, k = rep(0, length(t_off_seq)), alpha, beta, fit_gamma)
  1075. # Apply scaling and translation to match scVelo plotting
  1076. u_plot <- u_plot * fit_scaling + fit_u0_offset
  1077. s_plot <- s_plot + fit_s0_offset
  1079. df_off <- data.frame(s_plot = s_plot, u_plot = u_plot)
  1080. gg <- gg + geom_path(data = df_off, aes(x = s_plot, y = u_plot), lty = 2, col = colDyn, linewidth = lineSize)
  1082. # If a finite t0_off_real is provided, compute the transient branches:
  1083. # - ON branch from fit_t0_on to fit_t0_off (solid line)
  1084. # - OFF branch from fit_t0_off to fit_t0_off + 5000 (solid line)
  1085. # Else (t0_off_real == Inf) we reuse the previously computed branch(s).
  1086. if(fit_t0_off != Inf){
  1087. # compute initial conditions at the switching times
  1088. u0_off <- u0(t0_off = fit_t0_off, k = 0, alpha = alpha, beta = beta)
  1089. u0_on <- u0(t0_on = fit_t0_on, k = 2, alpha = alpha, beta = beta)
  1090. s0_off <- s0(t0_off = fit_t0_off, k = 0, alpha = alpha, beta = beta, gamma = fit_gamma)
  1091. s0_on <- s0(t0_on = fit_t0_on, k = 2, alpha = alpha, beta = beta, gamma = fit_gamma)
  1093. # ON transient: from fit_t0_on up to fit_t0_off (solid line)
  1094. t_on_trans <- seq(fit_t0_on, fit_t0_off, 0.01)
  1095. u_plot <- u(t_on_trans, t0_off = fit_t0_off, t0_on = fit_t0_on, u0_off = u0_off, u0_on = u0_on, k = rep(2, length(t_on_trans)), alpha, beta)
  1096. s_plot <- s(t_on_trans, t0_off = fit_t0_off, t0_on = fit_t0_on, u0_off = u0_off, u0_on = u0_on, s0_off = s0_off, s0_on = s0_on, k = rep(2, length(t_on_trans)), alpha, beta, fit_gamma)
  1097. # Apply scaling and translation to match scVelo plotting
  1098. u_plot <- u_plot * fit_scaling + fit_u0_offset
  1099. s_plot <- s_plot + fit_s0_offset
  1101. df_on_trans <- data.frame(s_plot = s_plot, u_plot = u_plot)
  1102. gg <- gg + geom_path(data = df_on_trans, aes(x = s_plot, y = u_plot), lty = 1, col = colDyn, linewidth = lineSize)
  1104. # OFF transient: from fit_t0_off onward (solid line)
  1105. t_off_trans <- seq(fit_t0_off, fit_t0_off + 5000, 0.01)
  1106. u_plot <- u(t_off_trans, t0_off = fit_t0_off, t0_on = fit_t0_on, u0_off = u0_off, u0_on = u0_on, k = rep(0, length(t_off_trans)), alpha, beta)
  1107. s_plot <- s(t_off_trans, t0_off = fit_t0_off, t0_on = fit_t0_on, u0_off = u0_off, u0_on = u0_on, s0_off = s0_off, s0_on = s0_on, k = rep(0, length(t_off_trans)), alpha, beta, fit_gamma)
  1108. # Apply scaling and translation to match scVelo plotting
  1109. u_plot = u_plot * fit_scaling + fit_u0_offset
  1110. s_plot = s_plot + fit_s0_offset
  1112. df_off_trans <- data.frame(s_plot = s_plot, u_plot = u_plot)
  1113. gg <- gg +
  1114. geom_path(data = df_off_trans, aes(x = s_plot, y = u_plot), lty = 1, col = colDyn, linewidth = lineSize)
  1115. }else{
  1116. # If infinite off-time: overlay the earlier computed branches.
  1117. gg <- gg +
  1118. geom_path(data = df_on, aes(x = s_plot, y = u_plot), lty = 1, col = colDyn, linewidth = lineSize) +
  1119. geom_path(data = df, aes(x = s_plot, y = u_plot), lty = 1, col = colDyn, linewidth = lineSize)
  1120. }
  1122. # If fit_t is provided (no NAs) add points (pos_s, pos_u)
  1123. # Colors are auto-chosen if colCell is NA.
  1124. n_cells <- length(fit_t)
  1125. pos_u <- rep(NA, n_cells)
  1126. pos_s <- rep(NA, n_cells)
  1127. k <- rep(0, n_cells)
  1128. k[which(fit_t < fit_t0_off)] <- 2
  1129. pos_u <- u(fit_t, t0_off = fit_t0_off, t0_on = 0, u0_off = u0_off, u0_on = u0_on, k = k, alpha = alpha, beta = beta)
  1130. pos_s <- s(fit_t, t0_off = fit_t0_off, t0_on = 0, u0_off = u0_off, u0_on = u0_on, s0_off = s0_off, s0_on = s0_on, k = k, alpha = alpha, beta = beta, gamma = fit_gamma)
  1131. # Apply scaling and translation to match scVelo plotting
  1132. pos_u <- pos_u * fit_scaling + fit_u0_offset
  1133. pos_s <- pos_s + fit_s0_offset
  1135. df <- data.frame(s = pos_s, u = pos_u)
  1137. # Choose palette if colCell not provided
  1138. subGrC <- unique(subGrLabels)
  1139. if(sum(is.na(colCell)) == 1){
  1140. colCell <- brewer.pal(length(subGrC), "Spectral")
  1141. }
  1142. color <- c()
  1143. for(i in 1:length(subGrC)){
  1144. x <- subGrC[i]
  1145. color <- c(color, rep(colCell[i], length(which(subGrLabels == x))))
  1146. }
  1147. if(nrow(df) == length(color)){
  1148. df$colCell <- color
  1149. }else{
  1150. df$colCell <- color[1:nrow(df)]
  1151. }
  1152. # Add points to the plot
  1153. gg <- gg + geom_point(data = df, aes(x = s, y = u, fill = colCell), fill = df$colCell, col = colDyn, size = sizePoint, shape = shapePoint)
  1155. return(gg)
  1156. }
  1158. # -------------------------------------------------------------------------
  1159. # Function: plot_GeneDynamic_withNotes
  1160. # Purpose:
  1161. # Draw the s-vs-u phase diagram for a single gene with annotation for Figure 1
  1162. #
  1163. # Arguments:
  1164. # t0_off : numeric scalar; time of OFF switch (use Inf if no switch before SS)
  1165. # u0_off : numeric scalar; u-coordinate of the switching point (u^omega)
  1166. # s0_off : numeric scalar; s-coordinate of the switching point (s^omega)
  1167. # t0_on : numeric scalar; time when ON phase starts (default 0)
  1168. # alpha : numeric vector length 2 with transcription rates: c(alpha_off, alpha_on)
  1169. # beta : numeric scalar for splicing rate
  1170. # gamma : numeric scalar fore degradation rate
  1171. # r : optional label used in title (default NA)
  1172. # g : gene identifier (used in title; may be NA)
  1173. #
  1174. # Returns:
  1175. # A ggplot object (s vs u) with annotations.
  1176. # -------------------------------------------------------------------------
  1177. plot_GeneDynamic_withNotes <- function(
  1178. t0_off,
  1179. u0_off,
  1180. s0_off,
  1181. t0_on = 0,
  1182. alpha,
  1183. beta,
  1184. gamma,
  1185. r = NA,
  1186. g
  1187. ){
  1189. # compute the coordinate of the steady states
  1190. u_SS_off <- alpha[1]/beta
  1191. u_SS_on <- alpha[2]/beta
  1192. s_SS_off <- alpha[1]/gamma
  1193. s_SS_on <- alpha[2]/gamma
  1195. # Compute the points lying on the on branch of the dynamic. We assume here the switching point is infinite, in order to plot the potential behavior up to the upper steady state.
  1196. u_plotON <- seq(u_SS_off, u_SS_on, 0.01)
  1197. s_plotON <- s_di_u(u_plotON, u_SS_on, s_SS_on, u_SS_off, s_SS_off, k = 2, alpha, beta, gamma)
  1199. # Compute the points lying on the off branch of the dynamic, describing the potential behavior from to the upper steady state.
  1200. u_plotOFF <- seq(u_SS_off, u_SS_on, 0.01)
  1201. s_plotOFF <- s_di_u(u_plotOFF, u_SS_on, s_SS_on, u_SS_off, s_SS_off, k = 0, alpha, beta, gamma)
  1203. df <- data.frame(uON_SS = as.vector(u_plotON), sON_SS = as.vector(s_plotON), uOFF_SS = as.vector(u_plotOFF), sOFF_SS = as.vector(s_plotOFF), gene = g)
  1204. # Plot on and off branches
  1205. pl <- ggplot(df) +
  1206. geom_path(aes(x = sON_SS, y = uON_SS, color = "red"), linetype = "dotted", linewidth = 3) +
  1207. geom_path(aes(x = sOFF_SS, y = uOFF_SS, color = "blue"), linetype = "dotted", linewidth = 3)
  1209. # If a finite t0_off_real is provided, compute the transient branches.
  1210. if(t0_off != Inf){
  1211. # ON transient: from the lower steady state to the switching point.
  1212. u_plotON_switch <- seq(u_SS_off, u0_off, 0.001)
  1213. s_plotON_switch <- s_di_u(u_plotON_switch, u0_off, s0_off, u_SS_off, s_SS_off, k = 2, alpha, beta, gamma)
  1215. # OFF transient: from the switching point to the lower steady state
  1216. u_plotOFF_switch <- seq(u_SS_off, u0_off, 0.001)
  1217. s_plotOFF_switch <- s_di_u(u_plotOFF_switch, u0_off, s0_off, u_SS_off, s_SS_off, k = 0, alpha, beta, gamma)
  1219. dfSwitchON <- data.frame(uON_switch = u_plotON_switch, sON_switch = s_plotON_switch)
  1220. dfSwitchOFF <- data.frame(uOFF_switch = u_plotOFF_switch, sOFF_switch = s_plotOFF_switch)
  1222. maxSoff <- which.max(s_plotOFF_switch)
  1224. # add the two switching branches
  1225. pl <- pl +
  1226. geom_path(data = dfSwitchON, aes(x = sON_switch, y = uON_switch, color = "red"), linewidth = 3) +
  1227. geom_path(data = dfSwitchOFF, aes(x = sOFF_switch, y = uOFF_switch, color = "blue"), linewidth = 3)
  1229. }else{
  1230. # If infinite off-time: overlay the earlier computed branches.
  1231. pl <- pl +
  1232. geom_path(aes(x = sON_SS, y = uON_SS, color = "red"), linewidth = 3) +
  1233. geom_path(aes(x = sOFF_SS, y = uOFF_SS, color = "blue"), linewidth = 3)
  1234. }
  1236. # Title for the plot
  1237. if(is.na(g)){
  1238. g <- "g"
  1239. }
  1240. main <- ""
  1241. if(!is.na(r)){
  1242. main <- paste(main, " for group r")
  1243. }
  1245. # Annotations: steady-state and switching point labels
  1246. pl <- pl +
  1247. theme(legend.position = "none") +
  1248. annotate(geom = "text", y = alpha[1]*0.88, x = (alpha[1]/gamma)*0.9, label = TeX("$SS^{off}$", output = "character"), size = 35, parse = TRUE, family = "serif") +
  1249. annotate(geom = "text", y = alpha[2]*1.05, x = (alpha[2]/gamma)*1.02, label = TeX("$SS^{on}$", output = "character"), size = 35, parse = TRUE, family = "serif") +
  1250. annotate(geom = "text", y = u0_off*1.05, x = s0_off*0.92,
  1251. label =TeX(r"($(s^{omega}, u^{omega})$)", output = "character"),
  1252. size = 35, parse = TRUE, family = "serif") +
  1253. annotate(geom = "text", x = min(df$sON_SS) + min(df$sON_SS)*0.45, y = min(df$uON_SS) + (max(df$uON_SS) - min(df$uON_SS))/2, label = "Induction", size = 35, angle = 62, family = "serif") +
  1254. annotate(geom = "text", x = max(df$sON_SS)*0.82, y = min(df$uON_SS) + (max(df$uON_SS) - min(df$uON_SS))/2 -0.1, label = "Repression", size = 35, angle = 56, family = "serif") +
  1255. theme(plot.title = element_text(family = "serif", size=70, hjust = 0.5)) +
  1256. labs(x = "", y = "")
  1258. # add arrows for direction of transient induction branch
  1259. ind_on_trans <- which.min(abs(dfSwitchON$uON_switch - (min(dfSwitchON$uON_switch) + (max(dfSwitchON$uON_switch) - min(dfSwitchON$uON_switch))/2)))
  1260. arr_on_u1 <- dfSwitchON$uON_switch[ind_on_trans]
  1261. arr_on_s1 <- dfSwitchON$sON_switch[which(dfSwitchON$uON_switch == arr_on_u1)[1]]
  1262. arr_on_u2 <- dfSwitchON$uON_switch[ind_on_trans + 1]
  1263. arr_on_s2 <- dfSwitchON$sON_switch[which(dfSwitchON$uON_switch == arr_on_u2)[1]]
  1265. pl <- pl +
  1266. geom_segment(aes(color = "red"), x = arr_on_s1, y = arr_on_u1, xend = arr_on_s2, yend = arr_on_u2, arrow = arrow( length = unit(0.3, "inches")), size = 3)
  1268. # add arrows for direction of transient repressive branch
  1269. ind_off_trans <- which.min(abs(dfSwitchOFF$uOFF_switch - (min(dfSwitchOFF$uOFF_switch) + (max(dfSwitchOFF$uOFF_switch) - min(dfSwitchOFF$uOFF_switch))/2)))
  1270. arr_off_u1 <- dfSwitchOFF$uOFF_switch[ind_off_trans]
  1271. arr_off_s1 <- dfSwitchOFF$sOFF_switch[which(dfSwitchOFF$uOFF_switch == arr_off_u1)[1]]
  1272. arr_off_u2 <- dfSwitchOFF$uOFF_switch[ind_off_trans + 1]
  1273. arr_off_s2 <- dfSwitchOFF$sOFF_switch[which(dfSwitchOFF$uOFF_switch == arr_off_u2)[1]]
  1275. pl <- pl + geom_segment(aes(color = "blue"), x = arr_off_s2, y = arr_off_u2, xend = arr_off_s1, yend = arr_off_u1, arrow = arrow( length = unit(0.3, "inches")), size = 3)
  1277. # add arrows for direction of not-transient repressive branch
  1278. ind_off <- which.min(abs(df$uOFF_SS - (min(df$uOFF_SS) + (max(df$uOFF_SS) - min(df$uOFF_SS))/2)))
  1279. arr_off_u1 <- df$uOFF_SS[ind_off]
  1280. arr_off_s1 <- df$sOFF_SS[which(df$uOFF_SS == arr_off_u1)[1]]
  1281. arr_off_u2 <- df$uOFF_SS[ind_off + 5]
  1282. arr_off_s2 <- df$sOFF[which(df$uOFF_SS == arr_off_u2)[1]]
  1284. pl <- pl + geom_segment(aes(colour = "blue"), x = arr_off_s2, y = arr_off_u2, xend = arr_off_s1, yend = arr_off_u1, arrow = arrow(length = unit(0.3, "inches")), size = 3)
  1286. pl <- pl +
  1287. coord_cartesian(xlim = c(min(rbind(df$sON_SS, df$sOFF_SS)) - 0.5, max(rbind(df$sON_SS, df$sOFF_SS)) + 0.5), ylim = c(min(rbind(df$uON_SS, df$uOFF_SS)) - 0.2, max(rbind(df$uON_SS, df$uOFF_SS)) + 0.2)) + theme(axis.ticks.y = element_blank(), axis.ticks.x = element_blank(), axis.text.y = element_blank(), axis.text.x = element_blank())
  1289. # draw cartesian axes on the side
  1290. pl <- pl + geom_segment(aes(x=3.7, y=0.85, xend=4.5, yend=0.85), arrow = arrow(length=unit(.5, 'cm')), color='black', linewidth=2) +
  1291. geom_segment(aes(x=3.7, y=0.846, xend=3.7, yend=1.35), arrow = arrow(length=unit(.5, 'cm')), color='black', linewidth=2) +
  1292. annotate(geom = "text", y = 0.78, x = 4.1, label = TeX("$s$", output = "character"), size = 30, parse = TRUE, family = "serif") +
  1293. annotate(geom = "text", y = 1.046, x = 3.55, label = TeX("$u$", output = "character"), size = 30, parse = TRUE, family = "serif")
  1295. return(pl)
  1296. }

functions.R at commit 73b6f0b, no license · at the source

Overview

Authors: Jiachen Li1, Zhe Wang1, Hong-Bin Shen2, Ye Yuan1,2
ORCID iDs: Jiachen Li, Ye Yuan
  1. State Key Laboratory of Biopharmaceutical Preparation and Delivery, Institute of Process Engineering, Chinese Academy of Sciences Beijing China
  2. Institute of Image Processing and Pattern Recognition, Shanghai Jiao Tong University, and Key Laboratory of System Control and Information Processing, Ministry of Education of China Shanghai China
Institutions: Chinese Academy of Sciences (China); Shanghai Jiao Tong University (China)
Journal: eLife, volume 14, article RP108950
Dates: published online 15 September 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.7554/elife.108950 · PMID 42742132 · PMCID PMC13577664 · OpenAlex W4416668627
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: genetics / omics (modality), mouse (organism), cellular / molecular (subfield)
Methods: Smoothing, state filtering, decompositions, Machine learning, Connectivity
Keywords: RNA velocity, cell trajectory, scRNA-seq analysis, Mouse
MeSH: Gene Expression Regulation*, RNA*, RNA Splicing*, Sequence Analysis, RNA*, Single-Cell Analysis*, Transcription, Genetic*, Animals, Mice, Single-Cell Gene Expression Analysis (* major topic)
Journal subjects: Computational and Systems Biology
Topic: Single-cell and spatial transcriptomics (Molecular Biology, Biochemistry, Genetics and Molecular Biology), according to OpenAlex
Funding: National Key Research and Development Program of China (2023YFF1204500); National Natural Science Foundation of China (62503452)
Citations: not cited yet (Europe PMC); 64 references in the paper

Abstract

RNA velocity approaches fit gene dynamics and infer cell fate by modeling the splicing process using single-cell RNA sequencing (scRNA-seq) data. However, due to the short time scale of splicing, high noise, and large complexity of data, existing RNA velocity methods often fail to precisely capture the complex velocity dynamics for individual genes and single cells, which makes their downstream analysis less reliable and less robust. We propose TSvelo, a comprehensive RNA velocity mathematics framework that can model the cascade of gene regulation, Transcription and Splicing using highly interpretable neural ordinary differential equations. TSvelo can precisely capture the transcription–unspliced–spliced 3D dynamics of all genes simultaneously, infer unified latent time shared by genes within a single cell, and be applied to multi-lineage datasets. Experiments on six scRNA-seq datasets, including two multi-lineage datasets, demonstrate TSvelo’s superiority.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repositories

Its files are read in the Code ↔ Paper reader above, with 8 matches between paragraphs and lines of code.

lijc0804/TSvelo

License: MIT
State: the link answers, verified on 26 September 2026
Evidence: files inventoried
Commit: 6e0ab8264e02aa511576278219fd7c1c65efec80, 8 September 2026
Languages: Python (7), Jupyter (2)
Size: 15 files, 9 scripts
Software Heritage: not archived
Found in: “Code availability statement”
Holds: README, license file, 2 notebooks
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: anndata (9 files), Matplotlib (9 files), NumPy (9 files), pandas (9 files), Scanpy (9 files), scVelo (4 files), SciPy (3 files), NetworkX (1 file), PyTorch (1 file), seaborn (1 file)
Availability: 1 check, the latest on 26 September 2026: the link answers
  • 26 September 2026: the link answers
11 files
At the source: github.com/lijc0804/TSvelo

elenasabbioni/bayvel_notebooks

License: none: the authors keep all their rights
State: the link answers, verified on 26 September 2026
Evidence: files inventoried
Commit: 73b6f0b05eaeddad1d155d1b1dce7b8f1aa88a4f, 26 October 2025
Languages: R (21), Jupyter (3), Julia (1)
Size: 144 files, 25 scripts
Software Heritage: not archived
Found in: “Data availability”
Holds: README, environment (Manifest.toml, Project.toml), 3 notebooks
Not found: license file, CITATION.cff, tests, continuous integration, documentation
Tools: ggplot2 (9 files), data.table (5 files), NumPy (3 files), pandas (3 files), Scanpy (3 files), SciPy (3 files), scVelo (3 files), anndata (1 file), cowplot (1 file), Distributions.jl (1 file), Matplotlib (1 file)
Availability: 1 check, the latest on 26 September 2026: the link answers
  • 26 September 2026: the link answers
26 files
At the source: github.com/elenasabbioni/BayVel_notebooks/tree/main/realdata/Pancreas/moments/output

Code availability statement

TSvelo is implemented in Python. The source code can be downloaded from the GitHub repository, https://github.com/lijc0804/TSvelo (copy archived at Li, 2026).

Reproduced under the paper's license (CC BY), from the paper cited above.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 34 scripts, each with its path and the digest of its content;
  • 8 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

Data availability

The pancreatic endocrinogenesis dataset comprises the single-cell RNA-seq (10X) data of pancreatic epithelial and Ngn3-Venus fusion cells sampled from mouse embryonic day 15.5, which could be loaded using scVelo's package scvelo.datasets.pancreas(). The gastrulation erythroid dataset, which is selected from the transcriptional profiles of mouse embryos39, which could be loaded using scVelo's package scvelo.datasets.pancreas(). 10x embryonic mouse brain dataset is provided at the 10x website at https://www.10xgenomics.com/resources/datasets/fresh-embryonic-e-18-mouse-brain-5-k-1-standard-1-0-0. The data preprocessed by Multivelo is utilized in this study, (https://multivelo.readthedocs.io/en/latest/MultiVelo_Fig2.html). The dentate gyrus neurogenesis data is available at http://pklab.med.harvard.edu/velocyto/DentateGyrus/DentateGyrus.loom. The LARRY dataset has been shared by pyrovelocity, which could be accessed at https://figshare.com/articles/dataset/larry_invitro_adata_sub_raw_h5ad/20780344. The raw data of Hindbrain (pons) of adolescent mice is from https://pklab.med.harvard.edu/ruslan/velocity/oligos/. The ENCODE TF-target database website: https://maayanlab.cloud/Harmonizome/dataset/ENCODE+Transcription+Factor+Targets. The ChEA TF–target database website: https://maayanlab.cloud/Harmonizome/dataset/CHEA+Transcription+Factor+Targets. The results of BayVel on the pancreas dataset are downloaded from its GitHub page at https://github.com/elenasabbioni/BayVel_notebooks/tree/main/real%20data/Pancreas/moments/output (https://github.com/elenasabbioni/BayVel_notebooks/tree/main/real data/Pancreas/moments/output) (Sabbioni, 2025).

The following previously published datasets were used:

QinQ 2022larry_invitro_adata_sub_raw.h5adfigshare10.6084/m9.figshare.20780344

TritschlerS 2019Comprehensive single cell mRNA profiling reveals a detailed roadmap for pancreatic endocrinogenesisNCBI Gene Expression OmnibusGSE13218810.1242/dev.17384931160421

Pijuan-SalaB GriffithsJ 2018Timecourse single-cell RNAseq of whole mouse embryos harvested between days 6.5 and 8.5 of developmentArrayExpressE-MTAB-6967

10x Genomics 2020Fresh Embryonic E18 Mouse Brain (5k)10x Genomicsfresh-embryonic-e-18-mouse-brain-5-k-1-standard-1-0-0

LinnarssonS 2016RNA-seq analysis of single cells of the oligodendrocyte lineage from nine distinct regions of the anterior-posterior and dorsal-ventral axis of the mouse juvenile central nervous systemNCBI Gene Expression OmnibusGSE75330

WeinrebC Rodriguez-FraticelliA CamargoF KleinAM 2019Lineage tracing on transcriptional landscapes links state to fate during differentiationNCBI Gene Expression OmnibusGSE14080210.1126/science.aaw3381PMC760807431974159

HochgernerH ZeiselA LönnerbergP LinnarssonS 2017Transcriptome analysis of single cells from the mouse dentate gyrusNCBI Gene Expression OmnibusGSE95753

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, pages, dates, 4 authors, 4 keywords, 9 MeSH terms, 2 funders, 60 references.

Cite

This paper

Li, J., Wang, Z., Shen, H.-B., & Yuan, Y. (2026). Comprehensive RNA velocity by modeling the cascade of gene regulation, transcription, and splicing from single-cell RNA sequencing data with TSvelo. eLife, 14, RP108950. https://doi.org/10.7554/elife.108950

BibTeX

@article{li2026comprehensive,
author = {Li, Jiachen and Wang, Zhe and Shen, Hong-Bin and Yuan, Ye},
title = {{Comprehensive RNA velocity by modeling the cascade of gene regulation, transcription, and splicing from single-cell RNA sequencing data with TSvelo}},
journal = {eLife},
year = {2026},
month = sep,
volume = {14},
pages = {RP108950},
publisher = {eLife Sciences Publications, Ltd},
issn = {2050-084X},
doi = {10.7554/elife.108950},
url = {https://doi.org/10.7554/elife.108950},
pmid = {42742132},
pmcid = {PMC13577664}
}

RIS

TY - JOUR
AU - Li, Jiachen
AU - Wang, Zhe
AU - Shen, Hong-Bin
AU - Yuan, Ye
TI - Comprehensive RNA velocity by modeling the cascade of gene regulation, transcription, and splicing from single-cell RNA sequencing data with TSvelo
T2 - eLife
J2 - Elife
PY - 2026
DA - 2026/09/15
VL - 14
SP - RP108950
SN - 2050-084X
PB - eLife Sciences Publications, Ltd
DO - 10.7554/elife.108950
UR - https://doi.org/10.7554/elife.108950
LA - en
ER -

CSL-JSON

{
"id": "10.7554/elife.108950",
"type": "article-journal",
"title": "Comprehensive RNA velocity by modeling the cascade of gene regulation, transcription, and splicing from single-cell RNA sequencing data with TSvelo",
"container-title": "eLife",
"author": [
{
"family": "Li",
"given": "Jiachen"
},
{
"family": "Wang",
"given": "Zhe"
},
{
"family": "Shen",
"given": "Hong-Bin"
},
{
"family": "Yuan",
"given": "Ye"
}
],
"container-title-short": "Elife",
"volume": "14",
"page": "RP108950",
"DOI": "10.7554/elife.108950",
"PMID": "42742132",
"PMCID": "PMC13577664",
"ISSN": "2050-084X",
"publisher": "eLife Sciences Publications, Ltd",
"URL": "https://doi.org/10.7554/elife.108950",
"language": "en",
"issued": {
"date-parts": [
[
2026,
9,
15
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1016/j.crmeth.2026.101342 [code]
Interpretable learning of temporal cellular dynamics from single-cell data.
Journal: Cell reports methods
In common: scVelo, anndata, Scanpy, 7 other tools, genetics / omics, cellular / molecular, 13 references
[2] doi:10.1093/bioinformatics/btag652 [code]
mmVelo: a deep generative model for estimating cell state-dependent dynamics across multiple modalities.
Journal: Bioinformatics (Oxford, England)
In common: scVelo, anndata, Scanpy, 6 other tools, genetics / omics, mouse, 10 references
[3] doi:10.1038/s41467-026-74000-4 [code]
ArchVelo: archetypal velocity modeling for single-cell multi-omic trajectories.
Journal: Nature communications
In common: scVelo, anndata, Scanpy, 5 other tools, genetics / omics, mouse, cellular / molecular, 10 references
[4] doi:10.1038/s41586-026-10629-x [code]
Whole-genome duplication shaped cell-type evolution in the vertebrate brain.
Journal: Nature
In common: Distributions.jl, scVelo, anndata, 9 other tools, genetics / omics, mouse, cellular / molecular, 2 references
[5] doi:10.1038/s41586-026-10490-y [code]
Lineage and organ signals sequentially build organ intrinsic nervous systems.
Journal: Nature
In common: scVelo, anndata, Scanpy, 9 other tools, mouse, 3 references
[6] doi:10.1016/j.xcrm.2026.102766 [code]
A longitudinal single-cell and spatial multiomic atlas of pediatric high-grade glioma.
Journal: Cell reports. Medicine
In common: scVelo, anndata, Scanpy, 9 other tools, genetics / omics, cellular / molecular, 2 references
[7] doi:10.7554/elife.106347 [code]
Esr1-dependent signaling and transcriptional maturation in the medial preoptic area of the hypothalamus shape the development of mating behavior during adolescence.
Journal: eLife
In common: anndata, Scanpy, NetworkX, 7 other tools, genetics / omics, mouse, 3 references
[8] doi:10.1038/s41467-026-71803-3 [code]
Charting the transition from in vitro gliogenesis to the in vivo maturation of human glial progenitor cells transplanted into the hypomyelinated mouse brain.
Journal: Nature communications
In common: anndata, Scanpy, cowplot, 8 other tools, genetics / omics, mouse, cellular / molecular, 2 references
[9] doi:10.1371/journal.pcbi.1014346 [code]
StPedf: Cell trajectory inference of spatial transcriptomics via spatial proximity embedding and spatial density-adaptive fusion.
Journal: PLoS computational biology
In common: anndata, Scanpy, NetworkX, 6 other tools, genetics / omics, 4 references
[10] doi:10.1016/j.celrep.2026.117073 [code]
Single-cell epigenomics uncovers heterochromatin instability and transcription factor dysfunction during mouse brain aging.
Journal: Cell reports
In common: anndata, Scanpy, NetworkX, 8 other tools, genetics / omics, mouse, cellular / molecular, 1 reference

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.