Informe práctica 4 - Inferencia de trayectorias con Monocle 3 1.0.0
Adam Casas
Compilado: 07 de julio del 2021
1 Introducción
En esta práctica usaremos los mismos datos que los de la práctica 1, o sea el dataset de SC3-seq (un tipo de scRNA-seq) del desarrollo embrionario en Macaca fascicularis (Nakamura et al., 2016). En esta ocasión, en lugar de usar TSCAN u ouija, usaremos Monocle sobre estos datos para inferir la trayectoria del desarrollo embrionario.
Las células cambian de un estado funcional a otro en respuesta a estímulos, durante patologías o a lo largo del desarrollo vital del individuo. Cada estado funcional expresa un conjunto de genes concretos, lo que dota a la célula en cuestión un repertorio de proteínas y metabolitos específicos para llevar a cabo la tarea para la cual esté diseñado el estado en cuestión.
No obstante, cuando una célula transiciona de un estado a otro, esta sufre una reconfiguración transcripcional gradual, con genes que se silencian y genes que se activan. Purificar células en transformación puede resultar difícil, lo que dificulta la caracterización de dichos estados transitorios. Por suerte, este problema no lo tenemos con scRNA-seq, dado que no es necesario purificar poblaciones celulares.
2 Esquema del protocolo de Monocle 3
El protocolo de Monocle consta de 5 pasos principales, detallados a continuación:
Protocolo de Monocle. Sacado de Trapnell et al. (2014)
2.1 Paso 1: Preprocesamiento y normalizado de los datos
Primero Monocle normaliza los valores de expresión génica para controlar la varianza técnica debida a la profundidad del secuenciado y la recuperación del ARN de la muestra.
2.2 Paso 2: Reducir la dimensionalidad de los datos
Luego se proyectan las células al espacio del PCA (por defecto se calculan las 50 1as componentes) para eliminar ruido y facilitar la computación de los siguientes pasos. Tras esto, y de manera optativa, puedes reducir la dimensionalidad más aún con t-SNE o UMAP. Sea como sea, después de la reducción de la dimensionalidad Monocle puede aplicar un clustering a las células, inferir la trayectoria, o hacer ambas.
2.3 Paso 3: Clustering y particionado de las células
Monocle 3 puede aprender múltiples trayectorias desconectadas, lo cual es importante dado que muchos experimentos capturan una comunidad de células que estén respondiendo a un estímulo o diferenciándose, y cada célula puede responder de una manera u otra.
Monocle 2 asume que todos tus datos son parte de una única trayectoria, por lo que para construir varias trayectorias individuales, se tendría que dividir manualmente cada grupo de células, generar datasets de manera acorde y ejecutar Monocle 2 sobre cada dataset. Por otro lado, Monocle 3 puede detectar que algunas células forman parte de un proceso biológico distinto y generar acordemente varias trayectorias en paralelo, sin tener que manipular y dividir el dataset como sería el caso con Monocle 2. Monocle 3 consigue esto particionando las células en supergrupos usando un método derivado del AGA o Approximate Graph Abstraction (Wolf et al., 2019). Células de supergrupos distintos no pueden ser parte de la misma trayectoria.
2.4 Paso 4: Generar el grafo principal
Versiones anteriores de Monocle 3 y Monocle 2 disponían de 3 maneras distintas de organizarlas las células en trayectorias, todas ellas basadas en el concepto de embebimiento en grafo reverso o “reversed graph embedding”. DDRTree era el método usado por Monocle 2 para generar trayectorias con forma de árbol, y este algoritmo fue mejorado en veriones anteriores de Monocle 3. En concreto, se optimizó su rendimiento, lo que le permitía procesar millones de células en minutos. SimplePPT genera también trayectorias en forma de árbol, pero a diferencia de DDRTree, no aplica reducciones de la dimensionalidad adicionales. L1Graph es un algoritmo de optimización avanzado capaz de generar trayectorias que incluyan bucles.
No obstante, el autor ha expresado recientemente que Monocle 3 1.0.0 en adelante incorporará sólo SimplePPT, puesto que al combinarlo con UMAP, es más rápido, robusto y fácil de interpretar que DDRTree.
“DDRTree is not likely to be included in Monocle 3 going forward, because it doesn’t scale well to the kind of datasets Monocle 3 is designed to process. Beyond a 100,000 cells DDRTree is impractically slow. However, SimplePPT and DDRTree are closely related techniques, and in our experience, UMAP + SimplePPT is faster, more robust, and gives easier to interpret trajectories than UMAP or PCA + DDRTree.”
– Cole Trapnell, 2019
Una vez que Monocle ha aprendido un grafo principal que se ajusta a los datos, se proyecta cada célula al susodicho grafo. A continuación el usuario selecciona uno o más puntos en el grafo que definan los puntos de partida de la trayectoria. Monocle define entonces el pseudotiempo como la distancia de cada célula al punto de inicio más cercano.
2.5 Paso 5: Análisis de expresión diferencial y visualización
Por último, detectamos biomarcadores específicos de cada clúster celular, genes con expresión cambiante a lo largo de la trayectoria y graficamos los resultados.
Monocle 3 provee una suite de tests estadísticos para detectar genes que difieren entre clusters y condiciones experimentales (covariables tales como tratamiento, tiempo, género…). Además, Monocle 3 incluye un novedoso test que usa el grafo principal directamente y permite detectar genes cuyas expresiones varían de manera compleja a lo largo de trayectorias con bucles y estructuras más complicadas.
3 El protocolo de Monocle 3 en la práctica
Antes de nada debemos asegurarnos de tener instalados tanto Seurat 3 como Monocle3 1.0.0. A continuación se proporciona un chunk de código para instalar dichos paquetes. Este chunk no se ejecutará en la compilación de este informe, dado que sendos paquetes ya se encuentran instalados en mi máquina.
# Instalación de Monocle 3 V1.0.0
BiocManager::install(c('BiocGenerics', 'DelayedArray', 'DelayedMatrixStats',
'limma', 'S4Vectors', 'SingleCellExperiment',
'SummarizedExperiment', 'batchelor', 'Matrix.utils'))
install.packages("devtools")
devtools::install_github('cole-trapnell-lab/leidenbase')
devtools::install_github('cole-trapnell-lab/monocle3')
# Si usas R 3.6, instala spatstat 1.64-1, ya que R 3.6 no es compatible con
# versiones posteriores de spatstat
install.packages("https://cran.r-project.org/src/contrib/Archive/spatstat/spatstat_1.64-1.tar.gz",
repo=NULL, type="source")## [1] "Seurat version 3.2.0 is installed"## [1] "Monocle 3 version 1.0.0 is installed as well"Hecho eso, procedemos a cargar las librerías de interés y cargar los datos:
# Cargamos librerías
library(monocle3)
library(reticulate)
library(Seurat)
library(ggplot2)
library(dplyr)
# Cargamos en memoria el objeto seurat de la práctica 1
data <- readRDS(file = "./Cosas accesorias Informe 1/nakamura_P4.rds")
head(data@meta.data)# Creamos un nuevo objeto CellDataSet a partir del objeto seurat. Nótese que
# tenemos que añadir la columna "gene_short_name" dado que varios métodos de
# Monocle 3 lo usan, como por ejemplo plot_cells(c3, genes = "ZNF692")
gene.annotation <- data.frame(gene_short_name = rownames(data@assays$RNA@meta.features),
row.names = rownames(data@assays$RNA@meta.features))
head(gene.annotation)# data[["RNA"]]@meta.features$gene_short_name <- rownames(data@assays$RNA@meta.features)
c3 <- new_cell_data_set(expression_data = data@assays$RNA@counts,
cell_metadata = data@meta.data,
gene_metadata = gene.annotation) # gene_metadata = data@assays$RNA@meta.featuresPara los objetos de tipo CellDataSet, podemos acceder a su información con los métodos pData(), fData() y exprs() de Monocle 3 (nótese que al crear un nuevo objeto CellDataSet, Monocle 3 ha calculado automáticamente el Size_Factor de las células, pues esa columna no la generamos en la práctica 1 y por tanto no ha aparecido al examinar el bolsillo @meta.data de nuestro objeto seurat):
# Vemos los metadatos de las células
head(pData(c3))## DataFrame with 6 rows and 10 columns
## orig.ident nCount_RNA nFeature_RNA RNA_snn_res.0.6 seurat_clusters
## <factor> <numeric> <integer> <factor> <factor>
## EXMC Nakamura 1550084 12492 6 6
## EXMC.1 Nakamura 1424376 11652 6 6
## EXMC.2 Nakamura 1539115 11988 6 6
## EXMC.3 Nakamura 1552198 11693 6 6
## EXMC.4 Nakamura 1907241 12143 6 6
## EXMC.5 Nakamura 2709088 12389 6 6
## Cell.Identity S.Score G2M.Score Phase
## <factor> <numeric> <numeric> <factor>
## EXMC EXMC -0.277548682038991 0.20165794549314 G2M
## EXMC.1 EXMC 0.102056042052932 -0.167133047794692 S
## EXMC.2 EXMC -0.206478768990954 0.0470766073806089 G2M
## EXMC.3 EXMC -0.150290022846147 0.267351913826893 G2M
## EXMC.4 EXMC 0.0821853205230935 -0.0541958056812388 S
## EXMC.5 EXMC 0.0861077841728103 -0.216778102457795 S
## Size_Factor
## <numeric>
## EXMC 0.886814042911187
## EXMC.1 0.814895605132151
## EXMC.2 0.880538600266341
## EXMC.3 0.888023477294558
## EXMC.4 1.09114609402844
## EXMC.5 1.5498884459695# Accedemos a la columna Cell.Identity para ver a qué tejido pertenecen las
# células
head(pData(c3)$Cell.Identity)## [1] EXMC EXMC EXMC EXMC EXMC EXMC
## 11 Levels: EXMC Gast Hypoblast ICM Post-paTE PostE-EPI PostL-EPI ... VEYE# Alternativamente:
# c3$Cell.Identity
# Ídem para los metadatos de los genes con fData:
colnames(fData(c3))## [1] "gene_short_name"head(fData(c3))## DataFrame with 6 rows and 1 column
## gene_short_name
## <factor>
## LOC102116043 LOC102116043
## PGBD2 PGBD2
## ZNF692 ZNF692
## ZNF672 ZNF672
## SH3BP5L SH3BP5L
## LYPD8 LYPD8# Visualizamos una parte de la matriz de conteos
exprs(c3)[1:5,1:5]## 5 x 5 sparse Matrix of class "dgCMatrix"
## EXMC EXMC.1 EXMC.2 EXMC.3 EXMC.4
## LOC102116043 . 3 . 2 6
## PGBD2 3 17 . 1 7
## ZNF692 14 51 8 8 7
## ZNF672 57 63 42 33 84
## SH3BP5L . 30 9 14 163.1 Paso 1: Preprocesamiento y normalizado de los datos
En el preprocesamiento y normalizado de datos se suele proceder a estimar el factor de tamaño de las muestras y la dispersión (=varianza) de la expresión de los genes. Monocle 3 realiza estas y otras operaciones con el paquete DelayedArray, lo que le permite escalar con millones de células. Dicho paquete divide las operaciones en bloques para prevenir el consumo excesivo de RAM del ordenador. Puedes controlar el tamaño de cada bloque y los mensajes que se devuelven por consola con los siguientes comandos:
Nota: Puesto que la estimación del factor de tamaño de las muestras se realiza automáticamente al generar un nuevo objeto
CellDataSet, podríamos incluir el siguiente chunk de código antes de la carga de datos, pero en aras de mantener un orden, lo incluiremos aquí. Si tienes experiencia trabajando con Monocle 3, puedes correr este chunk antes de generar el susodicho objeto.
# Selecciona TRUE para ver en consola el progreso de algunas operaciones de
# Monocle 3
DelayedArray:::set_verbose_block_processing(TRUE)## [1] FALSE# Tamaño de cada bloque, en bytes (es 1e8 bytes = 95.37MBs por defecto). A mayor
# tamaño del bloque, más rápidas serán algunos cálculos pero usará más memoria
# RAM, así que ajústalo con cautela.
DelayedArray::getAutoBlockSize()## [1] 1e+08DelayedArray::setAutoBlockSize(2e8)
# Alternativamente:
# getOption("DelayedArray.block.size")
# options(DelayedArray.block.size=1e9)Establecido el tamaño de bloque y los mensajes por consola, procedemos al normalizado de los datos con el comando preprocess_cds. Nótese que se recomienda usar el método PCA para datos de RNA-seq y LSI para datos de ATAC-seq.
# Si no se ha calculado previamente el size factor de cada muestra/célula, se
# puede calcular aquí con el comando:
# c3 <- estimate_size_factors(c3)
c3 <- preprocess_cds(c3, method = "PCA", num_dim = 20,
norm_method = "log")
plot_pc_variance_explained(c3)
3.2 Paso 2: Reducir la dimensionalidad de los datos
Ahora procedemos a computar la proyección de los datos en un espacio de dimensionalidad reducida (la proyección por defecto en Monocle 3 1.0.0 es UMAP, pero también podemos calcular el t-SNE, PCA, LSI o Aligned) con el comando reduce_dimension(). Tras calcular la proyección, podemos visualizarla con el comando plot_cells().
Nótese que si calculamos el UMAP, podemos acelerar su computación modificando los parámetros umap.fast_sgd y cores, pero hacer esto incurre en que cada vez que calculemos el UMAP, generaremos proyecciones ligeramente distintas, aunque sólo hayamos modificado uno de los 2 parámetros mentados. Esto puede ser aceptable según la situación, pero por motivos de reproducibilidad de este informe, optaremos por no modificar ninguno de esos 2 parámetros Por fortuna, el dataset de Nakamura et al. es relativamente pequeño, por lo que no notaremos mucho la diferencia.
c3 <- reduce_dimension(c3, preprocess_method = "PCA",
reduction_method = "UMAP",
umap.fast_sgd = F,
cores = 1)
plot_cells(c3, color_cells_by = "Cell.Identity",
group_label_size = 3, norm_method = "log") +
ggtitle("Tejidos embrionarios (UMAP)") +
theme(plot.title = element_text(hjust = 0.5))
# Biomarcador de early epiblast, expresado en tejidos Pre-EPI y PostE-EPI
plot_cells(c3, genes = "GBX2", color_cells_by = "Cell.Identity", group_label_size = 3,
norm_method = "log", cell_size = 1) +
ggtitle("Expresión de GBX2, biomarcador del epiblasto temprano (UMAP)") +
theme(plot.title = element_text(hjust = .5))
# Biomarcador del hipoblasto
plot_cells(c3, genes = "GATA6", color_cells_by = "Cell.Identity", group_label_size = 3,
norm_method = "log", cell_size = 1) +
ggtitle("Expresión de GATA6, biomarcador del hipoblasto (UMAP)") +
theme(plot.title = element_text(hjust = .5))
3.3 Paso 3: Partición de las células en supergrupos/particiones
En vez de forzar a todas las células en una única trayectoria, Monocle 3 permite inferir un conjunto de trayectorias, de manera que te permite estudiar varios procesos biológicos simultáneos con un mismo dataset. Por ejemplo, durante una infección (ya sea vírica, bacteriana, fúngica…), los diversos componentes del sistema inmune responden a un cóctel de estímulos (antígenos, citocinas, superantígenos…), por lo que cada célula inmune ejercerá un rol determinado. Por ende, dicho proceso se caracteriza mejor por un conjunto de trayectorias, en lugar de una única trayectoria.
Otra ventaja de inferir más de una trayectoria es que el algoritmo se vuelve resistente a la presencia de células extremas o outliers, las cuales normalmente confundirían a la trayectoria original. Al particionar las células en supergrupos o particiones, Monocle 3 agrupa entre sí dichas células extremas y las diferencia de otros supergrupos.
Monocle 3 incluye un test para la conectividad de comunidades celulares el cual divide las células en los ya mencionados supergrupos. Este test se deriva del concepto de participación en grafos abstractos (abstract graph participation), desarrollado por Wolf et al. (2019) para su algoritmo PAGA. Para realizar dicho clustering, basta con llamar a la función cluster_cells().
Maticemos que la función cluster_cells() genera clusters mediante el método de Leiden (Traag et al., 2019), mientras que los supergrupos, particiones o superclusters son generados mediante el ya mencionado algoritmo de Wolf et al..
c3 <- cluster_cells(c3, reduction_method = "UMAP",
cluster_method = "leiden",
verbose = T)## resolution_parameter quality modularity significance cluster_count selected
## 1e-04 15954.84 0.673237 7712.089 5 *
## cluster cell_count cell_fraction
## 1 206 0.516
## 2 54 0.135
## 3 53 0.133
## 4 52 0.130
## 5 34 0.085plot_cells(c3, color_cells_by = "cluster", group_label_size = 4) +
ggtitle("Clusters generados") +
theme(plot.title = element_text(hjust = .5))
plot_cells(c3, color_cells_by = "partition", group_label_size = 4) +
ggtitle("Particiones generadas") +
theme(plot.title = element_text(hjust = .5))
Una vez calculados los clusters y las particiones, podemos acceder a ellas mediante los siguientes comandos de Monocle 3:
head(clusters(c3))## EXMC EXMC.1 EXMC.2 EXMC.3 EXMC.4 EXMC.5
## 1 1 1 1 1 1
## Levels: 1 2 3 4 5head(partitions(c3))## EXMC EXMC.1 EXMC.2 EXMC.3 EXMC.4 EXMC.5
## 1 1 1 1 1 1
## Levels: 13.4 Paso 4: Generar el grafo principal
Como ya se comentó en la sección 2.4, emplearemos UMAP + SimplePPT para generar el grafo a partir del cual inferiremos la trayectoria(s). Cada supergrupo genera su propia trayectoria, y puesto que nuestros datos forman parte de un único supergrupo, esta vez generaremos una única trayectoria. Esto se consigue con el comando learn_graph():
c3 <- learn_graph(c3, use_partition = T,
close_loop = T, verbose = T)## Running louvain iteration 1 ...
## -Number of clusters: 11Ahora cada vez que llamemos a plot_cells(), veremos el grafo aprendido, y como ya comentamos en el párrafo anterior, en nuestro caso constará de una única trayectoria:
cell_type_color <- c("EXMC" = "blue",
"Gast" = "#46C7EF",
"PostL-EPI" = "magenta",
"Hypoblast" = "red",
"ICM" = "darkgreen",
"PreE-TE" = "#4EB859",
"Post-paTE" = "#EFAD1E",
"PostE-EPI" = "gold",
"Pre-EPI" = "gray",
"PreL-TE" = "black",
"VEYE" = "purple")
plot_cells(c3, color_cells_by = "Cell.Identity", cell_size = 2,
label_cell_groups = F, graph_label_size = 3,
label_branch_points = T, label_leaves = T) +
scale_color_manual(values = cell_type_color)
En el grafo superior, los números encapsulados en círculos negros son los puntos donde se ramifica el grafo/árbol/trayectoria (aquí son términos intercambiables), mientras que los círculos grises corresponden a las hojas del árbol, i.e. los posibles destinos que pueden tomar nuestras células al final de la trayectoria. A continuación desglosamos el grafo anterior para hacerlo más comprensible. Se observa que el árbol presenta 4 puntos de ramificación, y 6 destinos celulares (=hojas).
Se observa también que Monocle 3 no asigna automáticamente el nodo raíz del árbol. Para ello y con ayuda del comando order_cells() elegiremos en el gráfico el nodo raíz, a partir del cual Monocle 3 calculará el pseudotiempo y ordenará en él las células. Si no se proporcionan los nodos/células a esta función, se abrirá una ventana interactiva mediante la cual podremos elegir visualmente los nodos raíz. Este modo interactivo es incompatible con el compilado de informes, por lo que cuando ejecutes el código en tu máquina, podrás hacerlo en modo interactivo, pero en mi caso le pasaré a la función los nodos raíz en forma de argumento. Para ello, Visualizaremos el grafo, pero esta vez mostrando los nombres de los puntos principales.
Téngase en cuenta que a diferencia de Monocle 2, Monocle 3 1.0.0 permite elegir nodos internos como nodos raíz (en Monocle 2, sólo los nodos distales podían ser elegidos como nodos raíz, de manera similar a como pasaba con los mapas de difusión en la práctica 2). Además, como Monocle 3 puede generar más de una trayectoria, resulta coherente que se puedan elegir más de un nodo raíz.
plot_cells(c3, label_principal_points = T, cell_size = 2)
A priori sabemos que la masa interna de células (ICM o inner cell mass) del embrión dará lugar al hipoblasto y al epiblasto. Por tanto para ajustar el comienzo del pseudotiempo, tenemos que especificarle a Monocle 3 que dichas células son los nodos raíz del grafo. En este caso, el punto principal Y_67 se corresponde con la ICM, por lo que se lo pasaremos al comando order_cells() a través del argumento root_pr_nodes:
c3 <- order_cells(c3, reduction_method = "UMAP", root_pr_nodes = "Y_67")
plot_cells(c3, color_cells_by = "pseudotime", cell_size = 2,
label_leaves = F, label_branch_points = F,
label_roots = T, graph_label_size = 3)
En la página web de Monocle 3 se facilita una función auxiliar para facilitar la extracción de nodos raíz en el grafo, cortesía del autor de Monocle 3, Cole Trapnell. Dicha función se ha modificado ligeramente para adaptarla a nuestro caso, en el cual elegimos como nodo raíz aquellas células etiquetadas como ICM:
# Función auxiliar para identificar nodos raíz:
get_earliest_principal_node <- function(cds, cell_metadata, cell_identity){
# El argumento cell_metadata recibe el nombre de la columna en pData(c3)
# El argumento cell_identity recibe el valor a buscar en dicha columna
cell_ids <- which(colData(cds)[, cell_metadata] == cell_identity)
closest_vertex <-
cds@principal_graph_aux[["UMAP"]]$pr_graph_cell_proj_closest_vertex
closest_vertex <- as.matrix(closest_vertex[colnames(cds), ])
root_pr_nodes <-
igraph::V(principal_graph(cds)[["UMAP"]])$name[as.numeric(names
(which.max(table(closest_vertex[cell_ids,]))))]
root_pr_nodes
}
pData(c3)## DataFrame with 399 rows and 10 columns
## orig.ident nCount_RNA nFeature_RNA RNA_snn_res.0.6 seurat_clusters
## <factor> <numeric> <integer> <factor> <factor>
## EXMC Nakamura 1550084 12492 6 6
## EXMC.1 Nakamura 1424376 11652 6 6
## EXMC.2 Nakamura 1539115 11988 6 6
## EXMC.3 Nakamura 1552198 11693 6 6
## EXMC.4 Nakamura 1907241 12143 6 6
## ... ... ... ... ... ...
## Gast.45 Nakamura 1128788 10876 4 4
## Gast.46 Nakamura 1414743 12058 4 4
## Gast.47 Nakamura 1017712 11047 4 4
## Gast.48 Nakamura 2124849 16919 4 4
## Gast.49 Nakamura 2639447 12642 4 4
## Cell.Identity S.Score G2M.Score Phase
## <factor> <numeric> <numeric> <factor>
## EXMC EXMC -0.277548682038991 0.20165794549314 G2M
## EXMC.1 EXMC 0.102056042052932 -0.167133047794692 S
## EXMC.2 EXMC -0.206478768990954 0.0470766073806089 G2M
## EXMC.3 EXMC -0.150290022846147 0.267351913826893 G2M
## EXMC.4 EXMC 0.0821853205230935 -0.0541958056812388 S
## ... ... ... ... ...
## Gast.45 Gast -0.0265269285621769 -0.259579599141627 G1
## Gast.46 Gast -0.124061713175624 0.522947602788739 G2M
## Gast.47 Gast 0.333157471667443 -0.203806463700028 S
## Gast.48 Gast -0.0719841488390909 0.022089556620582 G2M
## Gast.49 Gast 0.159078376606547 -0.195129560788938 S
## Size_Factor
## <numeric>
## EXMC 0.886814042911187
## EXMC.1 0.814895605132151
## EXMC.2 0.880538600266341
## EXMC.3 0.888023477294558
## EXMC.4 1.09114609402844
## ... ...
## Gast.45 0.645787615296741
## Gast.46 0.809384497556456
## Gast.47 0.58224024842475
## Gast.48 1.21564117316597
## Gast.49 1.51004633627584# Ordenamos las células con el uso de la función auxiliar
c3 <- order_cells(c3,
root_pr_nodes = get_earliest_principal_node(c3, cell_metadata = "Cell.Identity", cell_identity = "ICM"))
# Graficamos
plot_cells(c3, color_cells_by = "pseudotime", cell_size = 2,
label_leaves = F, label_branch_points = F,
label_roots = T, graph_label_size = 3)
Por último, comentar que el pseudotiempo se almacena en el bolsillo c3@principal_graph_aux$UMAP$pseudotime:
head(c3@principal_graph_aux$UMAP$pseudotime)## EXMC EXMC.1 EXMC.2 EXMC.3 EXMC.4 EXMC.5
## 12.33988 12.25955 12.33161 12.29574 12.41787 12.141403.4.1 Visualización de la trayectoria en 3D
Monocle 3 también permite la visualización de las trayectorias en gráficos 3D interactivos. Para ello, durante la reducción de la dimensionalidad tenemos que pedirle a la función reduce_dimension() que guarde como máximo 3 dimensiones mediante el parámetro max_components:
# Protocolo resumido de Monocle 3 para generar gráficos en 3D
cds_3d <- reduce_dimension(c3, max_components = 3,
reduction_method = "UMAP", preprocess_method = "PCA")
cds_3d <- cluster_cells(cds_3d, reduction_method = "UMAP")
cds_3d <- learn_graph(cds_3d)
cds_3d <- order_cells(cds_3d, root_pr_nodes=get_earliest_principal_node(cds_3d, cell_metadata = "Cell.Identity", cell_identity = "ICM"))
# Generamos el gráfico 3D y lo visualizamos:
cds_3d_plot_obj <- plot_cells_3d(cds_3d, color_cells_by="pseudotime")
cds_3d_plot_obj
# También podemos usar el comando plot_cells_3d():
plot_cells_3d(cds_3d, color_cells_by = "Cell.Identity",
color_palette = cell_type_color)
# NOTA: Ejecuta este chunk en local, pues es muy pesado como para generarlo en
# un informe
3.5 Paso 5: Análisis de expresión diferencial
Monocle 3 permite afrontar este problema mediante dos funciones distintas:
Usando modelos de regresión mediante la función
fit_models()para evaluar si la expresión de un gen varía en función de covariables tales como tratamiento, género, edad, etc…Mediante el análisis de autocorrelación espacial del grafo usando el comando
graph_test()(función paralelizable en UNIX pero no en Windows) para encontrar genes que varían entre clusters o a lo largo de la trayectoria.
En esta práctica cubriremos el uso del comando graph_test(), pues permite detectar genes diferencialmente expresados tanto en clusters como a lo largo de la trayectoria (dependiendo de si le pasamos al argumento neighbor_graph el valor knn o principal_graph, respectivamente). Esta función se basa en un novedoso método denominado test I de Moran. La I de Moran es un estadístico que mide autocorrelaciones espaciales multidimensionales y multidireccionales. Esta medida indica si señales espaciales colindantes (en este caso, las células del grafo y sus vecinas) tienen valores similares (=niveles de expresión similares) u opuestos.
Si bien la correlación de Pearson y la I de Moran devuelven medidas entre -1 y 1, la interpretación de la I de Moran es ligeramente distinta: +1 significa que células adyacentes tienen una expresión génica idéntica (imagen inferior, panel derecho); 0 indica que no hay correlación alguna (panel central) y -1 significa que las células vecinas están inversamente correlacionadas (panel izquierdo). Cao et al. demostraron exitosamente en 2017 la aplicabilidad del test de la I de Moran a datasets de scRNA-seq para el análisis de expresión diferencial.
En esta práctica nos centraremos en el análisis de autocorrelación espacial, puesto que según busquemos variaciones entre clusters o a lo largo de la trayectoria, tendremos que cambiar ligeramente el código empleado. Nosotros mostraremos cómo buscar genes que varían a lo largo de la trayectoria:
# Expresión génica diferencial a lo largo de la trayectoria (test I de Moran).
# Este comando es paralelizable en UNIX, pero no en Windows
moran_test_principal_graph <- graph_test(c3, neighbor_graph = "principal_graph", cores = 2, verbose = F)
# Visualizamos el dataframe
head(moran_test_principal_graph)El comando graph_test nos devuelve un dataframe cuyas filas son los genes del dataset y las columnas son:
El p-valor y el q-valor (equivalente a un p-valor a posteriori bayesiano ajustado mediante el metodo pFDR para evitar falsos positivos durante tests múltiples) del gen
El estadístico de la I de Moran (cualquier valor) y la I de Moran (valores entre [-1,1]). Recordemos que valores del test I de Moran cercanos a 1 significa que la expresión del gen en cuestión se focaliza en una región concreta de la trayectoria (i.e. es específico de un clúster o estado celular).
Los metadatos de cada gen (recordemos que al generar el objeto
CellDataSet, uno de los posibles inputs eran dichos metadatos, en nuestro caso incluíamos sólo los nombres de los genes en la columnagene_short_name)El estado (
status).
Procedemos a quedarnos con los genes cuyo q-valor sea menor a 0.05 y los ordenamos en función de su significancia estadística:
# Nos quedamos con los descubrimientos con un q-valor menor a 0.05
moran_test_principal_graph <- subset(moran_test_principal_graph, q_value < 0.05)
# Ordenamos los genes según su q-valor en orden ascendente (de más significativo
# a menos)
moran_test_principal_graph <- arrange(moran_test_principal_graph, q_value)
# Visualizamos los 6 genes diferencialmente expresados a lo largo de la
# trayectoria más significativos
head(moran_test_principal_graph, n = 6L)# Nos quedamos con los nombres de dichos genes
dif_exp_genes <- rownames(moran_test_principal_graph)Una vez identificados los genes diferencialmente expresados a lo largo de la trayectoria, podemos visualizarlos sobre el propio grafo. A continuación mostramos la expresión génica de los 2 genes más significativos: LOC102145562 y HMGA2.
plot_cells(c3, genes = dif_exp_genes[c(1,2)], alpha = .7,
color_cells_by = "Cell.Identity", label_leaves = F,
label_branch_points = F, label_cell_groups = T,
cell_size = 2, label_roots = T, group_label_size = 3)
Además, Monocle 3 permite analizar módulos de genes co-expresados, de manera similar a como lo hace el paquete WGCNA:
# Agrupamos los genes diferencialmente expresados (y estadísticamente
# significativos) en módulos de genes co-expresados
gene_module_df <- find_gene_modules(c3[dif_exp_genes,])
# Dataframe con nombres de las células y el tipo de tejido al que pertenecen
cell_group_df <- tibble::tibble(cell = rownames(pData(c3)),
cell_group = pData(c3)$Cell.Identity)
# Agregamos expresión génica de los módulos en una matriz
agg_mat <- aggregate_gene_expression(c3, gene_group_df = gene_module_df,
cell_group_df = cell_group_df)
# Renombramos
rownames(agg_mat) <- stringr::str_c("Module ", row.names(agg_mat))
# Visualizamos el heatmap
library(pheatmap)
pheatmap(agg_mat, scale = "column", clustering_method = "ward.D2")
# Aislamos el módulo 1 y pasamos los nombres de los genes que lo componen a un
# archivo de texto para poder buscarlo en gProfiler
modulo_1 <- gene_module_df %>% filter(module == 1)
nombres_genes_modulo1 <- modulo_1$id
write(nombres_genes_modulo1, file = "./Cosas accesorias Informe 4/modulo1.txt")Podemos visualizar la información obtenida en el heatmap directamente en la trayectoria. Dicho de otra manera, podemos ver en qué destinos celulares se sobreexpresan módulos de genes:
plot_cells(c3,
genes = gene_module_df %>% filter(module %in% c(1)),
color_cells_by = "Cell.Identity", group_label_size = 3,
show_trajectory_graph = T, label_leaves = F,
label_branch_points = F, cell_size = 2) +
theme(plot.title = element_text(hjust = .5)) + ggtitle("Expresión de genes del módulo 1")
La anotación génica aportada por gProfiler para el módulo 1 implica transporte celular de proteínas (aunque no nos centraremos en este tema en este informe, pues eso ya se practicó en otra asignatura).
Al igual que hicimos en la práctica 2 con scatter y ouija, Monocle 3 1.0.0 te permite visualizar cómo cambia la expresión génica a lo largo del pseudotiempo y observar el orden en el que se activan/desactivan genes. Para realizar esto, tenemos que pasarle a la función plot_genes_in_pseudotime un objeto CellDataSet que incluya sólo los genes que queremos visualizar. En este caso, vamos a crear otro objeto CellDataSet a partir del original c3 que contenga los top 3 genes diferencialmente expresados del módulo 1:
subconjunto.modulo1 <- subset(c3, gene_short_name %in% c("MPEG1", "HMGA2", "LOC102145562"))
plot_genes_in_pseudotime(subconjunto.modulo1, color_cells_by = "pseudotime", min_expr = .5, horizontal_jitter = T)
En el gráfico se observa que MPEG1 y LOC102145562 muestran una tendencia a la baja (MPEG1 llega a silenciarse) tras la activación del gen HMGA2, lo cual podría darnos una pista sobre los patrones de regulación transcripcional en juego.
4 Bibliografía
Nakamura, T., Okamoto, I., Sasaki, K., Yabuta, Y., Iwatani, C., Tsuchiya, H., Seita, Y., Nakamura, S., Yamamoto, T., and Saitou, M. (2016). A developmental coordinate of pluripotency among mice, monkeys and humans. Nature, 537(7618):57–62.
Wolf, F.A., Hamey, F.K., Plass, M. et al. PAGA: graph abstraction reconciles clustering with trajectory inference through a topology preserving map of single cells. Genome Biol 20, 59 (2019). https://doi.org/10.1186/s13059-019-1663-x
Traag, V.A., Waltman, L. & van Eck, N.J. From Louvain to Leiden: guaranteeing well-connected communities. Sci Rep 9, 5233 (2019). https://doi.org/10.1038/s41598-019-41695-z
Cao J, Packer JS, Ramani V, Cusanovich DA, Huynh C, Daza R, Qiu X, Lee C, Furlan SN, Steemers FJ, Adey A, Waterston RH, Trapnell C, Shendure J. Comprehensive single-cell transcriptional profiling of a multicellular organism. Science. 2017 Aug 18;357(6352):661-667. doi: 10.1126/science.aam8940. PMID: 28818938; PMCID: PMC5894354.
5 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 19043)
##
## 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] stats4 parallel stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] pheatmap_1.0.12 dplyr_1.0.7
## [3] ggplot2_3.3.4 Seurat_3.2.0
## [5] reticulate_1.20 monocle3_1.0.0
## [7] SingleCellExperiment_1.8.0 SummarizedExperiment_1.16.1
## [9] DelayedArray_0.12.3 BiocParallel_1.20.1
## [11] matrixStats_0.59.0 GenomicRanges_1.38.0
## [13] GenomeInfoDb_1.22.1 IRanges_2.20.2
## [15] S4Vectors_0.24.4 Biobase_2.46.0
## [17] BiocGenerics_0.32.0
##
## loaded via a namespace (and not attached):
## [1] plyr_1.8.6 igraph_1.2.6 lazyeval_0.2.2
## [4] sp_1.4-5 splines_3.6.3 listenv_0.8.0
## [7] digest_0.6.27 htmltools_0.5.1.1 viridis_0.6.1
## [10] gdata_2.18.0 fansi_0.5.0 magrittr_2.0.1
## [13] tensor_1.5 cluster_2.1.2 ROCR_1.0-11
## [16] globals_0.14.0 gmodels_2.18.1 colorspace_2.0-1
## [19] ggrepel_0.9.1 xfun_0.24 crayon_1.4.1
## [22] RCurl_1.98-1.3 jsonlite_1.7.2 spatstat_1.64-1
## [25] spatstat.data_2.1-0 survival_3.2-11 zoo_1.8-9
## [28] ape_5.5 glue_1.4.2 polyclip_1.10-0
## [31] gtable_0.3.0 zlibbioc_1.32.0 XVector_0.26.0
## [34] leiden_0.3.8 future.apply_1.7.0 abind_1.4-5
## [37] scales_1.1.1 DBI_1.1.1 miniUI_0.1.1.1
## [40] Rcpp_1.0.6 viridisLite_0.4.0 xtable_1.8-4
## [43] spData_0.3.10 units_0.7-2 spdep_1.1-8
## [46] rsvd_1.0.3 proxy_0.4-26 htmlwidgets_1.5.3
## [49] httr_1.4.2 speedglm_0.3-3 RColorBrewer_1.1-2
## [52] ellipsis_0.3.2 ica_1.0-2 pkgconfig_2.0.3
## [55] farver_2.1.0 sass_0.4.0 uwot_0.1.10
## [58] deldir_0.2-10 utf8_1.2.1 tidyselect_1.1.1
## [61] labeling_0.4.2 rlang_0.4.11 reshape2_1.4.4
## [64] later_1.2.0 pbmcapply_1.5.0 munsell_0.5.0
## [67] tools_3.6.3 generics_0.1.0 ggridges_0.5.3
## [70] evaluate_0.14 stringr_1.4.0 fastmap_1.1.0
## [73] yaml_2.2.1 goftest_1.2-2 knitr_1.33
## [76] fitdistrplus_1.1-5 purrr_0.3.4 RANN_2.6.1
## [79] pbapply_1.4-3 future_1.21.0 nlme_3.1-152
## [82] mime_0.10 slam_0.1-48 grr_0.9.5
## [85] compiler_3.6.3 plotly_4.9.4.1 png_0.1-7
## [88] e1071_1.7-7 spatstat.utils_2.2-0 tibble_3.1.2
## [91] bslib_0.2.5.1 stringi_1.6.2 highr_0.9
## [94] RSpectra_0.16-0 lattice_0.20-44 Matrix_1.3-4
## [97] classInt_0.4-3 vctrs_0.3.8 LearnBayes_2.15.1
## [100] pillar_1.6.1 lifecycle_1.0.0 lmtest_0.9-38
## [103] jquerylib_0.1.4 RcppAnnoy_0.0.18 data.table_1.14.0
## [106] cowplot_1.1.1 bitops_1.0-7 irlba_2.3.3
## [109] Matrix.utils_0.9.8 raster_3.4-13 httpuv_1.6.1
## [112] patchwork_1.1.1 R6_2.5.0 promises_1.2.0.1
## [115] KernSmooth_2.23-20 gridExtra_2.3 parallelly_1.26.0
## [118] codetools_0.2-18 gtools_3.9.2 boot_1.3-28
## [121] MASS_7.3-54 assertthat_0.2.1 leidenbase_0.1.3
## [124] withr_2.4.2 sctransform_0.3.2 GenomeInfoDbData_1.2.2
## [127] expm_0.999-6 mgcv_1.8-36 grid_3.6.3
## [130] rpart_4.1-15 coda_0.19-4 class_7.3-19
## [133] tidyr_1.1.3 rmarkdown_2.9 DelayedMatrixStats_1.8.0
## [136] Rtsne_0.15 sf_1.0-0 shiny_1.6.0