Informe práctica 3 - Integración de datasets de distintas especies con Análisis de Correlación Canónica (CCA)
Adam Casas
Compilado: 17 de abril del 2021
1 Introducción
Integrar datasets en Seurat V2 era tedioso, pues había que ejecutar varias líneas de código para obtener los resultados deseados. La 3ª versión de Seurat simplificó enormemente el protocolo con la aparición de los comandos FindIntegrationAnchors() y IntegrateData, pues ejecutan automáticamente la retahíla de pasos necesarios en Seurat V2.
Dicho esto, usaremos para este informe 2 datasets de scRNA-seq provenientes del trabajo de Butler et al. (2018). Los datasets en cuestión consisten en tejido pancreático de humano y de ratón, respectivamente. Aplicaremos el análisis de correlación canónica (CCA) para integrar sendos datasets.
library(future)
# Usar "multicore" en Linux y "multisession" en Windows
parallel_plan <- plan(strategy = "multisession", workers = 4)
# Ampliamos el límite de RAM de cada trabajador de R (500 MBs por defecto) a
# 1 GBs
options(future.globals.maxSize = 1 * 1024^3)
# Comprobamos el parámetro "estrategy" y el nº de hilos de nuestro plan actual
plan()## multisession:
## - args: function (..., envir = parent.frame(), workers = 4)
## - tweaked: TRUE
## - call: plan(strategy = "multisession", workers = 4)2 Carga y preprocesado de los datasets humano y murino
Primero cargamos las matrices de conteos de los datasets de interés:
set.seed(1)
human.data <- read.table("./Archivos accesorios/Codigo y datos de las Practicas 1-4/data/pancreas_human.expressionMatrix.txt", sep = "\t")
mouse.data <- read.table("./Archivos accesorios/Codigo y datos de las Practicas 1-4/data/pancreas_mouse.expressionMatrix.txt", sep = "\t")A continuación creamos dos objetos seurat a partir de dichas matrices de conteos y añadimos a sus metadatos la columna species, donde se recoge la información sobre la especie de la que provienen los datos. Nótese que las matrices de conteo de Butler et al. ya tienen filtradas las células mal secuenciadas, por lo que podemos omitir el paso del filtrado de células.
Por tanto el siguiente paso a realizar sería el normalizado, escalado de datos y detección de genes diferencialmente expresados, pero dado que el comando FindIntegrationAnchors integra las funciones de escalado y detección de genes altamente variables, no es necesario hacerlo ahora, aunque sería una buena práctica hacerlo (por tanto dejaremos las funciones ScaleData y FindVariableFeatures comentadas en el siguiente chunk).
Téngase en cuenta que de acuerdo al paper de Stuart et al., 2019, las células ancla o anchors son pares de células de distintos datasets que presentan perfiles de expresión génica similares y por tanto son vecinas mutuas o mutual nearest neighbours.
library(Seurat)
# Seurat V3.2 no es compatible con la ultima version de spatstat (2.0-1). Puedes
# instalar la versión previa con el siguiente comando:
#
# devtools::install_version(package = "spatstat", version = "1.64-1")
# Objeto seurat del dataset humano
human <- CreateSeuratObject(counts = human.data)
human@meta.data$species <- "human"
human <- NormalizeData(human)
# human <- ScaleData(human)
# human <- FindVariableFeatures(human, selection.method = "vst")
# Objeto seurat del dataset murino
mouse <- CreateSeuratObject(counts = mouse.data)
mouse@meta.data$species <- "mouse"
mouse <- NormalizeData(mouse)
# mouse <- ScaleData(mouse)
# mouse <- FindVariableFeatures(mouse, selection.method = "vst")3 Integración con CCA
Con los datasets ya cargados, procedemos a integrarlos en un único objeto seurat denominado pancreas.integrated. Para integrar los datasets, usaremos el ya mencionado comando FindIntegrationAnchors, el cual acepta como input una lista nombrada de los objetos seurat a integrar. Nótese que en el parámetro normalization.method hay que indicar el algoritmo empleado para el normalizado de los datos en el paso previo (este puede ser el método LogNormalize o el método SCT).
La función FindIntegrationAnchors realiza una reducción de la dimensionalidad y busca en dicha proyección las anchors de los datasets que se le pasen como input. Las técnicas de reducción de la dimensionalidad que incluye este comando son el CCA y el RPCA ( Reciprocal PCA). De acuerdo al laboratorio Satija, el CCA es apto para integrar datasets de distintas especies, pero no es computacionalmente eficiente para datasets grandes, por lo que si se quiere integrar datasets grandes y de una misma especie, se recomienda usar en su lugar el RPCA.
# Lista nombrada de los objetos seurat a integrar
pancreas.list <- list(human, mouse)
names(pancreas.list) <- c("human", "mouse")
# Buscamos las anchors de ambos datasets y realizamos el CCA con 20 componentes
# canónicas
pancreas.anchors <- FindIntegrationAnchors(object.list = pancreas.list, scale = T,
normalization.method = "LogNormalize",
reduction = "cca",
dims = 1:20)
# Integramos ambos datasets en uno nuevo
pancreas.integrated <- IntegrateData(anchorset = pancreas.anchors,
new.assay.name = "integrated",
normalization.method = "LogNormalize",
dims = 1:20)Con los datasets ya integrados en un único objeto seurat, basta con escalar los datos y ya podemos emplear dicho objeto en nuestros protocolos de visualizado y clustering.
# Cargamos librería gráfica
library(ggplot2)
# Escalamos el dataset integrado para emplearlo en protocolos posteriores
pancreas.integrated <- ScaleData(pancreas.integrated, verbose = T)
# Con el dataset ya escalado, podemos calcular y visualizar las células en las
# distintas proyecciones disponibles en Seurat
pancreas.integrated <- RunPCA(pancreas.integrated, npcs = 20, verbose = T)
pancreas.integrated <- RunTSNE(pancreas.integrated, reduction = "pca", dims = 1:20)
pancreas.integrated <- RunUMAP(pancreas.integrated, reduction = "pca", dims = 1:20)
# Gráficos
DimPlot(pancreas.integrated, reduction = "pca", group.by = "species") +
ggtitle("PCA (20 componentes canónicas)") + theme(plot.title = element_text(hjust = 0.5))
DimPlot(pancreas.integrated, reduction = "tsne", group.by = "species") +
ggtitle("t-SNE (20 componentes canónicas)") + theme(plot.title = element_text(hjust = 0.5))
DimPlot(pancreas.integrated, reduction = "umap", group.by = "species") +
ggtitle("UMAP (20 componentes canónicas)") + theme(plot.title = element_text(hjust = 0.5))
4 Clustering
Ahora encontramos clusters de células como haríamos en un protocolo cualquiera de scRNA-seq. Dado que tenemos muchas células en el dataset integrado, y por motivos de eficiencia computacional, vamos a usar el método igraph en el comando FindClusters en lugar del método por defecto matrix:
# Encontramos clusters
pancreas.integrated <- FindNeighbors(pancreas.integrated, dims = 1:20)
pancreas.integrated <- FindClusters(pancreas.integrated, resolution = 0.25, method = "igraph", algorithm = 1)## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 10455
## Number of edges: 390529
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9529
## Number of communities: 11
## Elapsed time: 0 seconds# Visualizamos los clusters
DimPlot(pancreas.integrated, reduction = "tsne", label = T, group.by = "seurat_clusters") +
ggtitle("t-SNE (20 componentes canónicas)") + theme(plot.title = element_text(hjust = 0.5))
DimPlot(pancreas.integrated, reduction = "umap", label = T, group.by = "seurat_clusters") +
ggtitle("UMAP (20 componentes canónicas)") + theme(plot.title = element_text(hjust = 0.5))
5 Ejercicios
Ejercicio 3.1: Las células β de los islotes de Langerhans producen insulina. ¿Puedes identificar dichas células en el dataset?
Solución:
Pedro Madrigal nos recomienda aquí buscar aquellas células que expresan el gen arquitecto PDX1, el cual codifica para un factor de transcripción de los genes de la insulina y la somatostatina (entre otros) y que participa en la organogénesis del páncreas.
No obstante, si en su lugar buscamos el propio gen de la insulina, INS, podremos mejorar la resolución del gráfico (pues es específico de las células \(\beta\) y se expresa en mayor cantidad que PDX1, pues los genes que codifican para factores de transcripción suelen tener niveles bajos de expresión).
# Nota sobre `FeaturePlot`:
# Seurat V2 usa los parámetros `reduction.use`, `features.plot` y `cols.use`;
# Seurat V3 usa los parámetros `reduction`, `features` y `cols`, respectivamente
# Usamos el factor de transcripción PDX1. Al ser un TF, usamos un "min.cutoff"
# pequeño:
FeaturePlot(pancreas.integrated, reduction ='umap', features = "PDX1",
min.cutoff = "q10", cols = c("lightgrey", "blue"), pt.size = 0.5)
# Usamos la insulina. Al tener mayor expresión que un TF, usamos un "min.cutoff"
# mayor:
FeaturePlot(pancreas.integrated, reduction ='umap', features = "INS1",
min.cutoff = "q60", cols = c("lightgrey", "blue"), pt.size = 0.5,
split.by = "species")
# Adicionalmente, podemos usar el factor de transcripción expecífico de células
# beta MAFA:
FeaturePlot(pancreas.integrated, reduction ='umap', features = c("MAFA"),
min.cutoff = "q5", cols = c("lightgrey", "blue"), pt.size = 0.5)
En vista de los resultados obtenidos, parece ser que las células de los clusters 1 y 2 (grupo central derecho en la proyección del UMAP) corresponden a las células pancreáticas \(\beta\). Podemos validar estas conclusiones con los metadatos proporcionados por el estudio original:
# Puede haber ligeras discrepancias entre las proyecciones t-SNE y UMAP del
# paper y las nuestros, de ahí que surjan NAs
pancreas.metadata <- read.table("./Archivos accesorios/Codigo y datos de las Practicas 1-4/data/Supplementary_Table_PancreasCellData.tsv", sep = "\t", skip = 2, header = T, row.names = 1)
pancreas.integrated <- AddMetaData(pancreas.integrated, metadata = pancreas.metadata)
DimPlot(pancreas.integrated, group.by = "Cluster_ID", label = T, reduction = "tsne")
DimPlot(pancreas.integrated, group.by = "Cluster_ID", label = T,
reduction = "umap", pt.size = 0.3)
Adicionalmente, podemos inspeccionar la tabla de metadatos en busca de la correspondencia de cada clúster con la verdad experimental:
# Seleccionamos las células que pertenecen a los clusters 1 y 2
celulas.clusters.1.2 <- which(pancreas.integrated@meta.data[,6] == 1 | pancreas.integrated@meta.data[,6] == 2)
head(pancreas.integrated@meta.data[celulas.clusters.1.2, c(6,9)], n = 20)En vista de la verdad empírica, confirmamos que hemos identificado en nuestro análisis las células β de los islotes de Langerhans como los clusters 1 y 2.
Ejercicio 3.2: ¿Puedes identificar genes marcadores específicos de las células pancreáticas β?
Solución:
Ahora que sabemos que los clusters 1 y 2 se corresponden con dichas células, podemos usar la función FindMarkers para detectar sus biomarcadores específicos.
cluster.beta.markers.integrated <- FindMarkers(pancreas.integrated, ident.1 = c(1,2),
min.pct = 0.25, only.pos = T, assay = "integrated")
head(cluster.beta.markers.integrated, n = 10)Nótese que no todos los genes de interés se encuentran en el ensayo integrated, pues tal y como vimos en el ejercicio anterior, los genes INS1 y MAFA se localizaban en el ensayo RNA. En el siguiente chunk comprobamos que, efectivamente, dichos biomarcadores están ausentes en el ensayo integrated:
which(rownames(cluster.beta.markers.integrated) %in% c("INS1", "INS2", "MAFA"))## integer(0)Podemos solucionar esto cambiando el parámetro assay del comando en cuestión:
cluster.beta.markers.RNA <- FindMarkers(pancreas.integrated, ident.1 = c(1,2),
min.pct = 0.25, only.pos = T, assay = "RNA")
head(cluster.beta.markers.RNA, n = 10)Ahora debería detectar correctamente dichos biomarcadores:
# Confirmamos que detectamos genes diferencialmente expresados en células β
which(rownames(cluster.beta.markers.RNA) %in% c("INS1", "INS2", "MAFA"))## [1] 1 2 11cluster.beta.markers.RNA[c("INS1", "INS2", "MAFA"),]Ejercicio 3.3: ¿Puedes identificar genes con expresión diferencial entre las células pancreáticas α y β?
Solución:
Recordemos que las células α expresan el glucagón, mientras que las células β producen insulina, por lo que esperamos ver en la siguiente lista los genes INS1 y GCG. Podemos identificar las células α de la misma manera que hicimos en el ejercicio anterior, pero ahora que hemos cargado los metadatos del estudio original, podemos averiguar de un vistazo que el clúster 0 se corresponde con las células α.
FeaturePlot(pancreas.integrated, reduction ='umap',
features = "GCG",
min.cutoff = "q9", cols = c("lightgrey", "blue"),
pt.size = 0.5)
# "only.pos" devuelve sólo los genes sobreexpresados en "ident.1" respecto los
# niveles de expresión en "ident.2"
biomarcadores.diferenciales <- FindMarkers(pancreas.integrated, ident.1 = c(1,2), ident.2 = 0,
min.pct = 0.25, only.pos = F, assay = "RNA")
# Confirmamos que detectamos biomarcadores correspondientes con células α y β
biomarcadores.diferenciales[c("INS1", "GCG"), ]# Visualizamos
head(biomarcadores.diferenciales, n = 10)# REG1A se expresa en células acinares (páncreas exocrino)
# PPY se produce en las células PP de los islotes (anotadas en el paper como "gamma") y genera sensación de estar
# lleno (inhibe ingestión de alimentos)
# SST somatostatina, producida por células delta. Inhibe producción de glucagón
# e insulina
# GHRL codifica el péptido precursor o pre-proteína grelina-obestatina (se detecta en unas pocas células PP o gamma)
# VWF o Factor de Von Willebrand, glicoproteína transportadora del Factor VIII
# plasmático (hemostasis). Se expresa en casi todo el cuerpo, aunque aquí se detecta en las células endoteliales
# SOX10, detectado en algunas células "stellate" (¿neuronas o pancreatic stellate cell?). This protein acts as a nucleocytoplasmic shuttle protein and is important for neural crest and peripheral nervous system development.
FeaturePlot(pancreas.integrated, reduction ='umap',
features = c("REG1A", "PPY", "SST", "GHRL", "VWF", "SOX10"),
min.cutoff = "q9", cols = c("lightgrey", "blue"),
pt.size = 0.5,
combine = F)## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
6 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] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggplot2_3.3.3 Seurat_3.2.0 future_1.21.0
##
## loaded via a namespace (and not attached):
## [1] nlme_3.1-152 matrixStats_0.58.0 RcppAnnoy_0.0.18
## [4] RColorBrewer_1.1-2 httr_1.4.2 sctransform_0.3.2
## [7] tools_3.6.3 bslib_0.2.4 utf8_1.2.1
## [10] R6_2.5.0 irlba_2.3.3 rpart_4.1-15
## [13] KernSmooth_2.23-18 uwot_0.1.10 mgcv_1.8-34
## [16] lazyeval_0.2.2 colorspace_2.0-0 withr_2.4.1
## [19] gridExtra_2.3 tidyselect_1.1.0 compiler_3.6.3
## [22] plotly_4.9.3 labeling_0.4.2 sass_0.3.1
## [25] scales_1.1.1 spatstat.data_2.1-0 lmtest_0.9-38
## [28] ggridges_0.5.3 pbapply_1.4-3 goftest_1.2-2
## [31] spatstat_1.64-1 stringr_1.4.0 digest_0.6.27
## [34] spatstat.utils_2.1-0 rmarkdown_2.7 pkgconfig_2.0.3
## [37] htmltools_0.5.1.1 parallelly_1.24.0 limma_3.42.2
## [40] highr_0.8 fastmap_1.1.0 htmlwidgets_1.5.3
## [43] rlang_0.4.10 shiny_1.6.0 farver_2.1.0
## [46] jquerylib_0.1.3 generics_0.1.0 zoo_1.8-9
## [49] jsonlite_1.7.2 ica_1.0-2 dplyr_1.0.5
## [52] magrittr_2.0.1 patchwork_1.1.1 Matrix_1.3-2
## [55] Rcpp_1.0.6 munsell_0.5.0 fansi_0.4.2
## [58] abind_1.4-5 ape_5.4-1 reticulate_1.18
## [61] lifecycle_1.0.0 stringi_1.5.3 yaml_2.2.1
## [64] MASS_7.3-53.1 Rtsne_0.15 plyr_1.8.6
## [67] grid_3.6.3 parallel_3.6.3 listenv_0.8.0
## [70] promises_1.2.0.1 ggrepel_0.9.1 crayon_1.4.1
## [73] deldir_0.2-10 miniUI_0.1.1.1 lattice_0.20-41
## [76] cowplot_1.1.1 splines_3.6.3 tensor_1.5
## [79] knitr_1.31 pillar_1.5.1 igraph_1.2.6
## [82] future.apply_1.7.0 reshape2_1.4.4 codetools_0.2-18
## [85] leiden_0.3.7 glue_1.4.2 evaluate_0.14
## [88] data.table_1.14.0 vctrs_0.3.7 png_0.1-7
## [91] httpuv_1.5.5 polyclip_1.10-0 gtable_0.3.0
## [94] RANN_2.6.1 purrr_0.3.4 tidyr_1.1.3
## [97] xfun_0.22 rsvd_1.0.3 mime_0.10
## [100] xtable_1.8-4 RSpectra_0.16-0 later_1.1.0.1
## [103] survival_3.2-10 viridisLite_0.3.0 tibble_3.1.0
## [106] cluster_2.1.1 globals_0.14.0 fitdistrplus_1.1-3
## [109] ellipsis_0.3.1 ROCR_1.0-11