Análisis de datos scATAC-seq en Signac: cerebro murino
Adam Casas
Compilado: 14 de septiembre del 2021
1 Instalación
Instalamos Signac en R versión 4.1:
BiocManager::install("Rsamtools")
install.packages("Signac")2 Introducción
Signac (Stuart et al., 2020) es una extensión del famoso paquete Seurat para explorar, analizar e interpretar datasets de scATAC-seq. ATAC-seq (y su variante de célula única) es, por decirlo de alguna manera, la sucesora espiritural de técnicas de secuenciación como ChIP-seq y RRBS-Seq (un tipo de secuenciación por bisulfito). Si en dichas técnicas se estudia el estado epigenético de las células/tejidos de interés, cierto es que se debe tener un conocimiento previo de cómo funcionan los mecanismos epigenéticos del organismo en cuestión para poder sacarles el mayor rendimiento. ATAC-seq, por otro lado, evade dicha imposición al estudiar directamente la accesibilidad de la cromatina a la transposasa hiperactiva Tn5 de manera “agnóstica” (sin necesidad de conocer a priori los mecanismos epigenéticos subyacentes).
En este informe vamos a presentar Signac mediante el análisis de un dataset compuesto por 5000 células provenientes de la corteza visual primaria de cerebros de ratones adultos. Para ello debemos descargar primero los siguientes archivos, disponibles en la página web de 10X Genomics:
# El archivo con los datos crudos del ATAC-seq
download.file(
"http://cf.10xgenomics.com/samples/cell-atac/1.1.0/atac_v1_adult_brain_fresh_5k/atac_v1_adult_brain_fresh_5k_filtered_peak_bc_matrix.h5",
destfile = "./archivos_accesorios/atac_v1_adult_brain_fresh_5k_filtered_peak_bc_matrix.h5")
# Sus metadatos
download.file(
"http://cf.10xgenomics.com/samples/cell-atac/1.1.0/atac_v1_adult_brain_fresh_5k/atac_v1_adult_brain_fresh_5k_singlecell.csv",
destfile = "./archivos_accesorios/atac_v1_adult_brain_fresh_5k_singlecell.csv")
# El archivo de los fragmentos
download.file(
"http://cf.10xgenomics.com/samples/cell-atac/1.1.0/atac_v1_adult_brain_fresh_5k/atac_v1_adult_brain_fresh_5k_fragments.tsv.gz",
destfile = "./archivos_accesorios/atac_v1_adult_brain_fresh_5k_fragments.tsv.gz")
# El archivo con el índice de los fragmentos
download.file(
"http://cf.10xgenomics.com/samples/cell-atac/1.1.0/atac_v1_adult_brain_fresh_5k/atac_v1_adult_brain_fresh_5k_fragments.tsv.gz.tbi",
destfile = "./archivos_accesorios/atac_v1_adult_brain_fresh_5k_fragments.tsv.gz.tbi")Ahora cargamos las librerías que vamos a usar:
# Análisis
library("Signac")
library("Seurat")
# Genomas
library("GenomeInfoDb")
library("EnsDb.Mmusculus.v79")
# Gráficos
library("ggplot2")
library("patchwork")3 Preprocesado
Signac requiere dos archivos generados por el software CellRanger para poder preprocesar datos de ATAC-Seq (aunque al final del informe se muestra cómo actuar si los archivos de secuenciado se generaron con un software distinto):
Peak/Cell matrix o matriz de picos/célula. La matriz de picos/célula es similar a la matriz de expresión génica de dimensiones genes x células usada en scRNA-Seq. En ATAC-Seq la matriz es de regiones cromosómicas x células. Cada valor de la matriz es el nº de sitios de integración de Tn5 en cada célula que mapean en dichas regiones o peaks. A mayor nº de sitios de integración mapeados en un locus, más abierta está la cromatina en dicha región.
Fragment file o archivo de fragmentos. Representa la lista completa de fragmentos únicos (= lecturas) detectados en todas las células y es por tanto un archivo MUY pesado. Al ser tan pesado, se guarda en el disco duro en lugar de la RAM y a consecuencia manipularlo es un proceso lento. Su ventaja es que contiene todos los sitios de integración de la Tn5 en cada célula.
Ya que Signac es una extensión de Seurat, vamos a trabajar con un objeto seurat generado a partir de la matriz de picos/célula, sus metadatos creados en cellranger-atac y el archivo de fragmentos:
# Cargamos y visualizamos las 6 primeras filas de la matriz de ATAC_Seq en el
# objeto `matriz_picos_ATAC`:
matriz_picos_ATAC <- Read10X_h5(
"E:/Git/R/Signac MB/archivos_accesorios/atac_v1_adult_brain_fresh_5k_filtered_peak_bc_matrix.h5")
matriz_picos_ATAC[c(1:5), c(1:5)]## 5 x 5 sparse Matrix of class "dgCMatrix"
## AAACGAAAGTAATCAG-1 AAACGAACACGCTGTG-1 AAACGAATCCTGGGAC-1
## chr1:3094708-3095552 . . .
## chr1:3119414-3121782 2 . .
## chr1:3204809-3205178 . . .
## chr1:3217330-3217359 . . .
## chr1:3228123-3228463 2 . .
## AAACGAATCGGGAAAC-1 AAACTCGAGAAAGCAG-1
## chr1:3094708-3095552 . .
## chr1:3119414-3121782 . .
## chr1:3204809-3205178 . .
## chr1:3217330-3217359 . .
## chr1:3228123-3228463 . .# Cargamos y exploramos el archivo de metadatos
metadatos <- read.csv(
file = "E:/Git/R/Signac MB/archivos_accesorios/atac_v1_adult_brain_fresh_5k_singlecell.csv",
header = TRUE,
row.names = 1)
metadatos[c(1:5), c(1:5)]dim(metadatos)## [1] 499704 17summary(metadatos)## total duplicate chimeric unmapped
## Min. : 1 Min. : 0.0 Min. : 0.000 Min. : 0.0
## 1st Qu.: 1 1st Qu.: 0.0 1st Qu.: 0.000 1st Qu.: 0.0
## Median : 3 Median : 1.0 Median : 0.000 Median : 0.0
## Mean : 488 Mean : 248.7 Mean : 3.099 Mean : 9.9
## 3rd Qu.: 9 3rd Qu.: 2.0 3rd Qu.: 0.000 3rd Qu.: 2.0
## Max. :4032398 Max. :2070733.0 Max. :21036.000 Max. :880066.0
## lowmapq mitochondrial passed_filters cell_id
## Min. : 0.00 Min. : 0.00 Min. : 0.0 Length:499704
## 1st Qu.: 0.00 1st Qu.: 0.00 1st Qu.: 0.0 Class :character
## Median : 0.00 Median : 0.00 Median : 1.0 Mode :character
## Mean : 25.02 Mean : 8.34 Mean : 193.3
## 3rd Qu.: 1.00 3rd Qu.: 0.00 3rd Qu.: 4.0
## Max. :205322.00 Max. :58820.00 Max. :796421.0
## is__cell_barcode TSS_fragments DNase_sensitive_region_fragments
## Min. :0.00000 Min. : 0.00 Min. : 0.0
## 1st Qu.:0.00000 1st Qu.: 0.00 1st Qu.: 0.0
## Median :0.00000 Median : 0.00 Median : 1.0
## Mean :0.01068 Mean : 62.47 Mean : 112.6
## 3rd Qu.:0.00000 3rd Qu.: 1.00 3rd Qu.: 2.0
## Max. :1.00000 Max. :36136.00 Max. :71824.0
## enhancer_region_fragments promoter_region_fragments on_target_fragments
## Min. : 0.00 Min. : 0.00 Min. : 0.0
## 1st Qu.: 0.00 1st Qu.: 0.00 1st Qu.: 0.0
## Median : 0.00 Median : 0.00 Median : 1.0
## Mean : 56.27 Mean : 54.83 Mean : 138.6
## 3rd Qu.: 1.00 3rd Qu.: 1.00 3rd Qu.: 2.0
## Max. :38517.00 Max. :30549.00 Max. :89107.0
## blacklist_region_fragments peak_region_fragments peak_region_cutsites
## Min. : 0.000 Min. : 0.0 Min. : 0.0
## 1st Qu.: 0.000 1st Qu.: 0.0 1st Qu.: 0.0
## Median : 0.000 Median : 1.0 Median : 2.0
## Mean : 0.467 Mean : 129.2 Mean : 242.4
## 3rd Qu.: 0.000 3rd Qu.: 2.0 3rd Qu.: 4.0
## Max. :572.000 Max. :78021.0 Max. :144369.0# Creamos un objeto de tipo `ChromatinAssay`, necesario para generar el objeto
# de tipo `seurat`
ensayo_ATAC_Seq <- CreateChromatinAssay(
counts = matriz_picos_ATAC,
sep = c(":", "-"),
genome = "mm10",
fragments = "E:/Git/R/Signac MB/archivos_accesorios/atac_v1_adult_brain_fresh_5k_fragments.tsv.gz",
min.cells = 1)
# Finalmente creamos el objeto `seurat`
MB_Signac <- CreateSeuratObject(
counts = ensayo_ATAC_Seq,
assay = "peaks",
project = "ATAC",
meta.data = metadatos)
MB_Signac## An object of class Seurat
## 157203 features across 5337 samples within 1 assay
## Active assay: peaks (157203 features, 0 variable features)MB_Signac[["peaks"]]## ChromatinAssay data with 157203 features for 5337 cells
## Variable features: 0
## Genome: mm10
## Annotation present: FALSE
## Motifs present: FALSE
## Fragment files: 1Los datos de la secuenciación ATAC se guardan en un objeto de tipo ChromatinAssay, el cual se caracteriza por presentar bolsillos para almacenar información adicional sobre motifs, anotaciones de los genes e información del genoma. A partir de dicho objeto creamos el principal objeto seurat donde realizaremos los análisis.
Podemos usar la función granges() en un objeto de tipo seurat que tenga como experimento activo uno de tipo ChromatinAssay para observar las regiones cromosómicas (provenientes del archivo de fragmentos) asociadas a cada gen del objeto. Para más información sobre la clase ChromatinAssay, consulte la viñeta de interacción con objetos.
# Exploración del ensayo/experimento de ATAC-Seq
ensayo_ATAC_Seq[1:5, 1:5]## 5 x 5 sparse Matrix of class "dgCMatrix"
## AAACGAAAGTAATCAG-1 AAACGAACACGCTGTG-1 AAACGAATCCTGGGAC-1
## chr1-3094708-3095552 . . .
## chr1-3119414-3121782 2 . .
## chr1-3204809-3205178 . . .
## chr1-3217330-3217359 . . .
## chr1-3228123-3228463 2 . .
## AAACGAATCGGGAAAC-1 AAACTCGAGAAAGCAG-1
## chr1-3094708-3095552 . .
## chr1-3119414-3121782 . .
## chr1-3204809-3205178 . .
## chr1-3217330-3217359 . .
## chr1-3228123-3228463 . .# Las columnas son las células secuenciadas y las filas son regiones
# cromosómicas... 5337 células y 157.203 regiones cromosómicas estudiadas
dim(ensayo_ATAC_Seq) ## [1] 157203 5337granges(MB_Signac)## GRanges object with 157203 ranges and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## [1] chr1 3094708-3095552 *
## [2] chr1 3119414-3121782 *
## [3] chr1 3204809-3205178 *
## [4] chr1 3217330-3217359 *
## [5] chr1 3228123-3228463 *
## ... ... ... ...
## [157199] chrY 90761049-90762169 *
## [157200] chrY 90800417-90800831 *
## [157201] chrY 90804515-90805440 *
## [157202] chrY 90808545-90809111 *
## [157203] chrY 90810741-90810960 *
## -------
## seqinfo: 21 sequences from an unspecified genome; no seqlengthsTambién podemos añadir al objeto MB_Signac las anotaciones del genoma murino para que ciertos comandos puedan hacer uso de dicha información.
# Descargamos las anotaciones del genoma murino mm10 desde la base de datos
# Ensembl
annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Mmusculus.v79, verbose = F)
# Las anotaciones están en formato NCBI. Vamos a cambiarlas al formato UCSC y
# especificaremos que estamos trabajando con el genoma mm10
seqlevelsStyle(annotations) <- "UCSC"
genome(annotations) <- "mm10"
# Añadimos las anotaciones al objeto `seurat`
Annotation(MB_Signac) <- annotations
Annotation(MB_Signac)## GRanges object with 1763965 ranges and 5 metadata columns:
## seqnames ranges strand | tx_id
## <Rle> <IRanges> <Rle> | <character>
## ENSMUSE00001236884 chr3 3508030-3508332 + | ENSMUST00000108393
## ENSMUSE00000676606 chr3 3634150-3634347 + | ENSMUST00000108394
## ENSMUSE00001345708 chr3 3638059-3638230 + | ENSMUST00000108393
## ENSMUSE00001345708 chr3 3638059-3638230 + | ENSMUST00000108394
## ENSMUSE00000149313 chr3 3641223-3641317 + | ENSMUST00000108393
## ... ... ... ... . ...
## ENSMUST00000082414 chrM 10167-11544 + | ENSMUST00000082414
## ENSMUST00000082418 chrM 11742-13565 + | ENSMUST00000082418
## ENSMUST00000082419 chrM 13552-14070 - | ENSMUST00000082419
## ENSMUST00000082421 chrM 14145-15288 + | ENSMUST00000082421
## ENSMUST00000084013 chrM 9877-10173 + | ENSMUST00000084013
## gene_name gene_id gene_biotype type
## <character> <character> <character> <factor>
## ENSMUSE00001236884 Hnf4g ENSMUSG00000017688 protein_coding exon
## ENSMUSE00000676606 Hnf4g ENSMUSG00000017688 protein_coding exon
## ENSMUSE00001345708 Hnf4g ENSMUSG00000017688 protein_coding exon
## ENSMUSE00001345708 Hnf4g ENSMUSG00000017688 protein_coding exon
## ENSMUSE00000149313 Hnf4g ENSMUSG00000017688 protein_coding exon
## ... ... ... ... ...
## ENSMUST00000082414 mt-Nd4 ENSMUSG00000064363 protein_coding cds
## ENSMUST00000082418 mt-Nd5 ENSMUSG00000064367 protein_coding cds
## ENSMUST00000082419 mt-Nd6 ENSMUSG00000064368 protein_coding cds
## ENSMUST00000082421 mt-Cytb ENSMUSG00000064370 protein_coding cds
## ENSMUST00000084013 mt-Nd4l ENSMUSG00000065947 protein_coding cds
## -------
## seqinfo: 22 sequences from mm10 genome4 Métricas de calidad (QC)
Ahora podemos comenzar con el control de calidad. Los desarrolladores de Signac recomiendan prestar atención a las siguientes métricas:
- Nucleosome banding pattern: El histograma de los tamaños de los fragmentos de ADN debe exhibir un marcado patrón de bandeo de nucleosomas correspondiente a la longitud del ADN enrollado alrededor de un único nucleosoma. Esta métrica se calcula mediante la siguiente ecuación:
\[nucleosome\ signal = \frac{fragmentos\_mononucleosomales}{fragmentos\_sin\_nucleosomas}\]
Transcriptional start site (TSS) enrichment score. El proyecto ENCODE define el TSS enrichment score como el ratio de fragmentos centrados en el TSS vs fragmentos flanqueantes al TSS. Los experimentos de ATAC-Seq mal hechos suelen tener un TSS enrichment score bajo. Esta métrica se calcula para cada célula con el comando
TSSEnrichment()y se almacena en los metadatos en la columnaTSS.enrichment.Número total de fragmentos (= lecturas) en picos: Las células con muy pocas lecturas se han secuenciado con poca profundidad, mientras que las células con muchas lecturas podrían ser doublets (dos células secuenciadas como una), agregados de núcleos u otro tipo de artefactos. Por consiguiente se deben eliminar todas las células con valores extremos de lecturas en picos.
Ratio de lecturas en picos: Representa el porcentaje de fragmentos (i.e. regiones cromosómicas) que caen en picos. Las células con un ratio bajo (<15~20%) suelen ser células mal secuenciadas o artefactos técnicos (ruido) y por tanto se deberían eliminar. Tenga en cuenta que el umbral de rechazo de células depende del set de picos usado.
Ratio de lecturas en regiones genómicas ruidosas (blacklists). El proyecto ENCODE también ha publicado para el genoma humano (hg19 y GRCh38) y otros una lista negra de regiones cromosómicas que consiste en un listado de las lecturas normalmente asociadas a ruido técnico. Las células con un porcentaje elevado lecturas que mapeen a dichos loci suelen ser ruido y deben eliminarse. Las listas negras de ENCODE actualmente disponibles en
Signacson para genoma humano(hg19 y GRCh38), ratón (mm10), Drosophila (dm3), y C. elegans (ce10).
Téngase en cuenta que las métricas Total number of fragments in peaks, Fraction of fragments in peaks y Ratio reads in genomic blacklist regions pueden obtenerse del output de CellRanger de 10X Genomics. También es posible analizar datasets no provenientes de 10X Genomics gracias a las utilidades que incorpora Signac.
Calculamos la intensidad de la señal de los nucleosomas con el comando NucleosomeSignal():
MB_Signac <- NucleosomeSignal(object = MB_Signac)Ahora podemos ver la distribución de la longitud de las lecturas en células con una nucleosome signal < 4 (bueno) y > 4 (malo). Se aprecia que cada grupo de células tiene una distribución de fragmentos distinta. En las células bien secuenciadas (con un NS < 4), el primer pico de < 100pb se corresponde a los fragmentos sin nucleosomas (estas son lecturas de eucromatina, las cuales deberían mapear en los TSSs), el segundo pico se corresponde a lecturas provenientes de un nucleosoma o mononucleosoma (segundo pico entre 147pb y 147*2pb, flanquea regiones abiertas) y el tercer pico corresponde a lecturas dinucleosomas. Un buen ensayo de ATAC-Seq debe contener regiones libres dre nucleosomas, así como patrones de bandeos de nucleosomas, tal como se ve en nuestro caso en la gráfica NS < 4:
MB_Signac$nucleosome_group <- ifelse(MB_Signac$nucleosome_signal > 4, "NS > 4", "NS < 4")
FragmentHistogram(object = MB_Signac, group.by = "nucleosome_group", region = "chr1-1-10000000")
El enriquecimiento de eventos de integración de la Tn5 en los TSSs es también una importante métrica de control de calidad. El consorcio ENCODE define la puntuación de enriquecimiento de TSSs como el nº de sitios de integración de la Tn5 alrededor del sitio de comienzo de la transcripción normalizado por el nº de sitios de integración de la TN5 en regiones flanqueantes al TSS (hasta 2000pb aguas arriba y aguas abajo, respectivamente).
En Signac calculamos la puntuación de enriquecimiento de los TSSs con el comando TSSEnrichment(). Si usamos el parámetro fast = T, calcularemos sólo la puntuación de enriquecimiento de los TSSs y ahorraremos tiempo y memoria, pero no computaremos la matriz necesaria para visualizar dicho enriquecimiento mediante la función TSSPlot(). El significado del valor del enriquecimiento de TSSs depende del genoma que estemos analizando. De acuerdo al proyecto ENCODE, para el genoma mm10 valores < 10 son malos; valores ente 10 y 15 son aceptables y valores > 15 son ideales. En nuestro caso vemos que el enriquecimiento de nuestras células alrededor del TSS está entre 10 y 15, por lo que el experimento en cuestión es válido.
MB_Signac <- TSSEnrichment(MB_Signac, fast = F)
MB_Signac$high.tss <- ifelse(MB_Signac$TSS.enrichment > 2, "High", "Low")
TSSPlot(MB_Signac, group.by = "high.tss") + NoLegend()
Habiendo calculado ya el TSS enrichment score y el nucleosome signal, procedemos a calcular el ratio de lecturas que mapean en picos y el ratio de lecturas que caen en regiones prohibidas y mostramos todas estas métricas de calidad en gráficos de violín para analizar la calidad de nuestro experimento:
MB_Signac$pct_reads_in_peaks <- MB_Signac$peak_region_fragments / MB_Signac$passed_filters * 100
MB_Signac$blacklist_ratio <- MB_Signac$blacklist_region_fragments / MB_Signac$peak_region_fragments
VlnPlot(object = MB_Signac,
features = c("pct_reads_in_peaks", "peak_region_fragments",
"TSS.enrichment","blacklist_ratio", "nucleosome_signal"),
pt.size = 0.1,
ncol = 5)
Concluimos la sección de QC eliminando las células que no cumplan con los estándares de calidad. Para ello nos quedaremos con aquellas células que cumplan las siguientes condiciones:
Contener entre 3000 y 100000 lecturas que mapeen en picos
Un ratio de lecturas en picos > 40%
Un porcentaje de lecturas en regiones genómicas ruidosas < 0.025%
Una señal de nucleosoma < 4
Una puntuación de enriquecimiento de TSSs > 2
MB_Signac <- subset(
x = MB_Signac,
subset = peak_region_fragments > 3000 &
peak_region_fragments < 100000 &
pct_reads_in_peaks > 40 &
blacklist_ratio < 0.025 &
nucleosome_signal < 4 &
TSS.enrichment > 2)
MB_Signac## An object of class Seurat
## 157203 features across 3512 samples within 1 assay
## Active assay: peaks (157203 features, 0 variable features)4.1 Normalizado y reducción lineal de la dimensionalidad (LSI)
Para el normalizado, Signac realiza una técnica de dos pasos denominada TF-IDF o Term Frequency - Inverse Document Frequency, la cual por un lado normaliza las células (para corregir las posibles diferencias en la profundidad de secuenciado entre células) y por otro normaliza los picos (para dar más importancia a los picos poco comunes).
Si estuviésemos analizando un experimento de scRNA-Seq con Seurat, haríamos una selección de características (=genes) después del normalizado, pero en scATAC-Seq es más difícil realizar dicha selección debido al bajo rango dinámico de estos datos (i.e. bajo ratio señal/ruido). En este caso podemos usar el top x% de picos (=regiones cromosómicas) variables para la reducción de la dimensionalidad, o eliminar los picos que estén presentes en menos de x células (usando para ello el comando FindTopFeatures()). Nosotros usaremos el top 25% (min.cutoff = "q75").
Para la reducción de la dimensionalidad, realizamos una técnica matemática denominada SVD o Singular Value Decomposition sobre la matriz devuelta por el TD-IDF (el SVD se realiza también para generar PCAs) y usando sólo los picos seleccionados en el paso previo. Llamamos LSI o latent semantic indexing a la técnica que combina un paso de TF-IDF seguido de un paso de SVD, y fue usado exitósamente por primera vez en el análisis de scATAC-seq por Cusanovich et al., 2015.
MB_Signac <- RunTFIDF(MB_Signac) # Normalizado
MB_Signac <- FindTopFeatures(MB_Signac, min.cutoff = 'q75') # Selección de características
MB_Signac <- RunSVD(MB_Signac) # Reducción de la dimensionalidadLa primera componente del LSI suele capturar la profundidad del secuenciado (variación técnica) y no la varianza biológica. Si ese es el caso, esta componente de debe eliminar de los próximos análisis. Podemos observar la correlación entre cada componente del LSI y la profundidad de secuenciado con la función DepthCor():
DepthCor(MB_Signac)
Se aprecia que la primera componente del LSI está inversamente muy correlacionada (R ≈ -1) con el nº total de lecturas en células, por lo que omitiremos esta componente en los próximos análisis. Por último, graficamos el LSI, omitiendo la primera componente. El output del mismo recuerda al de un PCA en el sentido de que parece ser una reducción lineal de la dimensionalidad:
DimPlot(object = MB_Signac, label = TRUE, reduction = "lsi", dims = c(2,3)) + NoLegend()
5 Reducción no lineal de la dimensionalidad y clustering
Con las células embebidas en el LSI, podemos agruparlas mediante las técnicas ya aprendidas en los protocolos de Seurat (clustering basado en grafos) y posteriormente visualizarlas con un t-SNE o un UMAP. Este paso es por tanto exactamente igual que en anteriores workflows para scRNA-Seq ya mostrados, con la salvedad de que omitiremos la primera componente del LSI (dims = 2:30).
MB_Signac <- FindNeighbors(MB_Signac, reduction = "lsi", dims = 2:30)
MB_Signac <- FindClusters(MB_Signac, algorithm = 3, resolution = 1.2,
verbose = F)
MB_Signac <- RunUMAP(MB_Signac, reduction = "lsi", dims = 2:30)
DimPlot(MB_Signac, reduction = "umap", label = TRUE) + NoLegend()
6 Matriz de actividad génica
En el UMAP apreciamos la presencia de hasta 21 clusters, los cuales se corresponden con estirpes celulares presentes en el cerebro de ratones adultos. Por desgracia, anotar estos clusters no es tan fácil en ATAC-seq dado que a diferencia de RNA-seq (que se centra en genes), a día de hoy se sabe poco sobre el rol de las regiones genómicas no codificantes.
No obstante, podemos inferir la actividad de cada gen en el genoma a través de la accesibilidad de la cromatina asociada a dichos ORFs y crear así un nuevo ensayo para nuestro objeto seurat derivado del experimento de scATAC-seq. Una manera sencilla de lograrlo es sumar los fragmentos que caigan en la región promotora o en el ORF del gen en cuestión, aunque también es posible integrar en un worflow de Signac la herramienta Cicero, desarrollada por el autor de Monocle3 (otro paquete integrable en workflows de Seurat), Cole Trapnell.
Para generar una matriz de actividad génica extraemos las coordenadas de los genes, las extendemos 2kb aguas arriba (para incluir al promotor) y a continuación contamos para cada célula el nº de lecturas que mapean en dichas regiones usando la función FeatureMatrix(). Estos pasos se realizan automáticamente al invocar al comando GeneActivity(). El objetivo de esta actividad es inferir una matriz de conteos capaz de poder ser usada como input para un nuevo ensayo, como si hubiésemos secuenciado el transcriptoma de estas células:
# Generamos la matriz de actividad génica
matriz_act_gen <- GeneActivity(MB_Signac)
# Añadimos dicha matriz a nuestro objeto `seurat` en forma de un nuevo ensayo de
# scRNA-seq y normalizamos los datos
MB_Signac[["RNA"]] <- CreateAssayObject(counts = matriz_act_gen)
MB_Signac <- NormalizeData(object = MB_Signac, assay = "RNA",
normalization.method = "LogNormalize",
scale.factor = median(MB_Signac$nCount_RNA))Ahora somos capaces de visualizar la actividad de biomarcadores canónicos para guiar nuestra interpretación de los clusters de scATAC-seq, si bien este nuevo experimento putativo de “scRNA-seq” derivado de scATAC-seq será más ruidoso que uno canónico de scRNA-seq, pero nos servirá. Dicho ruido se debe a que al generar la matriz de actividad génica asumimos una correlación perfecta entre accesibilidad del promotor/ORF y expresión génica, lo cual a veces no es así.
DefaultAssay(MB_Signac) <- "RNA"
FeaturePlot(object = MB_Signac,
features = c("Sst", "Pvalb", "Gad2", "Neurod6", "Rorb", "Syt6"),
pt.size = 0.1, max.cutoff = 'q95', ncol = 3)
7 Integración multiómica/multimodal (scRNA-seq)
Para ayudarnos más aún con la intepretación de los clusters de scATAC-seq, podemos clasificar las células basándonos en un experimento real de scRNA-seq del mismo tipo de tejido. Para combinar datos de ATAC-seq y RNA-seq, Signac (y Seurat) usa los métodos de integración multiómicos y transferencia de etiquetas descritos en este paper. El objetivo es identificar patrones de correlación compartidos entre la matriz de actividad génica generada en el apartado anterior y un dataset anotado de scRNA-seq para descubrir tipos celulares compatidos entre ambas ómicas. Este procedimiento se realiza mediante un CCA (procedimiento descrito en otros workflows de Seurat y en la traducción de Seurat 2 a Seurat 3, 3ª práctica), y devuelve una puntuación de clasificación para cada célula de los clusters anotados en el experimento real de scRNA-seq.
Para el experimento real de scRNA-seq, puedes descargar los datos crudos en la página del Allen Institute y preprocesarlos con Seurat, o bien descargar el objeto seurat ya preprocesado. Para preprocesar los datos, se usó este script. Para este ensayo se secuenciaron cerebros de varios ratones con distintos genotipos (concretamente la corteza visual primaria).
# Cargamos el ensayo de scRNA-seq (objeto tipo `seurat`) del Instituto Allen
allen_rna <- readRDS("./archivos_accesorios/allen_brain.rds")
# Ampliamos el nº de genes diferencialmente expresados a 5000
allen_rna <- FindVariableFeatures(object = allen_rna, nfeatures = 5000)
# Detectamos células ancla en nuestros datasets mediante un CCA/RCPA
transfer.anchors <- FindTransferAnchors(reference = allen_rna, query = MB_Signac,
reduction = "rpca", dims = 1:40)
# Nota: ya que usamos datasets pesados, no divergentes (mirar UMAP) y de la
# misma especie, es computacionalmente más eficiente usar un RCPA en vez de un
# CCA
# Predecimos los tipos celulares presentes en el ensayo de scATAC-seq usando
# como input las células ancla y el dataset del Instituto Allen como referencia
predicted.labels <- TransferData(
anchorset = transfer.anchors,
refdata = allen_rna$subclass,
weight.reduction = MB_Signac[["lsi"]], # reducción de dim. del objeto `seurat` original
dims = 2:30)
# Añadimos a los metadatos del objeto principal `MB_Signac` las predicciones
# sobre los tipos celulares y sustituimos las etiquetas antiguas por estas
# nuevas
MB_Signac <- AddMetaData(object = MB_Signac, metadata = predicted.labels)
for(i in levels(MB_Signac)) {
cells_to_relabel <- WhichCells(MB_Signac, idents = i)
newid <- names(sort(table(MB_Signac$predicted.id[cells_to_relabel]),decreasing=TRUE))[1]
Idents(MB_Signac, cells = cells_to_relabel) <- newid
}plot1 <- DimPlot(allen_rna, group.by = "subclass", label = TRUE) + NoLegend() + ggtitle("scRNA-seq (Allen Inst., anotado)")
plot2 <- DimPlot(MB_Signac, group.by = "predicted.id", label = TRUE) + NoLegend() + ggtitle("scATAC-seq (predicción)")
plot3 <- DimPlot(MB_Signac, label = T, group.by = "ident") + NoLegend() + ggtitle("scATAC-seq (anot. Allen)")
plot1 + plot2 + plot3
Se aprecia que aunque las visualizaciones de scRNA-seq y scATAC-seq no son exactamente iguales, la agrupación de los clusters es consistente entre ambos ensayos (por ejemplo, las neuronas productoras de somatostatina Sst y las interneuronas productoras de parvalbúmina Pvalb forman clusters adyacentes en ambos experimentos).
8 Detección de biomarcadores de clusters (regiones cromosómicas diferencialmente accesibles)
Si en estudios de transcriptómica analizamos los genes diferencialmente transcribidos, es lógico estudiar en ATAC-seq aquellas regiones del genoma que son diferencialmente accesibles a la transposasa Tn5. Para estudiar la apertura diferencial de la cromatina se usan regresiones logísticas, tal como recomienda Ntranos et al. 2018, y se añade el nº total de lecturas en forma de variable latente para mitigar el efecto negativo en los resultados de tener librerías/muestras con distintas profundidades de secuenciado.
Para datos dispersos como es scATAC-seq, es necesario ajustar el parámetro min.pct del comando FindMarkers() con valores más bajos, ya que el valor que usa por defecto (0.1) está pensado para datos de scRNA-seq. En nuestro caso detectaremos regiones diferencialmente accesibles entre las neuronas intratelencefálicas de capas corticales 2/3 y las neuronas de capas 4-6 (siendo todas ellas neuronas excitatorias):
# Volvemos a trabajar con picos en vez de actividad génica
DefaultAssay(MB_Signac) <- "peaks"
# `L5 IT` está en tipos celulares predichos
Idents(MB_Signac) <- MB_Signac$predicted.id
da_peaks <- FindMarkers(
object = MB_Signac,
ident.1 = c("L2/3 IT"),
ident.2 = c("L4", "L5 IT", "L6 IT"), # IT = IntraTelencefálico; L = Layer
min.pct = 0.05, # 0.1 para scRNA-seq; 0.05 para scATAC-seq
test.use = "LR",
latent.vars = "peak_region_fragments"
)
head(da_peaks)Visualizamos los resultados del test de accesibilidad diferencial en un violinplot y sobre la proyección del UMAP:
plot1 <- VlnPlot(object = MB_Signac, features = rownames(da_peaks)[1],
pt.size = 0.1, idents = c("L4","L5 IT","L2/3 IT", "L6 IT")) +
ylab("Insertion/Aperture level") + NoLegend()
plot2 <- FeaturePlot(object = MB_Signac, features = rownames(da_peaks)[1],
pt.size = 0.1, max.cutoff = "q95")
plot1 | plot2
En vista de los resultados podemos ver que la región del cromosoma 4 comprendida entre las bases 86.523.678-86.525.285 está diferencialmente abierta (i.e. expresada) en las neuronas piramidales de capas corticales 2/3 respecto de las de capas 4-6. Dicha región abarca el gen Fam154a/Saxo1:
CoveragePlot(object = MB_Signac, region = "chr4-86523678-86525285",
extend.upstream = 1000, extend.downstream = 1000,
ncol = 1, idents = c("L2/3 IT", "L4", "L5 IT", "L6 IT"))
Puede ser difícil interpretar las coordenadas de los picos de integración de transposones (al fin y al cabo, no es común saberse de memoria las bases exactas que ocupa un gen). Por suerte, Signac cuenta con una utilidad para descubrir genes cercanos a (o que mapeen en) las regiones cromosómicas que le indiquemos. Esta utilidad es el comando ClosestFeature(). Con dicho comando podemos comprobar que, efectivamente, el gen más cercano a la región chr4-86523678-86525285 es Fam154a/Saxo1:
# Regiones cromosómicas abiertas con un log2FC > |0.25| en neuronas de capas
# L2/3 y L4-5
open_l23 <- rownames(da_peaks[da_peaks$avg_log2FC > 0.25, ])
open_l456 <- rownames(da_peaks[da_peaks$avg_log2FC < -0.25, ])
# Genes que mapean en dichas regiones
closest_l23 <- ClosestFeature(MB_Signac, open_l23)
closest_l456 <- ClosestFeature(MB_Signac, open_l456)
head(closest_l23)head(closest_l456)9 Visualizar regiones cromosómicas
Tal como veíamos 2 párrafos antes, podemos visualizar el nº de integraciones de la Tn5 en el genoma mediante la función CoveragePlot(). Esta función acepta loci tanto en formato coordenadas cromosómicas (i.e. chr4-86523678-86525285) como en formato nombre del gen (i.e. Neurod6). Cada tipo celular presente en el gráfico resultante muestra el nº medio de sitios de integración en todas las células del tipo en cuestión. También podemos graficar clusters o cualquier otro metadato, en lugar de los tipos celulares:
# Establecemos el orden en el gráfico de las células a mostrar:
levels(MB_Signac) <- c("L2/3 IT", "L4", "L5 IT", "L5 PT", "L6b", "L6 CT", "L6 IT", "NP", "Sst" ,"Pvalb" ,"Vip" ,"Lamp5" , "Oligo", "Astro", "Endo", "VLMC", "Macrophage", "Sncg", "Peri")
# Coordenadas cromosómicas
CoveragePlot(
object = MB_Signac,
region = "chr4-86523678-86525285",
extend.upstream = 1000,
extend.downstream = 1000,
ncol = 1)
# Nombres de genes
CoveragePlot(
object = MB_Signac,
region = c("Neurod6", "Gad2"),
extend.upstream = 1000,
extend.downstream = 1000,
ncol = 1,
idents = levels(MB_Signac)[1:5])
# Apertura media en todas las células, se aprecia la caja TATA del gen Sst
CoveragePlot(
object = MB_Signac,
region = "Sst",
extend.upstream = 1000,
extend.downstream = 1000,
ncol = 1,
group.by = "orig.ident")
También se puede generar mediante la función CoverageBrowser() una versión interactiva de estos gráficos en la cual se puede explorar el genoma y ajustar al vuelo los parámetros de graficado. El botón “Save plot” guarda el gráfico en un objeto ggplot y se carga en el ambiente al cerrar la página interactiva.
10 Addendum: Qué hacer si el dataset no se generó con CellRanger
Como ya se comentó al principio del informe, el programa CellRanger de 10X Genomics genera varias métricas de control de calidad por célula, así como también crea la matriz de picos/célula y el archivo del índice de los fragmentos/lecturas. El dataset empleado para este informe se generó con dicho software y se pudo emplear tal cual, pero en caso de generar el dataset con otro software, Signac cuenta con las herramientas necesarias para generar esos archivos.
10.1 Generar la matriz de picos/célula
Para generar la matriz de picos/célula se puede usar la función FeatureMatrix(), la cual devuelve como es de esperar una matriz de conteos de picos (filas) en celulas (columnas):
# peak_ranges debe ser un objeto de tipo `GRanges`, o sea un conjunto de
# coordenadas cromosómicas que abarque el conjunto de picos a cuantificar por
# célula
peak_matrix <- FeatureMatrix(
fragments = Fragments(objeto_seurat),
features = peak_ranges)10.2 Calcular el ratio de lecturas que mapean en picos
Dado un archivo de fragmentos, la función FRiP() es capaz de computar el porcentaje de lecturas que caen en picos por cada célula y la almacena en el hueco de metadatos del objeto de tipo seurat empleado como input:
# Contamos los fragmentos totales en nuestro ensayo
total_fragments <- CountFragments(
fragments = "E:/Git/R/Signac MB/archivos_accesorios/atac_v1_adult_brain_fresh_5k_fragments.tsv.gz")
objeto_seurat$fragments <- total_fragments[colnames(objeto_seurat), "frequency_count"]
objeto_seurat <- FRiP(
object = objeto_seurat,
assay = "peaks",
total.fragments = "fragments")10.3 Calcular el porcentaje de lecturas que mapean en regiones ruidosas (blacklists)
El paquete Signac incluye para conveniencia del usuario las coordenadas cromosómicas de regiones ruidosas para varios genomas (a fecha de redacción incluye las versiones del genoma humano hg19 y hg38, genoma murino versiones mm9 y mm10, genoma de Caenorhabditis elegans versiones ce10 y ce11 y genoma de Drosophila melanogaster versiones dm3 y dm6). Puedes emplear la función FractionCountsInRegion() y una blacklist para calcular el porcentaje de lecturas que caen en estas regiones problemáticas.
objeto_seurat$blacklist_fraction <- FractionCountsInRegion(
object = objeto_seurat,
assay = "peaks",
regions = blacklist_mm10) # blacklist_hg19 para datasets humanosEstas blacklists fueron provistas por el consorcio ENCODE, por lo que si las empleas, es recomendable citar el artículo correspondiente.
11 Bibliografía
Stuart et al. Multimodal single-cell chromatin analysis with Signac. bioRxiv (2020).
ENCODE portal (PMID: 29126249; PMCID: PMC5753278)
Amemiya, H.M., Kundaje, A. & Boyle, A.P. The ENCODE Blacklist: Identification of Problematic Regions of the Genome. Sci Rep 9, 9354 (2019). https://doi.org/10.1038/s41598-019-45839-z
Srivatsan, Sanjay & McFaline-Figueroa, José & Ramani, Vijay & Saunders, Lauren & Cao, Junyue & Packer, Jonathan & Pliner, Hannah & Jackson, Dana & Daza, Riza & Christiansen, Lena & Zhang, Fan & Steemers, Frank & Shendure, Jay & Trapnell, Cole. (2019). Massively multiplex chemical transcriptomics at single cell resolution. Science. 367. eaax6234. 10.1126/science.aax6234.
Hannah A. Pliner, Jonathan S. Packer, José L. McFaline-Figueroa, Darren A. Cusanovich, Riza M. Daza, Sanjay Srivatsan, Xiaojie Qiu, Dana Jackson, Anna Minkina, Andrew C. Adey, Frank J. Steemers, Jay Shendure, Cole Trapnell. Molecular Cell 71, 1–14 2018
Stuart T, Butler A, Hoffman P, Hafemeister C, Papalexi E, Mauck WM 3rd, Hao Y, Stoeckius M, Smibert P, Satija R. Comprehensive Integration of Single-Cell Data. Cell. 2019 Jun 13;177(7):1888-1902.e21. doi: 10.1016/j.cell.2019.05.031. Epub 2019 Jun 6. PMID: 31178118; PMCID: PMC6687398.
Billeh YN, Cai B, Gratiy SL, Dai K, Iyer R, Gouwens NW, Abbasi-Asl R, Jia X, Siegle JH, Olsen SR, Koch C, Mihalas S, Arkhipov A. Systematic Integration of Structural and Functional Data into Multi-scale Models of Mouse Primary Visual Cortex. Neuron. 2020 May 6;106(3):388-403.e18. doi: 10.1016/j.neuron.2020.01.040. Epub 2020 Mar 5. PMID: 32142648.
Ntranos, V., Yi, L., Melsted, P. et al. A discriminative learning approach to differential expression analysis for single-cell RNA-seq. Nat Methods 16, 163–166 (2019). https://doi.org/10.1038/s41592-018-0303-9
12 sessionInfo()
Click para mostrar
## R version 4.1.0 (2021-05-18)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19043)
##
## 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] stats4 parallel stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] patchwork_1.1.1 ggplot2_3.3.5
## [3] EnsDb.Mmusculus.v79_2.99.0 ensembldb_2.16.4
## [5] AnnotationFilter_1.16.0 GenomicFeatures_1.44.2
## [7] AnnotationDbi_1.54.1 Biobase_2.52.0
## [9] GenomicRanges_1.44.0 GenomeInfoDb_1.28.4
## [11] IRanges_2.26.0 S4Vectors_0.30.0
## [13] BiocGenerics_0.38.0 SeuratObject_4.0.2
## [15] Seurat_4.0.4 Signac_1.3.0.9001
## [17] future_1.22.1
##
## loaded via a namespace (and not attached):
## [1] utf8_1.2.2 reticulate_1.20
## [3] tidyselect_1.1.1 RSQLite_2.2.8
## [5] htmlwidgets_1.5.4 grid_4.1.0
## [7] docopt_0.7.1 BiocParallel_1.26.2
## [9] Rtsne_0.15 munsell_0.5.0
## [11] codetools_0.2-18 ica_1.0-2
## [13] miniUI_0.1.1.1 withr_2.4.2
## [15] colorspace_2.0-2 filelock_1.0.2
## [17] highr_0.9 knitr_1.34
## [19] rstudioapi_0.13 ROCR_1.0-11
## [21] tensor_1.5 listenv_0.8.0
## [23] labeling_0.4.2 MatrixGenerics_1.4.3
## [25] slam_0.1-48 GenomeInfoDbData_1.2.6
## [27] polyclip_1.10-0 bit64_4.0.5
## [29] farver_2.1.0 parallelly_1.28.1
## [31] vctrs_0.3.8 generics_0.1.0
## [33] xfun_0.25 biovizBase_1.40.0
## [35] BiocFileCache_2.0.0 lsa_0.73.2
## [37] ggseqlogo_0.1 R6_2.5.1
## [39] hdf5r_1.3.4 bitops_1.0-7
## [41] spatstat.utils_2.2-0 cachem_1.0.6
## [43] DelayedArray_0.18.0 assertthat_0.2.1
## [45] promises_1.2.0.1 BiocIO_1.2.0
## [47] scales_1.1.1 nnet_7.3-16
## [49] gtable_0.3.0 globals_0.14.0
## [51] goftest_1.2-2 rlang_0.4.11
## [53] RcppRoll_0.3.0 splines_4.1.0
## [55] rtracklayer_1.52.1 lazyeval_0.2.2
## [57] dichromat_2.0-0 checkmate_2.0.0
## [59] spatstat.geom_2.2-2 yaml_2.2.1
## [61] reshape2_1.4.4 abind_1.4-5
## [63] backports_1.2.1 httpuv_1.6.3
## [65] Hmisc_4.5-0 tools_4.1.0
## [67] ellipsis_0.3.2 spatstat.core_2.3-0
## [69] jquerylib_0.1.4 RColorBrewer_1.1-2
## [71] ggridges_0.5.3 Rcpp_1.0.7
## [73] plyr_1.8.6 base64enc_0.1-3
## [75] progress_1.2.2 zlibbioc_1.38.0
## [77] purrr_0.3.4 RCurl_1.98-1.4
## [79] prettyunits_1.1.1 rpart_4.1-15
## [81] deldir_0.2-10 pbapply_1.4-3
## [83] cowplot_1.1.1 zoo_1.8-9
## [85] SummarizedExperiment_1.22.0 ggrepel_0.9.1
## [87] cluster_2.1.2 magrittr_2.0.1
## [89] RSpectra_0.16-0 data.table_1.14.0
## [91] scattermore_0.7 lmtest_0.9-38
## [93] RANN_2.6.1 SnowballC_0.7.0
## [95] ProtGenerics_1.24.0 fitdistrplus_1.1-5
## [97] matrixStats_0.60.1 hms_1.1.0
## [99] mime_0.11 evaluate_0.14
## [101] xtable_1.8-4 XML_3.99-0.7
## [103] jpeg_0.1-9 sparsesvd_0.2
## [105] gridExtra_2.3 compiler_4.1.0
## [107] biomaRt_2.48.3 tibble_3.1.4
## [109] KernSmooth_2.23-20 crayon_1.4.1
## [111] htmltools_0.5.2 mgcv_1.8-36
## [113] later_1.3.0 Formula_1.2-4
## [115] tidyr_1.1.3 DBI_1.1.1
## [117] tweenr_1.0.2 dbplyr_2.1.1
## [119] MASS_7.3-54 rappdirs_0.3.3
## [121] Matrix_1.3-4 igraph_1.2.6
## [123] pkgconfig_2.0.3 GenomicAlignments_1.28.0
## [125] foreign_0.8-81 plotly_4.9.4.1
## [127] spatstat.sparse_2.0-0 xml2_1.3.2
## [129] bslib_0.3.0 XVector_0.32.0
## [131] VariantAnnotation_1.38.0 stringr_1.4.0
## [133] digest_0.6.27 sctransform_0.3.2
## [135] RcppAnnoy_0.0.19 spatstat.data_2.1-0
## [137] Biostrings_2.60.2 rmarkdown_2.10
## [139] leiden_0.3.9 fastmatch_1.1-3
## [141] htmlTable_2.2.1 uwot_0.1.10
## [143] restfulr_0.0.13 curl_4.3.2
## [145] shiny_1.6.0 Rsamtools_2.8.0
## [147] rjson_0.2.20 lifecycle_1.0.0
## [149] nlme_3.1-153 jsonlite_1.7.2
## [151] BSgenome_1.60.0 viridisLite_0.4.0
## [153] fansi_0.5.0 pillar_1.6.2
## [155] lattice_0.20-44 KEGGREST_1.32.0
## [157] fastmap_1.1.0 httr_1.4.2
## [159] survival_3.2-13 glue_1.4.2
## [161] qlcMatrix_0.9.7 png_0.1-7
## [163] bit_4.0.4 ggforce_0.3.3
## [165] stringi_1.7.4 sass_0.4.0
## [167] blob_1.2.2 latticeExtra_0.6-29
## [169] memoise_2.0.0 dplyr_1.0.7
## [171] irlba_2.3.3 future.apply_1.8.1