Informe práctica 1 - Análisis de datos scRNA-seq en Seurat V3
Adam Casas
Compilado: 27 de abril del 2021
1 Introducción y pre-requisitos
En este documento se muestran las prácticas y ejercicios de Análisis de datos scRNA-seq, traducidas de Seurat 2 a Seurat 3, ya que las prácticas usaban la versión 2 de Seurat, la cual está actualmente obsoleta y tiene dependencias obsoletas.
Las prácticas 1 y 2 se basan en el código modificado de las prácticas originales.
En la práctica 3 se integran en Seurat datos de scRNA-seq de páncreas humano y de ratón ( mus musculus) siguiendo la metodología de Butler et al., 2018.
En la práctica 4 aprenderemos a inferir trayectorias celulares a partir de datos de expresión single cell usando Monocle3.
Los archivos necesarios para realizar las prácticas se pueden encontrar en Dropbox o en GitHub.
# Instalación Seurat V3 en R 3.6
remotes::install_version("Seurat", version = "3.2")
# 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")Comenzamos cargando Seurat
# Cargamos librerías
library("Seurat")
library("ggplot2")
library("dplyr")Hemos cargado la siguiente versión de Seurat:
packageVersion("Seurat")## [1] '3.2.0'1.1 Matriz de conteos de UMIs
En la práctica 1 de Seurat usaremos datos de SC3-seq (un tipo de scRNA-seq con amplificación del extremo 3 prima del ARNm de la célula) de embriones del macaco cangrejero Macaca fascicularis (Nakamura et al., 2016). Los datos se encuentran en el archivo D3Ecounts.txt (tiene valores separados por tabuladores) y vienen en forma de matriz de conteos de UMIs en las que las filas son los genes y las columnas son las células. El dataset consta de 399 células y les aplicaremos reducción de la dimensionalidad, clustering e identificación de genes marcadores.
La tabla 2.1 muestra detalles sobre los tipos celulares encontrados en el dataset.
1.2 Preprocesado del dataset (ya realizado)
El dataset empleado se descargó del GEO en formato FastQ y contenía originalmente 474 células. Las células de tipo cyESC, ePGC e IPGC fueron filtradas y eliminadas del dataset, lo que generó el tamaño actual de 399 células. Tal como se hizo en el artículo de Nakamura et al., las lecturas de ARNm fueron preprocesadas con cutadapt 1.16 para eliminar el adaptador, la cola poli-A y los nucleótidos mal secuenciados (baja calidad) con cutadapt -colorspace -e 0.1 -q 20 -n 2 -O 1 -m 30 -a CTCGAGGGCGCGCCGGATCCATATACGCCTTGGCCGTACAGCAG -a “A{100}” -g CTCGAGGGCGCGCCGGATCCATATACGCCTTGGCCGTACAGCAG.
Posteriormente las lecturas se mapearon al genoma de Macaca_fascicularis_5.0 descargado de ftp://ftp.ncbi.nih.gov/genomes/Macaca_fascicularis/Assembled_chromosomes/seq/ (chr1-20, chrX y chrMT) y a las secuencias ERCC spike-in RNA obtenidas de https://tools.thermofisher.com/content/sfs/manuals/cms_095047.txt. Para el mapeado, se empleó el programa Bowtie V1.2.2 con los parámetros bowtie -q -C -S -k 1 -p 1.
Tras el mapeado, analizamos aquellas lecturas alineadas que tenían una calidad de mapeado o Mapping Quality >= 10. Para la cuantificación de genes se usó el programa featureCounts v1.6.0 con los parámetros -T 2 -t exon -g gene y el archivo de anotación génica macFas5_Annot_wo_pseudo_v1e_add10k.gtf (disponible gracias a Tomonori Nakamura y Mitinori Saitou).
2 Comenzamos el workflow
Primero leemos la matriz de conteos de UMIs, almacenada en el archivo D3Ecounts.txt, y observamos los 10 primeros genes y células.
# Leemos matriz de conteos UMIs (contiene valores separados por tabuladores)
matriz_umis_nakamura <- read.table("./Cosas accesorias Informe 1/D3Ecounts.txt", sep = "\t", header = T)
head(matriz_umis_nakamura[1:10, 1:10])Ahora formateamos la matriz para que quede más bonita y ordenada. Dado que cada fila corresponde a un gen, vamos a renombrar las filas con los nombres de los susodichos genes:
rownames(matriz_umis_nakamura) <- matriz_umis_nakamura$Geneid
matriz_umis_nakamura <- matriz_umis_nakamura[,2:ncol(matriz_umis_nakamura)]
head(matriz_umis_nakamura[1:10, 1:10])3 El objeto Seurat
El paquete Seurat no contiene una de las clases más usadas en Bioconductor, SingleCellExperiment. En su lugar usa su propio tipo de objeto, seurat. Todos los cálculos realizados en esta práctica se ejecutan sobre este objeto. Para comenzar el análisis primero tenemos que crear el objeto de tipo seurat a partir de los datos SC3-seq crudos ( i.e. no normalizados). Conservaremos todos los genes expresados en >= 3 células (parámetro min.cells) y todas las células que expresen como mínimo 200 genes (parámetro min.features). Esto se hace para evitar analizar células mal secuenciadas.
nakamura_seurat <- CreateSeuratObject(counts = matriz_umis_nakamura,
project = "Nakamura",
min.cells = 3,
min.features = 200)4 Control de calidad
Seurat te permite visualizar con simpleza métricas QC y filtrar células en base a cualquier criterio que imponga el usuario. Podemos visualizar conteos de moléculas de ARNm y genes y graficar su relación.
# Seurat 2 usa los nombres "nGene" y "nUMI";
# Seurat 3 usa los nombres "nFeature_RNA" y "nCount_RNA", respectivamente
VlnPlot(object = nakamura_seurat,
features = c("nFeature_RNA", "nCount_RNA"),
ncol = 2)
De los gráficos de violín se observa que la mayoría de células expresan entre 11000 y 14000 genes, y la mayoría de células presentan una expresión génica de entre 1 y 3 millones de moléculas de ARNm.
Los UMIs (Unique Molecular Identifiers) son los identificadores de los primers empleados en las tecnologías de scRNA-seq basadas en gotículas. Se usan para identificar lecturas originadas de una misma molécula de ARNm, de manera que evitamos contar dos veces la misma molécula de ARNm (ver imagen inferior, sacada de Bach et al., 2017).
# En Seurat 2, la función se llama "GenePlot()";
# mientras que en Seurat 3 se renombró a "FeatureScatter()"
FeatureScatter(object = nakamura_seurat,
feature1 = "nFeature_RNA",
feature2 = "nCount_RNA")
En el gráfico anterior se observa que hay una relación (aunque no muy buena, dado que R = 0.44) entre el nº de genes que expresa una célula y el nº de transcritos de ARN que contiene la susodicha, lo cual es lógico.
Podemos excluir del análisis aquellas células que presenten un valor anormalmente elevado de transcritos con el comando subset():
# En Seurat 2 la función se llama "FilterCells()";
# mientras que en Seurat 3 se renombró a "subset()"
nakamura_seurat <- subset(x = nakamura_seurat, subset = nCount_RNA < 1e+07)
max(nakamura_seurat@meta.data$nCount_RNA)## [1] 42999815 Normalización
Tras eliminar células indeseadas del dataset, el siguiente paso es normalizar los datos. Por defecto usamos LogNormalize, método ya explicado en el tutorial de 2700 PBMCs.
Nótese que siempre podemos utilizar el método sctransform en lugar de LogNormalize para obtener una mejor resolución de los datos biológicos en cuestión (véase el tutorial pertinente del laboratorio Satija)
nakamura_seurat <- NormalizeData(nakamura_seurat,
normalization.method = "LogNormalize",
scale.factor = 10000)6 Genes altamente variables (=diferencialmente expresados)
Seurat calcula los genes altamente variables y centra los posteriores análisis en ellos. Seurat 2 usaba el método mvp (mean.var.plot) del comando FindVariableGenes(), el cual calculaba la expresión media y la dispersión (desviación estándar) de cada gen, luego colocaba los genes en grupos (bins) y finalmente calculaba el z-score para la dispersión dentro de cada grupo. Esto permitía controlar la fuerte relación entre la varianza y la expresión media de cada gen. No obstante, el método por defecto en Seurat 3 es vst.
# En Seurat 2 la función se llama "FindVariableGenes()";
# mientras que en Seurat 3 se renombró a "FindVariableFeatures()"
nakamura_seurat <- FindVariableFeatures(nakamura_seurat,
selection.method = "mvp")
VariableFeaturePlot(nakamura_seurat)
# Identificamos los 10 genes más variables
top10 <- head(VariableFeatures(nakamura_seurat), n = 10L)
# Graficamos dichos genes
grafico_genes_variable <- VariableFeaturePlot(nakamura_seurat)
LabelPoints(grafico_genes_variable, points = top10, repel = TRUE)
7 Escalado de datos
Después de la normalización de los datos, procede escalarlos. 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). En este paso podemos controlar algunas variables confusoras o confounders mediante el parámetro vars.to.regress (no se hace en esta práctica 1).
nakamura_seurat <- ScaleData(nakamura_seurat)8 Reducción lineal de la dimensionalidad (PCA)
A continuación realizamos el PCA en los datos escalados. Por defecto el comando RunPCA() actúa sobre los genes variables detectados previamente (almacenados en nakamura_seurat[["RNA"]]@meta.features), aunque podemos configurar los genes sobre los que realizar el PCA mediante el parámetro features (los genes seleccionados necesitan estar escalados y tener varianza distinta de 0).
La razón por la que las técnicas de reducción de la dimensionalidad en Seurat usan como input sólo los genes variables es que con algunos tipos de datos (como las matrices de conteos de UMIs), y sobre todo si se controlan variables técnicas, dichas técnicas devuelven resultados similares tanto si se ejecutan sólo sobre los genes variables como si se ejecutan en muchos más genes o incluso en el transcriptoma completo (lo cual es más lento de computar).
# Seurat 2 usa los parámetros "pcs.print", "genes.print" y "pcs.compute";
# mientras que en Seurat 3 se renombraron a "ndims.print", "nfeatures.print" y
# "npcs", respectivamente.
nakamura_seurat <- RunPCA(object = nakamura_seurat,
ndims.print = 1:5,
nfeatures.print = 5,
npcs = 20)## PC_ 1
## Positive: ANXA2, BASP1, NLRP7, HAVCR1, DPP4
## Negative: NME4, SNRPN, MARCKSL1, RAB34, PSIP1
## PC_ 2
## Positive: LOC102131610, COL4A2, COL6A1, PITX2, EGFLAM
## Negative: IFITM1, POU5F1, NANOG, DNMT3B, TDGF1
## PC_ 3
## Positive: AKAP12, GPX2, HNF1B, IGF1, LOC102120281
## Negative: TSPAN15, GATA2, RAB25, GRHL1, CTSD
## PC_ 4
## Positive: LOC101867055, S100A16, NR6A1, CD59, S100A13
## Negative: STOM, MT2A, DPRX, DPPA2, CABLES1
## PC_ 5
## Positive: PYGL, GPR56, GJA5, ECM1, MRGPRX1
## Negative: RALGDS, ACSL6, SLCO2B1, MTUS2, SLC47A1Seurat ofrece varias maneras de visualizar las células y los genes que definen las componentes principales encontradas:
# Células
head(nakamura_seurat@reductions$pca@cell.embeddings[,1:10])## PC_1 PC_2 PC_3 PC_4 PC_5 PC_6 PC_7
## EXMC -11.309357 35.86274 -2.825483 -15.74335 -4.642842 3.949886 -0.9967165
## EXMC.1 -8.234667 24.52330 -1.134669 -11.11603 -3.015602 3.265595 -1.0060007
## EXMC.2 -10.545690 33.11574 -2.031776 -15.54442 -5.653357 4.811282 -0.8679818
## EXMC.3 -8.888462 30.56810 -2.654018 -14.50364 -3.798891 4.155416 0.2685116
## EXMC.4 -9.432198 32.25424 -2.674476 -14.85510 -3.907804 4.562688 0.6795201
## EXMC.5 -10.183890 35.52257 -3.026999 -15.54798 -3.810246 5.239302 -0.7745654
## PC_8 PC_9 PC_10
## EXMC -1.571345 4.166363 3.066461
## EXMC.1 -2.585341 5.840654 7.735386
## EXMC.2 -1.978976 5.505307 3.062014
## EXMC.3 -1.590057 4.485082 5.467402
## EXMC.4 -2.088337 4.320722 5.283886
## EXMC.5 -2.148220 5.215708 4.916689# Genes
head(nakamura_seurat@reductions$pca@feature.loadings[,1:10])## PC_1 PC_2 PC_3 PC_4 PC_5
## CGA -0.006349136 0.01570605 -0.01923307 0.010420608 0.066699043
## COL3A1 -0.016023672 0.06161703 -0.01257933 -0.046679852 -0.015084367
## DEFB1 0.009147953 0.01510659 -0.01575525 0.023197315 0.065322521
## FN1 -0.006931077 0.05915392 0.04344355 0.009592335 -0.005997288
## COL1A1 -0.031923608 0.04676786 -0.01560889 -0.036962789 -0.013478514
## NID2 0.018207297 0.04342587 0.06068622 0.015042142 -0.011123386
## PC_6 PC_7 PC_8 PC_9 PC_10
## CGA -0.004501473 0.026770424 0.007568193 -0.021481023 -0.03392886
## COL3A1 0.017411778 0.005091450 -0.002892462 -0.029324464 -0.03127701
## DEFB1 0.015490525 0.014030564 0.007763182 -0.001943517 -0.02762655
## FN1 0.009496249 -0.025080270 -0.004245887 0.013319604 -0.05923133
## COL1A1 0.029597770 0.008912462 0.004546361 -0.017396086 -0.04695004
## NID2 0.044726697 0.010612259 0.022470025 -0.015478820 -0.01808552# Seurat 2 usa la función "VizPCA()";
# mientras que en Seurat 3 se renombró a "VizDimLoadings()"
VizDimLoadings(nakamura_seurat, dims = 1:2, nfeatures = 30, reduction = "pca")
# Nota: También puedes usar "PCAPlot()" tal que así:
# PCAPlot(nakamura_seurat)
DimPlot(nakamura_seurat, reduction = "pca", dims = c(1,2))
Concretamente, los comandos PCHeatmap() y DimHeatmap() (realizan exactamente la misma función) permiten explorar fácilmente las fuentes de varianza en un dataset, y pueden ser útiles a la hora de decidir la verdadera dimensionalidad del dataset, o sea, decidir qué componentes principales conservar para los próximos análisis. Tanto las células como los genes están ordenados en función de sus loadings del PCA. El parámetro cells permite visualizar un número reducido de células (útil para graficar datasets grandes), e incluye las células con los loadings de ambos extremos (i.e. incluye en el heatmap tanto las células que más infraexpresan un gen concreto como aquellas que más lo sobreexpresan). Colores violáceos indican infraexpresión, mientras que colores amarillos indican sobreexpresión.
# PCHeatmap(object = nakamura_seurat,
# dims = 1:6,
# cells = 500,
# balanced = T,
# fast = F)
DimHeatmap(object = nakamura_seurat,
dims = 1:6,
cells = 500,
balanced = T,
fast = F)
Nótese que el parámetro fast (TRUE por defecto) permite graficar rápidamente los heatmaps, pero elimina la leyenda y previene que podamos modificar los gráficos resultantes (por ejemplo, no podremos cambiar el título de los heatmaps).
9 Componentes principales significativas (determinar la verdadera dimensionalidad del dataset)
Para superar el ruido técnico presente en cualquier dataset de scRNA-seq, Seurat agrupa las células en función de sus loadings del PCA, siendo cada componente principal un “metagén” que combina información sobre un set de genes correlacionados. Por tanto, determinar el nº correcto de componentes principales que conservar para futuros análisis es un paso muy importante, y se puede averiguar con varios métodos, entre ellos, el método del Jackstraw.
El comando JackStraw() permuta aleatoriamente un subconjunto del dataset original (1% de células por defecto) y ejecuta de nuevo el PCA y repite este proceso para generar una distribución nula de puntuaciones de genes. Las componentes principales significativas serán aquellas que estén enriquecidas con genes que presenten un p-valor pequeño.
# Nota: El Jackstraw puede ser lento en datasets medianos y grandes. En esos
# casos se puede optar por usar el método del codo, mostrado a continuación
nakamura_seurat <- JackStraw(nakamura_seurat, num.replicate = 100, dims = 9)
nakamura_seurat <- ScoreJackStraw(nakamura_seurat, dims = 1:9, do.plot = F)Las funciones ScoreJackStraw (parámetro do.plot) y JackStrawPlot permiten visualizar qué componentes principales son estadísticamente significativas. Serán significativas aquellas componentes que se sitúen por encima de la línea discontinua (distribución uniforme). En el gráfico a continuación se observa que podríamos conservar las primeras 9 componentes principales.
JackStrawPlot(nakamura_seurat, dims = 1:9)
El método del codo, ya explicado en bioestadística y en otros informes, se puede realizar con el siguiente comando:
# Seurat 2 usa la función "PCElbowPlot()",
# mientras que en Seurat 3 se renombró a "ElbowPlot()"
ElbowPlot(nakamura_seurat)
10 Clustering basado en grafos
El método de clustering que incorpora Seurat se basa en grafos y ya viene explicado en otro informe.
La función FindClusters() contiene un parámetro de resolución que establece la granularidad del clustering. A mayor valor de resolución/granularidad, más clusters se generan. Tanto 10x Genomics como el Dr. Pedro Madrigal recomiendan establecer una resolución de entre 0,4 ~ 1,2 para datasets de 3000 células. Para datasets mayores, la resolución debe aumentar proporcionalmente. Los clusters con sus respectivas células se guardan en nakamura_seurat@active.ident, nakamura_seurat$RNA_snn_res.0.6 y nakamura_seurat$seurat_clusters (nakamura_seurat@ident en Seurat V2), y se puede acceder a ellos con el comando Idents(nakamura_seurat).
Nótese que una vez realizado el clustering, las visualizaciones en PCA, tSNE y UMAP pasarán a incluir información sobre a qué clúster pertenece cada célula representada en las mismas.
# Seurat 2 prescinde del comando "FindNeighbors()" y en su lugar agrupa todos
# sus parámetros en la función "FindCluster()"
nakamura_seurat <- FindNeighbors(object = nakamura_seurat,
reduction = "pca",
dims = 1:9,
k.param = 11)## Computing nearest neighbor graph## Computing SNNnakamura_seurat <- FindClusters(object = nakamura_seurat,
resolution = 0.6)## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 399
## Number of edges: 4692
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9194
## Number of communities: 11
## Elapsed time: 0 seconds# 20 primeras células y sus clusters
head(Idents(nakamura_seurat), n = 20L)## EXMC EXMC.1 EXMC.2 EXMC.3 EXMC.4 EXMC.5 EXMC.6 EXMC.7 EXMC.8 EXMC.9
## 6 6 6 6 6 6 6 6 6 6
## EXMC.10 EXMC.11 EXMC.12 EXMC.13 EXMC.14 EXMC.15 EXMC.16 EXMC.17 EXMC.18 EXMC.19
## 6 6 6 6 6 6 6 6 6 6
## Levels: 0 1 2 3 4 5 6 7 8 9 10# Nº de células por clúster
#
# table(nakamura_seurat$seurat_clusters)
# table(nakamura_seurat$RNA_snn_res.0.6)
table(nakamura_seurat@active.ident)##
## 0 1 2 3 4 5 6 7 8 9 10
## 55 54 52 49 44 34 33 32 21 14 11Una utilidad muy práctica de Seurat es la posibilidad de obtener los parámetros con los que configuramos los comandos previamente ejecutados sobre el objeto seurat nakamura_seurat. Con el comando Command(nakamura_seurat) podemos ver una lista de dichas funciones, y con la estructura nakamura_seurat$comando_usado podemos obtener los parámetros con los que ejecutamos previamente esa función:
# Comandos ejecutados
Command(nakamura_seurat)## [1] "NormalizeData.RNA" "FindVariableFeatures.RNA"
## [3] "ScaleData.RNA" "RunPCA.RNA"
## [5] "JackStraw.RNA.pca" "ScoreJackStraw"
## [7] "FindNeighbors.RNA.pca" "FindClusters"# Parámetros del comando "FindClusters()"
nakamura_seurat$FindClusters## Command: FindClusters(object = nakamura_seurat, resolution = 0.6)
## Time: 2021-04-27 15:27:28
## graph.name : RNA_snn
## modularity.fxn : 1
## resolution : 0.6
## method : matrix
## algorithm : 1
## n.start : 10
## n.iter : 10
## random.seed : 0
## group.singletons : TRUE
## verbose : TRUESeurat permite usar también gráficos tSNE y UMAP para visualizar el clustering. tSNE usa como input las dimensiones significativas del PCA (por ejemplo, las 9 que detallábamos en el paso del JackStraw), aunque también podemos pedirle que use como input un set de genes escalados concretos mediante el argumento features (genes.use en Seurat V2).
nakamura_seurat <- RunTSNE(nakamura_seurat, reduction = "pca", dims = 1:9, tsne.method = "Rtsne", verbose = F)
# Nota: También puedes usar "TSNEPlot()", tal que así:
# TSNEPlot(nakamura_seurat)
DimPlot(nakamura_seurat, reduction = "tsne", dims = c(1,2))
Ídem para UMAP:
nakamura_seurat <- RunUMAP(nakamura_seurat, reduction = "pca",
dims = 1:9, umap.method = "uwot",
metric = "cosine", min.dist = 0.3)## 15:27:30 UMAP embedding parameters a = 0.9922 b = 1.112## 15:27:30 Read 399 rows and found 9 numeric columns## 15:27:30 Using Annoy for neighbor search, n_neighbors = 30## 15:27:30 Building Annoy index with metric = cosine, n_trees = 50## 0% 10 20 30 40 50 60 70 80 90 100%## [----|----|----|----|----|----|----|----|----|----|## **************************************************|
## 15:27:30 Writing NN index file to temp file C:\Users\Adam_\AppData\Local\Temp\Rtmp8IQjZH\file2a285ccfa5b
## 15:27:30 Searching Annoy index using 1 thread, search_k = 3000
## 15:27:30 Annoy recall = 100%
## 15:27:30 Commencing smooth kNN distance calibration using 1 thread
## 15:27:30 Initializing from normalized Laplacian + noise
## 15:27:30 Commencing optimization for 500 epochs, with 13218 positive edges
## 15:27:31 Optimization finished# Nota: También puedes usar "UMAPPlot()", tal que así:
# UMAPPlot(nakamura_seurat)
DimPlot(nakamura_seurat, reduction = "umap", dims = c(1,2), label = T)
11 Genes marcadores (biomarcadores de clusters/genes diferencialmente expresados)
Seurat te permite encontrar marcadores que definen clusters mediante expresión diferencial. Por defecto, identifica marcadores positivos y negativos de un único clúster, comparado con el resto de células. Puedes testar clusters entre sí o un grupo de clusters vs el resto de células. Por ejemplo, para encontrar los biomarcadores del clúster 2 ejecutamos este código:
cluster2.markers <- FindMarkers(nakamura_seurat,
ident.1 = 2,
min.pct = 0.25)
head(cluster2.markers, n = 5)Para ver los genes diferencialmente expresados entre el clúster 2 y el 4 ( o sea, comparar los patrones de expresión génica del clúster 2 vs los patrones del clúster 4):
cluster2_4.markers <- FindMarkers(nakamura_seurat, ident.1 = 2,
ident.2 = 4, min.pct = 0.25)
head(cluster2_4.markers, n = 5)Para ver los genes diferencialmente expresados del clúster 1 y 5 vs el resto del dataset:
cluster1_5.markers <- FindMarkers(nakamura_seurat, ident.1 = c(1,5),
min.pct = 0.25)
head(cluster1_5.markers, n = 5)Podemos visualizar los biomarcadores de los clusters con:
VlnPlot(nakamura_seurat, features = "BMP6") +
theme(axis.text.x = element_text(angle = 0, hjust = .5))
FeaturePlot(nakamura_seurat, features = "IL6R", reduction = "umap",
pt.size = 2)
El comando FindAllMarkers() automatiza este proceso y encuentra los biomarcadores de todos los clusters (seleccionamos solo los genes sobreexpresados de cada cluster con el parámetro only.pos = T):
# Nota: En Seurat V2 el parámetro "logfc.threshold" se llama "thresh.use"
markers <- FindAllMarkers(nakamura_seurat, min.pct = 0.25,
logfc.threshold = 0.25, only.pos = T)## Calculating cluster 0## Calculating cluster 1## Calculating cluster 2## Calculating cluster 3## Calculating cluster 4## Calculating cluster 5## Calculating cluster 6## Calculating cluster 7## Calculating cluster 8## Calculating cluster 9## Calculating cluster 10head(markers)El comando DoHeatmap() genera un heatmap de expresión génica para las células y genes seleccionados. En este caso graficamos los top 5 biomarcadores de cada clúster:
top5 <- markers %>% group_by(cluster) %>% slice_max(order_by = avg_logFC, n = 5)
DoHeatmap(nakamura_seurat, features = top5$gene)
Ahora podemos asignar las etiquetas dadas a las células (disponibles en el archivo accesorio colnames.txt):
archivo_etiquetas_celulas <- t(read.table("./Cosas accesorias Informe 1/colnames.txt"))
# Descartamos el primer elemento del archivo dado que no es un nombre de célula
NewIdent <- as.factor(archivo_etiquetas_celulas[2:length(archivo_etiquetas_celulas)])
names(NewIdent) <- names(nakamura_seurat@active.ident)
nakamura_seurat@active.ident <- NewIdent
nakamura_seurat <- AddMetaData(nakamura_seurat, metadata = NewIdent, col.name = "Cell.Identity")
head(nakamura_seurat@meta.data)# UMAP con etiquetas
DimPlot(nakamura_seurat , group.by='Cell.Identity', reduction = "umap", pt.size=2, label = T) +
ggtitle("UMAP") + theme(plot.title = element_text(hjust = 0.5))
# t-SNE con etiquetas
DimPlot(nakamura_seurat , group.by='Cell.Identity', reduction = "tsne", pt.size=2, label = T) +
ggtitle("t-SNE") + theme(plot.title = element_text(hjust = 0.5))
12 Inferencia del ciclo celular
Para la inferencia del ciclo celular, asiganmos a cada célula una puntuación basada en sus niveles de expresión de genes marcadores de las fases S y G2/M. Los niveles de expresión de los biomarcadores de una fase celular deberían estar inversamente correlacionados con los de la otra fase (están en extremos opuestos del ciclo celular), y las células que no expresen ninguno de estos biomarcadores se encontrarán probablemente en la fase G1 (no se están duplicando).
En Seurat, asiganmos esta puntuación a las células con la función CellCycleScoring(), la cual almacena las puntuación de las fases S y G2/M en nakamura_seurat@meta.data, junto con la fase celular predicha de cada célula. Si estableces el parámetro set.ident = TRUE, la identidad de las células pasará a ser su fase celular y la identidad previa se guardará en la columna old.ident.
Los biomarcadores de estos puntos de control en Macaca fascicularis se encuentran en el archivo de texto regev_lab_cell_cycle_genes.txt. Los primeros 43 genes corresponden a los biomarcadores de la fase S, y los 54 genes restantes corresponden a los de las fases G2/M.
# Leemos el archivo de texto con los biomarcadores y separamos los genes de cada
# fase
cellcycle.genes <- readLines(con = "./Cosas accesorias Informe 1/regev_lab_cell_cycle_genes.txt")
s.genes <- cellcycle.genes[1:43]
g2.m.genes <- cellcycle.genes[44:97]
# Predecimos la fase celular de nuestras células
nakamura_seurat <- CellCycleScoring(nakamura_seurat, s.features = s.genes, g2m.features = g2.m.genes, set.ident = F)
DimPlot(nakamura_seurat, reduction = "umap", group.by = "Phase", pt.size = 2) +
ggtitle("UMAP") + theme(plot.title = element_text(hjust= 0.5))
Mediante un ridgeplot podemos comprobar visualmente la expresión diferencial de los marcadores PCNA, MCM6, TOP2A y MKI67 en función de la fase en la que se haye la célula:
# PCNA y MCM6 son marcadores de la fase S;
# TOP2A y MKI67 son marcadores de la fase G2/M
RidgePlot(nakamura_seurat, features = c("PCNA", "MCM6", "TOP2A", "MKI67"), group.by = "Phase", ncol = 2)## Picking joint bandwidth of 0.129## Picking joint bandwidth of 0.126## Picking joint bandwidth of 0.14## Picking joint bandwidth of 0.114
Por último procedemos a guardar el trabajo realizado hasta ahora, pues lo necesitaremos en la práctica 4
saveRDS(object = nakamura_seurat, file = "./Cosas accesorias Informe 1/nakamura_P4.rds")13 Bibliografía
Butler, A., Hoffman, P., Smibert, P. et al. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat Biotechnol 36, 411–420 (2018). https://doi.org/10.1038/nbt.4096
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
14 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] dplyr_1.0.5 ggplot2_3.3.3 Seurat_3.2.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-35
## [16] lazyeval_0.2.2 colorspace_2.0-0 withr_2.4.2
## [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.9 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.4 generics_0.1.0 zoo_1.8-9
## [49] jsonlite_1.7.2 ica_1.0-2 magrittr_2.0.1
## [52] patchwork_1.1.1 Matrix_1.3-2 Rcpp_1.0.6
## [55] munsell_0.5.0 fansi_0.4.2 abind_1.4-5
## [58] ape_5.5 reticulate_1.19 lifecycle_1.0.0
## [61] stringi_1.5.3 yaml_2.2.1 MASS_7.3-53.1
## [64] Rtsne_0.15 plyr_1.8.6 grid_3.6.3
## [67] parallel_3.6.3 listenv_0.8.0 promises_1.2.0.1
## [70] ggrepel_0.9.1 crayon_1.4.1 deldir_0.2-10
## [73] miniUI_0.1.1.1 lattice_0.20-41 cowplot_1.1.1
## [76] splines_3.6.3 tensor_1.5 knitr_1.33
## [79] pillar_1.6.0 igraph_1.2.6 future.apply_1.7.0
## [82] reshape2_1.4.4 codetools_0.2-18 leiden_0.3.7
## [85] glue_1.4.2 evaluate_0.14 data.table_1.14.0
## [88] vctrs_0.3.7 png_0.1-7 httpuv_1.6.0
## [91] polyclip_1.10-0 gtable_0.3.0 RANN_2.6.1
## [94] purrr_0.3.4 tidyr_1.1.3 future_1.21.0
## [97] xfun_0.22 rsvd_1.0.3 mime_0.10
## [100] xtable_1.8-4 RSpectra_0.16-0 later_1.2.0
## [103] survival_3.2-11 viridisLite_0.4.0 tibble_3.1.1
## [106] cluster_2.1.2 globals_0.14.0 fitdistrplus_1.1-3
## [109] ellipsis_0.3.1 ROCR_1.0-11