ANÁLISIS DE DATOS scRNA-Seq CON BIOCONDUCTOR
Cristina Gil Martínez
Febrero, 2020
NGS y RNA-Seq
El término Next-Generation Sequencing o NGS es un término utilizado para hacer referencia a técnicas modernas de secuenciación masiva de DNA o RNA, es decir, de determinación de la secuencia de nucleótidos de dichos materiales genéticos. Una de dichas técnicas es la conocida como RNA-Seq o “secuenciación de RNA”, que se puede realizar sobre poblaciones de células o sobre células individuales (scRNA-Seq).
En el proceso de expresión genética, determinadas secuencias de DNA se transcriben a moléculas de RNA mensajero o mRNA que posteriormente son traducidas a proteínas:
Todas las células de un organismo contienen el mismo material genético, pero cada tipo celular expresa solo una parte específica y diferente a otros tipos celulares. En función del tipo celular, contexto o proceso biológico llevado a cabo, se produce una expresión diferencial del conjunto de genes, es decir, genes que se expresan en unas células y no en otras.
En este contexto, el RNA-Seq se puede aplicar para revelar la presencia y cantidad de mRNA en una muestra biológica. Se trata de una técnica de análisis de expresión diferencial, es decir, nos permite analizar cambios en el transcriptoma (conjunto de todas las moléculas de RNA presentes en una célula) o expresión de genes. De forma general, esta técnica se basa en fragmentar los transcritos de mRNA y revertirlos a DNA complementario (cDNA), los cuales se secuencian con la ayuda de adaptadores de secuencia, utilizando tecnologías como por ejemplo Illumina o SOLiD sequencing. Estas secuencias se alinean posteriormente con un genoma de referencia para determinar de qué genes proceden y el nivel de expresión por el conteo del número de secuencias detectadas:
MATRIZ DE CONTEO
La técnica RNA-Seq nos proporciona el número de secuencias alineadas con los genes o zonas de interés de un genoma de referencia, lo que nos da información sobre el nivel de expresión de cada gen en las células analizadas. Así, con un conjunto de \(\small{m}\) experimentos obtendremos una matriz de conteo con cada columna asociado a los datos de una muestra de células o célula individual (dependiendo del tipo de RNA-Seq aplicado), y cada fila asociada a los datos de un gen específico \(\small{g}\):
El método puede ser propenso a errores, por lo que suelen añadirse controles para corregirlos, como los RNA spike-ins: moléculas de RNA de secuencia conocida diseñados para hibridar o unirse a sondas (fragmentos pequeños de DNA de secuencia también conocida) presentes en las muestras. Los spike-ins se añaden con una concentración determinada a las muestras a secuenciar para la calibración de las medidas obtenidas de hibridación mRNA-cDNA. En este escenario diferenciaríamos el RNA endógeno proveniente de las células, y el RNA proveniente de los spike-ins. En ausencia de spike-ins también puede hacerse uso de la proporción de secuencias asociadas a genes mitocondriales.
Normalización
El conteo del número de secuencias obtenidas de cada gen no nos aporta directamente el nivel de expresión de los mismos, ya que este número puede verse afectado por factores como:
Tamaño de la genoteca o profundidad de secuenciación: para cada muestra, el total de lecturas varía.
Composición del RNA: para cada muestra obtenemos datos de abundancia relativa de cada gen sin disponer de la información de la cantidad total de RNA por célula. Si no se realizan las correcciones apropiadas, pueden aparecer genes infraregulados donde otros aparecen con altos grados de expresión.
Longitud de los genes: los genes más cortos pueden presentar menos lecturas que los más largos.
Estos factores hacen necesario normalizar los datos de conteos para minimizar el ruido técnico introducido durante el proceso de secuenciación y así poder comparar correctamente las distintas muestras o experimentos. Al normalizar, pues, buscamos convertir a una misma escala los datos de todas las muestras (ya que una genoteca de mayor tamaño tiene más probabilidad de tener genes de expresión diferencial respecto a otra).
Uno de los métodos de normalización comunes es la transformación \(\small{log_2}\), que corrige y equilibra los conteos debido a las diferencias en tamaños de las genotecas. Esto se consigue dividiendo primero cada valor de conteo \(\small{y_{gm}}\) entre un factor de tamaño de la genoteca (tamaño normalizado de la genoteca de la célula correspondiente, cuya media conjunta para todas las muestras sea 1) que representa una estimación del bias relativo de dicha muestra, por lo que la división de cada conteo por este valor pretende eliminar dicho bias. Por último se calcula el \(\small{log_2}\) del valor obtenido tras la división. Este método interpreta cualquier diferencia en la cantidad total de RNA entre muestras como parte del bias y lo elimina. En R puede aplicarse con la función librarySizeFactors().
En caso de que el contenido total de RNA de cada muestra sea de interés (por ejemplo en estudios de actividad del ciclo celular), se puede realizar una normalización en base un factor de contenido de spike-ins, también con media 1 en el conjunto de muestras, y calculados sobre la cantidad de spike-ins, y no sobre el tamaño de la genoteca de cada muestra. Esta normalización asegura que el efecto del RNA total sobre la expresión entre muestras no se elimine. En R puede aplicarse con la función computeSpikeFactors().
La elección de un tipo u otro de normalización, no siendo estos dos métodos los únicos disponibles, dependerá de la hipótesis biológica de partida.
BIOCONDUCTOR
Bioconductor es un paquete desarrollado para proporcionar herramientas para el análisis de datos ómicos en R. Para instalar Bioconductor y paquetes asociados adicionales, usar:
# Instalar Bioconductor version 3.10:
install.packages("BiocManager")
BiocManager::install(version = "3.10")
# Instalar paquetes específicos de Bioconductor
BiocManager::install("paquete")
# Actualizar paquetes de Bioconductor
BiocManager::install()OBJETO DE CLASE SingleCellExperiment
El objeto de clase SingleCellExperiment presenta un esquema como el mostrado en la siguiente imagen:
Fuente: https://osca.bioconductor.org/data-infrastructure.html
Este tipo de objeto dispone de varias entidades, cada una con un formato propio para los datos que almacenan:
assays: almacena los datos principales del análisis en formato matriz dentro de una lista. Se almacenan las características (ej.: genes) por filas y las muestras (células) por columnas.rawData: DataFrame con información sobre metadatos de las características (genes) correspondientes a las filas deassays.colData: incluye metadatos de las muestras (células) de cada experimento, en formato DataFrame con células por filas e información de las mismas por columnas.reducedDims: almacena el resultado, en formato matriz o lista de matrices, de técnicas de reducción de dimensionalidad sobre los datos principales enassays.
A continuación se muestra un ejemplo simple de creación de un objeto de clase SingleCellExperiment con datos artificiales para cada una de las entidades comentadas.
library(SingleCellExperiment)assays
Para completar esta entidad, el procedimiento constaría en importar los datos RNA-Seq para crear el objeto SingleCellExperiment. Para este ejemplo crearemos los datos artificialmente a partir de una distribución Poisson:
# Creamos datos de conteo de secuencias de 10 genes en 3 experimentcon con 3
# lineas celulares
set.seed(565)
matriz_conteo <- matrix(c(rpois(10, 15), rpois(10, 15), rpois(10, 35)),
byrow = FALSE,
nrow = 10,
ncol = 3,
dimnames = list(paste0("gen_", 1:10),
paste0("cel_", c("a", "b", "c"))))
# Creamos el objeto "sce", almacenando la matriz asignada a un nombre (counts)
# en una lista
sce <- SingleCellExperiment(assays = list(counts = matriz_conteo))Si a la hora de inspeccionar el contenido queremos acceder a una matriz en particular, haremos uso de la función assay():
assay(sce, "counts")## cel_a cel_b cel_c
## gen_1 16 10 32
## gen_2 16 17 28
## gen_3 17 18 33
## gen_4 12 14 43
## gen_5 21 12 32
## gen_6 8 20 33
## gen_7 21 23 37
## gen_8 15 20 28
## gen_9 16 22 27
## gen_10 14 26 32La entidad assays puede contener varias versiones de los datos principales, lo cual es útil si además de los datos originales, queremos mantener también una versión normalizada o transformada de los mismos. Como ejemplo vamos a almacenar una versión normalizada aplicando la transformación \(\small{log_2}\). En concreto, toma cada valor de conteo y lo divide por el tamaño escalado (con media 1 sobre todas las muestras de células) de la genoteca correspondiente (valor obtenido por la función librarySizeFactors()), sumando un pseudoconteo (normalmente +1) para evitar valores no definidos en 0, y aplicando finalmente el \(\small{log_2}\):
library(scater)
# Factores de normalizacion por el que se dividira cada conteo antes de sumarle
# 1 y aplicar el log2
librarySizeFactors(sce)## cel_a cel_b cel_c
## 0.7058824 0.8235294 1.4705882summary(librarySizeFactors(sce)) # (media 1 sobre todas las muestras)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.7059 0.7647 0.8235 1.0000 1.1471 1.4706# Transformacion log2
sce <- logNormCounts(sce)
assays(sce)## List of length 2
## names(2): counts logcounts# Datos transformados
assay(sce, "logcounts")## cel_a cel_b cel_c
## gen_1 4.564785 3.716207 4.508429
## gen_2 4.564785 4.435819 4.324811
## gen_3 4.648657 4.514573 4.550901
## gen_4 4.169925 4.169925 4.918386
## gen_5 4.942515 3.960829 4.508429
## gen_6 3.624491 4.660251 4.550901
## gen_7 4.942515 4.854423 4.709291
## gen_8 4.475733 4.660251 4.324811
## gen_9 4.564785 4.792558 4.275007
## gen_10 4.380822 5.025535 4.508429Por ejemplo, el conteo del gen1 para la célula a (\(\small{y_{1a}}\)) se ha obtenido de la siguiente forma:
colData
Siguiendo el ejemplo anterior, vamos a añadir metadatos para cada una de las células de cada experimento. Supongamos que la célula a pertenece a la tanda o bloque (batch) 1, y las células b y c a una tanda o bloque 2:
# Creamos los metadatos
metadatos <- data.frame(batch = c(1, 2, 2),
row.names = paste0("cel_", c("a", "b", "c")))
# Asignamos los metadatos al objeto SingleCellExperiment
colData(sce) <- DataFrame(metadatos)
sce## class: SingleCellExperiment
## dim: 10 3
## metadata(0):
## assays(2): counts logcounts
## rownames(10): gen_1 gen_2 ... gen_9 gen_10
## rowData names(0):
## colnames(3): cel_a cel_b cel_c
## colData names(1): batch
## reducedDimNames(0):
## spikeNames(0):
## altExpNames(0):La función addPerCellQC() del paquete scater nos permite añadir varias métricas de control sobre cada célula/experimento, como el tamaño de la genoteca o número total de secuencias de genes por célula (sum), número de genes por célula con conteos por encima del límite de detección, por defecto 0 (detected), etc:
sce <- addPerCellQC(sce)
colData(sce)## DataFrame with 3 rows and 8 columns
## batch sum detected percent_top_50 percent_top_100
## <numeric> <integer> <integer> <numeric> <numeric>
## cel_a 1 156 10 100 100
## cel_b 2 182 10 100 100
## cel_c 2 325 10 100 100
## percent_top_200 percent_top_500 total
## <numeric> <numeric> <integer>
## cel_a 100 100 156
## cel_b 100 100 182
## cel_c 100 100 325También tenemos la posibilidad de filtrar subgrupos de datos, por ejemplo los experimentos de la tanda 2 (células b y c):
sce[ , sce$batch == 2]## class: SingleCellExperiment
## dim: 10 2
## metadata(0):
## assays(2): counts logcounts
## rownames(10): gen_1 gen_2 ... gen_9 gen_10
## rowData names(0):
## colnames(2): cel_b cel_c
## colData names(8): batch sum ... percent_top_500 total
## reducedDimNames(0):
## spikeNames(0):
## altExpNames(0):rowData
Con la función addPerFeatureQC() podemos añadir campos con información sobre cada gen en la entidad rowData del objeto:
sce <- addPerFeatureQC(sce)
rowData(sce)## DataFrame with 10 rows and 2 columns
## mean detected
## <numeric> <numeric>
## gen_1 19.3333333333333 100
## gen_2 20.3333333333333 100
## gen_3 22.6666666666667 100
## gen_4 23 100
## gen_5 21.6666666666667 100
## gen_6 20.3333333333333 100
## gen_7 27 100
## gen_8 21 100
## gen_9 21.6666666666667 100
## gen_10 24 100reducedDims
Sobre esta entidad podemos almacenar los datos reducidos dimensionalmente con técnicas como PCA, t-SNE y UMAP. Por ejemplo, para obtener las matrices con el resultado de estas tres técnicas, podemos aplicar las funciones wrapper runPCA(), runTSNE() y runUMAP() del paquete scater:
# PCA: devuelve los vectores de loadings y el porcentaje de varianza explicada
sce <- runPCA(sce)
reducedDim(sce, "PCA")## PC1 PC2
## cel_a -0.8486629 -0.4276706
## cel_b 0.9629685 -0.2915349
## cel_c -0.1143056 0.7192055
## attr(,"percentVar")
## [1] 67.89737 32.10263# t-SNE
sce <- runTSNE(sce, perplexity = 0.15)## Perplexity should be lower than K!reducedDim(sce, "TSNE")## [,1] [,2]
## cel_a 5479.042 -1551.352
## cel_b -1394.043 5520.570
## cel_c -4084.999 -3969.218# UMAP
sce <- runUMAP(sce, n_neighbors = 2)
reducedDim(sce, "UMAP")## [,1] [,2]
## cel_a -0.2927511 -0.6327727
## cel_b 0.5187562 0.4168413
## cel_c -0.2260051 0.2159314
## attr(,"scaled:center")
## [1] 4.795692 11.756641# Objeto completado
sce## class: SingleCellExperiment
## dim: 10 3
## metadata(0):
## assays(2): counts logcounts
## rownames(10): gen_1 gen_2 ... gen_9 gen_10
## rowData names(2): mean detected
## colnames(3): cel_a cel_b cel_c
## colData names(8): batch sum ... percent_top_500 total
## reducedDimNames(3): PCA TSNE UMAP
## spikeNames(0):
## altExpNames(0):Nota: en el apartado metadata del objeto sce podemos almacenar otros metadatos de tipo variado de los que dispongamos.
altExpNames
Entidad donde puede almacenarse datos sobre los controles spike-ins (ERCC, etc.).
CASO DE ESTUDIO
En este ejemplo se muestra el análisis de datos scRNA-Seq sobre la línea celular mieloide de ratón 416b (Lun et al. 2017), donde intervienen dos fenotipos: wild type y expresión inducida del oncogen CBFB-MYH11. El experimento se ha llevado a cabo en dos bloques (20160113 y 20160325) y se ha hecho uso de dos sets de spike-ins (ERCC y SIRV).
Estos datos podemos encontrarlos en el paquete scRNAseq:
library(scRNAseq)
# Datos scRNA-Seq de la linea celular de raton 416b
sce_416b <- LunSpikeInData(which = "416b")
sce_416b## class: SingleCellExperiment
## dim: 46604 192
## metadata(0):
## assays(1): counts
## rownames(46604): ENSMUSG00000102693 ENSMUSG00000064842 ...
## ENSMUSG00000095742 CBFB-MYH11-mcherry
## rowData names(1): Length
## colnames(192): SLX-9555.N701_S502.C89V9ANXX.s_1.r_1
## SLX-9555.N701_S503.C89V9ANXX.s_1.r_1 ...
## SLX-11312.N712_S508.H5H5YBBXX.s_8.r_1
## SLX-11312.N712_S517.H5H5YBBXX.s_8.r_1
## colData names(9): Source Name cell line ... spike-in addition block
## reducedDimNames(0):
## spikeNames(0):
## altExpNames(2): ERCC SIRVEl set contiene 46604 datos de conteos (assays(1): counts) sobre 192 células, todas provenientes de la misma línea celular. En rownames se recogen los códigos Ensembl de cada gen, en colnames se almacena el código de cada célula, y en colData names información diversa sobre cada una (tipo, nombre, etc.), además del bloque experimental al que pertenecen, que convertiremos a factor para los análisis posteriores:
# Convertimos a factor los datos de "block" en colData para analisis posteriores
sce_416b$block <- factor(sce_416b$block)A partir de aquí, los pasos a seguir en el análisis serán:
1. Calcular métricas de control de calidad (número de conteos y características por célula, proporción de lecturas mitocondriales, etc.) para eliminar células con baja calidad.
2. Convertir los conteos en valores normalizados de expresión, para poder realizar comparaciones entre células. Se aplica una transformación a los datos, típicamente una transformación log para ajustar la relación media-varianza.
3. Selección de características para mantener los genes con mayor variabilidad, y reducir el ruido aportado por genes irrelevantes.
4. Aplicar técnicas de reducción de la dimensionalidad, con los que compactar los datos y obtener visualizaciones.
5. Obtener agrupaciones (clustering) en función de la similitud de los valores normalizados de expresión para obtener posibles estados biológicos diferenciados.
Anotación
Antes de continuar, vamos a proceder con la anotación genómica, es decir, la identificación de los genes expresados en cada muestra analizada y la obtención de información disponible de dichos genes en bases de datos públicas. Para ello utilizaremos el AnnotationHub:
library(AnnotationHub)
# Recuperamos de la base de datos la informacion referente al raton domestico
hub <- subset(AnnotationHub(), species == "Mus musculus")
hub## AnnotationHub with 2978 records
## # snapshotDate(): 2019-10-29
## # $dataprovider: Ensembl, Haemcode, UCSC, BioMart, Gencode, KEGG, ftp://ftp....
## # $species: Mus musculus
## # $rdataclass: GRanges, TwoBitFile, BigWigFile, ChainFile, Rle, EnsDb, list,...
## # additional mcols(): taxonomyid, genome, description,
## # coordinate_1_based, maintainer, rdatadateadded, preparerclass, tags,
## # rdatapath, sourceurl, sourcetype
## # retrieve records with, e.g., 'object[["AH6057"]]'
##
## title
## AH6057 | Contigs
## AH6058 | Assembly
## AH6059 | Gap
## AH6060 | GRC Incident
## AH6061 | UCSC Genes
## ... ...
## AH79286 | Mus_musculus_pwkphj.PWK_PhJ_v1.99.chr.gtf
## AH79287 | Mus_musculus_pwkphj.PWK_PhJ_v1.99.gtf
## AH79288 | Mus_musculus_wsbeij.WSB_EiJ_v1.99.abinitio.gtf
## AH79289 | Mus_musculus_wsbeij.WSB_EiJ_v1.99.chr.gtf
## AH79290 | Mus_musculus_wsbeij.WSB_EiJ_v1.99.gtf# Descargamos el archivo AH73905 con anotaciones sobre el raton domestico
# (se carga desde la cache)
hub <- AnnotationHub()[["AH73905"]]
hub## EnsDb for Ensembl:
## |Backend: SQLite
## |Db type: EnsDb
## |Type of Gene ID: Ensembl Gene ID
## |Supporting package: ensembldb
## |Db created by: ensembldb package from Bioconductor
## |script_version: 0.3.4
## |Creation time: Sun Jul 7 08:07:59 2019
## |ensembl_version: 97
## |ensembl_host: localhost
## |Organism: Mus musculus
## |taxonomy_id: 10090
## |genome_build: GRCm38
## |DBSCHEMAVERSION: 2.1
## | No. of genes: 56393.
## | No. of transcripts: 144404.
## |Protein data available.# Informacion que podemos recuperar del hub
columns(hub)## [1] "DESCRIPTION" "ENTREZID" "EXONID"
## [4] "EXONIDX" "EXONSEQEND" "EXONSEQSTART"
## [7] "GENEBIOTYPE" "GENEID" "GENEIDVERSION"
## [10] "GENENAME" "GENESEQEND" "GENESEQSTART"
## [13] "INTERPROACCESSION" "ISCIRCULAR" "PROTDOMEND"
## [16] "PROTDOMSTART" "PROTEINDOMAINID" "PROTEINDOMAINSOURCE"
## [19] "PROTEINID" "PROTEINSEQUENCE" "SEQCOORDSYSTEM"
## [22] "SEQLENGTH" "SEQNAME" "SEQSTRAND"
## [25] "SYMBOL" "TXBIOTYPE" "TXCDSSEQEND"
## [28] "TXCDSSEQSTART" "TXID" "TXIDVERSION"
## [31] "TXNAME" "TXSEQEND" "TXSEQSTART"
## [34] "TXSUPPORTLEVEL" "UNIPROTDB" "UNIPROTID"
## [37] "UNIPROTMAPPINGTYPE"En la matriz de conteos del objeto sce, los genes aparecen con sus identificadores Ensembl (ENSMUSGxxxxxxxxxxx). Utilizaremos la base de datos disponible para cambiar estos identificadores por el nombre de los genes, para facilitar la interpretación y representación de los resultados del análisis:
# Identificadores Ensembl de los genes en los datos originales
head(rownames(sce_416b))## [1] "ENSMUSG00000102693" "ENSMUSG00000064842" "ENSMUSG00000051951"
## [4] "ENSMUSG00000102851" "ENSMUSG00000103377" "ENSMUSG00000104017"# Utilizamos la base de datos descargada para mapear los identificadores de cada
# gen (GENEID) del objeto sce y recuperar su localizacion cromosomica ("SEQNAME")
# y el nombre ("SYMBOL")
map_seqname <- mapIds(hub,
keys = rownames(sce_416b),
keytype = "GENEID",
column = "SEQNAME")
# Numero de genes por localizacion cromosomica
table(map_seqname)## map_seqname
## 1 10 11 12 13
## 3410 1598 2969 1506 1472
## 14 15 16 17 18
## 1771 1255 1148 1692 879
## 19 2 3 4 5
## 1026 3700 2876 2840 3235
## 6 7 8 9 CHR_MG65_PATCH
## 3071 3421 1762 2143 1
## GL456210.1 GL456211.1 GL456212.1 GL456216.1 GL456219.1
## 5 7 2 2 2
## GL456221.1 GL456233.1 GL456239.1 GL456350.1 GL456354.1
## 7 5 1 6 5
## GL456372.1 GL456381.1 GL456385.1 JH584292.1 JH584293.1
## 1 1 2 1 9
## JH584294.1 JH584295.1 JH584296.1 JH584297.1 JH584298.1
## 6 1 3 3 3
## JH584299.1 JH584303.1 JH584304.1 MT X
## 12 1 1 37 2581
## Y
## 1562map_symbol <- mapIds(hub,
keys = rownames(sce_416b),
keytype = "GENEID",
column = "SYMBOL")
# Nombre de cada gen correspondiente a su identificador Ensembl
head(map_symbol)## ENSMUSG00000102693 ENSMUSG00000064842 ENSMUSG00000051951 ENSMUSG00000102851
## "4933401J01Rik" "Gm26206" "Xkr4" "Gm18956"
## ENSMUSG00000103377 ENSMUSG00000104017
## "Gm37180" "Gm37363"# Añadimos los identificadores Ensembl, nombre y tipo de gen en la seccion
# rowData del objeto sce
rowData(sce_416b)$ENSEMBL <- rownames(sce_416b)
rowData(sce_416b)$SYMBOL <- map_symbol
rowData(sce_416b)$SEQNAME <- map_seqname
# Para una identificacion mas facil, renombramos las filas del objeto sce
# (identificadores de los genes) con sus nombres (SYMBOL), manteniendo el
# identificador original (ENSEMBL) en caso de que el nombre no este disponible.
# Se asegura que cada nombre sea unico
rownames(sce_416b) <- uniquifyFeatureNames(rowData(sce_416b)$ENSEMBL,
rowData(sce_416b)$SYMBOL)
head(rownames(sce_416b))## [1] "4933401J01Rik" "Gm26206" "Xkr4" "Gm18956"
## [5] "Gm37180" "Gm37363"Control de calidad
Las genotecas de baja calidad en los datos RNA-Seq, es decir, células con recuentos totales bajos pueden surgir por diversas causas, como el daño celular, o ineficiencia en la transcripción inversa o amplificación de los fragmentos genéticos. La presencia de estas genotecas de baja calidad pueden suponer un inconveniente ya que pueden alterar el proceso de normalización de los datos y contribuir a proporcionar resultados engañosos. Por tanto, antes de comenzar el análisis es conveniente detectar y eliminar la información procedente de estos casos. Síntomas indicativos de células o muestras de baja calidad serían:
Altas proporciones de secuencias relacionadas con los spike-ins en relación al conteo total de secuencias de todos los genes para cada muestra celular, teniendo en cuenta que los spike-ins se añaden a concentración conocida e igual en todas las muestras al inicio del proceso.
Altas proporciones de secuencias mitocondriales (cuya salida de la célula dañada es menos probable que la salida de las presentes en el citoplasma).
# Identificamos las secuencias mitocondriales (SEQNAME = MT)
mitocondrial <- which(rowData(sce_416b)$SEQNAME == "MT")
# Calculo de las metricas de calidad por celula indicando el subconjunto de
# secuencias que son mitocondriales
metricas <- perCellQCMetrics(sce_416b, subsets = list(Mt = mitocondrial))
metricas[1:3, ]## DataFrame with 3 rows and 16 columns
## sum detected percent_top_50 percent_top_100 percent_top_200
## <integer> <integer> <numeric> <numeric> <numeric>
## 1 865936 7618 26.7218362557972 32.2773276546997 39.7208338722492
## 2 1076277 7521 29.4043262097025 35.0354044544295 42.258080401235
## 3 1180138 8306 27.3453613052033 32.4769645583822 39.3295529844815
## percent_top_500 subsets_Mt_sum subsets_Mt_detected subsets_Mt_percent
## <numeric> <integer> <integer> <numeric>
## 1 52.9037942757894 78790 20 9.09882485541657
## 2 55.7454075484285 98613 20 9.16241822504801
## 3 51.9336721637639 100341 19 8.50248021841514
## altexps_ERCC_sum altexps_ERCC_detected altexps_ERCC_percent altexps_SIRV_sum
## <integer> <integer> <numeric> <integer>
## 1 65278 39 6.80658407035354 27828
## 2 74748 40 6.28029958040595 39173
## 3 60878 42 4.78949297995239 30058
## altexps_SIRV_detected altexps_SIRV_percent total
## <integer> <numeric> <integer>
## 1 7 2.9016456005055 959042
## 2 7 3.29130111124368 1190198
## 3 7 2.36477183861837 1271074Las métricas muestran el número de secuencias por cada gen (sum), número de genes por célula con conteos superiores a 0 (detected), el porcentaje de secuencias mitocondriales (subsets_Mt_percent), y el porcentaje de secuencias procedentes de spike-ins (altexps_ERCC_percent, altexps_SIRV_percent), entre otros. Destacar que dichas métricas son independientes del estado biológico de cada muestra/célula.
Para identificar células de baja calidad, aplicaremos filtros sobre las métricas obtenidas. Por ejemplo, podemos considerar que la calidad es baja si el tamaño de la genoteca es menor a 90000 lecturas, si se han expresado menos de 4000 genes, si la proporción de secuencias mitocondriales está por encima del 13% o si la de spike-ins está por encima del 10%:
# Criterios de celulas de baja calidad definidos por el usuario
cc_genoteca <- metricas$sum < 9e4
cc.ngenes <- metricas$detected < 4e3
cc_mito <- metricas$subsets_Mt_percent > 13
cc_spike <- metricas$altexps_ERCC_percent > 10
criterios_eliminacion <- cc_genoteca | cc.ngenes | cc_mito | cc_spike
# El numero de celulas eliminadas por los filtros anteriores seria
data.frame(tamanoGenoteca = sum(cc_genoteca),
nGenes = sum(cc.ngenes),
propMito = sum(cc_mito),
propSpikes = sum(cc_spike),
total = sum(criterios_eliminacion))## tamanoGenoteca nGenes propMito propSpikes total
## 1 3 0 2 19 22Con los filtros comentados (que no deben escogerse aleatoriamente), se eliminarían un total de 22 celulas de las 192 disponibles.
En caso de no disponer de una experiencia suficiente para establecer estos filtros, se puede hacer uso de la función quickPerCellQC(), la cual eliminará las células de baja calidad que considera como outliers en base a la mediana de la desviación absoluta (MDA) del valor mediano de cada métrica de calidad. Un outlier será aquel que sobrepase 3 MDAs por debajo de la mediana:
# Se tiene en cuenta el bloque (batch) de cada experimento para identificar
# outliers
eliminacion <- quickPerCellQC(metricas,
percent_subsets = c("subsets_Mt_percent",
"altexps_ERCC_percent"),
batch = sce_416b$block)
colSums(as.matrix(eliminacion))## low_lib_size low_n_features high_subsets_Mt_percent
## 5 0 2
## high_altexps_ERCC_percent discard
## 2 7Como vemos, con este último método serían menos las células eliminadas.
Gráficos de diagnóstico
Graficar las métricas de calidad puede ayudarnos a interpretar su distribución y detectar posibles problemas. En los casos más ideales, muestran una distribución normal, lo que justifica el uso del umbral de 3 MDAs en la detección de outliers.
# Creamos una copia del objeto sce original
sce_original <- sce_416b
# Añadimos a la copia las metricas calculadas y la informacion de que celulas
# se eliminarian y el bloque al que pertenecen cada una (todo esto se almacena
# en colData)
colData(sce_original) <- cbind(colData(sce_original), metricas)
sce_original$discard <- eliminacion$discard
sce_original$block <- factor(sce_original$block)
library(gridExtra)
grid.arrange(
plotColData(sce_original, x = "block", y = "sum", colour_by = "discard") +
scale_y_log10() +
ggtitle("Tamaño genoteca") +
theme(plot.title = element_text(hjust = 0.5)),
plotColData(sce_original, x = "block", y = "detected", colour_by = "discard") +
scale_y_log10() +
ggtitle("Genes detectados") +
theme(plot.title = element_text(hjust = 0.5)),
plotColData(sce_original, x = "block", y = "subsets_Mt_percent",
colour_by = "discard") +
scale_y_log10() +
ggtitle("% sec mitocondriales") +
theme(plot.title = element_text(hjust = 0.5)),
plotColData(sce_original, x = "block", y = "altexps_ERCC_percent",
colour_by = "discard") +
scale_y_log10() +
ggtitle("% sec spike-in") +
theme(plot.title = element_text(hjust = 0.5)),
nrow = 2, ncol = 2
)
En los gráficos anteriores, cada punto representa una célula en función al bloque al que pertenece, coloreada en función de si sería excluida o no.
En este ejemplo optamos por filtrar los datos eliminando las celulas detectadas por la función quickPerCellQC():
# Eliminamos de los datos principales los datos de celulas de baja calidad
sce_416b <- sce_416b[, !eliminacion$discard]Normalización
Normalizamos los datos para asegurarnos en lo posible que las diferencias de expresión de genes entre células se deba a factores biológicos y no a factores técnicos. Esta normalización es independiente de la estructura en bloques (batch) de los experimentos.
Podemos antes inspeccionar los 30 genes más expresados o con más conteos:
plotHighestExprs(sce_416b, exprs_values = "counts", n = 30)
Aplicaremos la transformación \(\small{log_2}\) explicada en secciones anteriores:
library(scran)
# Aplicamos la normalizacion logaritmica
sce_416b <- computeSumFactors(sce_416b)
sce_416b <- logNormCounts(sce_416b)Podemos ver la distribución de conteos sin normalizar y tras la normalización:
par(mfrow = c(1, 2))
hist(assay(sce_416b, "counts"), xlab = "Conteos", main = "")
hist(assay(sce_416b, "logcounts"), xlab = "Conteos normalizados", main = "")
Selección de genes
En este paso pretendemos seleccionar los genes que contienen la información útil acerca de procesos biológicos diferenciales (esto es, los genes que presentan más variabilidad entre células) y eliminar los que proporcionan ruido o un nivel base de expresión poco interesante.
Existen varios métodos para cuantificar la variación por gen y seleccionar un conjunto óptimo con alta variabilidad. Uno de ellos consiste en cuantificar la varianza de los valores de expresión o conteos normalizados logarítmicamente. Los que presenten mayor varianza biológica serán los que más contribuyan a la distancia eclídea entre células.
Este método de selección requiere modelar la relación media-varianza. La varianza de un gen depende más de su abundancia que de su heterogeneidad biológica subyacente. Para tener en cuenta este efecto, utilizamos la función modelGeneVar () para modelar la tendencia de la varianza con respecto a la abundancia en todos los genes. Con esta función podemos ajustar una línea de tendencia correspondiente a la variación “técnica” para cada gen (tech) en función de la abundancia, bajo el supuesto que la mayoría de genes expresan un nivel bajo inicial de variación que no es biológicamente interesante. Sin embargo, como en este estudio se han usado spike-ins, utilizaremos la versión modelGeneVarWithSpikes(), que ajusta la línea de tendencia sobre la expresión o conteos de los spike-ins, y no sobre los genes endógenos. La premisa aquí es que los spike-ins no deben verse afectados por la variación biológica, ya que son transcritos añadidos por el experimentador en una cantidad conocida, por lo que el valor ajustado de la tendencia sobre ellos debería representar una mejor estimación del componente de variación técnica para cada gen.
# Media y varianza de cada gen en base a los spike-ins, y teniendo en cuenta el
# bloqueo del diseño para eliminar la variabilidad debida al mismo. Los valores
# se calculan para cada bloque por separado, calculandose tambien una media
# conjunta.
# Utilizaremos los datos de los spike-ins ERCC
varGenesS <- modelGeneVarWithSpikes(sce_416b, "ERCC", block = sce_416b$block)
head(varGenesS)## DataFrame with 6 rows and 7 columns
## mean total tech
## <numeric> <numeric> <numeric>
## 4933401J01Rik 0 0 0
## Gm26206 0 0 0
## Xkr4 0 0 0
## Gm18956 0 0 0
## Gm37180 0.00790909509636549 0.0116350040552636 0.0377206732031878
## Gm37363 0 0 0
## bio p.value FDR
## <numeric> <numeric> <numeric>
## 4933401J01Rik 0 NA NA
## Gm26206 0 NA NA
## Xkr4 0 NA NA
## Gm18956 0 NA NA
## Gm37180 -0.0260856691479243 0.999748220046482 0.999998947604598
## Gm37363 0 NA NA
## per.block
## <DataFrame>
## 4933401J01Rik 0:0:0:...:0:0:0:...
## Gm26206 0:0:0:...:0:0:0:...
## Xkr4 0:0:0:...:0:0:0:...
## Gm18956 0:0:0:...:0:0:0:...
## Gm37180 0.015818190192731:0.0232700081105271:0.0754413464063756:...:0:0:0:...
## Gm37363 0:0:0:...:0:0:0:...Del resultado obtenemos la media y varianza total de los valores normalizados de cada gen (mean y total), el componente técnico y biológico de la varianza (tech y bio), y el p-valor original y ajustado (p.value y FDR) del test que evalúa la hipótesis nula bio <= 0. El componente biológico de la varianza de cada gen se calcula como el residuo o diferencia de su valor (total) con respecto al valor ajustado tech.
# Grafico media vs varianza de todos los genes para cada bloque experimental
par(mfrow = c(1, 2))
for (i in names(varGenesS$per.block)) {
bloque <- varGenesS$per.block[[i]]
plot(bloque$mean,
bloque$total,
main = i,
cex = 0.3,
col = "darkgrey",
xlab = "Media log(expresión)",
ylab = "Varianza log(expresión)")
curva <- metadata(bloque)
points(curva$mean, curva$var, col = "red2", pch = 16, cex = 0.5)
curve(curva$trend(x), col = "deepskyblue3", add = TRUE, lwd = 3)
}
Como podemos observar, no se da una diferencia notable entre bloques, lo cual indica que la replicación del proceso en cada uno se ha llevado a cabo adecuadamente. Así pues, podemos ver que la variación ajustada en base a los spike-ins (representados en rojo) aumenta linealmente con un aumento de la media de la abundancia normalizada, llegando un punto a partir del cual disminuye. El resto de puntos representados en gris corresponden a cada gen.
El p-valor se calcula asumiendo que las estimaciones de las varianzas se distribuyen de forma normal en torno a la línea ajustada de tendencia.
El componente biológico de la variación de cada gen es el que nos interesa para utilizarlo como métrica en la selección de los genes más variables. Podemos hacer una inspección rápida de qué genes tienen el mayor componente biológico de variación:
# Genes con mayor variación biologica
head(varGenesS[order(varGenesS$bio, decreasing = TRUE),
c("bio", "p.value", "FDR")])## DataFrame with 6 rows and 3 columns
## bio p.value FDR
## <numeric> <numeric> <numeric>
## Lyz2 12.2777343943255 0 0
## Ccl9 11.814342172529 0 0
## Top2a 11.2734488957067 3.89854931473513e-137 8.43397982206835e-135
## Cd200r3 11.2719357730414 1.17783154981232e-54 7.00721434772093e-53
## Ccnb2 10.5591416318493 1.20380054696742e-151 2.98404600168581e-149
## Hbb-bt 10.5323022856577 2.52638584279935e-49 1.34197330136375e-47Vemos que los seis genes con mayor componente de variabilidad biológica son: Lyz2, Ccl9, Top2a, Cd200r3, Ccnb2 y Hbb-bt. Además, con la función plotExpression() podemos obtener una visualización de la distribución de los valores de expresión normalizados logarítmicamente de dichos genes para cada célula, en función del fenotipo del estudio:
#(acortamos el nombre de los fenotipos para que encajen en el grafico)
sce_416b$phenotype <- factor(sce_416b$phenotype,
labels = c("induced CBFB-MYH11", "wild type"))
plotExpression(sce_416b,
c("Lyz2", "Ccl9", "Top2a", "Cd200r3", "Ccnb2", "Hbb-bt"),
x = "phenotype",
xlab = "",
colour_by = "phenotype")
Una vez que hemos cuantificado la variación por gen, el siguiente paso es seleccionar el subconjunto con mayor variabilidad biológica para continuar con análisis posteriores. Mantener un subconjunto mayor reducirá el riesgo de descartar señales biológicas interesantes, a costa de aumentar el ruido de genes irrelevantes. Por tanto, la selección del límite es importante. Algunas estrategias que nos pueden guiar, todas con sus ventajas e inconvenientes, son:
Seleccionar los genes con los mayores valores de variabilidad biológica. En el caso de contar con células más heterogéneas, es conveniente escoger un valor mayor que en el caso de trabajar con células más similares. Ventajas: el usuario puede controlar el número de genes mantenidos. Desventajas: los genes con mayor variabilidad pueden desplazar a otros con información interesante.
Seleccionar los genes con valores de FDR más significativos. Ventajas: mantiene la proporción de falsos positivos por debajo del límite establecido. Desventajas: resultado menos predecible, puede mantener genes biológicamente irrelevantes, o incluso muy pocos o demasiados genes (depende de la tasa de error tipo II del test y la severidad de las correcciones múltiples).
Seleccionar los genes con varianzas por encima de la linea de tendencia, es decir, con valores positivos de varianza biológica. Puede suponer una ventaja si trabajamos con poblaciones de células muy heterogéneas. Inconveniente: puede mantener más ruido.
Como ejemplo, en este caso mantendremos el top 10% de los genes con mayor variabilidad biológica (HVG), teniendo en cuenta que se asume que la mayoría de genes no presenta hetereogeneidad de expresión entre células (genes con funcionalidades comunes, etc.).
Nota: puede estudiarse el resultado de aplicar distintos valores, no hay un único valor óptimo.
# Creamos una lista con el 10% de los genes con mayor variabilidad biologica
genes_selec <- getTopHVGs(varGenesS, prop = 0.1)
str(genes_selec)## chr [1:1067] "Lyz2" "Ccl9" "Top2a" "Cd200r3" "Ccnb2" "Hbb-bt" "Cenpa" ...En el caso de haber querido aplicar el filtro sobre el FDR, usaríamos el argumento fdr.threshold de la función getTopHVGs(), y en caso de querer utilizar el método de mantener los genes con variabilidad biológica positivos, el argumento de la función a usar sería var.threshold = 0.
Corrección del batch
Eliminaremos el efecto de la posible covariable debida al efecto sistemático del batch o bloque del experimento (puede suponer uso de diferentes equipos, lotes de productos, personal, etc.), aunque por como se ha llevado a cabo la secuenciación en este experimento de ejemplo se espera que este efecto no sea notable. La función removeBatchEffect() ajusta un modelo lineal a los datos logarítmicos de expresión, incluyendo los bloques y tratamientos, y elimina el componente debido al efecto bloque:
library(limma)
# Incorporamos al objeto sce la correcion del efecto batch sobre los logcounts
assay(sce_416b, "bcorrected") <- removeBatchEffect(
logcounts(sce_416b),
design = model.matrix(~ sce_416b$phenotype), # condiciones a preservar
batch = sce_416b$block)Reducción de la dimensionalidad
Cada gen representa una dimensión de los datos. Si contáramos con datos de expresión de solo dos genes, podríamos obtener un gráfico bidimensional, siendo cada dimensión un gen, y cada punto del gráfico una célula. Sin embargo, los análisis de RNA cuentan con miles de genes, por lo que se hace necesario aplicar métodos de reducción de dimensionalidad para poder obtener tales representaciones, manteniendo toda la variabilidad posible en los datos. Esto es posible ya que muchos genes pueden estar correlacionados si se ven afectados por el mismo proceso biológico, por lo que se puede comprimir la información de múltiples de ellos en una sola dimensión.
En esta sección se muestra la aplicación del PCA, t-SNE y UMAP.
PCA
El análisis PCA lo llevamos a cabo sobre los datos normalizados logarítmicamente, con la función runPCA(). Por defecto se calculan las primeras 50 componentes principales (si no se especifica el argumento ncomponents), almacenándose en el apartado reducedDims del objeto sce. Además, lo aplicaremos sobre el subconjunto de genes seleccionados en el paso anterior y los datos corregidos por el bloque.
# Añadimos el PCA al objeto sce
sce_416b <- runPCA(sce_416b,
exprs_values="bcorrected",
subset_row = genes_selec,
ncomponents = 20)
str(reducedDim(sce_416b, "PCA"))## num [1:185, 1:20] -34.47 -19.25 -55.98 4.73 50.42 ...
## - attr(*, "dimnames")=List of 2
## ..$ : chr [1:185] "SLX-9555.N701_S502.C89V9ANXX.s_1.r_1" "SLX-9555.N701_S503.C89V9ANXX.s_1.r_1" "SLX-9555.N701_S504.C89V9ANXX.s_1.r_1" "SLX-9555.N701_S505.C89V9ANXX.s_1.r_1" ...
## ..$ : chr [1:20] "PC1" "PC2" "PC3" "PC4" ...
## - attr(*, "percentVar")= num [1:20] 22.42 2.96 1.38 1.27 1.17 ...# Grafico PCA: celulas coloreadas en funcion del fenotipo (contenido en el
# apartado colData del objeto sce)
plotReducedDim(sce_416b, dimred = "PCA", colour_by = "phenotype") +
ggtitle("PCA")
En este caso podemos ver cierta agrupación en función del fenotipo, con tan solo las dos primeras componentes que explican el 25% de la variabilidad total en los datos.
El PCA es una técnica de reducción lineal, y puede no dar resultados muy satisfactorios para visualizaciones de poblaciones celulares complejas. Otras técnicas disponibles de reducción no lineales son t-SNE y UMAP.
t-SNE
Una de las técnicas más populares en análisis RNA-Seq es t-SNE (T-distributed Stochastic Neighbor Embedding). Al igual que anteriormente, su cálculo se realiza sobre los valores normalizados logarítmicamente, y el resultado se almacena en el apartado reducedDims del objeto sce.
Una de las principales desventajas de t-SNE es que es más costoso computacionalmente comparado con otros métodos de visualización. Esto puede mitigarse indicando el argumento dimred = "PCA" en la función runtTSNE(), para realizar los cálculos sobre las componentes principales, aprovechando la compactación de datos y la eliminación de ruido proporcionado por ellas.
Tener en cuenta que esta técnica necesita que el usuario especifique parámetros específicos como el perplexity, el cual es aconsejable probar con distintos valores:
set.seed(395) # (importante establecer una semilla)
# t-SNE sobre las componentes principales
sce_416b <- runTSNE(sce_416b, dimred = "PCA", perplexity = 9)
# Grafico t-SNE
plotTSNE(sce_416b, colour_by = "phenotype") +
ggtitle("t-SNE (perplexity = 9)")
Tener en cuenta que con t-SNE no podemos interpretar las posiciones de los puntos como una relación directa de distancia clusters, ya que la técnica tiende a dispersar los clusters densos y comprimir los más dispersos, de forma que tampoco podemos interpretar el tamaño de los clusters como una medida de heterogeneidad de las subpoblaciones celulares que podrían estar representando.
UMAP
La técnica de reducción no lineal de la dimensionalidad UMAP (Uniform Manifold Approximation and Projection) tiende a generar clusters más compactos y a preservar más la estructura global de los datos en comparación con t-SNE. También es más rápida computacionalmente que esta última. Al igual que el t-SNE, cuenta con hiperparámetros que deben ser ajustados por el usuario. Entre estos, los que más efecto tienen son el número de vecinos (n_neighbors) y la distancia mínima entre puntos (min_dist) de la función runUMAP():
set.seed(596)
# UMAP
sce_416b <- runUMAP(sce_416b, n_neighbors = 15)
plotReducedDim(sce_416b,
dimred = "UMAP",
colour_by = "phenotype") + ggtitle("UMAP (n_neigbors = 15)")
Clustering
El clustering es una técnica de aprendizaje no supervisado que suele aplicarse en los análisis RNA-Seq para definir grupos de células con perfiles de expresión génica similares. Se trata de una técnia de exploración, con la cual podemos aplicar distintos algoritmos de agrupamiento y ajustar los parámtros respectivos para obtener distintas perspectivas de los datos. Algunas técnicas disponibles son:
Métodos basados en grafos
k-means clustering
Clustering jerárquico
Para este ejemplo aplicaremos el clustering jerárquico, ya que contamos con un grupo de datos no demasiado grande, teniendo en cuenta que puede ser lento de computar. Con esta técnica podemos generar dendogramas, que podemos “cortar” a diferentes alturas para definir subgrupos celulares con distinta granularidad.
Un procedimiento que puede aplicarse es realizar el clustering calculando la matriz de distancias célula-célula sobre las primeras componentes principales para luego visualizar el resultado sobre un gráfico t-SNE, de forma que este pueda ser utilizado para la comparación de los grupos obtenidos y si ciertos clusters pueden ser subdivididos:
library(dendextend)
# Distancia euclidea celula-celula
distancias_416b <- dist(reducedDim(sce_416b, "PCA"))
# Clustering jerarquico (hclust) sobre la matriz de distancias con el metodo de
# Ward. Dendograma
dendograma_416b <- as.dendrogram(hclust(distancias_416b, method = "ward.D2"))Podemos cortar el dendograma resultante a una determinada altura para obtener un número concreto de clústers:
# Corte del dendograma para obtener 4 clusters
clusters <- cutree(dendograma_416b, k = 4)
# Numero de celulas asignadas a cada cluster
table(clusters)## clusters
## 1 2 3 4
## 67 25 14 79# Grafico dendograma coloreado en funcion de los clusters obtenidos
labels_colors(dendograma_416b) <- clusters[order.dendrogram(dendograma_416b)]
plot(dendograma_416b, main = "Clustering jerárquico Ward")
Por último visualizamos los clústers sobre el gráfico t-SNE:
# Almacenamos en el objeto sce el cluster al que pertenece cada celula en base
# al algoritmo aplicado
sce_416b$cluster <- factor(clusters)
plotReducedDim(sce_416b, "TSNE", colour_by = "cluster")
Evaluación
Un paso importante es evaluar la calidad del clústering obtenido. Existen varios métodos, pero como ejemplo mostraremos el uso del índice silueta, el cual mide la validez y consistencia de las particiones. Concretamente, el índice asigna un valor a cada elemento en un rango [-1, 1]. Un valor próximo a -1 indica que el elemento posiblemente debería asignarse a otro clúster vecino. Valores próximos a 1 indican una correcta agrupación. Para calcular este índice aplicamos la función silhoutte() sobre el dendograma con 4 clústers, utilizano la matriz de distancias calculada anteriormente:
library(cluster)
# Calculo del indice Silueta para cada celula
siluetas <- silhouette(clusters, dist = distancias_416b)
plot(siluetas)
Cada línea representa el índice de una célula. Para cada clúster se indica el número de células y la media del índice silueta.
Un paso extra sería aplicar técnicas de bootstrapping para evaluar la consistencia y reproducibilidad del clústering, y las consecuencias de pequeños cambios en los datos originales.
Unión de resultados
Bioconductor nos ofrece también la función findMarkers() para encontrar marcadores o genes candidatos para grupos o clústers de células, evaluando la expresión diferencial entre pares de grupos mediante un test estadístico (t-test para valores escalados y normalizados logarítmicamente, Wilcoxon si el logaritmo no se ha aplicado, o binomial para conteos sin transformar).
# Genes marcadores candidatos, indicando el cluster asignado a cada celula
marcadores <- findMarkers(sce_416b,
groups = clusters,
test.type = "t",
block = sce_416b$block)[["1"]]
head(marcadores, 5)## DataFrame with 5 rows and 6 columns
## Top p.value FDR logFC.2
## <integer> <numeric> <numeric> <numeric>
## Spc25 1 1.46689125602357e-42 2.13634375299134e-39 5.81595068364625
## Ccna2 1 1.15687595017324e-83 5.39150467818736e-79 4.82456997326014
## Dctd 1 3.23303183109781e-26 7.84751122169182e-24 1.49996232468023
## Ckap2l 2 4.73204325678472e-27 1.28966166046313e-24 4.86140012070698
## Foxs1 2 1.03228582565941e-23 1.83620796255844e-21 -0.86370478292967
## logFC.3 logFC.4
## <numeric> <numeric>
## Spc25 0.158082606995869 6.29426356962192
## Ccna2 -0.183553794514964 7.28554766199255
## Dctd 5.62942347573833 3.90309323042851
## Ckap2l -1.66809301794946 5.0951567721544
## Foxs1 4.88357796761942 1.71440944757383marcadores_top10 <- rownames(marcadores)[marcadores$Top <= 10]
# Heatmap con informacion conjunta
plotHeatmap(sce_416b,
features = marcadores_top10,
order_columns_by = "cluster",
colour_columns_by = c("cluster", "block", "phenotype"),
center = TRUE,
symmetric = TRUE,
zlim = c(-5, 5))
En esta representación podemos ver qué genes (de los más variables entre grupos) presentan mayor o menor expresión en función del clúster y el fenotipo celular.
Podríamos llevar a cabo lo mismo centrándonos en un determinado tipo de genes de interés, por ejemplo los relacionados con las ciclinas, una familia de proteínas que controlan el ciclo celular y sus patrones de expresión en las fases del ciclo están bien caracterizados. Conociendo el patrón de los nombres de estos genes comienza podemos buscarlos en la matriz de expresión y recuperar sus datos de expresión:
genes_ciclinas <- grep("^Ccn[abde][0-9]$", rowData(sce_416b)$SYMBOL)
genes_ciclinas <- rownames(sce_416b)[genes_ciclinas]
plotHeatmap(sce_416b,
features = genes_ciclinas,
order_columns_by = "cluster",
colour_columns_by = c("cluster", "block", "phenotype"),
cluster_rows = FALSE)
BIBLIOGRAFÍA
Lun, A. T. L., F. J. Calero-Nieto, L. Haim-Vilmovsky, B. Gottgens, and J. C. Marioni. 2017. “Assessing the reliability of spike-in normalization for analyses of single-cell RNA sequencing data.” Genome Res. 27 (11):1795–1806.
Shield: 
This work by Cristina Gil Martínez is licensed under a Creative Commons Attribution 4.0 International License.
