Tutorial eficiencia Seurat V3 MCA 250k
Adam Casas
Compilado: 15 de enero del 2021
Adaptado de la guía publicada por el laboratorio Satija.
1 Introducción al tutorial
En este tutorial vamos a mostrar un workflow para manejar datasets muy grandes, haciendo hincapié en el uso eficiente de la memoria y el tiempo de CPU de nuestro sistema. Por tanto, no analizaremos el sentido biológico de los clusters de células resultantes del workflow, aunque si te sientes preparado, te invitamos a que los analices por tu cuenta. Los análisis que realicemos en este tutorial se harán en la memoria RAM, aunque Seurat ha añadido recientemente soporte para realizar el análisis desde el disco duro mediante el framework loom. En esta viñeta puedes ver este mismo workflow, pero usando loomR.
2 Workflow: Creamos el objeto Seurat mca
El dataset MCA (Mouse Cell Atlas) fue generado por Xiaoping y compañeros, 2018, contiene 250000 células y está disponible en este enlace. No obstante, para agilizar el workflow partimos directamente de un archivo .rds que contiene la matriz de expresión (=matriz de conteos UMI) y los metadatos ya combinados (archivo disponible en este enlace)
##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(patchwork)
library(ggplot2)
matriz.conteo.UMIs <- readRDS(file = "./MCA_merged_mat.rds")
metadatos <- read.csv(file = "./MCA_All-batch-removed-assignments.csv", row.names = 1)Vamos a analizar las 242000 células a las que se les asignó un clúster y fueron anotadas en el estudio original del que provienen los datos. Por tanto no hace falta hacer el control de calidad (cosa que sí hicimos en el tutorial de las 2700 PBMCs).
Nótese que en la matriz de conteos hay 405192 células (=columnas), pero en los metadatos hay menos células, ya que no todas fueron asignadas a un clúster y anotadas. Por tanto vamos a quedarnos sólo con aquellas células que hayan sido aisgnadas a un clúster con el comando subset del paquete Seurat.
mca <- CreateSeuratObject(counts = matriz.conteo.UMIs, project = "Mouse-Cell-Atlas",
meta.data = metadatos, assay = "RNA")
mca <- subset(x = mca, cells = which(!is.na(mca$ClusterID)))
mca## An object of class Seurat
## 39855 features across 242533 samples within 1 assay
## Active assay: RNA (39855 features, 0 variable features)2.1 Preprocesado de datos
Los siguientes pasos se benefician de ser paralelizados en cuanto a tiempo de CPU se refiere, pero al estar trabajando con un dataset muy pesado y teniendo en cuenta que los requerimientos de RAM escalan con el nº de hilos que usemos, se vuelve necesario tener un sistema con mucha memoria RAM para poder paralelizar este workflow. Dado que mi ordenador de trabajo cuenta sólo con 16 GB, no usaré el siguiente código, pero siéntete libre de usarlo si cuentas con más RAM.
library(future)
plan(strategy = "multisession", workers = 4)
# Establecemos el límite de tamaño que puede usar `future` en 2 GBs dado que de
# lo contrario superamos el límite por defecto de 500 MBs al trabajar con objetos
# tan grandes
options(future.globals.maxSize = 2 * 1024^3)2.1.1 Normalización
Normalizamos los datos con el método estándar LogNormalize, tal como hicimos en el tutorial de 2700 PBMCs.
mca <- NormalizeData(mca, assay = "RNA", normalization.method = "LogNormalize", scale.factor = 10000)2.1.2 Detección de genes diferencialmente expresados
El comando FindVariableGenes calcula la varianza y la media de cada gen del dataset y guarda los genes más variables en mca[["RNA]]@meta.features). Para datasets grandes que usen matrices de conteos de UMIs, seleccionar genes muy variables basándose solo en el ratio varianza-media (VMR en inglés) es una estrategia eficiente y robusta. En este caso seleccionamos los 1000 genes con mayor varianza para posteriores análisis.
2.1.3 Escalado y control de contaminación mitocondrial
Ahora calculamos la expresión de genes mitocondriales de las células con el comando PercentageFeatureSet, y con la función ScaleData escalamos los genes diferencialmente expresados (opción por defecto) y controlamos la contaminación mitocondrial para evitar que esta actúe como un factor de confusión o confounder.
mca[["percent.mt"]] <- PercentageFeatureSet(mca, pattern = "^mt-")
mca <- ScaleData(mca, vars.to.regress = "percent.mt")## Regressing out percent.mt## Centering and scaling data matrixNota: Para datasets grandes, procura paralelizar el workflow de Seurat con el paquete
future(Si tu capacidad RAM te lo permite) y escalar sólo los genes variables. Como ya se comentó en el tutorial de 2700 PBMCs, sólo los genes diferencialmente expresados se usan como input del PCA, por lo que no nos reporta ningún beneficio escalar el resto de genes.
2.2 Reducción de la dimensionalidad (PCA)
## PC_ 1
## Positive: Prm1, Prm2, Fabp9, Meig1, Ldhc
## Negative: Lyz2, S100a9, S100a6, Igkc, Ngp
## PC_ 2
## Positive: Reg3b, Reg3g, Defa24, Lypd8, Defa30
## Negative: Sparc, Col2a1, Col9a2, Col9a1, Col9a3
## PC_ 3
## Positive: Col2a1, Col9a2, Col9a1, Col9a3, Cnmd
## Negative: Mal, Psca, Gkn1, Tff1, Gkn2
## PC_ 4
## Positive: Psca, Tff1, Gkn1, Gkn2, Sptssb
## Negative: Plp1, Mobp, Mag, Ermn, Cldn11
## PC_ 5
## Positive: Lyz2, Camp, Chil3, Elane, Ngp
## Negative: Plp1, Mobp, Mag, Ermn, Cldn11Para elegir la dimensionalidad del dataset, se podría realizar un Jackstraw, pero este dataset en concreto es muy grande, por lo que tal método queda descartado. En su lugar, podemos recurrir al elbow plot.
Podríamos quedarnos con todas las componentes principales que tengan una desviación estándar >=1 (línea roja), lo que nos dejaría con unas 60 componentes principales. No obstante, el laboratorio Satija decidió quedarse con 75, y como almacenar más componentes de las necesarias no afecta negativamente a los análisis, elegimos quedarnos con dichas 75 dimensiones.
Por último, podemos visualizar nuestras células representadas en el espacio del PCA (omitimos la leyenda debido al gargantuesco número de clusters que hay):
2.3 Clustering
Ahora agrupamos las células mediante la técnica de clustering basadas en grafos de Seurat, ya descrita en el tutorial de 2700 PBMCs.
Para el comando FindNeighbors usamos tantas dimensiones como las que hayamos conservado en el PCA (75), y el parámetro nn.eps vamos a establecerlo en \(1.5\). El valor por defecto de nn.eps es 0 (busca vecinos exactos), pero para datasets grandes, se aconseja incrementar este valor (busca vecinos aproximados), ya que así conseguimos acelerar nuestro código.
Respecto al comando FindClusters, podemos elegir el nº de veces que se ejecuta el clustering con el parámetro n.start, ya que dicho parámetro determina el número de comienzos aleatorios del clustering. Al introducir esta variación al comienzo del algoritmo, generamos varios agrupamientos para nuestro dataset, y podemos elegir de todos ellos aquel clustering que genere el grafo con mayor modularidad, o sea, que genere clusters con individuos muy similares entre sí, pero muy distintos de aquellos presentes en otros clusters. El valor por defecto de n.start es 10, pero puedes reducirlo para acelerar el clustering. Por otro lado, se aconseja usar un valor de resolution > 1.2 para datasets de más de 3000 células (ya que de lo contrario generaremos muy pocos clusters en comparación al tamaño del dataset).
mca <- FindNeighbors(mca, reduction = "pca", dims = 1:75, assay = "RNA", nn.eps = 0.5) # `nn.eps = 0` sería más lento ## Computing nearest neighbor graph## Computing SNNmca <- FindClusters(mca, resolution = 3, n.start = 10) # resolution > 1.2 ya que el dataset tiene > 3000 células## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 242533
## Number of edges: 11061154
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9356
## Number of communities: 178
## Elapsed time: 94 seconds2.4 Visualización de clusters en tSNE y UMAP
Ahora visualizamos nuestro dataset en espacios dimensionales reducidos no lineales con RunTSNE y RunUMAP, respectivamente. Usamos como input de dichos comandos las componentes principales que conservamos en el paso del PCA.
Cabe destacar que si bien UMAP es una técnica de reducción de la dimensionalidad más eficiente y rápida que tSNE, también es cierto que ahora Seurat soporta el uso de FIt-SNE, una implementación del laboratorio Kluger de t-SNE acelerada por la transformada rápida de Fourier. Esta implementación es más rápida que el método que usa RunTSNE por defecto, Rtsne, pero para poder aprovecharlo, necesitas instalar el software requerido. Afortunadamente, el repositorio incluye indicaciones claras para instalarlo en función de tu sistema operativo. Una vez instalado, puedes seleccionar la implementación a través del parámetro tsne.method. Además, el método FIt-SNE permite ser paralelizado sin necesidad de usar el paquete future. Aumentar el valor del parámetro max_iter del método FIt-SNE puede ayudar a mejorar la resolución de los clusters.
En cuanto a UMAP se refiere, el comando RunUMAP en Seurat V3.2.2 usa por defecto el método nativo de R uwot, por lo que no es necesario instalar paquetes adicionales. No obstante podría ser útil tomarse las molestias de instalar las dependencias necesarias para poder usar el método de Python umap-learn debido a que dicho wrapper permite el paralelizado de la computación en CUDA, aumentando drásticamente la velocidad de ejecución de nuestro código. Dado que mi máquina de trabajo no tiene una GPU compatible con CUDA, mostraré el código pero no lo ejecutaré. Nótese que al usar el método umap-learn, debemos cambiar también la métrica de cosine (valor por defecto) a correlation. Para más detalles sobre cómo ejecutar código paralelizable en CUDA en R, puedes revisar los apartados 1.1 y 1.2 de mi informe al respecto. Finalmente, puedes ajustar los parámetros n.neighbors y min.dist de la función RunUMAP para mejorar la resolución de los clusters en datasets grandes.
# FIt-SNE
# Los parámetros `nthreads` y `max_iter` son argumentos específicos del método
# `FIt-SNE`. Para más información sobre los argumentos que acepta `FIt-SNE`,
# examine la función `fftRtsne` del archivo fast_tsne.R
mca <- RunTSNE(mca, reduction = "pca", dims = 1:75, tsne.method = "FIt-SNE", nthreads = 4,
max_iter = 2000)## Using rsvd() to compute the top PCs for initialization.# UMAP
mca <- RunUMAP(mca, reduction = "pca", dims = 1:75, umap.method = "uwot", metric = "cosine",
min.dist = 0.75)
## 23:10:25 UMAP embedding parameters a = 0.2734 b = 1.622
## 23:10:26 Read 242533 rows and found 75 numeric columns
## 23:10:26 Using Annoy for neighbor search, n_neighbors = 30
## 23:10:26 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 23:11:13 Writing NN index file to temp file C:\Users\Adam_\AppData\Local\Temp\Rtmpq8wquT\file19103d326e5b
## 23:11:13 Searching Annoy index using 1 thread, search_k = 3000
## 23:12:29 Annoy recall = 100%
## 23:12:30 Commencing smooth kNN distance calibration using 1 thread
## 23:12:39 Initializing from normalized Laplacian + noise
## 23:15:50 Commencing optimization for 200 epochs, with 11349282 positive edges
## 23:19:20 Optimization finished# UMAP en CUDA (No se ejecuta en este informe)
# Necesario para usar el método 'umap-learn', independientemente de si usas CUDA
# o no:
reticulate::py_install(packages = "umap-learn")
mca <- RunUMAP(mca, reduction = "pca", dims = 1:75, umap.method = "umap-learn", metric = "correlation",
min.dist = 0.75)# Gráfico FIt-SNE crudo
grafico_FItsne <- DimPlot(mca, reduction = "tsne", label = F) + ggtitle("FIt-SNE") +
NoLegend()
grafico_FItsne
# Gráfico UMAP crudo
grafico_umap <- DimPlot(mca, reduction = "umap", label = F) + ggtitle("UMAP") + NoLegend()
grafico_umap
# Expresión génica sobre FIt-SNE
expresion_genica_FItsne <- FeaturePlot(mca, features = c("S100a9", "Sftpc"), dims = c(1,
2), reduction = "tsne", combine = F)
expresion_genica_FItsne <- lapply(X = expresion_genica_FItsne, FUN = function(x) AugmentPlot(x +
DarkTheme())) # Aplicamos fondo oscuro y vectorizamos el gráfico
# expresion_genica_FItsne # Para ver los graficos por separado
wrap_plots(expresion_genica_FItsne) # Combinamos ambos gráficos en uno solo
# Expresión génica sobre UMAP
expresion_genica_umap <- FeaturePlot(mca, features = c("S100a9", "Sftpc"), dims = c(1,
2), reduction = "umap", combine = F)
expresion_genica_umap <- lapply(X = expresion_genica_umap, FUN = function(x) AugmentPlot(x +
DarkTheme())) # Aplicamos fondo oscuro y vectorizamos el gráfico
# expresion_genica_umap # Para ver los graficos por separado
wrap_plots(expresion_genica_umap) # Combinamos ambos gráficos en uno solo
Finalmente procedemos a guardar nuestro trabajo.
2.5 Identificación tipos celulares
Si bien el archivo de origen con el que hemos trabajado ya tenía las células anotadas en función de su tejido origen, desconocemos sus tipos celulares. Por ello, podemos usar el paquete clustifyr para anotar nuestras células con mayor detalle.
El paquete clustifyr no incluye el material de referencia de nuestro dataset, pero incluye las URLs de los materiales de referencia almacenados en el github de su desarrollador. Podemos ver dichos enlaces con el siguiente comando:
# Puedes encontrar más información sobre los materiales de referencia en el
# github de clustifyrdata: https://github.com/rnabioco/clustifyrdata
library(clustifyr)
downrefsUna vez localizado el archivo de referencia de nuestro dataset, copiamos su URL y lo descargamos con download.file. Acto seguido lo cargamos en nuestro Environment.
# download.file('https://github.com/rnabioco/clustifyrdata/raw/master/data/ref_MCA.rda',
# destfile = './ref_MCA.rda')
# Cargamos el material de referencia para el dataset Mouse Cell Atlas
load(file = "ref_MCA.rda")Ahora podemos proseguir con el workflow de clustifyr.
# Calculamos el coeficiente de correlación de Spearman entre los clusters y los
# tipos celulares disponibles
correlaciones <- clustify(input = mca[["RNA"]]@data, metadata = mca@meta.data, cluster_col = "seurat_clusters",
ref_mat = ref_MCA, query_genes = mca@assays$RNA@var.features, obj_out = T, seurat_out = T)A continuación usamos el comando cor_to_call para obtener el tipo celular más probable de cada clúster y ordenamos alfabéticamente el dataframe resultante en función del nombre de los clusters.
# Mostramos los tipos celulares más probables de nuestros clusters
asignaciones <- cor_to_call(correlaciones, metadata = mca@meta.data, cluster_col = "seurat_clusters")
head(asignaciones) # Desordenado# Ordenamos los clusters del objeto `asignaciones` por orden alfabético
# install.packages('remotes') remotes::install_github('jmw86069/jamba')
library("jamba")
asignaciones <- mixedSortDF(asignaciones, byCols = 1) # mixedSortDF evita 1, 11, 12... 19, 2, 20, 21...
head(asignaciones) # OrdenadoActo segido asignamos a los clusters el tipo de célula que representan y los visualizamos sobre el UMAP (también podemos visualizarlos sobre el FIt-SNE de así desearlo).
# Hacemos un vector con el tipo celular de los clusters ordenados
tipo.celular.clusters <- asignaciones$type
# Ordenamos los clusters en el objeto mca@meta.data para que esté ordenado de la
# misma manera que el objeto tipo.celular.clusters
mca@meta.data <- mixedSortDF(mca@meta.data, byCols = 10) # columna 10 = seurat_clusters
# Asignamos a los clusters los tipos celulares predichos por `clustifyr`
names(tipo.celular.clusters) <- levels(mca)
mca <- RenameIdents(mca, tipo.celular.clusters)
grafico_umap_clustify <- DimPlot(mca, reduction = "umap", label = T, repel = T, label.box = T,
label.size = 2, pt.size = 0.5) + theme(legend.text = element_text(size = 10))
# Visualizamos el gráfico
grafico_umap_clustify <- AugmentPlot(plot = grafico_umap_clustify, dpi = 600, width = 15,
height = 15)
grafico_umap_clustify # Las etiquetas podrian verse mejor...
Tiempo de ejecución e información del sistema
## Este informe ejecutó todos los comandos en modo secuencial excepto por el comando RunTSNE con la implementación FIt-SNE (4 hilos) debido a la insuficiente capacidad RAM del ordenador empleado para el dataset en cuestión, y tardó en compilarse 55.21193 mins .
##
## El SO empleado fue Windows 10 x64 compilación 19042, y el hardware se compone de un ordenador de sobremesa con un procesador Ryzen 9 3900X @ 3.8GHz, 16GB de RAM y una GPU R9 Fury (no compatible con CUDA). ¿Puedes mejorarlo?