La obtención de transcriptomas de célula única es una metodología reciente y en continuo desarrollo, con un enorme potencial para describir diferentes sistemas biológicos. Entre otras aplicaciones, permite la detección de poblaciones celulares no abundantes y complejas, explorar la regulación de la expresión génica en diferentes condiciones, y en las diferentes etapas del desarrollo celular.
Aunque todas las células de nuestro cuerpo comparten genotipos casi idénticos, la información del transcriptoma en cualquier célula refleja la actividad de sólo un subconjunto de genes. Además, debido a que los diversos tipos de células en nuestro cuerpo expresan un transcriptoma único, la secuenciación clásica de bulk RNAseq puede proporcionar únicamente la señal de expresión promedio para un conjunto de células. La evidencia creciente sugiere que la expresión génica es heterogénea, incluso en tipos de células similares, y esto refleja la composición del tipo celular y también puede desencadenar decisiones sobre el destino celular.
Más de 20 protocolos diferentes han sido desarrollados desde el primer experimento de single cell RNAseq, y para la elección del protocolo adecuado es clave tener en cuenta la pregunta biológica que se quiere responder. Aquí trabajaremos con uno de estos protocolos, pero es importante tener en cuenta que es una técnica en continuo desarrollo y con muchas variantes.
En el práctico re-analizaremos parte de un set de datos publicados en 2019 (https://www.nature.com/articles/s41593-019-0491-3). En este artículo se realizó un análisis transcriptómico de single cell de cerebros de 8 ratones jóvenes y 8 ratones viejos utilizando la tecnología Chromium desarrollada por 10x Genomics. El análisis comparativo entre cerebros jovenes y envejecidos permitió identificar firmas genéticas que variaron de manera coordinada entre tipos celulares y conjuntos de genes que se regularon de manera específica para cada tipo de célula, incluso en direcciones opuestas. Estos datos revelaron que el envejecimiento, en lugar de inducir un programa universal, impulsó un curso transcripcional distinto en cada población celular, y destacaron procesos moleculares clave, incluyendo la biogénesis de ribosomas, subyacentes al envejecimiento cerebral.
En la actualidad, existen varias estrategias experimentales que permiten realizar single cell RNAseq, y las mismas se pueden subdividir en dos grandes grupos: aquellas que obtienen secuencias de transcriptos completos (Smartseq2, MATQ-seq, SUPeR-seq), y aquellas que obtienen secuencias de los extremos 5´ o 3´ de los transcriptos (Drop-seq, Seq-Well, Chromium, DroNC-seq, STRT-seq). La estrategia óptima dependerá de la pregunta biológica que se quiera responder y de los fondos disponibles para realizar el experimento, ya que cada estrategia tiene ventajas y desventajas.
Todas las estrategias tienen como principales objetivos: i) aislar células y ii) generar librerías de secuenciación que permitan identificar la célula de la que provienen. Los principales pasos para la generación de librerías de single cell RNAseq incluyen: lisis celular, transcripción reversa a ADNc de simple hebra, síntesis de la segunda hebra de ADNc y amplificación del ADNc. En los pasos de transcripción reversa, y de amplificación del ADNc, se suelen introducir barcodes (secuencias conocidas de unos pocos pb) únicos para cada célula. Luego de secuenciar, las lecturas con un mismo barcode serán asignadas a una misma célula, obteniendo los transcriptomas individuales. Además, algunas estrategias permiten agregar un barcode único para cada molécula (UMIs) en el paso de la transcripción reversa, lo cual permite luego corregir por los sesgos de la PCR.
Una vez obtenidas las lecturas del secuenciador, los primeros pasos son:
Análisis de calidad de las secuencias
Alineamiento frente al genoma de referencia asignando barcodes celulares
Corrección por sesgos de PCR (UMIs), en caso de ser posible.
Por razones de tiempo y consumo computacional nos saltearemos estos pasos. Los interesados pueden observar los detalles de estas actividades en los tutoriales disponibles en los siguiente links:
Como resultado de lo anterior, se obtiene una matriz de conteos donde las filas corresponden a los genes y las columnas a las distintas células analizadas. En muchas ocasiones los investigadores suelen publicar, además de los datos de secuenciación, las matrices de conteos, por lo que podemos partir de dichas matrices para las actividades del práctico.
Las matrices de conteos que vamos a utilizar en el presente práctico fueron obtenidas por los autores del artículo siguiendo el pipeline de CellRanger, y las descargamos del siguiente link:
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE129788
Utilizaremos los datos de scRNA-seq de 3 ratones jóvenes (réplicas) y 3 ratones envejecidos (réplicas).
Para comenzar el análisis, conéctese al servidor y diríjase a su carpeta personal. Dentro de su carpeta cree una nueva carpeta para realizar las tareas asociadas a este práctico.
mkdir Single_Cell_Analysis
cd Single_Cell_Analysis/ Realizaremos los análisis utilizando R:
RDentro de R, comenzamos cargando los paquetes necesarios para el análisis.
library(stringr)
library(ggplot2)
library(Seurat)
library(tidyverse)
library(cowplot)
library(Matrix.utils)
library(edgeR)
library(dplyr)
library(magrittr)
library(Matrix)
library(purrr)
library(reshape2)
library(S4Vectors)
library(tibble)
library(SingleCellExperiment)
library(pheatmap)
library(apeglm)
library(png)
library(DESeq2)
library(RColorBrewer)Luego, cargamos las matrices de conteos:
Y1 = read.table("/media/Disco10/curso_genomica/estudiantes/2024/Material_Practico_scRNAseq/GSM3722100_YX1L_10X.txt")
Y2= read.table("/media/Disco10/curso_genomica/estudiantes/2024/Material_Practico_scRNAseq/GSM3722101_YX2L_10X.txt")
Y3= read.table("/media/Disco10/curso_genomica/estudiantes/2024/Material_Practico_scRNAseq/GSM3722102_YX3R_10X.txt")
O1= read.table("/media/Disco10/curso_genomica/estudiantes/2024/Material_Practico_scRNAseq/GSM3722108_OX1X_10X.txt")
O2= read.table("/media/Disco10/curso_genomica/estudiantes/2024/Material_Practico_scRNAseq/GSM3722109_OX2X_10X.txt")
O3= read.table("/media/Disco10/curso_genomica/estudiantes/2024/Material_Practico_scRNAseq/GSM3722114_OX7X_10X.txt")Exploremos brevemente estas matrices. Con el comando View() se puede abrir las tablas que cargamos, ingresando el nombre del objeto dentro de los paréntesis.
View(Y1)También se puede observar las primeras líneas de la tabla con el comando head():
head(Y1)¿Qué información contiene esta matriz?
¿Por qué cree que hay tantos 0s en la matriz?
A continuación utilizaremos el paquete Seurat, que nos permite trabajar con las matrices de conteos, filtrarlas, y analizarlas para obtener resultados:
Seurat_Young_1 = CreateSeuratObject(Y1)
Seurat_Young_2 = CreateSeuratObject(Y2)
Seurat_Young_3 = CreateSeuratObject(Y3)
Seurat_Old_1 = CreateSeuratObject(O1)
Seurat_Old_2 = CreateSeuratObject(O2)
Seurat_Old_3 = CreateSeuratObject(O3)Y fusionaremos los objetos Seurat de las muestras en uno solo:
merged_seurat <- merge(Seurat_Young_1, y = c(Seurat_Young_2, Seurat_Young_3, Seurat_Old_1, Seurat_Old_2, Seurat_Old_3), add.cell.id = c("young", "young", "young", "old", "old", "old"))Además del control de calidad de las lecturas que se realiza normalmente en experimentos de RNAseq, los experimentos de single cell RNAseq deben controlar la calidad por célula, ya que mantener células de baja calidad nos influirá en los análisis posteriores. Para ello se pueden utilizar varias estrategias.
La proporción de lecturas provenientes de genes mitocondriales es un buen indicador de la viabilidad celular. Usualmente, una alta proporción de genes mitocondriales indica células dañadas. Para calcular dicha proporción necesitamos previamente definir los genes mitocondriales, y luego contar cuantas lecturas asociadas a dichos genes hay por cada célula.
El genoma mitocondrial del ratón contiene 37 genes: 13 genes que codifican proteínas involucradas en la cadena de transporte de electrones y la fosforilación oxidativa, 22 genes que codifican ARN de transferencia (ARNt) y 2 genes que codifican ARN ribosomales (ARNr).
Procederemos a calcular la proporción de lecturas provenientes de los 13 genes mitocondriales codificantes para proteínas, fácilmente identificables ya que sus nombres comienzan con “mt”.
Generaremos una tabla de metadatos que utilizaremos para el filtro de calidad:
merged_seurat$mitoRatio <- PercentageFeatureSet(object = merged_seurat, pattern = "^mt-")
merged_seurat$mitoRatio <- merged_seurat@meta.data$mitoRatio / 100Es relevante calcular cuantos genes se detectan por UMI alineado, para ello podremos calcular lo siguiente:
merged_seurat$log10GenesPerUMI <- log10(merged_seurat$nFeature_RNA) / log10(merged_seurat$nCount_RNA)
metadata <- merged_seurat@meta.data
metadata$cells <- rownames(metadata)
metadata <- metadata %>%
dplyr::rename(seq_folder = orig.ident,
nUMI = nCount_RNA,
nGene = nFeature_RNA)
metadata$sample <- NA
metadata$sample[which(str_detect(metadata$cells, "^young_"))] <- "young"
metadata$sample[which(str_detect(metadata$cells, "^old_"))] <- "old"
merged_seurat@meta.data <- metadataCon los comandos ejecutados hasta el momento tenemos una tabla de metadatos con la información necesaria para analizar la calidad de las células y realizar los filtros necesarios.
Visualizaremos la calidad de las células realizando gráficos con los datos previamente calculados.
En primer lugar, analizaremos la cantidad de células secuenciadas para cada condición:
metadata %>%
ggplot(aes(x=sample, fill=sample)) +
geom_bar() +
theme_classic() +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
theme(plot.title = element_text(hjust=0.5, face="bold")) +
ggtitle("NCells")Para visualizar la cantidad de UMIs detectados por célula:
metadata %>%
ggplot(aes(color=sample, x=nUMI, fill= sample)) +
geom_density(alpha = 0.2) +
scale_x_log10() +
theme_classic() +
ylab("Cell density") +
geom_vline(xintercept = 1000)Para visualizar la cantidad de genes detectados por célula:
metadata %>%
ggplot(aes(color=sample, x=nGene, fill= sample)) +
geom_density(alpha = 0.2) +
theme_classic() +
scale_x_log10() +
geom_vline(xintercept = 500)metadata %>%
ggplot(aes(x=sample, y=log10(nGene), fill=sample)) +
geom_boxplot() +
theme_classic() +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
theme(plot.title = element_text(hjust=0.5, face="bold")) +
ggtitle("NCells vs NGenes")¿Cuántos genes se detectan por célula?
Para visualizar el ratio de genes mitocondriales detectados:
metadata %>%
ggplot(aes(color=sample, x=mitoRatio, fill=sample)) +
geom_density(alpha = 0.2) +
scale_x_log10() +
theme_classic() +
geom_vline(xintercept = 0.2)Al graficar el número de genes en función del nUMI, podemos ver la relación entre el número de lecturas alineadas y el número de genes detectados. Las células de baja calidad suelen tener bajo número de lecturas y bajo número de genes detectados.
metadata %>%
ggplot(aes(x=nUMI, y=nGene, color=mitoRatio)) +
geom_point() +
scale_colour_gradient(low = "gray90", high = "black") +
stat_smooth(method=lm) +
scale_x_log10() +
scale_y_log10() +
theme_classic() +
geom_vline(xintercept = 1000) +
geom_hline(yintercept = 500) +
facet_wrap(~sample)¿Considera que es necesario realizar un filtrado?
También podemos analizar según la complejidad de las células presentes. En este caso, esperamos una complejidad similar ya que las células provienen del mismo tejido. Complejidad similar implica un número similar de genes detectado por lecturas mapeadas. Outliers en complejidad podrían ser detectado mediante el siguiente gráfico:
metadata %>%
ggplot(aes(x=log10GenesPerUMI, color = sample, fill=sample)) +
geom_density(alpha = 0.2) +
theme_classic() +
geom_vline(xintercept = 0.8)¿Qué opina de sus células? ¿Deberíamos filtrar por complejidad?
Considerando los gráficos y análisis anteriores, procederemos a filtrar en primer lugar según la cantidad de UMIs alineados por célula y la cantidad de genes detectados por célula. La idea es quitar aquellas células con baja cantidad de UMIs y pocos genes detectados. Discuta los parámetros numéricos a utilizar en clase.
filtered_seurat <- subset(x = merged_seurat,
subset= (nUMI >= 1000) &
(nGene >= 500) &
(log10GenesPerUMI > 0.80) &
(mitoRatio < 0.10))¿Cuántas células tiene el nuevo archivo filtrado? ¿Cuántas tiene el archivo sin filtrar? Además del filtrado de células, es usual realizar filtrado a nivel de genes. Este tipo de datos suele tener muchos genes sin lecturas asignadas en cada célula. Por lo tanto, usualmente se filtran aquellos genes que no se detectan para evitar ruido en posteriores análisis. Por lo tanto, eliminaremos todos los genes que no presenten al menos 10 células con al menos una lectura.
filtered_seurat[["RNA"]] <- JoinLayers(filtered_seurat[["RNA"]])
counts <- GetAssayData(object = filtered_seurat, assay = "RNA")
nonzero <- counts > 0
keep_genes <- Matrix::rowSums(nonzero) >= 10
filtered_counts <- counts[keep_genes, ]
filtered_seurat <- CreateSeuratObject(filtered_counts, meta.data = filtered_seurat@meta.data)¿Cuántos genes tiene el nuevo archivo?
Limpiemos el espacio de R:
remove(merged_seurat, metadata, O1, O2, O3, Y1, Y2, Y3, Seurat_Old_1, Seurat_Old_2, Seurat_Old_3, Seurat_Young_1, Seurat_Young_2, Seurat_Young_3)CHECKPOINT Nº1
Si se le desconecta la PC del servidor, o se le cierra R por algún motivo, puede cargar el archivo conteniendo los conteos y la metadata con el siguiente comando dentro de R:
load("/media/Disco10/curso_genomica/estudiantes/2024/Material_Practico_scRNAseq/CheckPoint1.RData"
)Luego de estos pasos, suele ser recomendable volver a analizar la calidad de los datos para confirmar que todo esté correcto y continuar los análisis posteriores.
El primer paso para analizar nuestra matriz de conteos será el de normalización. Para ello emplearemos dos métodos en paralelo, el método NormalizeData que consiste en una normalización básica dónde se escala los conteos por muestra y se logaritmiza, y el método SCTransform. El método SCTransform de Seurat es más complejo, y se basa en ajustar los datos a una binomial negativa regularizada, y ajusta la varianza basado en el pooleo de genes de abundancia similar. Sin embargo, contamos con matrices de conteos de dos condiciones diferentes, con tres réplicas cada una, por lo que además de normalizar deberemos integrar los datos de modo de descartar posibles efectos de batch. Para ello correremos los siguientes comandos:
# Separar objetos Seurat por grupo
split_seurat <- SplitObject(filtered_seurat, split.by = "sample")
split_seurat <- split_seurat[c("young", "old")]
# Normalizar datos y realizar an?lisis de integraci?n
for (i in 1:length(split_seurat)) {
split_seurat[[i]] <- NormalizeData(split_seurat[[i]], verbose = TRUE)
split_seurat[[i]] <- SCTransform(split_seurat[[i]], vars.to.regress = c("mitoRatio"))
}
integ_features <- SelectIntegrationFeatures(object.list = split_seurat,
nfeatures = 3000)
split_seurat <- PrepSCTIntegration(object.list = split_seurat,
anchor.features = integ_features)
integ_anchors <- FindIntegrationAnchors(object.list = split_seurat,
normalization.method = "SCT",
anchor.features = integ_features)
seurat_integrated <- IntegrateData(anchorset = integ_anchors,
normalization.method = "SCT")Por defecto, SCTransform calcula los 3000 genes más variables, los cuales utilizará para parte de los análisis posteriores. Puede ser recomendable aumentar el número de genes a utilizar. Además, esta función provoca que por defecto el Seurat utilice el assay SCT para posteriores análisis. Puede visualizar los assays con el siguiente comando y apretando la tecla Tab: seurat_integrated@assays$
CHECKPOINT Nº2
Si se le desconecta la PC del servidor, o se le cierra R por algún motivo, puede cargar el archivo conteniendo los conteos y la metadata con el siguiente comando dentro de R:
load("../../Material_Practico_scRNAseq/CheckPoint2.RData"
)Luego de realizada la normalización procedemos al análisis de reducción dimensional que nos permite visualizar y estudiar la variabilidad de nuestros datos en un gráfico de dos dimensiones. Para grandes set de datos, el método de PCA no logra una buena separación de las variables, por lo que se suelen realizar análisis no lineales como el t-SNE o el UMAP. De todas maneras, ambos métodos funcionan basados en una reducción de dimensiones previa realizada mediante PCA. Con los siguientes comandos realizaremos las reducciones por PCA y luego por UMAP, para obtener el gráfico de reducción dimensional.
# Realizar análisis de PCA
seurat_integrated <- RunPCA(object = seurat_integrated)
# Visualizar resultados del PCA
PCAPlot(seurat_integrated,
split.by = "sample")
# Realizar análisis de UMAP
seurat_integrated <- RunUMAP(seurat_integrated,
dims = 1:50,
reduction = "pca")
DimPlot(seurat_integrated)
# UMAP separado por muestra
DimPlot(seurat_integrated,
split.by = "sample") Limpiemos el espacio de R:
remove(split_seurat, nonzero, integ_anchors, filtered_seurat, filtered_counts, counts)A continuación, nos interesa saber qué tipos celulares hay en nuestra muestra, para lo que vamos a realizar el clustering. Es importante conocer la muestra que estamos analizando, para saber qué tipos celulares esperamos encontrar (células en diferenciación, células con alta complejidad, células con alta expresión de genes mitocondriales, etc). Un paso clave en el clustering es determinar cuántos componentes principales (PC) vamos a utilizar para el mismo. Es importante seleccionar una cantidad de PC tal que se capte la mayor variabilidad entre tipos celulares posible. El método SCTransform logra estimar la varianza de mejor manera que los métodos anteriores, generando que la selección de PC no sea un paso tan crítico. En teoría, al normalizar con SCTransform y utilizar un mayor número de PC significa considerar una mayor proporción de variabilidad. Sin embargo, a mayor número de PCs, mayor consumo computacional. Para el presente práctico probaremos utilizando 50 PCs. Por defecto, el Seurat se basa en una estructura de grafos (K Neareast Neighbors graph) relacionando a células con patrones similares de expresión génica. Utilizaremos la función FindClusters para realizar el clustering basado en grafos. El parámetro resolution setea la granularidad de cada cluster, y debe ser optimizado para cada experimento individualmente. Valores altos de resolución dan lugar a mayor número de clusters.
# Identificar clusters
seurat_integrated <- FindNeighbors(object = seurat_integrated,
dims = 1:50)
seurat_integrated <- FindClusters(object = seurat_integrated,
resolution = c(0.4, 0.6, 0.8, 1.0, 1.4))Podemos observar que los resultados son almacenados en la tabla de metadatos, y que podemos acceder a ellos de la siguiente forma:
seurat_integrated@meta.data$integrated_snn_resEn primer lugar analizaremos los clusters generados a una resolución intermedia (0.8):
Idents(object = seurat_integrated) <- "integrated_snn_res.0.8"
DimPlot(seurat_integrated,
reduction = "umap",
label = TRUE,
label.size = 6)Pruebe con otras resoluciones definiendo de nuevo las identidades y realizando el gráfico. ¿Qué observa?
Proseguiremos en el práctico con el valor de resolución 0.8 (si probó varios valores de resolución, termine ejecutando de nuevo el comando de identidad para resolución 0.8.
Idents(object = seurat_integrated) <- "integrated_snn_res.0.8"Tenga en cuenta que el número de clusters puede ser reajustado nuevamente si observamos incoherencias en análisis posteriores, por lo que nos será útil ya haber calculado clusters a diferentes resoluciones. A continuación procederemos a un nuevo control de calidad, esta vez de los clusters, analizando si existen propiedades que nos hagan dudar de los mismos. Primero, analicemos la cantidad de células por cluster.
n_cells <- FetchData(seurat_integrated,
vars = c("ident", "orig.ident")) %>%
dplyr::count(ident, orig.ident) %>%
tidyr::spread(ident, n)
View(n_cells)
DimPlot(seurat_integrated,
label = TRUE,
split.by = "sample") + NoLegend()También podemos visualizar la distribución de diferentes métricas de calidad en los clusters:
metrics <- c("nUMI", "nGene", "mitoRatio")
FeaturePlot(seurat_integrated,
reduction = "umap",
features = metrics,
pt.size = 0.4,
order = TRUE,
min.cutoff = 'q10',
label = TRUE)¿Qué opina de las métricas? ¿Se encuentran parejas a lo largo de los clusters?
La función FeaturePlot nos permite analizar la expresión de ciertos genes. Como ya conocemos la muestra, podemos esperar la presencia de ciertos tipos celulares y estudiar la expresión de marcadores conocidos. Para visualizar los niveles de expresión es recomendable estudiar los valores normalizados con el primer método.
DefaultAssay(seurat_integrated) <- "RNA"Dado que estamos trabajando con datos de cerebro, esperamos encontrar principalmente neuronas astrocitos, oligodendrocitos y células endoteliales. Por lo tanto resulta importante tener en cuenta los marcadores conocidos de cada tipo celular para poder identificar los clusters detectados. Se sabe que un marcador de neuronas maduras es Syt1, y que un marcador de Oligodendrocitos es Cldn11. Para más marcadores puede buscar en la literatura o dirigirse al artículo que estamos analizando.
# Visualizar expresión de genes espec?ficos
FeaturePlot(seurat_integrated,
reduction = "umap",
features = c("Cldn11"),
order = TRUE,
min.cutoff = 'q10',
label = TRUE)
FeaturePlot(seurat_integrated,
reduction = "umap",
features = c("Syt1"),
order = TRUE,
min.cutoff = 'q10',
label = TRUE)¿Qué funciones tienen Cldn11 y Syt1?
Basado en la lista de marcadores utilizados por los autores del artículo generamos un vector conteniendo algunos marcadores conocidos:
# Identificar genes marcadores
Known_Markers = c("Pdgfra",
"Cldn11",
"Npy",
"Thbs4",
"Cd44",
"Gja1",
"Cdk1",
"Sox11",
"Syt1",
"Baiap3",
"Sspo",
"Rax",
"Ttr",
"Cldn5",
"Kcnj8",
"Acta2",
"Alas2",
"Slc6a13",
"Slc47a1",
"Tmem119",
"Plac8",
"Pf4",
"Cd209a",
"S100a9")Podemos visualizar los niveles de expresión junto al procentaje de células que expresan estos genes en el siguiente gráfico:
# Visualizar genes marcadores
DotPlot(seurat_integrated, features=Known_Markers)Es importante notar que, si detectamos marcadores indicando que tenemos más de un tipo celular por cluster, sería recomendable fijar un valor más alto de resolución. En último lugar, la función FindAllMarkers permite realizar un análisis de expresión diferencial, comparando un cluster contra todo el resto de clusters. Cada célula es tomada como una réplica, y se realiza un test estadístico que permite determinar fold change y p-valor. Por defecto, el Seurat realiza el Wilcoxon Rank Sum Test. A la función se le fijan parámetros claves, como:
logfc.threshold - Logaritmo en base 2 del Fold Change límite
min.pct - Solo testear genes que tengan este valor mínimo de fracción de células que lo expresan en el grupo (default es 0.1)
seurat_integrated[["RNA"]] <- JoinLayers(seurat_integrated[["RNA"]])
markers <- FindAllMarkers(object = seurat_integrated,
only.pos = TRUE,
logfc.threshold = 1)
CHECKPOINT Nº3
Si se le desconecta la PC del servidor, o se le cierra R por algún motivo, puede cargar el archivo conteniendo los conteos y la metadata con el siguiente comando dentro de R:
load("../../Material_Practico_scRNAseq/CheckPoint2.RData"
)NOTA: Dado que cada célula se trata como una réplica, esto dará como resultado valores estadísticos de p extremadamente pequeños. Un gen puede tener un p-valor tan bajo como < 1E-50, aunque eso no implique que se trate de un gen marcador confiable.
Al estar comparando un cluster contra el resto, este análisis de expresión diferencial es muy útil para detectar marcadores de cluster. Se sugiere buscar como marcadores a genes que posean grandes diferencias en el porcentaje de células que lo expresan, y con un Fold Change grande.
Además de la búsqueda de marcadores, resulta interesante analizar la expresión diferencial entre dos clusters de interés, de modo de analizar diferencias biológicas entre dos subtipos neuronales. Ejemplo: Markers1_2 <- FindMarkers(seurat_integrated, ident.1 = 1, ident.2 = 2) Una vez que tengamos a los clusters identificados podemos renombrarlos, incluso será útil nombrar clusters de células que hayamos detectado como dañadas o que queramos descartar: filtered_seurat <- RenameIdents(object = seurat_integrated, “7” = “Cell_Name”)
DimPlot(seurat_integrated, reduction = “umap”, label = TRUE, label.size = 3)
También podríamos eliminar a los clusters de células que hayamos identificado como células dañadas, o no confiables.
seurat_subset_labeled <- subset(seurat_integrated, idents = “Nº”, invert = TRUE)
Luego de analizados los marcadores por clusters tendremos clasificados los tipos celulares, y se puede estudiar en mayor profundidad las diferencias entre cada subtipo celular.
En este set de datos se tiene células de ratones adultos y de ratones viejos, con el fin de estudiar cómo se altera la expresión génica en el proceso de envejecimiento. En single cell RNAseq, la forma más aceptada de hacer análisis de expresión diferencial entre condiciones consiste en hacer análisis denominados pseudobulkRNAseq. Básicamente, el mismo consiste en sumar los UMIs de una misma célula y réplica biológica de modo de construir una nueva tabla de conteos (ejemplo: sumo todas las lecturas de un mismo ratón que provienen de células que clasificamos previamente como neuronas). Luego, con estas tablas de pseudobulk RNAseq, uno puede realizar análisis de expresión diferencial con herramientas clásicas de bulk RNAseq como DESeq2 o edgeR, para preguntarse específicamente cómo varía la expresión génica de un tipo celular en particular entre condiciones.
En el presente práctico, nos planteamos estudiar la expresión diferencial en neuronas. Para ello, en primer lugar realizamos un subset del objeto seurat de modo de quedarnos únicamente con las neuronas:
neurons <- subset(seurat_integrated, idents = c(2,11,12,17,18,20))
DimPlot(neurons)Luego procederemos a generar nuestra tabla de pseudobulk RNAseq:
# Crear objeto SingleCellExperiment
neurons[["RNA"]] <- JoinLayers(neurons[["RNA"]])
counts <- GetAssayData(object = neurons, assay = "RNA")
metadata = neurons@meta.data
sce <- SingleCellExperiment(assays = list(counts = counts),
colData = metadata)
groups <- colData(sce)[, c("orig.ident", "seq_folder")]
kids <- purrr::set_names(unique(sce$orig.ident))
nk <- length(kids)
sids <- purrr::set_names(unique(sce$seq_folder))
ns <- length(sids)
n_cells <- as.numeric(table(sce$seq_folder))
m <- match(sids, sce$seq_folder)
ei <- data.frame(colData(sce)[m, ],
n_cells, row.names = NULL) %>%
select(-"orig.ident")
groups <- colData(sce)[, c("orig.ident", "seq_folder")]
# Agregar conteo de genes a grupos
pb <- aggregate.Matrix(t(counts(sce)),
groupings = groups, fun = "sum")
pb = round(pb)
pb[1:6, 1:6]
pb = t(pb)
metadata = neurons@meta.dataGeneramos una tabla que indique al DESeq2 las condiciones asociadas a nuestra tabla de conteos:
colData = cbind(row.names = colnames(pb), c(rep("Old", 3), rep("Young",3)))
colData = as.data.frame(colData, row.names = colData[,1])
colData$row.names = NULL
colnames(colData) = "Age"Ingresamos la tabla de DESeq2:
dds <- DESeqDataSetFromMatrix(pb,
colData = colData,
design = ~ Age)
rld <- rlog(dds, blind=TRUE)Visualizamos el PCA para ver como agrupan nuestras réplicas biológicas:
DESeq2::plotPCA(rld, intgroup = "Age")Realizamos la normalización y ajuste con DESeq2:
dds <- DESeq(dds)Generamos una tabla de genes expresados diferencialmente entre ratones jóvenes y viejos:
contrast <- c("Age", levels(colData$Age)[2], levels(colData$Age)[1])
res <- results(dds,
alpha = 0.05)
resultados = as.data.frame(res)Luego, con la tabla de expresión diferencial se podría continuar de manera similar a como ya vieron en prácticos anteriores (análisis de enriquecimiento funcional, vías metabólicas, etc).
Ejemplo con el VlnPlot:
tiff(filename="./plotprueba.tiff", width=4400, height=4400, units="px", pointsize=8, compression="lzw", res=1000,)
VlnPlot(object = filtered_seurat, features = "SLC1A3") dev.off();Ejemplo con tabla de marcadores de cluster:
write.csv(markers, "markers.csv")