Informe práctica 2 - Métodos para inferencia de la trayectoria en procesos biológicos
Adam Casas
Compilado: 17 de abril del 2021
1 Introducción al pseudotiempo
El pseudotiempo es el ordenamiento basado en la similaridad de las células a lo largo de la trayectoria de un proceso de desarrollo contínuo.
El pseudotiempo es una medida de cuánto progreso ha hecho una célula en un proceso biológico como por ejemplo la diferenciación celular. En muchos procesos biológicos, las células no progresan en perfecta sincronía. En estudios de expresión de single-cell sobre procesos tales como la diferenciación celular, las células estudiadas pueden mostrar diversos grados de desarrollo. Es decir, en una población de células capturadas exactamente a la vez, algunas células pueden estar más desarrolladas, mientras que otras pueden ni haber empezado el proceso en cuestión. Esta asincronía dificulta seriamente el estudio de la secuencia de cambios regulatorios que ocurren cuando las células transicionan de un estado a otro. Seguirle la pista a la expresión en todas las células capturadas a la vez produce un sentido muy comprimido de la cinética de un gen, y la variabilidad aparente de la expresión de ese gen será muy elevada. Al ordenar cada célula de acuerdo a su progreso a lo largo de una trayectoria aprendida, el análisis del pseudotiempo alivia los problemas derivados de la susodicha asincronía. En lugar de rastrear los cambios de la expresión génica en función del tiempo, diversas herramientas rastrean dichos cambios en función del progreso a lo largo de la trayectoria, la cual denominamos “pseudotiempo”. El pseudotiempo es una unidad abstracta de progreso: en realidad es la distancia entre una célula dada y el comienzo de la trayectoria, medida a lo largo del camino más corto. La longitud total de la trayectoria se define en términos de la cantidad total de cambios transcripcionales que sufre una célula mientras se mueve del estado inicial al final.
library("SingleCellExperiment")
library("TSCAN")
library("destiny")
library("SLICER")
library("ouija")
library("scater")
library("ggplot2")
library("ggthemes")
library("ggbeeswarm")
library("corrplot")
# library("M3Drop")
# library("monocle)
set.seed(1)En muchas situaciones, uno estudia un proceso en el cual las células cambian continuamente. Esto incluye, por ejemplo, muchos procesos de diferenciación que tienen lugar durante el desarrollo del organismo: tras un estímulo, las células cambian de un tipo a otro. Idealmente, nos interesaría monitorizar los niveles de expresión de una célula dada a lo largo del tiempo. Desafortunadamente, tal monitorización no es posible en scRNA-seq dado que la célula se lisa cuando el ARN se extrae.
En su lugar, debemos muestrear en múltiples puntos del tiempo y así obtener “fotografías” de los perfiles de expresión génicos. Puesto que algunas células se diferenciarán más rápido que otras, cada fotografía puede contener células en diversas etapas del desarrollo. Usamos métodos estadísticos para ordenar las células a lo largo de una o más trayectorias las cuales representan las trayectorias del desarrollo subyacentes, este ordenamiento se conoce como “pseudotiempo”.
En esta práctica nos centraremos en 4 herramientas para ordenar las células según su grado de desarrollo en el pseudotiempo: TSCAN, SLICER, destiny y Ouija. Para ilustrar su uso, emplearemos el dataset del desarrollo embrionario en ratón de Deng et al., 2014. Dicho dataset cuenta con 268 células provenientes de la fase temprana del desarrollo embrional y extraídas en 10 puntos de tiempo. En este caso no hace falta alinear las células a la trayectoria del pseudotiempo dado que están ya etiquetadas con información sobre la trayectoria embrionaria. Por tanto, las etiquetas son la ground truth o verdad experimental sobre la cual podemos comparar el rendimiento obtenido por los susodichos métodos.
Una revisión más o menos reciente de Cannoodt et al. publicada en 2016 nos aporta un resumen detallado de los distintos algoritmos empleados para la inferencia de trayectoria en scRNA-seq. Dicha revisión cubre:
Monocle - Paquete de Bioconductor
TSCAN - Paquete de Bioconductor
SLICER - Paquete de R, disponible sólo a través de GitHub
Waterfall - Paquete de R (se incluye un script en el suplemento del artículo original de Shin et al, 2015 y en los archivos asociados a este informe)
Mpath - Paquete de R disponible en GitHub (comprimido con formato
tar.gz). La función está documentada, pero el workflow no, ni tiene una vignette.Wishbone - Paquete de Python
SCUBA - Paquete de Matlab
Wanderlust - Paquete de Matlab
SCOUP - Programa en C++ por línea de comandos
Las figuras 1 y 2 del paper de Cannoodt et al. resumen algunas de las características de estos métodos.
Figura 1. Descripción de los diversos métodos de inferencia de trayectoria para datos de scRNA-seq (Figura 2 de Cannoodt et al., 2016).
Figura 2. Caracterización de los métodos de inferencia de trayectoria para scRNA-seq (Figura 3 de Cannoodt et al., 2016).
2 Visualización del dataset en PCA
Al igual que el dataset de Nakamura et al. usado en la práctica 1, el dataset de Deng et al. contiene células embrionarias en los primeros estadios de Carnegie, solo que en este caso las células son de ratón, y no de macaco cangrejero. El dataset de esta práctica se encuentra en el archivo deng-reads.rds.
# Leemos los datos
deng_SCE <- readRDS(file = "./Archivos accesorios/Codigo y datos de las Practicas 1-4/data/deng-reads.rds")
# Reordenamos los niveles, ya que zigoto < 2 células < 4 células, etc...
deng_SCE$cell_type2 <- factor(deng_SCE$cell_type2,
levels = c("zy", "early2cell", "mid2cell", "late2cell",
"4cell", "8cell", "16cell", "earlyblast", "midblast",
"lateblast"))
# Generamos y formateamos una matriz de conteos crudos (i.e. nº de transcritos
# por gen y célula)
raw_deng <- counts(deng_SCE)
# Establecemos el nombre de las células en el objeto `cellLabels`
cellLabels <- deng_SCE$cell_type2
colnames(raw_deng) <- cellLabelsTal como muestra el siguiente gráfico, el PCA muestra bastante bien la estructura de los datos y los separa correctamente en su mayoría, exceptuando aquellas células que corresponden al estadio de blastocisto (early/mid/lateblast), que están mezcladas entre sí. Se aprecia cómo las células están ordenadas en el espacio del PCA en sentido horario, desde la etapa de zigoto (zy) hasta la de blastocisto tardío (lateblast).
# Corremos PCA en un objeto de tipo SingleCellExperiment
deng_SCE <- runPCA(deng_SCE)
# Graficamos el PCA
BiocGenerics::plotPCA(deng_SCE, colour_by = "cell_type2") +
ggtitle("PCA de células embrionarias de ratón") +
theme(plot.title = element_text(hjust = 0.5))
Aquí el PCA nos proporciona un útil punto de referencia para comprobar diversos métodos de inferencia de trayectoria. A diferencia del t-SNE y el UMAP, el PCA podría revelarnos mejor el tiempo de desarrollo, a expensas de producir un clustering peor. Un algoritmo muy sencillo de inferencia podría ser simplemente utilizar las coordenadas en el espacio del PCA de la primera componente principal.
# Guardamos en el bolsillo `PC1` del objeto `SCE` el pseudotiempo calculado para
# las células con el PCA
deng_SCE$PC1 <- reducedDim(x = deng_SCE, type = "PCA")[,1]
# Graficamos el pseudotiempo según las coordenadas de la primera componente
# principal
ggplot(as.data.frame(colData(deng_SCE)),
aes(x = PC1, y = cell_type2, colour = cell_type2)) +
geom_quasirandom(groupOnX = F) + theme_classic() +
theme(plot.title = element_text(hjust = 0.5)) +
ggtitle("Células ordenadas por la 1ªCP") +
xlab("Primera componente principal") + ylab("Etapa de desarrollo") +
scale_colour_calc()
Aplicando dicho algoritmo básico, y tal y como muestra el gráfico anterior, vemos que a la primera dimensión le cuesta un poco diferenciar algunos tipos celulares, pero en general ordena bien las células según su estapa de desarrollo. Nótese que aunque los resultados estén al revés, esto se puede solucionar fácilmente invirtiendo el eje X.
# Invertimos el pseudotiempo (1ª componente principal)
deng_SCE$PC1 <- deng_SCE$PC1 * (-1)
# Graficamos la trayectoria final inferida (corregida)
ggplot(as.data.frame(colData(deng_SCE)),
aes(x = PC1, y = cell_type2, colour = cell_type2)) +
geom_quasirandom(groupOnX = F) + theme_classic() +
theme(plot.title = element_text(hjust = 0.5)) +
ggtitle("Células ordenadas por la 1ªCP (orden correcto)") +
xlab("Primera componente principal") + ylab("Etapa de desarrollo") +
scale_colour_calc()
En las siguientes secciones veremos qué rendimiento desempeñan los algoritmos previamente citados ( i.e. TSCAN y sucedáneos).
3 TSCAN
TSCAN combina el clustering con el análisis de pseudotiempo. El clustering de las células que ejerce el comando exprmclust se basa en mclust (Scrucca et al., 2016), el cual modela varias distribuciones normales. Tras el agrupamiento, se construye un árbol del menor tamaño posible ( minimum-spanning-tree o MST) que conecta los clusters entre sí. La rama principal del árbol será la que conecte el mayor nº de clusters, y es la rama que se usa para determinar el pseudotiempo.
Mostramos a continuación el workflow:
# El comando `preprocess` usa como input la matriz de conteos cruda y devuelve
# la matriz de conteos logarítmica a la vez que filtra genes con poca varianza
# y/o expresados en pocas células. Uso la forma `TSCAN::preprocess` para evitar
# confusión con el comando homónimo del paquete `caret`.
procdeng <- TSCAN::preprocess(raw_deng)
# Cambiamos el nombre de las columnas de tipo celular a índice numérico.
colnames(procdeng) <- 1:ncol(raw_deng)
# Agrupamos las células con `exprmclust`.
dengclust <- exprmclust(procdeng, clusternum = 10)
# Graficamos el clustering en el espacio del PCA generado por `exprmclust`.
plotmclust(dengclust, show_tree = T, show_cell_names = F)
# Computamos el pseudotiempo calculado por TSCAN mediante el comando
# `TSCANorder`. Con el parámetro `orderonly = F`, nos devuelve un dataframe con
# el nombre de la célula, el estadio celular en el que se encuentra y su
# posición en el pseudotiempo, respectivamente.
dengorderTSCAN <- TSCANorder(dengclust, orderonly = F)
# Inicializamos el bolsillo vacío "pseudotime_order_tscan" en el objeto `deng_SCE` y
# lo rellenamos con el pseudotiempo que asigna TSCAN a cada célula.
deng_SCE$pseudotime_order_tscan <- NA
deng_SCE$pseudotime_order_tscan[as.numeric(dengorderTSCAN$sample_name)] <- dengorderTSCAN$PseudotimeNota: Con el siguiente código podemos comprobar qué cluster asigna TSCAN a cada estadio celular. En este caso, el clúster nº 10 corresponde a células en fase “late2cell”.
cellLabels[dengclust$clusterid == 10]## [1] late2cell late2cell late2cell late2cell late2cell late2cell late2cell
## [8] late2cell late2cell late2cell
## 10 Levels: zy early2cell mid2cell late2cell 4cell 8cell 16cell ... lateblastAhora visualizamos la trayectoria inferida por TSCAN y el pseudotiempo asignado a cada estadio celular:
# Graficamos el pseudotiempo final de TSCAN
ggplot(as.data.frame(colData(deng_SCE)),
aes(x = pseudotime_order_tscan, y = cell_type2, colour = cell_type2)) +
geom_quasirandom(groupOnX = F) + theme_classic() +
theme(plot.title = element_text(hjust = 0.5)) +
ggtitle("Células ordenadas según TSCAN") +
xlab("Pseudotiempo") + ylab("Etapa de desarrollo") +
scale_colour_calc()
Al igual que el PCA, TSCAN nos devuelve la trayectoria al revés. Vamos a invertirla.
# Invertimos la trayectoria para ponerla en el sentido correcto
deng_SCE$pseudotime_order_tscan <- -deng_SCE$pseudotime_order_tscan
# Graficamos la trayectoria final inferida (corregida)
ggplot(as.data.frame(colData(deng_SCE)),
aes(x = pseudotime_order_tscan, y = cell_type2, colour = cell_type2)) +
geom_quasirandom(groupOnX = F) + theme_classic() +
theme(plot.title = element_text(hjust = 0.5)) +
xlab("Pseudotiempo") + ylab("Etapa de desarrollo") +
scale_colour_calc() + labs(title = "Células ordenadas según TSCAN (orden correcto)")
Observamos que la trayectoria inferida por TSCAN es similar a la obtenida en el PCA, y por tanto no podemos asegurar que TSCAN funcione mejor que el ordenamiento que hicimos mediante las coordenadas de la 1ª componente principal. Ello se debe a que TSCAN es un método de inferencia de trayectoria basado en el PCA.
Ejercicio 2.1: Compara los resultados obtenidos al variar el parámetro clusternum (nº de clusters) del comando exprmclust
Solución: en este dataset, TSCAN infiere bien en su mayoría la trayectoria del desarrollo embrionario del ratón si configuramos el parámetro clusternum entre 3 y 8.
for (i in c(2:9)) {
# clustering buscando 2-9 clusters
dengclust <- exprmclust(procdeng, clusternum = i)
# Computamos el pseudotiempo
dengorderTSCAN <- TSCANorder(dengclust, orderonly = F)
# Añadimos al objeto `SCE` el pseudotiempo calculado por `TSCAN`
deng_SCE$exercise_tscan[as.numeric(dengorderTSCAN$sample_name)] <- dengorderTSCAN$Pseudotime
# Graficamos la trayectoria final
print(ggplot(as.data.frame(colData(deng_SCE)),
aes(x = exercise_tscan, y = cell_type2, colour = cell_type2)) +
geom_quasirandom(groupOnX = F) + theme_classic() +
theme(plot.title = element_text(hjust = 0.5)) +
ggtitle(paste("Trayectoria inferida con", i, "clusters"))) +
xlab(paste("Pseudotiempo")) + ylab("Etapa de desarrollo") +
scale_colour_calc()
# Tienes que usar explícitamente print(ggplot()) si usas gráficos de ggplot2
# dentro de un bucle for
}
4 Mapas de difusión: destiny
Los mapas de difusión fueron inventados por Ronald Coifman y Stephane Lagon en 2006, y la idea principal tras dichos mapas es asumir que los datos son muestras de un proceso de difusión. Este método infiere la estructura de los datos en un espacio de baja dimensionalidad al estimar los eigenvalues y los eigenvectors (valores y vectores propios) del operador de difusión relacionado con los datos.
10 años después, Angerer et al. aplicaron exitosamente el concepto de los mapas de difusión al análisis de datos scRNA-seq para crear así un paquete de R llamado destiny. Para más información sobre el paquete y sus viñetas, véase la página correspondiente del ICB.
Si en la sección del PCA comentábamos que el pseudotiempo era las coordenadas de la 1ª componente principal, en destiny el pseudotiempo pasa a ser las coordenadas de la 1ª componente de difusión. Los mapas de difusión son no lineales.
# Primero creamos la matriz de conteos en forma logarítmica con el comando
# `logcounts`. Recuerda que usa la pseudocuenta de 1, o sea, log2(x+1), para
# definir el logaritmo cuando x = 0.
log_deng <- logcounts(deng_SCE)
# Formateamos bien los nombres de las células (columnas) en la matriz de conteos
# logarítmica.
colnames(log_deng) <- cellLabels
# Calculamos las componentes de difusión. Hay que transponer la matriz de
# conteos logarítmica.
dm <- DiffusionMap(t(log_deng))
# plot.DiffusionMap(dm)
# Creamos un dataframe _ad interim_ con las 2 primeras componentes de difusión y
# el tipo ce célula.
tmp <- data.frame(DC1 = eigenvectors(dm)[,1],
DC2 = eigenvectors(dm)[,2],
etapa_desarrollo = cellLabels) # = deng_SCE$cell_type2# Graficamos las dos primeras componentes del mapa de difusión
ggplot(tmp, aes(x = DC1, y = DC2, color = etapa_desarrollo)) + geom_point() +
theme_classic() + theme(plot.title = element_text(hjust = 0.5)) +
xlab("Componente de difusión 1") + ylab("Componente de difusión 2") +
scale_color_calc() + ggtitle("Mapa de difusión de destiny")
# Ordenamos las células según sus coordenadas en la primera componente de
# difusión y guardamos el resultado en el objeto `SCE`
deng_SCE$pseudotime_diffusionmap <- rank(eigenvectors(dm)[,1]) # = rank(dm@eigenvectors[,1])
# Graficamos la trayectoria final inferida
ggplot(as.data.frame(colData(deng_SCE)), aes(x = pseudotime_diffusionmap, y = cell_type2, colour = cell_type2)) +
geom_quasirandom(groupOnX = FALSE) + scale_color_calc() + theme_classic() +
theme(plot.title = element_text(hjust = 0.5)) + xlab("Pseudotiempo del mapa de difusión (1ª componente)") +
ylab("Etapa de desarrollo") + ggtitle("Células ordenadas según destiny")
# Invertimos el pseudotiempo
deng_SCE$pseudotime_diffusionmap <- -deng_SCE$pseudotime_diffusionmap
# Graficamos la trayectoria final inferida (corregida)
ggplot(as.data.frame(colData(deng_SCE)), aes(x = pseudotime_diffusionmap, y = cell_type2, colour = cell_type2)) +
geom_quasirandom(groupOnX = FALSE) + scale_color_calc() + theme_classic() +
theme(plot.title = element_text(hjust = 0.5)) + xlab("Pseudotiempo del mapa de difusión (1ª componente)") +
ylab("Etapa de desarrollo") + ggtitle("Células ordenadas según destiny (orden correcto)")
Al igual que TSCAN y el PCA, usar la 1ª componente de difusión a modo de pseudotiempo nos ayuda a ordenar las células tempranas, pero no puede distinguir las células en etapas más avanzadas del desarrollo (embrión con 8 o más células).
5 SLICER
SLICER es un algoritmo para construir trayectorias que describen cambios en la expresión génica durante un proceso biológico secuencial (al igual que Monocle y TSCAN). SLICER está diseñado para capturar cambios de expresión génica no lineales, seleccionar automáticamente aquellos genes que participen en el proceso biológico y detectar branches/caminos alternativos y loops en la trayectoria (Welch et al., 2016). El paquete de SLICER está disponible en GitHub y se puede instalar mediante las utilidades del paquete devtools. SLICER asume que la mayoría de la variación observada en el dataset proviene de un proceso biológico secuencial con capacidad de ramificarse (tomar caminos paralelos), por lo que se debe tener en cuenta si dicha suposición es acertada para el dataset en cuestión.
Para comenzar el workflow, emplearemos la función select_genes de SLICER para seleccionar automáticamente los genes que se usarán para construir la trayectoria celular. Dicha función usa la “neighbourhood variance” para identificar los genes que varían suavemente (en lugar de aquellos que fluctúan aleatoriamente) en las células estudiadas. A continuación, determinaremos mediante el comando select_k el nº de vecinos cercanos o nearest neighbours que devuelva un embedding o embebimiento (construcción/proyección en un espacio de dimensionalidad reducida) lo más parecido a una trayectoria. Finalmente estimaremos el embebimiento localmente lineal (LLE) de las células mediante el comando lle del paquete homónimo.
Téngase en cuenta que:
Una buena práctica para determinar el nº óptimo de
kpara el LLE es probar distintos valores deky examinar visualmente dichas las proyecciones.En algunos casos (sobretodo cuando el dataset cuenta con muchos estados celulares), usar un LLE con más de 2 dimensiones puede mejorar los resultados.
# Para seleccionar los genes en `SLICER`, usamos la matriz de conteos
# logarítmica traspuesta. Ídem para seleccionar el parámetro `k` del LLE
slicer_genes <- select_genes(t(log_deng))
k <- select_k(t(log_deng[slicer_genes,]), kmin = 30, kmax = 60)## finding neighbours
## calculating weights
## computing coordinates
## finding neighbours
## calculating weights
## computing coordinates
## finding neighbours
## calculating weights
## computing coordinates
## finding neighbours
## calculating weights
## computing coordinates
## finding neighbours
## calculating weights
## computing coordinates
## finding neighbours
## calculating weights
## computing coordinates
## finding neighbours
## calculating weights
## computing coordinates# La librería `lle` permite usar el comando homónimo para calcular el LLE o
# Locally Linear Embedding.
library("lle")
slicer_traj_lle <- lle(t(log_deng[slicer_genes,]), m = 2, k)$Y # Importante usar esta notación, si no, no funciona el comando `conn_knn_graph()`## finding neighbours
## calculating weights
## computing coordinates# Añadimos los resultados de la técnica de reducción de la dimensionalidad al
# objeto `SCE`
reducedDim(x = deng_SCE, type = "LLE") <- slicer_traj_lle
# Con el siguiente comando comprobamos que el LLE esté disponible para la
# función `reducedDim()`
reducedDimNames(deng_SCE)## [1] "PCA" "LLE"# Graficamos el espacio de dimensionalidad reducida LLE
plotReducedDim(deng_SCE, dimred = "LLE", colour_by = "cell_type2", ncomponents = 2) +
xlab("Componente de LLE 1") + ylab("Componente de LLE 2") +
ggtitle("Embebimiento localmente lineal (LLE) de las células - SLICER") +
theme(plot.title = element_text(hjust = 0.5))
Ahora que hemos calculado el LLE, podemos construir un grafo basado en k-nearest neighbours (k-NN) que sea completamente conectado ( fully connected, como se suele decir en el ámbito de los perceptrones multicapa) mediante el comando conn_knn_graph. En el grafo a continuación se representan las células como círculos naranjas con sus IDs en azul. Mostramos a continuación el grafo obtenido con 10 vecinos. El parámetro k del grafo se suele establecer en \(5\), aunque a veces variar ligeramente su valor puede mejorar los resultados.
set.seed(43765)
# El input del grafo es el LLE (de k = 60 vecinos)
slicer_traj_graph <- conn_knn_graph(slicer_traj_lle, k = 5)
# Graficamos
plot(slicer_traj_graph, main = "Grafo k-NN completamente conectado de SLICER")
A partir del grafo podemos identificar las células distales (de los extremos), las cuales serán candidatas para el comienzo y el fin de la trayectoria. Esto se lleva a cabo a través de la función find_extreme_cells.
# Encontramos las células distales del grafo
ends <- find_extreme_cells(traj_graph = slicer_traj_graph, embedding = slicer_traj_lle)
start <- ends[1]Tras definir la célula inicial y la final, ambas contenidas en el objeto ends, procedemos a ordenar las células a lo largo de la trayectoria. Para ello, primero ordenamos las células con el comando cell_order y acto seguido detectamos las ramas de la trayectoria con assign_branches (este comando es muy sensible a la elección de los parámetros k del LLE y del grafo, en el sentido de que puede devolver un error fácilmente. Para más información, véase la solución del ejercicio 2.2).
# Ordenamos las células en el pseudotiempo con la función `cell_order()` del
# paquete `SLICER`
pseudotime_order_slicer <- cell_order(slicer_traj_graph, start)
# Detectamos ramas/caminos alternativos
branches <- assign_branches(slicer_traj_graph, start)
# creamos dataframe con las ramas, los estadios celulares y el pseudotiempo
pseudotime_slicer <- data.frame(TimePoint = cellLabels,
pseudotime = NA,
State = branches)
# Añadimos a la columna pseudotime del dataframe anterior el orden de las
# células en el pseudotiempo
pseudotime_slicer$pseudotime[pseudotime_order_slicer] <- 1:length(pseudotime_order_slicer)Hecho esto, ahora podemos comparar la trayectoria inferida por SLICER con las etapas del desarrollo ya conocidas. SLICER no devuelve un pseudotiempo per se, sino simplemente el orden final de las células.
# Añadimos al objeto `SCE` el pseudotiempo del dataframe anterior para comparar
# los resultados con los métodos anteriores y la verdad experimental
deng_SCE$pseudotime_slicer <- pseudotime_slicer$pseudotime
# Graficamos el resultado
ggplot(as.data.frame(colData(deng_SCE)), aes(x = pseudotime_slicer, y = cell_type2, colour = cell_type2)) +
geom_quasirandom(groupOnX = F) + xlab(" \"Pseudotiempo\" (orden celular) ") + ylab("etapa de desarrollo") +
theme_classic() + ggtitle("Células ordenadas según SLICER") + theme(plot.title = element_text(hjust = 0.5))
Al igual que los métodos ya probados hasta ahora, SLICER ordena relativamente bien las primeras etapas del desarrollo, aunque ha invertido el orden de las etapas 16cell y 8cell (e ídem para las etapas zy y early2cell). Pese a ello, logra separar un poco mejor las blastómeras tardías del resto de blastómeras.
Ejercicio 2.2: ¿Cómo cambian los resultados del LLE en función del valor de k? ¿Y qué pasa al cambiar el nº de vecinos cercanos en la función conn_knn_graph()?
Solución:
Manteniendo el parámetro
k = 60del LLE y variando el parámetrokdel comandoconn_knn_graph(), vemos que a mayor nº de vecinos cercanos, la trayectoria se perfila cada vez más en el grafo, si bien el ordenamiento final de las células se mantiene igual. Nótese que para valores muy bajos de vecinos (k = 1~3), el comandoassign_branchesdevuelve el errorError in 1:num_branches : Argumento NA/NaNManteniendo el parámetro
k = 10del grafo y variando el parámetrokdel LLE, vemos que para valores extremos dek, a SLICER le cuesta ordenar las células. SLICER ordena exitosamente las células tempranas con unk = 30, mientras que conk = 55intenta separar las blastómeras (aunque sin mucho éxito). Nótese que para valores muy bajos de vecinos (k < 10), el comandoassign_branchesdevuelve el errorError in if (mean_dist < min_dist) { : valor ausente donde TRUE/FALSE es necesario(Si se observa la función con F2, se aprecia que se omitieron algunos casos en la sección de condicionales).
Grafo con k = 1; LLE con k = 60
Grafo con k = 5; LLE con k = 60
Grafo con k = 30; LLE con k = 60
Orden celular con LLE con k = 10
Orden celular con LLE con k = 30
Orden celular con LLE con k = 50
Orden celular con LLE con k = 55
Orden celular con LLE con k = 100
Orden celular con LLE con k = 150
6 ouija
ouija (viñeta online: http://kieranrcampbell.github.io/ouija/) infiere la trayectoria de manera distinta a como lo han hecho los paquetes previamente presentados. Mientras que el PCA, TSCAN, destiny y SLICER eran métodos no supervisados y no tenían en cuenta las etiquetas asignadas por el investigador, ouija es un algoritmo probabilístico capaz de inferir el pseudotiempo a partir de un pequeño panel de genes marcadores. En concreto, este algoritmo:
Infiere el pseudotiempo a partir de un nº reducido de biomarcadores, llevando la atención del investigador a dichos genes.
Proporciona estimaciones de los parámetros (incluye incertidumbre), lo que permite interpretar patrones de regulación génica.
Incluye un test para hipótesis Bayesiana para encontrar a lo largo de la trayectoria genes regulados antes que otros.
Identifica estados metaestables a lo largo de la trayectoria, i.e. tipos celulares.
En ouija podemos modelar los genes de varias maneras:
Que se sobre o infraexpresen de manera monotónica ( i.e. su expresión a lo largo del tiempo sea exclusivamente al alza o a la baja, respectivamente). Esto se conoce como comportamiento switch-like.
Que se sobre o infraexpresen de manera aguda ( i.e. la expresión del gen alcanza brevemente y de manera rápida un máximo/mínimo). Esto se conoce como comportamiento transitorio.
Por defecto, ouija asume que todos los genes usan el modelo monotónico, aunque el autor asegura que no pasa nada, pues modelar un gen transitorio como uno monotónico tiene un efecto mínimo en los resultados finales.
Proporcionaremos a ouija los siguientes biomarcadores, junto con la etapa del desarrollo en la que se sobreexpresan:
Etapa inicial: Dazl, Rnf17, Sycp3, Nanog, Pou5f1, Fgf8, Egfr, Bmp5, Bmp15
Etapa intermedia: Zscan4b, Foxa1, Prdm14, Sox21
Etapa tardía: Creb3, Gpx4, Krt8, Elf5, Eomes, Cdx2, Tdgf1, Gdf3
En nuestro caso particular con el dataset de Deng podemos hacer un poco de trampas y comprobar que nuestros genes marcadores seleccionados son realmente biomarcadores de las distintas etapas del desarrollo aquí estudiadas.
# Genes de la etapa inicial(excepto Nanog). Comienzan activos y se desactivan a
# medida que progresa el pseudotiempo (de ahí lo de "down").
ouija_markers_down <- c("Dazl", "Rnf17", "Sycp3", "Pou5f1", "Fgf8", "Egfr", "Bmp5", "Bmp15")
# Genes de la etapa intermedia. Comienzan inactivos, se activan un rato y
# vuelven a desactivarse hacia el final de la trayectoria. Presentan una
# activación transitoria
ouija_markers_transient <- c("Zscan4b", "Foxa1", "Prdm14", "Sox21")
# Genes de la etapa tardía. Comienzan inactivos y se activan a medida que
# progresa el pseudotiempo (de ahí lo de "up").
ouija_markers_up <- c("Creb3", "Gpx4", "Krt8", "Elf5", "Eomes", "Cdx2", "Tdgf1", "Gdf3") # Mirar gen Sox2
# Mirar differential activity of the CARM1 histone methyltransferase has been
# shown to determine cell fate decision in 4-cell mouse embryos
ouija_markers <- c(ouija_markers_up, ouija_markers_down, ouija_markers_transient)
plotExpression(deng_SCE, features = ouija_markers, x = "cell_type2", colour_by = "cell_type2") +
theme(axis.text.x = element_text(angle = 60, hjust = 1))
Para computar el pseudotiempo, basta con llamar a ouija y pasarle los modelos de genes esperados. Como ya se dijo en el párrafo anterior, el paquete modela por defecto los genes de manera monotónica (cosa que haremos a continuación). El input que espera ouija es cualquiera de las 3 opciones siguientes:
Una matriz de célula(fila) x gen(columna) con valores de expresión génica no negativos.
Una matriz de conteos logarítmica, producto de aplicar el comando
logcounts()a un objeto de tipoSingleCellExperiment.Un objeto de tipo
ExpressionSet.
Podemos añadir información conocida a priori respecto a la sobre o infraexpresión de los genes implicados en el proceso biológico estudiado, además de incluir conocimiento previo sobre el tipo de modelo que sigue cada gen (monotónico/transitorio).
Podemos contruir el modelo de ouija con:
El método de Monte Carlo Hamiltoniano o Monte Carlo Híbrido ( Hamiltonian Monte Carlo o HMC). Consiste en inferencia completa con el método de Monte Carlo basado en Cadenas de Markov ( Monte Carlo Markov Chain o MCMC) en la cual se usa la información del gradiente de la log-posterior (logaritmo de la probabilidad a posteriori) para guiar el paseo aleatorio en el espacio del parámetro.
El método de Automatic Differentiation Variational Bayes, ADVI o simplemente VI. Consiste en inferencia aproximada en la cual se minimiza la divergencia Kullback-Leibler o KL (una medida de distancia asimétrica y distinta de la euclídea; no es lo mismo medir la distancia forwards KL que la distancia backwards KL) a una distribución aproximada.
En general, el método de HMC infere de manera más precisa la trayectoria y devuelve una varianza posterior aproximadamente correcta para todos los parámetros. No obstante, el ADVI es significativamente más rápido que HMC y aunque pueda infraestimar la varianza a posteriori, el autor de ouija comenta que en ocasiones ADVI puede devolver resultados igual de buenos que los obtenidos con HMC, en cuanto a descubrimiento de pseudotiempos posteriores se refiere.
Para ayudar al modelo de ouija, vamos a proporcionarle información a priori sobre el comportamiento de nuestros genes sobre e infraexpresados. Al seleccionar la switch strength en +-10 (-10 para genes infraexpresados y +10 para genes sobreexpresados) y la prior strength standard deviation en 0.5, le estamos indicando al modelo que estamos muy seguros sobre el comportamiento esperado de estos genes a lo largo del proceso biológico estudiado (desarrollo embrionario en este caso).
# Tipo de respuesta de cada gen
response_type <- c(rep("switch", length(ouija_markers_down) + length(ouija_markers_up)),
rep("transient", length(ouija_markers_transient)))
# c(rep("switch", length(ouija_markers_down) + length(ouija_markers_up))) equivale a
# c(rep("switch", length(ouija_markers_down)), rep("switch, length(ouija_markers_up)))
# Certeza del comportamiento de nuestros biomarcadores
switch_strengths <- c(rep(-10, length(ouija_markers_down)),
rep(10, length(ouija_markers_up)))
# Desviación estándar de dicha certeza
switch_strength_sd <- c(rep(0.5, length(ouija_markers_down)),
rep(0.5, length(ouija_markers_up)))
# El comando `ouija()` recomienda usar la siguiente línea de código para
# paralelizar en local su computación:
options(mc.cores = parallel::detectCores())
oui_vb <- ouija(x = deng_SCE[ouija_markers,],
response_type = response_type,
switch_strengths = switch_strengths,
switch_strength_sd = switch_strength_sd,
inference_type = "vb", # "vb" = ADVI
single_cell_experiment_assay = "logcounts")## Chain 1: ------------------------------------------------------------
## Chain 1: EXPERIMENTAL ALGORITHM:
## Chain 1: This procedure has not been thoroughly tested and may be unstable
## Chain 1: or buggy. The interface is subject to change.
## Chain 1: ------------------------------------------------------------
## Chain 1:
## Chain 1:
## Chain 1:
## Chain 1: Gradient evaluation took 0.006 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 60 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Begin eta adaptation.
## Chain 1: Iteration: 1 / 250 [ 0%] (Adaptation)
## Chain 1: Iteration: 50 / 250 [ 20%] (Adaptation)
## Chain 1: Iteration: 100 / 250 [ 40%] (Adaptation)
## Chain 1: Iteration: 150 / 250 [ 60%] (Adaptation)
## Chain 1: Iteration: 200 / 250 [ 80%] (Adaptation)
## Chain 1: Success! Found best value [eta = 1] earlier than expected.
## Chain 1:
## Chain 1: Begin stochastic gradient ascent.
## Chain 1: iter ELBO delta_ELBO_mean delta_ELBO_med notes
## Chain 1: 100 -3774.104 1.000 1.000
## Chain 1: 200 -3717.426 0.508 1.000
## Chain 1: 300 -3691.324 0.341 0.015
## Chain 1: 400 -3676.025 0.257 0.015
## Chain 1: 500 -3577.174 0.211 0.015
## Chain 1: 600 -3524.273 0.178 0.015
## Chain 1: 700 -3529.497 0.153 0.015
## Chain 1: 800 -3522.582 0.134 0.015
## Chain 1: 900 -3526.697 0.119 0.007 MEDIAN ELBO CONVERGED
## Chain 1:
## Chain 1: Drawing a sample of size 1000 from the approximate posterior...
## Chain 1: COMPLETED.print(oui_vb)## A Ouija fit with 268 cells and 20 marker genes
## Inference type: Variational Bayes
## (Gene behaviour) Switch/transient: 16 / 4Podemos graficar la expresión génica a lo largo del pseudotiempo junto con las estimaciones de la función media-MAP ( maximum a posteriori) usando el comando plot_expression() del paquete ouija.
plot_expression(oui_vb) + theme_classic()
Nótese que genes tempranos como Zscan4b (expresado en los estadios de 2 y 4 células) aparecen sobreexpresados al final del pseudotiempo, lo que nos indica que el gráfico tiene el eje X invertido. Vamos a solucionarlo:
plot_expression(oui_vb) + theme_classic() +
scale_x_reverse() + xlab("Ouija pseudotime (corrected)")
Ahora podemos observar en el gráfico la expresión de cada gen a lo largo del proceso de desarrollo embrionario. Por ejemplo, el gen Elf5 comienza a expresarse en las estapas finales (blastocisto).
También podemos visualizar cuándo suceden en la trayectoria los comportamientos de regulación génica, ya sea en forma de switch time o el peak time (para genes monotónicos o transitorios) usando los comandos plot_switch_times y plot_transient_times.
# Graficamos los switch times de los genes monotónicos (switch-like)
plot_switch_times(oui_vb)
# Graficamos los peak times de los genes transitorios
plot_peak_times(oui_vb)
Se aprecia que igual que en el gráfico anterior, la trayectoria está invertida. Además, el gráfico que devuelve plot_switch_times() indica la regulación al revés (El gen Fgf8 se desactiva al comienzo de la trayectoria aunque en el gráfico indique que se activa, ídem para el gen Egfr). Procedemos a corregir ambos gráficos inviertiendo el eje X de ambos gráficos y la leyenda del gráfico plot_switch_times, teniendo en cuenta que:
Los comandos
plot_switch_timesyplot_peak_timesusan internamente el comandocoord_flip(), por lo que para invertir el eje x, tenemos que usar el comandoscale_y_reverse()en lugar descale_x_reverse().Para invertir los valores de la regulación del gráfico de
plot_switch_times, multiplicaremos por \(-1\) los valores contenidos en el objetooui_vb$fit@sim$samples, pues la funciónplot_switch_timesobtiene dichos valores de regulación de ahí.
Compendio de soluciones para corregir los valores de regulación génica de plot_switch_times
Mientras que para invertir el eje x de nuestros gráficos, basta con llamar a la función scale_y_reverse (y no scale_x_reverse dado que las funciones de graficado plot_switch_times y plot_peak_times usan internamente coord_flip), para invertir los valores de la regulación génica no es tan fácil.
La función de graficado plot_switch_times obtiene los valores de regulación génica de los objetos k.1-k.16 de la lista oui_vb$fit@sim$samples[[1]] (i.e. oui_vb$fit@sim$samples[[1]]$k.1, oui_vb$fit@sim$samples[[1]]$k.2…), ubicada dentro del objeto oui_vb. A continuación ofrecemos 3 soluciones:
################################################################
############### Solución 1 #########################
################################################################
# A partir de la función original `plot_switch_times()`, definimos una función
# de graficado con las modificaciones necesarias ya incluidas (requiere
# librerías `viridis` para el color, `coda` para el comando `mcmc` y `tibble`
# para el comando ´data_frame´)
plot_switch_times2 <- function (oui) {
Gene <- t0_mean <- lower <- upper <- NULL
vpal <- viridis_pal()(8)
k_trace <- extract(oui$fit, "k")$k
kmean <- colMeans(k_trace)
t0 <- extract(oui$fit, "t0")$t0
t0_means <- colMeans(t0)
t0_interval <- coda::HPDinterval(mcmc(t0))
t0_df <- data_frame(t0_mean = t0_means, lower = t0_interval[,
1], upper = t0_interval[, 2], kmean = -kmean) # Multiplicamos kmean por -1
t0_df$Gene <- colnames(oui$Y[, oui$response_type == "switch"])
t0_df$Gene <- factor(t0_df$Gene, t0_df$Gene[order(t0_means)])
ggplot(t0_df, aes(x = Gene, y = t0_mean, fill = kmean)) +
geom_errorbar(aes(ymin = lower, ymax = upper), color = "grey60",
width = 0.5, alpha = 0.5) + coord_flip() + geom_point(color = "grey50",
shape = 21, size = 3) + ylab("Switch point") + scale_fill_gradient2(name = "Regulation",
low = vpal[1], high = vpal[5]) + scale_color_gradient2(name = "Regulation",
low = vpal[1], high = vpal[5]) + theme(legend.position = "top") + scale_y_reverse() +
theme_classic() # Añadimos `scale_y_reverse()` y `theme_classic()`
}
library(viridis)
library(coda)
library(tibble)
# Graficamos los switch times de los genes monotónicos (switch-like)
plot_switch_times2(oui_vb)
################################################################
############### Solución 2 #########################
################################################################
# Solucion manual. Modificamos manualmente cada uno de los 16 objetos de interés
# de la lista `oui_vb$fit@sim$samples[[1]]`
oui_vb$fit@sim$samples[[1]]$k.1 <- -oui_vb$fit@sim$samples[[1]]$k.1
oui_vb$fit@sim$samples[[1]]$k.2 <- -oui_vb$fit@sim$samples[[1]]$k.2
oui_vb$fit@sim$samples[[1]]$k.3 <- -oui_vb$fit@sim$samples[[1]]$k.3
oui_vb$fit@sim$samples[[1]]$k.4 <- -oui_vb$fit@sim$samples[[1]]$k.4
oui_vb$fit@sim$samples[[1]]$k.5 <- -oui_vb$fit@sim$samples[[1]]$k.5
oui_vb$fit@sim$samples[[1]]$k.6 <- -oui_vb$fit@sim$samples[[1]]$k.6
oui_vb$fit@sim$samples[[1]]$k.7 <- -oui_vb$fit@sim$samples[[1]]$k.7
oui_vb$fit@sim$samples[[1]]$k.8 <- -oui_vb$fit@sim$samples[[1]]$k.8
oui_vb$fit@sim$samples[[1]]$k.9 <- -oui_vb$fit@sim$samples[[1]]$k.9
oui_vb$fit@sim$samples[[1]]$k.10 <- -oui_vb$fit@sim$samples[[1]]$k.10
oui_vb$fit@sim$samples[[1]]$k.11 <- -oui_vb$fit@sim$samples[[1]]$k.11
oui_vb$fit@sim$samples[[1]]$k.12 <- -oui_vb$fit@sim$samples[[1]]$k.12
oui_vb$fit@sim$samples[[1]]$k.13 <- -oui_vb$fit@sim$samples[[1]]$k.13
oui_vb$fit@sim$samples[[1]]$k.14 <- -oui_vb$fit@sim$samples[[1]]$k.14
oui_vb$fit@sim$samples[[1]]$k.15 <- -oui_vb$fit@sim$samples[[1]]$k.15
oui_vb$fit@sim$samples[[1]]$k.16 <- -oui_vb$fit@sim$samples[[1]]$k.16
# Graficamos los switch times de los genes monotónicos (switch-like)
plot_switch_times(oui = oui_vb) + theme_classic() + scale_y_reverse() +
scale_x_discrete(limits = rev)
# Opcionalmente, podemos usar `scale_x_discrete(limits = rev)` para invertir el
# orden del eje x, que muestra factores (ergo no podemos usar `scale_x_reverse`)
################################################################
############### Solución 3 #########################
################################################################
# Solución con bucle for. Automatizamos con un bucle `for` los pasos realizados en
# la solución 2 (los bucles for no son eficientes en R)
# Con `grep` obtenemos los índices de los k.1-k.16 y los visualizamos con
# `print()`.
indices <- grep("k.", names(oui_vb$fit@sim$samples[[1]]))
# Ahora los multiplicamos por -1 para obtener la regulación génica en el sentido
# correcto
for (i in indices) {
oui_vb$fit@sim$samples[[1]][[i]] <- -oui_vb$fit@sim$samples[[1]][[i]]
}
# Graficamos los switch times de los genes monotónicos (switch-like)
plot_switch_times(oui = oui_vb) + theme_classic() + scale_y_reverse()Finalmente, optamos por una 4ª solución: usar lapply, pues es mucho más eficiente que usar bucles for en R (ver siguiente chunk)
# Con `grep` (RegEx) obtenemos los índices de los k.1-k.16 y los visualizamos con
# `print()`
indices <- grep("k.", names(oui_vb$fit@sim$samples[[1]])); print(c(min(indices), max(indices)))## [1] 22 37# Multiplicamos sus valores por -1 para obtener la regulación génica en el
# sentido correcto
oui_vb$fit@sim$samples[[1]][c(22:37)] <- lapply(oui_vb$fit@sim$samples[[1]][c(22:37)], function(x) x*-1)
# Graficamos los switch times de los genes monotónicos (switch-like)
plot_switch_times(oui = oui_vb) + theme_classic() +
theme(legend.position ="top") + scale_y_reverse() +
ylab("Switch point (corrected)")
# Graficamos los peak times de los genes transitorios
plot_peak_times(oui_vb) + theme_classic() +
scale_y_reverse(limits = c(1,0)) + ylab("Peak time (corrected)")
Ahora tanto la trayectoria como los valores de regulación se muestran correctamente, i.e. Zscan4b muestra un pico de expresión al comienzo de la trayectoria; Egfr comienza a infraexpresarse desde el comienzo de la trayectoria; Elf5 comienza a sobreexpresarse en etapas tardías, etc…
Adicionalmente, con ouija podemos detectar estados celulares metaestables usando matrices de consistencia. Los estados celulares metaestables son un sinónimo de tipos celulares, ya que la idea detrás del nombre es que a medida que las células se diferencian durante el proceso biológico, pueden mantener un estado estable en un tramo de la trayectoria, lo que contaría como un tipo celular. Esto puede ser útil para detectar nuevos tipos celulares en procesos de diferenciación más o menos conocidos, como por ejemplo en la diferenciación de las células troncales hematopoyéticas (Ceredig et al., 2009).
Para identificar estos tipos celulares metaestables, primero tenemos que computar la matriz de consistencia con el comando consistency_matrix, y acto seguido agrupar las células con cluster_consistency. Tras ello, podemos computar el pseudotiempo con el comando map_pseudotime y visualizar cómo se suceden dichos estados celulares a lo largo del pseudotiempo. Nótese que comprobé posteriormente que ouija también nos devuelve el pseudotiempo al revés, así que lo corregiremos aquí para ahorrarnos futuros quebraderos de cabeza.
# Visualizamos la matriz de consistencia (opcional)
plot_consistency(oui_vb)
# Guardamos la matriz de consistencia en el objeto `cmo` y agrupamos las células
cmo <- consistency_matrix(oui_vb)
cell_classifications <- cluster_consistency(cmat = cmo, n_clusters = 2:9)
# Hemos detectado 5 estados celulares metaestables
levels(as.factor(cell_classifications))## [1] "1" "2" "3" "4"# Computamos el pseudotiempo en el orden correcto, lo añadimos a un dataframe
# con los estados celulares metaestables y graficamos dicho dataframe
map_pst <- map_pseudotime(oui_vb) * (-1)
ouija_pseudotime <- data.frame(map_pst, cell_classifications)
ggplot(ouija_pseudotime, aes(x = map_pst, y = cell_classifications)) +
geom_point() + xlab("Pseudotiempo MAP") + ylab("Estado celular metaestable") +
theme_classic() + theme(plot.title = element_text(hjust = 0.5)) +
ggtitle("Estados celulares metaestables a lo largo del pseudotiempo (orden correcto)")
Ahora que hemos identificado 4 estados celulares metaestables, procedemos a identificarlos mediante su comparación con la verdad experimental:
# Añadimos el pseudotiempo y los estados celulares metaestables computados en el
# chunk anterior al objeto `SCE` para poder graficarlos
deng_SCE$pseudotime_ouija <- ouija_pseudotime$map_pst
deng_SCE$ouija_cell_class <- ouija_pseudotime$cell_classifications
# Graficamos. Nótese que invertimos el pseudotiempo (`pseudotime_ouija`) en el
# chunk anterior
ggplot(as.data.frame(colData(deng_SCE)),
aes(x = as.factor(ouija_cell_class), y = pseudotime_ouija, colour = cell_type2)) +
geom_boxplot() + coord_flip() + scale_color_tableau() +
theme_classic() + theme(plot.title = element_text(size = 10, hjust = 0.35)) +
ggtitle("Estados celulares metaestables y su correspondencia con la verdad experimental (orden correcto)") +
xlab("Estado celular metaestable") + ylab("Pseudotiempo MAP")
En vista del gráfico, anotaremos esos 4 estados metaestables como “zygote/2 cell”, “4/8/16 cell”, “early/mid blast” y “late blast”.
Finalmente procedemos a visualizar la trayectoria inferida por ouija:
# Graficamos la trayectoria final inferida
ggplot(as.data.frame(colData(deng_SCE)),
aes(x = pseudotime_ouija, y = cell_type2, colour = cell_type2)) +
geom_quasirandom(groupOnX = F) + theme_classic() +
scale_color_calc() + xlab("Pseudotiempo MAP") +
ylab("estado de desarrollo") + ggtitle("Ordenamiento celular según ouija (orden correcto)") +
theme(plot.title = element_text(hjust = 0.5))
ouija ha ordenado bien en general los tipos celulares, aunque cabe destacar que el algoritmo puede ser sensible a la elección de biomarcadores y la información a priori proporcionada al mismo.
Ejercicio 2.3: ¿Cómo cambian los resultados al elegir otros biomarcadores o cambiar las probabilidades a priori?
Solución:
Efectivamente, el algoritmo es sensible a las probabilidades a priori de los genes. Resulta muy beneficial para el algoritmo tener muy claro el funcionamiento de los biomarcadores que le proporciones y elegir en consecuencia unos valores de probabilidades a prioris elevados (estás muy seguro de su funcionamiento) y de signo correcto (si un gen tiende a activarse, ponle una probabilidad a priori > 0, y viceversa para genes que se desactivan). De no elegir correctamente los biomarcadores y sus probabilidades a priori,
ouijano podrá inferir correctamente la trayectoria, ni podrá afirmar con precisión el orden en el que se activan/desactivan algunos genes.Respecto a la selección de otros biomarcadores, parece ser que
ouijapuede seguir infiriendo correctamente la trayectoria al seleccionar otros biomarcadores como los de la familia Hox (genes arquitecto, orquestran el desarrollo embrionario).
Un análisis típico en inferencia de trayectorias es averiguar el orden de regulación de los genes. Por ejemplo, ¿se activa el gen A antes que el gen B? ¿Alcanza el gen C su máximo de expresión antes de que el gen D se desactive?
ouija responde a estas preguntas en términos de un test a pares de hipótesis Bayesiana de si la diferencia en los tiempos de regulación, ya sea switch time o peak time, es significativamente distinta de 0. Esto se comprueba con el comando gene_regulation.
gene_regs <- gene_regulation(oui_vb)
head(gene_regs)Ejercicio 2.4: ¿Qué conclusiones puedes derivar de la tabla de regulación génica?
Solución: La tabla anterior nos confirma mediante el susodicho test a pares la información previamente obtenida en los gráficos de expresión génica a lo largo del pseudotiempo (comandos plot_expression,plot_switch_times y plot_peak_times).
Tanto en el gráfico superior como en el inferior, observamos que Zscan4b alcanza su máximo de expresión antes que Sox21, pero, ¿es esa diferencia significativa o se debe al azar)? Según la tabla, esa diferencia es efectivamente significativa. Podemos afirmar con un 95% de confianza que Zscan4b alcanza su máximo de expresión antes que lo hace Sox21. Un resultado así podría indicarnos que la activación de Sox21 podría depender del producto génico de Zscan4b, y el retardo observado entre ambos picos de expresión se correspondería con el tiempo que tarda en transcribirse Zscan4b, traducirse en el citoplasma, desplazarse al núcleo e inducir una activación detectable de Sox21.
Si comparamos dos genes monotónicos de la etapa temprana, Egfr y Dazl, el primero se desactiva/infraexpresa antes que el segundo. Un resultado de este tipo podría indicar que el silenciamiento de Dazl depende de que Egfr se desactive primero ( i.e. ¿inhibición por producto génico?).
Por otro lado, en el gráfico parece que Egfr y Bmp5 se activan a la vez. Tras consultar la tabla, lo confirmamos con respaldo estadístico. Un resultado de este tipo podría indicarnos que dichos genes se coexpresan y que participan en la misma etapa de desarrollo.
7 Comparación de métodos
¿Cómo se comparan las trayectorias inferidas por el PCA, TSCAN, destiny, SLICER y ouija? Para responder a esta pregunta, analizaremos su matriz de correlación. Dado que en los apartados previos hemos ido ajustando el pseudotiempo de los diversos métodos para que todos tuviesen el mismo sentido (i.e. el pseudotiempo progresa de izquierda a derecha, estando el zigoto en el comienzo y el blastocisto al final del mismo), no es necesario realizar más modificaciones.
# Obtenemos con `grep` (RegEx) los pseudotiempos inferidos por `TSCAN`,
# `destiny`, `SLICER` y `ouija` y los agrupamos en un dataframe
df_pseudotime <- as.data.frame(
colData(deng_SCE)[,grep("pseudotime", colnames(colData(deng_SCE)))]
)
# Añadimos también la trayectoria inferida con el PCA
df_pseudotime$PC1 <- deng_SCE$PC1
# Renombramos las columnas con RegEx
colnames(df_pseudotime) <- gsub("pseudotime_", "", colnames(df_pseudotime))
# Graficamos la matriz de correlación
corrplot.mixed(cor(df_pseudotime, use = "na.or.complete"),
order = "hclust", tl.col = "black",
main = "Matriz de correlación para los resultados del pseudotiempo",
mar = c(0,0,2,0))
Aquí vemos que ouija, TSCAN y SLICER dan trayectorias fuertemente correlacionadas con la 1ª CP (R > 0.9). La trayectoria computada por destiny está menos correlacionada con el resto de métodos (0.65 < R < 0.8).
Ejercicio 2.5: Compara los métodos en mayor profundidad. ¿Dónde y cómo difieren?
TSCANes un método basado en MST o Minimum Spanning Tree, y es conceptualmente muy parecido aMonocle. Además, hace uso de PCA para reducir la dimensionalidad de los datos y Mclust para agrupar las células en un nº óptimo de clusters, determinado automáticamente usando el criterio de información bayesiano (= BIC o criterio de Schwarz).TSCANcalcula un MST que pasa por el centroide de todos los clusters y de ahí infiere la trayectoria como el camino más largo a través del árbol (Monocle también busca la ruta más larga en un MST). Por último, las células son proyectadas al punto más cercano de la trayectoria, creando el orden final de las mismas.TSCANes además uno de los pocos métodos completamente no supervisado, ya que al igual queMonocle, determina automáticamente los clusters de partida y de fin de la trayectoria.TSCANes fácil de usar dado que determina todos sus parámetros automáticamente. Si bien la trayectoria inferida porTSCANes más fiable que la deMonocle(puesto queTSCANagrupa primero las células en clusters antes de calcular el camino más largo en el MST), dicha trayectoria sigue sin ser robusta, dado que ligeros cambios en la localización de los centroides puede provocar grandes cambios en la susodicha trayectoria inferida.SLICERes un método de inferencia de trayectoria basado en grafos k-NN y capaz de inferir trayectorias ramificadas a la vez que requiere de un escaso input por parte del usuario. Usa una envolutra alfa o alpha-hull (serie de rectas que forman en su conjunto curvas en el espacio euclídeo y están asociadas a un grupo de instancias/puntos) para estimar el parámetrokdel LLE. Dicho LLE se usa para computar el grafo k-NN y encontrar las células distales, candidatas a ser el punto de partida o fin de la trayectoria. El usuario debe especificar cuál de estas células será el punto de partida, y a partir de ahíSLICERcalcula la ruta más corta al resto de las células. Los puntos de ramificación de la trayectoria se detectan por vía de la entropía geodésica (zonas con una elevada entropía indican puntos de ramificación).SLICERcalcula automáticamente sus propios parámetros y muestra visualmente las células distales, para permitir al usuario seleccionar la célula de partida óptima. A día de hoy es la única técnica que ofrece una trayectoria ramificada con un sesgo mínimo en los resultados, y escala bien con miles de células y genes.Monoclees un método para detectar trayectorias ramificadas. Primero detecta genes diferencialmente expresados en las poblaciones celulares de interés y luego aplica un ICA sobre dichos genes para reducir la dimensionalidad de los datos (esto se hace así porque el ICA no escala bien con muchos genes). Después se calcula un MST sobre las células y el camino más largo a través del mismo (nótese que al calcular el MST sobre células individuales y no sobre clusters como haceTSCAN, este algoritmo es sensible al ruido). Cada célula se proyecta al punto más cercano de la trayectoria inferida. Al ejecutar el ICA sólo sobre genes diferencialmente expresados, se introduce un sesgo en los resultados.
Nótese que los métodos actuales de inferencia de trayectoria raramente son capaces de modelar caminos alternativos que convergan en un mismo destino celular. También asumen que siempre hay una trayectoria en los datos, aunque no siempre sea el caso, por lo que en esas situaciones se recomienda emplear métodos de clustering y/o visualización alternativos tales como SPADE, viSNE o FlowSOM.
Para finalizar, conviene tener en cuenta que las trayectorias generadas por estos métodos no dejan de ser hipótesis, las cuales deben ser validadas experimentalmente.
8 Expresión génica temporal: scater
La mayoría de los paquetes vistos en esta práctica cuentan con herramientas para visualizar la expresión génica a lo largo del pseudotiempo. No obstante, SLICER carece de estas herramientas, por lo que es conveniente suplir esa carencia con las utilidades del paquete scater. Seguir indivualmente el rastro de ciertos genes resulta útil para identificar aquellos que tengan un rol importante en el proceso de diferenciación celular estudiado. Ilustraremos el procedimiento usando los genes Rhoa, Hoxa1 y Carm1.
Dado que hemos añadido el pseudotiempo computado por cada algoritmo estudiado al bolsillo colData del objeto deng_SCE, podemos usar scater para investigar de manera visual la relación entre expresión génica, poblaciones celulares y pseudotiempo en nuestro dataset. Para ver en más detalle las utilidades que incorpora scater, véase su viñeta online.
Rhoa:
De acuerdo a la información disponible en GeneCards y MGI, el gen Rhoa:
Se localiza en el cromosoma 9 en ratón
Codifica para una GTPasa de la familia Rho, la cual actúa como interruptor molecular en cascadas de señalización
Promueve la reorganización del citoesqueleto, alterando así la forma, motilidad y adhesión celular. Además, participa en la citocinesis
Se expresa en el desarrollo embrionario
Visualizamos la expresión génica de Rhoa a lo largo de las trayectorias generadas:
PCA:
plotExpression(deng_SCE, features = "Rhoa", x = "PC1", colour_by = "cell_type2", show_smooth = T, show_violin = F)
TSCAN:
plotExpression(deng_SCE, "Rhoa", x = "pseudotime_order_tscan", colour_by = "cell_type2", show_smooth = T, show_violin = F)
Mapa de difusión:
plotExpression(deng_SCE, "Rhoa", x = "pseudotime_diffusionmap", colour_by = "cell_type2", show_smooth = T, show_violin = F)
SLICER:
plotExpression(deng_SCE, "Rhoa", x = "pseudotime_slicer", colour_by = "cell_type2", show_smooth = T, show_violin = F)
Ouija:
plotExpression(deng_SCE, "Rhoa", x = "pseudotime_ouija", colour_by = "cell_type2", show_smooth = T, show_violin = F)
Hoxa1:
De acuerdo a la información disponible en GeneCards y MGI, el gen Hoxa1:
En ratón, su locus se halla en el cromosoma 6
Codifica para un factor transcripción que se une a motivos ( motifs) de ADN
Regula procesos de desarrollo tales como la formación del tronco encefálico, oído interior y exterior, nervio VI ( abducens) y el sistema cardiovascular
El producto génico podría estar involucrado en la generación del tejido anterior-posterior (delante-arriba) del metencéfalo, el cual se diferenciará en cerebelo, puente de Valorio y los nervios craneales V (trigeminal o trigémino), VI ( abducens), VII (facial) y VIII (auditivo).
Codifica para dos transcritos de ARNm y produce 2 isoformas (aunque sólo una contiene la región del homeodominio)
Se expresa en las fases tempranas del desarrollo embrionario
Visualizamos la expresión génica de Hoxa1 a lo largo de las trayectorias generadas:
PCA:
plotExpression(deng_SCE, features = "Hoxa1", x = "PC1", colour_by = "cell_type2", show_smooth = T, show_violin = F)
TSCAN:
plotExpression(deng_SCE, "Hoxa1", x = "pseudotime_order_tscan", colour_by = "cell_type2", show_smooth = T, show_violin = F)
Mapa de difusión:
plotExpression(deng_SCE, "Hoxa1", x = "pseudotime_diffusionmap", colour_by = "cell_type2", show_smooth = T, show_violin = F)
SLICER:
plotExpression(deng_SCE, "Hoxa1", x = "pseudotime_slicer", colour_by = "cell_type2", show_smooth = T, show_violin = F)
Ouija:
plotExpression(deng_SCE, "Hoxa1", x = "pseudotime_ouija", colour_by = "cell_type2", show_smooth = T, show_violin = F)
Carm1:
De acuerdo a la información disponible en Genecards y MGI, el gen Carm1:
Se localiza en el cromosoma 9 en ratón
Codifica para una enzima que metila argininas de histonas y proteínas asociadas a la cromatina
Participa en procesos de regulación de la expresión génica a través del remodelado de cromatina y la regulación de la transcripción, del splicing alternativo de pre-ARNm y de la estabilidad del ARNm. También silencia los genes de la familia CREB
Se expresa en las fases tempranas del desarrollo embrionario
Visualizamos la expresión génica de Carm1 a lo largo de las trayectorias generadas:
PCA:
plotExpression(deng_SCE, features = "Carm1", x = "PC1", colour_by = "cell_type2", show_smooth = T, show_violin = F)
TSCAN:
plotExpression(deng_SCE, "Carm1", x = "pseudotime_order_tscan", colour_by = "cell_type2", show_smooth = T, show_violin = F)
Mapa de difusión:
plotExpression(deng_SCE, "Carm1", x = "pseudotime_diffusionmap", colour_by = "cell_type2", show_smooth = T, show_violin = F)
SLICER:
plotExpression(deng_SCE, "Carm1", x = "pseudotime_slicer", colour_by = "cell_type2", show_smooth = T, show_violin = F)
Ouija:
plotExpression(deng_SCE, "Carm1", x = "pseudotime_ouija", colour_by = "cell_type2", show_smooth = T, show_violin = F)
Ejercicio 2.6?: Para nuestro dataset, ¿qué algoritmos de los aquí estudiados riden mejor que el PCA a la hora de inferir el pseudotiempo?
Solución:
Nos encontramos ante una respuesta compleja, y es que al igual que sucede en Machine Learning, no hay un algoritmo que dé el mejor resultado en todos los casos, sino que cada uno de los algoritmos estudiados destaca en algo.
Para etapas tempranas del desarrollo (de “zy” a “4 cell”),
destinyes el método que mejor separa las etapas, y en su orden correcto.Para las etapas “8 cell” y “16 cell”,
SLICERes el algoritmo que mejor las separa, pero les asigna un orden incorrecto, pues segúnSLICER, la etapa de 16 células se sucede antes que la etapa de 8 células. Por otro lado,TSCANordena correctamente ambas etapas, y las separa relativamente bien (exceptuando unas pocas células de la etapa “8 células”)Para las 3 etapas del blastocisto, el PCA es el método que mejor las separa del resto de etapas, aunque si nos centramos en esas 3 etapas, le cuesta discernir cuándo una célula es del “early blast”, “mid blast” o “late blast”. En ese aspecto,
SLICERlogra distinguir con mejores resultados las 3 etapas de blastocisto entre sí, si bien no hay ningún algoritmo que produjese una separación limpia entre esas 3 etapas.Finalmente, comentar que si bien
SLICERno es perfecto, parece producir los clusters más compactos en general.
A continuación se recopilan los resultados obtenidos por los métodos aquí estudiados:
PCA:
TSCAN:
ouija:
SLICER:
destiny:
9 sessionInfo()
Click para mostrar
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19042)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=Spanish_Spain.1252 LC_CTYPE=Spanish_Spain.1252
## [3] LC_MONETARY=Spanish_Spain.1252 LC_NUMERIC=C
## [5] LC_TIME=Spanish_Spain.1252
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] rstan_2.21.2 StanHeaders_2.21.0-7
## [3] lle_1.1 snowfall_1.84-6.1
## [5] snow_0.4-3 MASS_7.3-53.1
## [7] scatterplot3d_0.3-41 corrplot_0.84
## [9] ggbeeswarm_0.6.0 ggthemes_4.2.4
## [11] scater_1.14.0 ggplot2_3.3.3
## [13] ouija_0.99.1 Rcpp_1.0.6
## [15] SLICER_0.2.0 destiny_3.0.1
## [17] TSCAN_1.24.0 SingleCellExperiment_1.8.0
## [19] SummarizedExperiment_1.16.1 DelayedArray_0.12.3
## [21] BiocParallel_1.20.1 matrixStats_0.58.0
## [23] Biobase_2.46.0 GenomicRanges_1.38.0
## [25] GenomeInfoDb_1.22.1 IRanges_2.20.2
## [27] S4Vectors_0.24.4 BiocGenerics_0.32.0
##
## loaded via a namespace (and not attached):
## [1] utf8_1.2.1 R.utils_2.10.1 tidyselect_1.1.0
## [4] grid_3.6.3 combinat_0.0-8 ranger_0.12.1
## [7] munsell_0.5.0 codetools_0.2-18 withr_2.4.1
## [10] colorspace_2.0-0 fastICA_1.2-2 highr_0.8
## [13] knitr_1.31 robustbase_0.93-7 vcd_1.4-8
## [16] tensor_1.5 VIM_6.1.0 TTR_0.24.2
## [19] labeling_0.4.2 splancs_2.01-40 GenomeInfoDbData_1.2.2
## [22] polyclip_1.10-0 farver_2.1.0 coda_0.19-4
## [25] vctrs_0.3.7 generics_0.1.0 xfun_0.22
## [28] R6_2.5.0 rsvd_1.0.3 RcppEigen_0.3.3.9.1
## [31] bitops_1.0-6 spatstat.utils_2.1-0 promises_1.2.0.1
## [34] scales_1.1.1 nnet_7.3-15 beeswarm_0.3.1
## [37] gtable_0.3.0 processx_3.5.0 goftest_1.2-2
## [40] rlang_0.4.10 splines_3.6.3 hexbin_1.28.2
## [43] inline_0.3.17 yaml_2.2.1 abind_1.4-5
## [46] httpuv_1.5.5 tensorA_0.36.2 tools_3.6.3
## [49] cubature_2.0.4.1 ellipsis_0.3.1 gplots_3.1.1
## [52] jquerylib_0.1.3 proxy_0.4-25 MCMCglmm_2.32
## [55] plyr_1.8.6 zlibbioc_1.32.0 purrr_0.3.4
## [58] RCurl_1.98-1.3 ps_1.6.0 prettyunits_1.1.1
## [61] rpart_4.1-15 alphahull_2.2 deldir_0.2-10
## [64] viridis_0.5.1 cowplot_1.1.1 zoo_1.8-9
## [67] haven_2.3.1 magrittr_2.0.1 data.table_1.14.0
## [70] RSpectra_0.16-0 openxlsx_4.2.3 lmtest_0.9-38
## [73] pcaMethods_1.78.0 hms_1.0.0 mime_0.10
## [76] evaluate_0.14 xtable_1.8-4 smoother_1.1
## [79] rio_0.5.26 mclust_5.4.7 readxl_1.3.1
## [82] gridExtra_2.3 compiler_3.6.3 tibble_3.1.0
## [85] KernSmooth_2.23-18 V8_3.4.0 crayon_1.4.1
## [88] R.oo_1.24.0 htmltools_0.5.1.1 mgcv_1.8-34
## [91] corpcor_1.6.9 later_1.1.0.1 tidyr_1.1.3
## [94] RcppParallel_5.0.3 boot_1.3-27 Matrix_1.3-2
## [97] car_3.0-10 cli_2.4.0 sgeostat_1.0-27
## [100] R.methodsS3_1.8.1 igraph_1.2.6 forcats_0.5.1
## [103] pkgconfig_2.0.3 foreign_0.8-75 laeken_0.5.1
## [106] sp_1.4-5 vipor_0.4.5 bslib_0.2.4
## [109] XVector_0.26.0 stringr_1.4.0 callr_3.6.0
## [112] digest_0.6.27 spatstat.data_2.1-0 rmarkdown_2.7
## [115] cellranger_1.1.0 DelayedMatrixStats_1.8.0 curl_4.3
## [118] shiny_1.6.0 gtools_3.8.2 ggplot.multistats_1.0.0
## [121] lifecycle_1.0.0 nlme_3.1-152 jsonlite_1.7.2
## [124] carData_3.0-4 BiocNeighbors_1.4.2 viridisLite_0.3.0
## [127] fansi_0.4.2 pillar_1.5.1 lattice_0.20-41
## [130] loo_2.4.1 fastmap_1.1.0 DEoptimR_1.0-8
## [133] pkgbuild_1.2.0 glue_1.4.2 xts_0.12.1
## [136] zip_2.1.1 spatstat_1.64-1 class_7.3-18
## [139] stringi_1.5.3 sass_0.3.1 RcppHNSW_0.3.0
## [142] BiocSingular_1.2.2 caTools_1.18.2 dplyr_1.0.5
## [145] knn.covertree_1.0 irlba_2.3.3 e1071_1.7-6
## [148] ape_5.4-1 tripack_1.3-9.1