El pseudotiempo proporciona orden, pero no dirección

En S32, reconstruimos el orden de las células mediante el pseudotiempo. Sin embargo, existía una limitación importante: no podíamos determinar si una célula se estaba diferenciando o permaneciendo en su estado actual. El pseudotiempo es una medida estática de distancia, no una dirección dinámica.

En 2018, La Manno y sus colaboradores propusieron una solución sorprendentemente elegante a este problema con el concepto de velocidad del ARN. La idea central es la siguiente: en una célula coexisten simultáneamente el ARNm inmaduro no empalmado y el ARNm maduro completamente empalmado, y la proporción entre ambos indica si el gen se está activando o desactivando. Si un gen se acaba de activar, habrá una gran cantidad de ARNm no empalmado y poca cantidad de ARNm empalmado. Por el contrario, si un gen se está desactivando, la nueva transcripción se detiene, por lo que la cantidad de ARNm no empalmado disminuye y solo queda el ARNm empalmado, que se degrada.

En esta sección, derivaremos matemáticamente cómo esta proporción genera un vector de velocidad para cada célula individual y visualizaremos directamente el campo vectorial utilizando scVelo.

Dinámica del empalme: un par de ecuaciones diferenciales

Consideremos un solo gen. Si denotamos la tasa de transcripción como $\alpha$, la tasa de empalme como $\beta$ y la tasa de degradación como $\gamma$, la cantidad de ARNm no empalmado $u$ y la cantidad de ARNm empalmado $s$ siguen las siguientes ecuaciones diferenciales:

En el estado estacionario, donde $du/dt = 0$ y $ds/dt = 0$, obtenemos $u_{\text{ss}} = \alpha / \beta$ y $s_{\text{ss}} = \alpha / \gamma$. Reorganizando esto, la relación en el estado estacionario es:

Esto representa una línea que pasa por el origen en el plano $(s, u)$. La pendiente es $\gamma / \beta$, y esta línea es lo que llamamos la "línea de estado estacionario".

Por encima y por debajo de la línea de estado estacionario

Si una célula se encuentra exactamente sobre esta línea, significa que la transcripción y la degradación están en equilibrio para ese gen. Sin embargo, las células reales se desvían de esta línea.

• Por encima de la línea (más ARNm no empalmado de lo esperado): El gen está en una fase de inducción: la transcripción acaba de comenzar y el ARNm no empalmado se está acumulando. Dado que $ds/dt > 0$, la cantidad de ARNm empalmado aumentará en el futuro.
• Debajo de la línea (menos unspliced de lo esperado): El gen está en una fase de represión: la transcripción se ha detenido, interrumpiendo el suministro de transcritos no empalmados (unspliced), mientras que los transcritos empalmados (spliced) se están degradando. Dado que $ds/dt < 0$, los niveles de transcritos empalmados (spliced) disminuirán.

El "residuo desde la posición actual hasta el estado estacionario" de cada gen constituye su velocidad. Al sintetizar las velocidades de miles de genes, se obtiene un vector de velocidad en el espacio génico para una sola célula. Al proyectar este vector en las coordenadas UMAP, obtenemos flechas.

Velocyto y scVelo: Dos generaciones

Velocyto (1.ª generación, 2018)

Es el método original de La Manno et al. Separa las lecturas de intrones (unspliced) y de exones (spliced) a partir de archivos BAM, ajusta una línea a los diagramas de dispersión $(s, u)$ por gen y utiliza el residuo como velocidad. Sin embargo, asume que todas las células comparten la misma relación $\gamma/\beta$ (supuesto de estado estacionario), por lo que no captura estados transitorios.

scVelo (2.ª generación, 2020)

Bergen et al. abandonaron el supuesto de estado estacionario e introdujeron un modelo dinámico. Para cada gen, estiman simultáneamente $\alpha$, $\beta$, $\gamma$ y el tiempo latente mediante el algoritmo EM, permitiendo que la tasa de transcripción $\alpha$ varíe según el estado celular. Esto permite modelar explícitamente las transiciones de fase de inducción/represión y proporciona una estimación direccional más precisa.

Práctica: Dibujar flechas con datos de diferenciación endocrina pancreática

Utilizamos los datos del páncreas del tutorial oficial de scVelo (Bastidas-Ponce et al. 2019). El proceso de diferenciación desde células precursoras (Ductal/Ngn3) hacia células alfa/beta/delta/epsilón es ideal para verificar la direccionalidad de la velocidad.

import scvelo as scv

# Carga de datos (conjunto de datos pancreas integrado)
adata = scv.datasets.pancreas()

# Preprocesamiento: filtrado + normalización + log + HVG
scv.pp.filter_and_normalize(adata, min_shared_counts=20, n_top_genes=2000)

# Cálculo de momentos (momentos 1° y 2° basados en vecinos kNN)
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)

# Ajuste del modelo dinámico (alpha, beta, gamma por gen + tiempo latente)
scv.tl.recover_dynamics(adata)
scv.tl.velocity(adata, mode="dynamical")

# Gráfico de velocidad + proyección UMAP
scv.tl.velocity_graph(adata)
scv.pl.velocity_embedding_stream(adata, basis="umap", color="clusters")

velocity_embedding_stream dibuja flujos de flechas que simulan el flujo de agua sobre el espacio UMAP. Se considera que el proceso es exitoso si se observa un flujo que comienza en las células precursoras ductales, pasa por Ngn3 y luego se ramifica hacia los distintos tipos de células endocrinas.

Retratos de fase de genes individuales

# Diagrama de fase (s, u) para genes específicos — discriminación de inducción/represión
scv.pl.velocity(adata, var_names=["Ins2", "Gcg"], color="clusters")

En el diagrama de fase de Ins2 (insulina), los grupos de células beta se concentran en el estado de inducción (parte superior, a lo largo de la línea recta), mientras que en Gcg (glucagón), las células alfa estarían en el estado de inducción. Esta única figura proporciona una respuesta intuitiva a la pregunta: "¿Este gen se está expresando actualmente en esta población celular?".

Tiempo potencial frente a pseudotiempo

scv.tl.latent_time(adata)
scv.pl.scatter(adata, color="latent_time", cmap="gnuplot")

El tiempo latente de scVelo corresponde al eje temporal absoluto estimado por el modelo dinámico, por lo que suele ser más fácil de interpretar biológicamente que el pseudotiempo de difusión de S32. Sin embargo, esto también representa un orden relativo y no un tiempo real (horas/días).

¿De dónde provienen los conteos de unspliced y spliced?

El supuesto fundamental del análisis de velocidad es que se deben poder separar las moléculas unspliced y spliced. Existen dos vías principales:

• Velocyto CLI: Clasifica como unspliced las lecturas mapeadas en regiones intrónicas a partir de archivos BAM (salida de 10x CellRanger). El resultado se genera en el archivo .loom.
• kallisto bustools (kb-python): Distingue directamente el estado de splicing a partir de archivos FASTQ. Funciona sin CellRanger y es más rápido.

# Cuantificación simultánea spliced/unspliced con kb-python
kb count -i index.idx -g t2g.txt -x 10xv3 -o output \
   --workflow lamanno \
   R1.fastq.gz R2.fastq.gz

--workflow lamanno genera automáticamente las referencias de intrones y produce tres matrices: spliced, unspliced e ambiguous. Al insertar estas matrices en las capas de AnnData (adata.layers["spliced"], adata.layers["unspliced"]), se obtiene la entrada para scVelo.

Correspondencia de conceptos

• Derivada: La clave de la velocidad es $ds/dt$, es decir, la tasa de cambio temporal del mRNA spliced. Estimar la derivada de una función continua a partir de datos discretos es la versión bioinformática de la diferenciación numérica.
• Campo vectorial: Cuando se asigna un vector de velocidad a cada célula, se forma un campo vectorial en el espacio génico (o UMAP 2D). Es la misma visualización que los diagramas de flujo en dinámica de fluidos.
• Punto de equilibrio de la EDO: La línea de estado estacionario es el punto de equilibrio del sistema de ecuaciones diferenciales ordinarias (EDO). La inducción/represión son respuestas transitorias alejadas del punto de equilibrio, con una estructura similar al análisis de respuesta transitoria en la teoría de control.

Problemas comunes

• Recuento insuficiente de transcritos no empalmados (unspliced) — Las bibliotecas 3' (10x Chromium 3') tienen baja cobertura de intrones, por lo que los transcritos no empalmados son escasos. Si el filtrado es demasiado estricto, no quedan genes para estimar la velocidad. Pruebe a reducir min_shared_counts o considere usar bibliotecas 5'.
• Esperar una velocidad en todos los genes — Para que la velocidad tenga significado, debe existir un cambio entre inducción y represión. Los genes de mantenimiento celular (housekeeping) siempre activos se sitúan sobre la línea de estado estacionario, por lo que su velocidad es cercana a 0. Esto no es un defecto, sino el estado normal.
• Confundir el efecto de lote con la velocidad — Las diferencias técnicas entre lotes pueden afectar la proporción spliced/unspliced. Debe aplicarse previamente la corrección de lote de S31.
• Sacar conclusiones solo por la dirección de las flechas — Los diagramas de flujo son visualmente impactantes, pero es necesario verificar siempre la confianza de la velocidad (scv.tl.velocity_confidence). Las flechas en áreas con baja confianza son ruido.

Para profundizar más

El texto ha sido reconstruido directamente por BPD. Para una mayor comprensión, utilice los siguientes materiales de referencia:

• Artículo original (RNA velocity): La Manno et al. (2018), RNA velocity of single cells, Nature 560:494. Es la fuente original del concepto de velocidad.
• Artículo original (scVelo): Bergen et al. (2020), Generalizing RNA velocity to transient cell states through dynamical modeling, Nature Biotechnology 38:1408. Es la base del modelado dinámico.
• Libro web gratuito de referencia: Single-cell best practices (sc-best-practices.org) — El capítulo sobre velocidad de ARN aborda en detalle las comparaciones metodológicas y el diagnóstico.
• Desarrollo posterior (UniTVelo): Gao et al. (2022), UniTVelo: temporally unified RNA velocity reinforces single-cell trajectory inference, Nature Communications 13:6586. Mejoró la precisión mediante la estimación unificada del tiempo.
• Desarrollo posterior (veloVI): Gayoso et al. (2024), Deep generative modeling of transcriptional dynamics for RNA velocity analysis in single cells, Nature Methods 21:50. Es un enfoque basado en la inferencia variacional.

La velocidad de ARN asignó una flecha a cada célula. Sin embargo, las células no se diferencian solas; determinan su destino intercambiando señales de ligando-receptor con las células vecinas. En la siguiente sección S34, reconstruiremos esta red de comunicación intercelular utilizando CellChat y CellPhoneDB.