Tutorial clustering Seurat V3 2700 PBMCs
Adam Casas
Compilado: 26 de abril del 2021
Adaptado de la guía gratuita publicada por 10X Genomics.
1 Introducción a Seurat V3
Seurat es un paquete de R diseñado para trabajar con datos de single-cell RNA-seq (o secuenciación de ARN de célula única, en español). Con este paquete podrás realizar control de calidad, exploración y análisis de los ya mencionados datos. Seurat permite al usuario identificar e interpretar fuentes de heterogeneidad en mediciones de transcriptómica de célula única, e integrar diversos tipos de datos de célula única.
Si usas este paquete en tus investigaciones, sus desarrolladores recomiendan citar las siguientes fuentes:
Seurat hace especial énfasis en las visualizaciones sencillas, interpretables y agradables a la vista, y han sido diseñadas para su fácil manejo tanto para investigadores computacionales como experimentales ( i.e. investigadores in-silico e in-vitro).
Este paquete ha sido desarrollado y es actualmente mantenido por Andrew Butler, Paul Hoffman, Tim Stuart, Christoph Hafemeister, y Shiwei Zheng, todos ellos investigadores del laboratorio Satija, en colaboración con otros tantos contribuyentes tales como Jeff Farrell y Karthik Shekhar, por mencionar algunos.
2 Workflow en Seurat para clustering de datos scRNA-seq
Para introducirnos en el workflow de Seurat para el clustering de datos, analizaremos el dataset de 2700 células mononucleares de sangre periférica (PBMCs en inglés), publicado por 10x Genomics. En este dataset se secuenciaron mediante técnicas de single-cell RNA-seq 2700 células en la plataforma NextSeq 500 de Illumina.
Nota: Las PBMCs corresponden a linfocitos y monocitos.
En este tutorial se trabajará el control de calidad o QC, filtrado de datos, cálculo de genes con elevada varianza, reducción de dimensionalidad, clustering basado en grafos e identificación de marcadores de clusters.
2.1 Creamos el objeto Seurat pbmc
Para comenzar el tutorial, primero necesitamos descargar las lecturas crudas de la plataforma NextSeq 500, las cuales están disponibles en mi repositorio de GitHub en forma de archivo comprimido, y ubicarlas en la misma carpeta que este script.
Empezamos cargando las librerías a emplear y leyendo los datos crudos. El comando Read10X lee el output del pipeline “Cell Ranger” de 10X Genomics y devuelve una matriz de conteo de Identificadores Moleculares Únicos (o UMI en inglés). Los valores en esta matriz representan el nº de moléculas para cada característica ( i.e. gen; fila) que son detectados en cada célula (columna).
A continuación usamos la matriz de conteos para crear un objeto de tipo Seurat. Este objeto sirve a modo de contenedor tanto para los datos (la matriz de conteos, por ejemplo) como para los análisis (resultados del PCA o un clustering) de un dataset de célula única. Para una descripción técnica de la estructura de los objetos de tipo Seurat, consulta la Wiki que han habilitado sus desarrolladores en GitHub.
# Nota: El paquete `future` permite paralelizar los comandos `NormalizeData`,
# `ScaleData`, `JackStraw`, `FindMarkers`, `RunUMAP` y quizás 'FindClusters'
# (solo si hacemos clustering con resolutions mayores; nosotros usamos
# `resolution = 0.5`)
library(future)
plan(strategy = "multisession", workers = 4)
# Nota: la estrategia 'multiprocess' de la guía original está obsoleta. En su
# lugar se usa la estrategia 'multicore' (no soportada en Windows) o
# 'multisession' (soportada en Windows)
# Comprobamos el parámetro 'estrategy' y el nº de hilos de nuestro plan actual
# con el comando `plan`
plan()## multisession:
## - args: function (..., envir = parent.frame(), workers = 4)
## - tweaked: TRUE
## - call: plan(strategy = "multisession", workers = 4)# Librerías
library(dplyr)
library(Seurat)
library(patchwork)
library(curl)
# La carpeta de trabajo será la ubicación de este script
setwd("./")
# Descargamos y descomprimimos las lecturas crudas:
url_lecturas <- "https://github.com/gloknar/R-Utilities/raw/master/Tutorial%20Seurat%202700%20PBMCs/pbmc3k_filtered_gene_bc_matrices.tar.gz"
curl_download(url_lecturas, destfile = "pbmc3k_filtered_gene_bc_matrices.tar.gz")
untar("pbmc3k_filtered_gene_bc_matrices.tar.gz", exdir = "./")
# Cargamos la matriz de conteos de nuestras 2.700 PBMCs y visualizamos los 10
# primeros genes y células
datos.pbmc <- Read10X(data.dir = "./filtered_gene_bc_matrices/hg19/")
head(datos.pbmc[1:10, 1:10])## 6 x 10 sparse Matrix of class "dgCMatrix"
##
## MIR1302-10 . . . . . . . . . .
## FAM138A . . . . . . . . . .
## OR4F5 . . . . . . . . . .
## RP11-34P13.7 . . . . . . . . . .
## RP11-34P13.8 . . . . . . . . . .
## AL627309.1 . . . . . . . . . .# Inicializamos el objeto Seurat con los datos crudos (no normalizados)
pbmc <- CreateSeuratObject(counts = datos.pbmc, project = "2700-PBMCs", assay = "RNA",
min.cells = 3, min.features = 200)
pbmc## An object of class Seurat
## 13714 features across 2700 samples within 1 assay
## Active assay: RNA (13714 features, 0 variable features)¿Qué apariencia tienen los datos en la matriz de conteo?
La previamente mencionada matriz de conteos de UMIs se guardó en el objeto datos.pbmc. Recordemos que cada fila de esta matriz corresponde a un gen, y cada columna corresponde al nº de transcritos detectados en cada célula, de manera que la primera columna corresponde a los transcritos detectados en la célula nº1, la segunda columna corresponde a los transcritos de la célula nº2, y así.
# Veamos la expresión de 4 genes concretos en las primeras 30 células
datos.pbmc[c("SERPINC1", "CD3D", "TCL1A", "MS4A1"), 1:30]## 4 x 30 sparse Matrix of class "dgCMatrix"
##
## SERPINC1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
## CD3D 4 . 10 . . 1 2 3 1 . . 2 7 1 . . 1 3 . 2 3 . . . . . 3 4 1 5
## TCL1A . . . . . . . . 1 . . . . . . . . . . . . 1 . . . . . . . .
## MS4A1 . 6 . . . . . . 1 1 1 . . . . . . . . . 36 1 2 . . 2 . . . .# En caso de no encontrar los genes de interés, Seurat V3 incluye el siguiente
# comando para buscar sus sinónimos:
GeneSymbolThesarus(symbols = c("SERPINC1", "CD3D", "TCL1A", "MS4A1"))## MS4A1
## "MS4A2"Los puntos (“.”) de la matriz representan ceros (no se detectaron transcritos de ARNm de ese gen, en esa célula). Así, vemos que en la primera célula sólo se detectaron 4 moléculas de ARNm del gen CD3D, mientras que en la 9ª célula se encontraron transcritos de CD3D, TCL1A y MS4A1.
Dado que la mayoría de valores en una matriz de scRNA-seq suelen ser 0, Seurat usa una “sparse-matrix” o matriz dispersa siempre que sea posible (como es en este caso). Ello disminuye la cantidad de memoria usada y acelera el procesado de los datos.
matriz_normal <- object.size(as.matrix(datos.pbmc))
matriz_normal## 709591472 bytesmatriz_dispersa <- object.size(datos.pbmc)
cat("Al guardar los datos en formato matriz dispersa en vez de matriz normal, conseguimos reducir el tamaño de los datos en un factor de",
matriz_normal[1]/matriz_dispersa[1])## Al guardar los datos en formato matriz dispersa en vez de matriz normal, conseguimos reducir el tamaño de los datos en un factor de 23.728042.2 Control de calidad y selección de células para posteriores análisis
Los pasos mostrados a continuación abarcan el workflow estándar para el preprocesado de datos de scRNA-seq. Estos incluyen la selección y filtrado de células en función del control de calidad, normalizado y escalado de datos y la detección de transcritos con gran varianza.
Seurat te permite explorar fácilmente las métricas del control de calidad y filtrar células según el criterio establecido por el usuario. Algunas de las métricas de calidad más usadas por la comunidad son:
Nº de genes expresados en cada célula
Células mal secuenciadas o gotículas vacías suelen expresar muy pocos genes
Doublets y multiplets (gotículas con más de una célula) suele expresar un nº anormalmente alto de genes
Nº total de transcritos detectados en una célula (lo cual está muy correlacionado con nº de genes expresados)
Porcentaje de transcritos que mapean en el genoma mitocondrial
Las células mal secuenciadas y las moribundas suelen mostrar muchos transcritos provenientes de genes mitocondriales
Para calcular las métricas de calidad mitocondriales, se emplea el comando
PercentageFeatureSet, el cual calcula el porcentaje de transcritos provenientes de un set de genes (en este caso, genes mitocondriales)Los genes cuyo nombre comienzan por “MT-” son genes mitocondriales
Nota: Durante el aislamiento de células individuales (antes de ser secuenciadas), estas son embebidas en gotículas de aceite de manera que, idealmente, en cada gotícula se halle una única célula (ver imagen inferior, sacada de Bach et al., 2017).
# Los operadores [[ ]] permiten añadir la columna 'percent.mt' a los metadatos
# del objeto Seurat 'pbmc', guardados en el bolsillo @meta.data (pbmc@meta.data).
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")¿Dónde se almacenan las métricas del control de calidad en Seurat?
El nº de genes expresados (
nFeature_RNA) y moléculas de ARNm detectadas (nCount_RNA) se calculan automáticamente al llamar a la funciónCreateSeuratObject- En concreto, se almacenan en los metadatos del objeto Seurat
# Mostramos las métricas del control de calidad de las 5 primeras células
head(pbmc@meta.data, 5)En el ejemplo a continuación, visualizamos las métricas de calidad y las empleamos para filtrar las células
Vamos a descartar las células que expresen más de 2500 genes y las que expresen menos de 200
También filtraremos las células cuyos transcritos mitocondriales supongan >=5% del total de transcritos presentes en ellas
# Visualizamos las métricas de calidad en un gráfico de violín
VlnPlot(pbmc, features = c("nFeature_RNA", "percent.mt"), ncol = 3)
# Filtramos las células según los criterios de calidad previamente mencionados
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt <
5)
# Visualizamos las células retenidas por el filtrado
VlnPlot(pbmc, features = c("nFeature_RNA", "percent.mt"), ncol = 3)
2.3 Normalización de datos
Tras filtrar las células indeseadas de nuestro conjunto de datos, el siguiente paso es normalizar los niveles de expresión génica mediante el comando NormalizeData. NormalizeData usa por defecto un factor de escala de 104 y el método LogNormalize, el cual normaliza los niveles de expresión génica de cada célula de la siguiente manera: Para un gen concreto, divide su nº de transcritos entre el nº total de transcritos presentes en esa célula, luego se multiplica por el factor de escala especificado, y al resultado de esa multiplicación (llamémoslo “x”) se le aplica la transformación \(ln(1+x)\) (se usa ln(1+x) en lugar de ln(x) para evitar perder precisión cuando x es un número flotante cercano a 0).
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)Los niveles de expresión génica normalizados de nuestras células se almacenan en pbmc[["RNA"]]@data (también puedes escribirlo como pbmc@assays$RNA@data).
# Nº de transcritos tras el normalizado, se muestran las 3 primeras células
pbmc[["RNA"]]@data[968:974, 1:3]## 7 x 3 sparse Matrix of class "dgCMatrix"
## AAACATACAACCAC-1 AAACATTGAGCTAC-1 AAACATTGATCAGC-1
## GON4L . . 1.429744
## SYT11 . . .
## RIT1 . . .
## KIAA0907 . 1.111715 .
## ARHGEF2 . . 1.429744
## RP11-336K24.12 . . .
## SSR2 1.635873 1.111715 1.4297442.4 Identificación de genes con gran varianza (selección de genes)
A continuación calculamos un subconjunto de genes cuya expresión génica presente una gran varianza (o sea, que se expresen mucho en algunas células, y que apenas se expresen en otras). 10X Genomics y otros tantos investigadores han descubierto que centrarse en el análisis de estos genes ayuda a resaltar señales biológicas en datasets de célula única.
El método vst del comando FindVariableFeatures modela la relación media-varianza de los genes para detectar aquellos con mayor varianza, y por defecto devuelve 2.000 genes. La implementación de dicho comando en Seurat se debe al trabajo de 10X Genomics.
# Tras usar `FindVariableFeatures`, los genes variables se guardan en
# pbmc[['RNA]]@meta.features
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
# Identificamos los 10 genes más variables con el comando `VariableFeatures`
top10 <- head(VariableFeatures(pbmc), 10)
# Graficamos los genes más variables con sus nombres
grafico_genes_variable <- VariableFeaturePlot(pbmc)
LabelPoints(grafico_genes_variable, points = top10, repel = TRUE)
Nota: Fíjate en que el gráfico generado se asemeja, en diseño y concepto, a los gráficos MA empleados en transcriptómica.
2.5 Escalado de datos
Después de la normalización de los datos, los escalamos. Este es un paso estándar y previo a la aplicación de técnicas de reducción de dimensionalidad (también llamadas en su conjunto como “extracción de características” en el ámbito de Machine Learning).
El comando ScaleData:
Cambia la expresión de cada gen de manera que su expresión media en todas las células sea 0
Escala la expresión de cada gen, de manera que su varianza a lo largo de todas las células sea 1
- Este paso previene que genes con expresión constitutiva tengan más peso en los próximos análisis que genes con escasa expresión (por ejemplo, aquellos que codifican para factores de transcripción)
Los datos escalados se guardan en
pbmc[["RNA"]]@scale.dataSólo nos interesa escalar los genes variables, ya que los usaremos como input del PCA. Por defecto,
ScaleDataescala sólo los genes variables identificados en el paso anterior- Si aún así quieres escalar todos los genes, debes pasarle a la función el nombre de todos los genes a través del parámetro
features(ver código comentado)
- Si aún así quieres escalar todos los genes, debes pasarle a la función el nombre de todos los genes a través del parámetro
# Si quieres escalar todos los genes, usa:
# nombres_genes <- rownames(pbmc) pbmc <- ScaleData(pbmc, features =
# nombres_genes)
# Para escalar sólo los genes variables que usaremos en el PCA (valor por
# defecto):
pbmc <- ScaleData(pbmc)Nota: Si escalas sólo los 2000 genes variables, el comando
ScaleDatatardará menos tiempo en ejecutarse que si escalas todos los genes disponibles. No obstante, para ver heatmaps no distorsionados, es necesario escalar todos los genes. Nosotros sólo escalamos esos 2000 genes variables, así que tenlo en cuenta cuando uses el comandoDoHeatmap(ver apartado 2.10 Encontrar genes diferencialmente expresados (biomarcadores de clusters)). El siguiente paso, donde usamos el comandoDimheatmap, no se verá afectado dado que dicha función sólo trabaja con los genes variables usados como input del PCA. En resumen: Escala todos los genes si usasDoHeatmap; escala solo los genes variables si usasDimheatmap.
¿Cómo puedo eliminar fuentes de varianza indeseadas?
Puedes eliminar fuentes de varianza indeseables para el experimento, tales como la contaminación mitocondrial o la varianza asociada al ciclo celular usando el parámetro vars.to.regress del comando ScaleData. Por ejemplo:
pbmc <- ScaleData(pbmc, vars.to.regress = "percent.mt")No obstante, el laboratorio Satija recomienda encarecidamente a los usuarios avanzados que deseen aprovechar esta funcionalidad usar la nueva función de normalización SCTransform. Dicho comando está descrito en este artículo del laboratorio Satija, y además también han publicado una viñeta donde usan tal función en Seurat V3. Al igual que ScaleData, la función SCTransform también incluye el parámetro vars.to.regress.
2.6 Reducción de la dimensionalidad del dataset (PCA)
Seguidamente realizamos un PCA a los datos escalados. Por defecto, solo se le aplica el PCA a los genes variables previamente seleccionados, pero puedes añadir más genes al análisis si quieres mediante el argumento features de la función RunPCA. El resultado se almacena en pbmc[["pca"]].
pbmc <- RunPCA(pbmc)Podemos visualizar el PCA con el comando DimPlot y seleccionando en el parámetro reduction la técnica empleada, que en este caso es el pca.
DimPlot(pbmc, reduction = "pca", dims = c(1, 2))
El comando DimHeatmap en particular ofrece una fácil exploración de las fuentes primarias de heterogeneidad en el dataset en cuestión, y nos puede ayudar a determinar qué componentes principales deberíamos analizar. Tanto los genes como las células son ordenados en función de sus loadings del PCA. Podemos seleccionar el nº de células a analizar en cada heatmap, de manera que el heatmap no tarde demasiado en dibujarse en pantalla cuando trabajemos con datasets grandes.
DimHeatmap(pbmc, dims = 1:6, cells = 500)
Nótese que se podría aplicar un análisis de enriquecimiento génico a los resultados de estos heatmaps, de manera que podamos identificar las facetas biológicas que describen cada componente principal.
2.7 Determinar la dimensionalidad del dataset
Una vez realizado el PCA, y sabiendo que cada componente principal explica una faceta distinta de las poblaciones celulares estudiadas, nos puede surgir la siguiente pregunta: ¿con cuántas componentes principales debería quedarme? ¿10, 20, 50?
Existen varias maneras o criterios a seguir a la hora de elegir el nº adecuado de dimensiones que estudiar. No obstante, en este tutorial nos centraremos en el método implementado por 10X Genomics en Macosko et al., 2015. En dicha publicación aplican exitosamente un test de remuestreo inspirado en la técnica JackStraw (Chung y Storey, 2015).
Dicho test consiste en permutar aleatoriamente un subconjunto del dataset (el 1% por defecto), volver a ejecutar el PCA (con lo que se genera una distribución nula de los loadings de los genes en el PCA) y repetir el proceso para generar una distribución nula de puntuaciones de genes. Se identifican como dimensiones significativas aquellas enriquecidas con muchos genes cuyos p-valores sean bajos ( i.e. aquellas CPs por encima de la línea discontinua que representa la distribución uniforme).
# Evaluamos la velocidad obtenida por paralelizar el comando `JackStraw` con el
# paquete `future`.
comparacion_tiempo <- data.frame(funcion = character(), tiempo = numeric(), strategy = character())
# Mononúcleo
plan("sequential")
start <- Sys.time()
pbmc <- JackStraw(pbmc, dims = 5)
end <- Sys.time()
comparacion_tiempo <- rbind(comparacion_tiempo, data.frame(funcion = "JackStraw",
tiempo = as.numeric(end - start, units = "secs"), strategy = "sequential"))
# Multinúcleo (4 hilos)
plan("multisession", workers = 4)
start <- Sys.time()
pbmc <- JackStraw(pbmc, dims = 5)
end <- Sys.time()
comparacion_tiempo <- rbind(comparacion_tiempo, data.frame(funcion = "JackStraw",
tiempo = as.numeric(end - start, units = "secs"), strategy = "multisession"))
# Nota: El JackStraw con 20 dimensiones puede tardar en completarse
pbmc <- JackStraw(pbmc, dims = 20)
pbmc <- ScoreJackStraw(pbmc, dims = 1:20)La función JackStrawPlot nos permite comparar de manera visual la distribución de los p-valores de cada componente principal con una distribución uniforme (línea discontínua). Las componentes principales significativas muestran un gran enriquecimiento de genes con p-valores bajos y se situarán muy por encima de la línea discontínua. En este caso se observa que el p-valor aumenta bastante (¡dos órdenes de magnitud!) entre la decimotercera y la decimocuarta componente principal.
# Visualizamos las componentes principales significativas
JackStrawPlot(pbmc, dims = 1:20)
Si bien el método JackStraw es efectivo, también es cierto que es computacionalmente costoso, sobre todo en datasets muy grandes, en cuyo caso uno puede optar por usar el Scree test/método del codo. En este caso, nos quedaríamos con las 8 primeras dimensiones, ya que la pendiente se acerca a 0 a partir de la octava componente principal:
# Usamos el comando `geom_vline` en vez de `abline` dado que la función
# 'ElbowPlot' devuelve un gráfico de ggplot2
ElbowPlot(pbmc) + geom_vline(xintercept = 8, color = "red", linetype = "dashed")
Identificar la verdadera dimensionalidad del dataset puede resultar ser una tarea laboriosa e incierta, y existen varios métodos para responder a esta cuestión. En nuestro caso, los métodos probados ofrecieron un nº parecido de dimensiones a estudiar. No obstante, los resultados de los próximos análisis no suelen variar mucho si cogemos más dimensiones de las necesarias (¡pero se ven muy perjudicados si usamos un nº insuficiente de dimensiones!), por lo que 10X Genomics recomienda tender hacia la parte superior del espectro. En nuestro caso, podríamos elegir conservar entre 8 y 14 componentes principales, por lo que el criterio de 10X Genomics nos puede encaminar hacia la elección de 12 componentes principales, por ejemplo.
2.8 Agrupar las células (Clustering basado en grafos)
Seurat V3 aplica un clustering basado en grafos, desarrollado a partir de las estrategias empleadas en el artículo previamente mencionado de Macosko et al., 2015 e inspirado en otros artículos donde se emplea clustering basado en grafos a datos de scRNA-seq SNN-Cliq, Xu y Su, Bioinformatics, 2015 y CyTOF PhenoGraph, Levine et al., Cell, 2015. En resumidas cuentas, estos métodos insertan las células en una estructura de forma de grafo (por ejemplo, un grafo basado en K-nearest neighbours (KNN) ), con los enlaces dibujados entre células con patrones de expresión génica similares, y luego intentan particionar este grafo en “comunidades” o clústers densamente interconectados ( quasi-cliques).
Al igual que el método de clustering PhenoGraph, primero construimos un grafo de KNN basado en la distancia euclídea en el espacio del PCA, y refinamos los pesos de los enlaces entre dos células, basándonos en el solape compartido entre sus vecindades locales (Similaridad/Distancia de Jaccard). Este paso se realiza con el comando FindNeighbors, y toma como input la dimensionalidad del dataset previamente definida (este comando usa por defecto las 10 primeras componentes principales).
Para agrupar las células en clusters, aplicamos seguidamente técnicas de optimizacion de la modularidad tales como el método de Louvain o SLM (SLM, Blondel et al., Journal of Statistical Mechanics), para agrupar células iterativamente con el objetivo de optimizar la función estándar de modularidad. El comando FindClusters implementa este procedimiento, y contiene el parámetro resolution(resolución) que establece la “granularidad” de los clusters, y a mayores valores de resolution, más clusters generamos. 10X Genomics recomienda usar una resolution de entre 0,4 y 1,2 para datasets de single-cell de 3000 células. Para datasets con más de 3000, se recomienda usar valores más altos. En nuestro caso usamos resolution = 0.5 dado que nuestro dataset es de menos de 3000 células polimorfonucleadas de sangre periférica.
El parámetro algorithm del comando FindClusters toma valores entre 1 y 4, y sirve para elegir el algoritmo de optimización de la modularidad: el método original de Louvain es el valor por defecto (algorithm = 1), mientras que el algoritmo SLM corresponde a algorithm = 3.
Los clusters con sus respectivas células se guardan en pbmc@active.ident, pbmc$RNA_snn_res.0.6 y pbmc$seurat_clusters. Alternativamente, podemos acceder a dicha información con el comando Idents(pbmc).
pbmc <- FindNeighbors(pbmc, dims = 1:12)
pbmc <- FindClusters(pbmc, resolution = 0.5, algorithm = 1)## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 2638
## Number of edges: 100911
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8759
## Number of communities: 9
## Elapsed time: 0 seconds# Observamos las IDs de los clusters a los que pertenecen las 5 primeras células.
head(Idents(pbmc), 5)## AAACATACAACCAC-1 AAACATTGAGCTAC-1 AAACATTGATCAGC-1 AAACCGTGCTTCCG-1
## 4 3 1 2
## AAACCGTGTATGCG-1
## 6
## Levels: 0 1 2 3 4 5 6 7 82.9 Reducción de dimensionalidad con métodos no lineales(UMAP/tSNE)
Seurat ofrece varias técnicas de reducción de dimensionalidad tales como t-SNE y UMAP para visualizar y explorar el dataset.
El comando RunUMAP permite ejecutar un UMAP en Seurat. 10X Genomics recomienda tanto para t-SNE como UMAP buscar el mismo nº de dimensiones que las establecidas en el paso del PCA (8 en nuestro caso).
# Si no tienes insalado el paquete `UMAP`, puedes instalarlo con:
# reticulate::py_install(packages = 'umap-learn')
pbmc <- RunUMAP(pbmc, dims = 1:8)# Puedes usar el parámetro 'label` del comando `DimPlot` para visualizar mejor la
# correspondencia a cada clúster
DimPlot(pbmc, reduction = "umap", label = T) # label = 'FALSE' por defecto
2.10 Encontrar genes diferencialmente expresados (biomarcadores de clusters)
Seurat puede ayudarte a encontrar biomarcadores (i.e. genes diferencialmente expresados) que definen clusters. Por defecto, identifica marcadores positivos (sobreexpresados) y negativos (infraexpresados) de un único cluster (especificado en el parámetro ident.1), comparado con el resto de células. El comando FindAllMarkers automatiza este proceso para todos los clusters, pero también puedes testar grupos de clusters entre ellos, o testar un clúster contra el resto de células.
Por defecto, el argumento min.pct hace que FindAllMarkers testeé solo aquellos genes expresados como mínimo en el 10% de las células de cualquiera de los grupos celulares estudiados, y el parámetro logfc.threshold selecciona para testar aquellos genes cuya diferencia de expresión entre dos grupos celulares sea de 0.25-fold (en escala logarítmica). Puedes poner ambos argumentos a 0, pero el comando irá mucho más lento (ya que se testarán muchos más genes, la mayoría de los cuales probablemente sean solo ruido). Otra opción para acelerar el comando es ajustar el valor del parámetro max.cells.per.ident, el cual establece el nº máximo de células a calcular por clúster (por defecto usa todas las células disponibles).
Si ajustamos estos 3 parámetros para reducir el tiempo de computación necesario, puede que perdamos un poco de potencia estadística, pero el ahorro en tiempo probablemente será significativo, y los genes más diferencialmente expresados (i.e. fold > 2) seguramente seguirán siendo detectados.
FindMarkers busca los marcadores del clúster(s) que le digas al comando, mientras que FindAllMarkers busca los marcadores de todos los clusters.
# Encontramos todos los marcadores del clúster 1
cluster1.markers <- FindMarkers(pbmc, ident.1 = 1, min.pct = 0.25)
head(cluster1.markers, n = 5)# Encontramos todos los marcadores que diferencian el clúster 5 de los clusters 0
# y 3
cluster5.markers <- FindMarkers(pbmc, ident.1 = 5, ident.2 = c(0, 3), min.pct = 0.25)
head(cluster5.markers, n = 5)# Encontramos los marcadores para cada clúster comparado con el resto de células
# y reportamos solo los marcadores positivos (`only.pos = T`).
pbmc.markers <- FindAllMarkers(pbmc, min.pct = 0.25, logfc.threshold = 0.25, only.pos = T)
pbmc.markers %>%
group_by(cluster) %>%
top_n(n = 2, wt = avg_log2FC)Seurat tiene varios tests para evaluar la expresión diferencial de los marcadores y puedes seleccionar el test a usar con el parámetro test.use del comando FindMarkers (mira la viñeta de expresión diferencial del laboratorio Satija para más detalles). Por ejemplo, el test de la curva ROC devuelve el poder de clasificación de cada marcador (yendo desde 0 (aleatorio) hasta 1 (perfecto) ).
cluster1.markers <- FindMarkers(pbmc, ident.1 = 0, logfc.threshold = 0.25, test.use = "roc",
only.pos = TRUE)Los comandos de Seurat más comunes para visualizar la expresión de biomarcadores de clusters son VlnPlot y FeaturePlot (visualiza la expresión génica sobre un gráfico UMAP, t-SNE o PCA). 10X Genomics recomienda también visualizar tu dataset con la ayuda de los comandos RidgePlot, CellScatter y DotPlot.
VlnPlot(pbmc, features = c("MS4A1", "CD79A"))
# También puedes graficar el nº de transcritos detectados. Ponemos el parámetro
# 'log = T' dado que se detectaron muchos transcritos. El parámetro `slot` hace
# referencia a los bolsillos ubicados en pbmc@assays$RNA@ , i.e.
# pbmc@assays$RNA@counts
VlnPlot(pbmc, features = c("NKG7", "PF4"), slot = "counts", log = T)
# Featureplot
FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A",
"LYZ", "PPBP", "CD8A"))
Ver gráficos adicionales
# Ridgeplot. Una visualización alternativa al gráfico de violín
RidgePlot(pbmc, features = "NKG7")
# Cellscatter. Compara los niveles de expresión génica entre dos células y
# muestra el coeficiente de correlación de Pearson
CellScatter(pbmc, cell1 = "AACGCCCTCGGGAA-1", cell2 = "AACGCCCTGGCATT-1")
# Dotplot de la expresión génica entre clusters (no confundir con el dotplot
# usado en alineamiento de secuencias). El tamaño del punto codifica el
# porcentaje de células del clúster que expresan el gen, mientras que el color
# codifica su nivel de expresión medio.
DotPlot(pbmc, features = c("LYZ", "CD3E"))
El comando DoHeatmap genera un heatmap de expresión génica para las células y genes seleccionados. En este caso graficamos los top 10 biomarcadores de cada clúster (el siguiente código es una modificación del encontrado en la sección 2.6 Reducción de la dimensionalidad del dataset (PCA)):
top10 <- pbmc.markers %>%
group_by(cluster) %>%
top_n(n = 10, wt = avg_log2FC)
DoHeatmap(pbmc, features = top10$gene) # + NoLegend() para eliminar la leyenda
2.11 Asignar el tipo de célula a cada clúster
En el caso de este dataset, podemos usar biomarcadores ya conocidos para asignar el tipo de célula que representa cada clúster:
| ID del clúster | Biomarcadores | Tipo de célula |
|---|---|---|
| 0 | IL7R, CCR7 | Linfocito T CD4+ de memoria |
| 1 | CD14, LYZ | Monocito CD14+ |
| 2 | IL7R, S100A4 | Linfocito T CD4+ naïve |
| 3 | MS4A1 | Linfocito B |
| 4 | CD8A | Linfocito T CD8+ |
| 5 | FCGR3A, MS4A7 | Monocito FCGR3A+ |
| 6 | GNLY, NKG7 | Natural Killer (NK) |
| 7 | FCER1A, CST3 | Célula dendrítica (DC) |
| 8 | PPBP | Megacariocito |
new.cluster.ids <- c("Naïve CD4 T", "Memory CD4 T", "CD14+ Mono", "B Linfo", "CD8 T",
"FCGR3A+ Mono", "NK", "DC", "Megakaryocyte")
names(new.cluster.ids) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, new.cluster.ids)
grafico_umap <- DimPlot(pbmc, reduction = "umap", label = T, pt.size = 0.5) + xlab("UMAP 1") +
ylab("UMAP 2") + theme(legend.text = element_text(size = 10))
grafico_umap
ggsave(filename = "./2700pbmc-umap.png", plot = grafico_umap)Una vez acabado el protocolo, podemos guardar nuestro objeto Seurat pbmc en un archivo .rds para así poder retomar nuestro trabajo en el futuro sin tener que volver a correr los algoritmos computacionalmente pesados del script. Guardar el objeto Seurat también nos permite compartirlo fácilmente con nuestros colaboradores:
saveRDS(pbmc, file = "./2700pbmc_final.rds")¿Cómo podemos clasificar los tipos de células representados en los clusters?
En este tutorial, 10X Genomics nos provee con la clasificación de los tipos celulares representados en cada clúster, pero cuando trabajemos con datasets de single cell, ¿cómo podemos asignar a nuestros clusters el tipo de célula que representan?
Rui Fu et al., 2019 han desarrollado clustifyr, un paquete para R diseñado específicamente para solventar este problema con ayuda de datos de referencia externos. Este paquete nos aporta funciones para anotar automáticamente clusters (o células, si así lo deseamos).
Con el comando clustify podemos calcular el coeficiente de correlación de Spearman (valor por defecto) de cada tipo celular posible para cada uno de nuestros clusters.
library(clustifyr)
# Calculamos el coeficiente de correlación de Spearman entre nuestros clusters y
# los tipos celulares disponibles
correlaciones <- clustify(input = pbmc@assays$RNA@data, # = pbmc[["RNA"]]@data
metadata = pbmc@meta.data,
cluster_col = "seurat_clusters",
ref_mat = cbmc_ref,
query_genes = pbmc@assays$RNA@var.features,
obj_out = T,
seurat_out = T)
correlaciones## B CD14+ Mono CD16+ Mono CD34+ CD4 T CD8 T DC
## 0 0.6536956 0.5287061 0.5683694 0.6534759 0.8973157 0.8774829 0.5440597
## 1 0.6251937 0.5361372 0.5856654 0.6278916 0.8585097 0.8301475 0.5555125
## 2 0.5280145 0.9144364 0.8916168 0.5106973 0.4632686 0.4377198 0.7644301
## 3 0.9095536 0.5350035 0.5911780 0.6249599 0.6403603 0.6259067 0.6041443
## 4 0.5815854 0.4865133 0.5487701 0.5674970 0.7701184 0.7732337 0.5214476
## 5 0.5370668 0.8439737 0.9298160 0.5322852 0.4853877 0.4642787 0.7241151
## 6 0.5014896 0.4436219 0.5129164 0.5379510 0.6725417 0.6874447 0.4835965
## 7 0.6088123 0.7415535 0.7520484 0.5694074 0.4951613 0.4789794 0.8493399
## 8 0.1456751 0.2571662 0.2867421 0.2455947 0.1422081 0.1208273 0.1898078
## Eryth Memory CD4 T Mk Naive CD4 T NK pDCs
## 0 0.5522939 0.7250839 0.2660267 0.7133310 0.6782599 0.4732272
## 1 0.5488699 0.7127095 0.2598159 0.6799978 0.6921106 0.4663640
## 2 0.4559577 0.4551942 0.4226791 0.4561654 0.5907618 0.4744751
## 3 0.4816215 0.4792129 0.2233366 0.5402612 0.5469465 0.5828493
## 4 0.4933390 0.6427027 0.2989328 0.6084035 0.8307007 0.4524823
## 5 0.4605121 0.4607598 0.3940008 0.4568271 0.6036479 0.4784231
## 6 0.4633664 0.5630126 0.2812534 0.5466574 0.8922740 0.4215019
## 7 0.4470039 0.4764929 0.2948557 0.4475783 0.5643991 0.6825388
## 8 0.3068445 0.1457308 0.7319575 0.1609963 0.1396979 0.2152614Ahora usamos el comando cor_to_call para obtener el tipo celular más probable de cada clúster, o sea, aquel tipo celular que muestra el mayor coeficiente de correlación para el clúster en cuestión.
# Mostramos los tipos celulares más probables de nuestros clusters
asignaciones <- cor_to_call(correlaciones, metadata = pbmc@meta.data, cluster_col = "seurat_clusters")
asignaciones <- asignaciones %>%
arrange(seurat_clusters)
asignacionesAhora podemos asignar a los clusters el tipo de célula que representan.
# Hacemos un vector con el nombre de los clusters, ordenados de 0 al 8
nuevas.ids.clusters <- asignaciones$type
# Asignamos a los clusters los tipos celulares predichos por `clustifyr`
names(nuevas.ids.clusters) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, nuevas.ids.clusters)
grafico_umap_clustify <- DimPlot(pbmc, reduction = "umap", label = T, pt.size = 0.5) +
xlab("UMAP 1") + ylab("UMAP 2") + theme(legend.text = element_text(size = 10))
grafico_umap_clustify
Ahora comparamos las anotaciones de 10X Genomics con las predichas por clustifyr. Vemos que las anotaciones de 10X Genomics son más detalladas, no obstante clustifyr no lo ha hecho nada mal. Podría servirnos para una primera toma de contacto y orientarnos para futuros análisis.
Nótese que clustifyr ha cometido un error al clasificar el clúster 4. Este clúster representa linfocitos T CD8+, aunque la anotación automática lo fusionó con el clúster 6 (linfocitos natural killer). Si prestamos atención a los coeficientes de correlación de Spearman, podemos entender el por qué del error: vemos que para el clúster 4, la segunda célula más probable es el linfocito T CD8+ con un r = 0.7732, mientras que el linfocito natural killer tiene un coeficiente ligeramente superior (r = 0.8307).
grafico_umap_clustify <- grafico_umap_clustify + ggtitle("Anotación automática clustifyr")
grafico_umap <- grafico_umap + ggtitle("Anotación manual 10X Genomics")
grafico_umap_clustify/grafico_umap # Usamos operador '/' del paquete `patchwork`
3 sessionInfo()
Click para mostrar
sessionInfo()## R version 4.0.3 (2020-10-10)
## 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] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] clustifyr_1.2.0 cowplot_1.1.1 ggplot2_3.3.3 curl_4.3
## [5] patchwork_1.1.1 SeuratObject_4.0.0 Seurat_4.0.1 dplyr_1.0.5
## [9] future_1.21.0
##
## loaded via a namespace (and not attached):
## [1] fastmatch_1.1-0 plyr_1.8.6
## [3] igraph_1.2.6 lazyeval_0.2.2
## [5] splines_4.0.3 entropy_1.3.0
## [7] BiocParallel_1.24.1 listenv_0.8.0
## [9] scattermore_0.7 GenomeInfoDb_1.26.7
## [11] digest_0.6.27 htmltools_0.5.1.1
## [13] fansi_0.4.2 magrittr_2.0.1
## [15] tensor_1.5 cluster_2.1.2
## [17] ROCR_1.0-11 limma_3.46.0
## [19] globals_0.14.0 matrixStats_0.58.0
## [21] spatstat.sparse_2.0-0 colorspace_2.0-0
## [23] ggrepel_0.9.1 xfun_0.22
## [25] crayon_1.4.1 RCurl_1.98-1.3
## [27] jsonlite_1.7.2 spatstat.data_2.1-0
## [29] survival_3.2-11 zoo_1.8-9
## [31] glue_1.4.2 polyclip_1.10-0
## [33] gtable_0.3.0 zlibbioc_1.36.0
## [35] XVector_0.30.0 leiden_0.3.7
## [37] DelayedArray_0.16.3 future.apply_1.7.0
## [39] SingleCellExperiment_1.12.0 BiocGenerics_0.36.1
## [41] abind_1.4-5 scales_1.1.1
## [43] DBI_1.1.1 miniUI_0.1.1.1
## [45] Rcpp_1.0.6 viridisLite_0.4.0
## [47] xtable_1.8-4 reticulate_1.19
## [49] spatstat.core_2.1-2 stats4_4.0.3
## [51] htmlwidgets_1.5.3 httr_1.4.2
## [53] fgsea_1.16.0 RColorBrewer_1.1-2
## [55] ellipsis_0.3.1 ica_1.0-2
## [57] pkgconfig_2.0.3 farver_2.1.0
## [59] sass_0.3.1 uwot_0.1.10
## [61] deldir_0.2-10 utf8_1.2.1
## [63] tidyselect_1.1.0 labeling_0.4.2
## [65] rlang_0.4.10 reshape2_1.4.4
## [67] later_1.2.0 munsell_0.5.0
## [69] tools_4.0.3 generics_0.1.0
## [71] ggridges_0.5.3 evaluate_0.14
## [73] stringr_1.4.0 fastmap_1.1.0
## [75] yaml_2.2.1 goftest_1.2-2
## [77] knitr_1.33 fitdistrplus_1.1-3
## [79] purrr_0.3.4 RANN_2.6.1
## [81] pbapply_1.4-3 nlme_3.1-152
## [83] mime_0.10 formatR_1.9
## [85] compiler_4.0.3 plotly_4.9.3
## [87] png_0.1-7 spatstat.utils_2.1-0
## [89] tibble_3.1.1 bslib_0.2.4
## [91] stringi_1.5.3 highr_0.9
## [93] lattice_0.20-41 Matrix_1.3-2
## [95] vctrs_0.3.7 pillar_1.6.0
## [97] lifecycle_1.0.0 spatstat.geom_2.1-0
## [99] lmtest_0.9-38 jquerylib_0.1.3
## [101] RcppAnnoy_0.0.18 data.table_1.14.0
## [103] bitops_1.0-7 irlba_2.3.3
## [105] httpuv_1.6.0 GenomicRanges_1.42.0
## [107] R6_2.5.0 promises_1.2.0.1
## [109] KernSmooth_2.23-18 gridExtra_2.3
## [111] IRanges_2.24.1 parallelly_1.24.0
## [113] codetools_0.2-18 MASS_7.3-53.1
## [115] assertthat_0.2.1 SummarizedExperiment_1.20.0
## [117] withr_2.4.2 sctransform_0.3.2
## [119] S4Vectors_0.28.1 GenomeInfoDbData_1.2.4
## [121] mgcv_1.8-35 parallel_4.0.3
## [123] grid_4.0.3 rpart_4.1-15
## [125] tidyr_1.1.3 rmarkdown_2.7
## [127] MatrixGenerics_1.2.1 Rtsne_0.15
## [129] Biobase_2.50.0 shiny_1.6.0