11-anotacion_clusters.Rmd

# Anotación de clusters de células

Instructora: [**Yalbi I. Balderas-Martínez**](http://Yalbibalderas.github.io/).

## Diapositivas de Peter Hickey

Ver las diapositivas originales [aquí](https://docs.google.com/presentation/d/1CW-7FjL34kIKVEoZ2L_RVdXoEKQVzRbEuFUDVGJqnXI/edit)

## Motivación

* Ahora estamos a punto de obtener la interpretación biológica de los resultados 
* Esta es la tarea más retadora en los análisis de datos scRNA-seq


* 👉 La obtención de clústeres es más o menos directa
* 🤔 ¿Cuál es el estado biológico que está representado por cada uno de los clústeres?

* 👉 Necesitamos hacer un puente entre el _gap_ del dataset actual y el conocimiento biológico a priori (no siempre está disponible en una forma consistente y cualitativa)

* 🤔 ¿Qué es un tipo celular?
* 🔬 "Lo sabré cuando lo vea"
* 💻 "No"

Aplicaremos varios métodos computacionales que explotan la información _a priori_ para asignar el significado a un dataset no caracterizado de scRNA-seq.

Algunas fuentes de información _a priori_

- Conjuntos de genes curados (e.g. Gene Ontology)
- Perfiles de expresión de bases de datos publicadas de referencia
- Los datos raros que tú hayas escondido en tu cerebro
- Google

## Dataset ilustrativo: PBMC4k 10X sin filtrar

```{r, warning=FALSE, message=FALSE}
# Usemos datos de pbmc4k
library(BiocFileCache)
bfc <- BiocFileCache()
raw.path <- bfcrpath(bfc, file.path(
    "http://cf.10xgenomics.com/samples",
    "cell-exp/2.1.0/pbmc4k/pbmc4k_raw_gene_bc_matrices.tar.gz"
))
untar(raw.path, exdir = file.path(tempdir(), "pbmc4k"))

library(DropletUtils)
library(Matrix)
fname <- file.path(tempdir(), "pbmc4k/raw_gene_bc_matrices/GRCh38")
sce.pbmc <- read10xCounts(fname, col.names = TRUE)
```


Dataset "Células mononucleares humanas de sangre periférica" de 10X Genomics

Descripción [aquí](https://osca.bioconductor.org/unfiltered-human-pbmcs-10x-genomics.html) ^[Zheng, G. X. Y. et al. Massively parallel digital transcriptional profiling of single cells. Nat. Commun. 8, 14049 (2017).]

### Anotación

```{r, warning=FALSE, message=FALSE}
# Anotación de los genes
library(scater)
rownames(sce.pbmc) <- uniquifyFeatureNames(
    rowData(sce.pbmc)$ID, rowData(sce.pbmc)$Symbol
)
library(EnsDb.Hsapiens.v86)
location <- mapIds(EnsDb.Hsapiens.v86,
    keys = rowData(sce.pbmc)$ID,
    column = "SEQNAME", keytype = "GENEID"
)

# Detección de _droplets_ con células
set.seed(100)
e.out <- emptyDrops(counts(sce.pbmc))
sce.pbmc <- sce.pbmc[, which(e.out$FDR <= 0.001)]
```

### Control de calidad

```{r, warning=FALSE, message=FALSE}
# Control de calidad
stats <- perCellQCMetrics(sce.pbmc,
    subsets = list(Mito = which(location == "MT"))
)
high.mito <- isOutlier(stats$subsets_Mito_percent,
    type = "higher"
)
sce.pbmc <- sce.pbmc[, !high.mito]

# Normalización de los datos
library(scran)
set.seed(1000)
clusters <- quickCluster(sce.pbmc)
sce.pbmc <- computeSumFactors(sce.pbmc, cluster = clusters)
sce.pbmc <- logNormCounts(sce.pbmc)
```

### Genes variables

```{r, warning=FALSE, message=FALSE}
## Identificación de genes altamente variables
set.seed(1001)
dec.pbmc <- modelGeneVarByPoisson(sce.pbmc)
top.pbmc <- getTopHVGs(dec.pbmc, prop = 0.1)
```

### Reducción de dimensiones

```{r, warning=FALSE, message=FALSE}
## Reducción de dimensiones
set.seed(10000)
sce.pbmc <- denoisePCA(sce.pbmc,
    subset.row = top.pbmc,
    technical = dec.pbmc
)
set.seed(100000)
sce.pbmc <- runTSNE(sce.pbmc, dimred = "PCA")
set.seed(1000000)
sce.pbmc <- runUMAP(sce.pbmc, dimred = "PCA")
```

### Clustering

```{r, warning=FALSE, message=FALSE}
# clustering
g <- buildSNNGraph(sce.pbmc, k = 10, use.dimred = "PCA")
clust <- igraph::cluster_walktrap(g)$membership
sce.pbmc$cluster <- factor(clust)
```

## Asignando las etiquetas celulares a partir de los datos de referencia

### Visión general

* 👉 Un enfoque directo es comparar los perfiles de expresión single-cell con datasets previamente anotados
* 👉 Las etiquetas pueden entonces ser asignadas a cada célula en nuestro dataset no caracterizado de prueba basado en la muestra de referencia más similar, por dar alguna definición de "similar"


* 👉 Cualquier dataset de expresión génica etiquetado (microarreglos, RNA-seq bulk, scRNA-seq) puede ser usado como una referencia
* ⚠️ Sin embargo, su confiabilidad depende enormemente en la calidad de los datos originales y la experiencia de los autores originales quienes asignaron las etiquetas en primer lugar


* 👉 Asignar las etiquetas a un dataset de "prueba" a partir de un dataset de "entrenamiento" (referencia), es un problema estándar en estadística / _machine learning_
* 👉 Usaremos el método [SingleR (Aran et al. 2019)](https://bioconductor.org/packages/release/bioc/html/SingleR.html)



## SingleR

* 🤓 Asigna las etiquetas a las células basado en las muestras de referencia con las correlaciones de rangos más altas de Spearman
* 🤓 Para reducir el ruido, identifica genes marcadores entre pares de etiquetas (en la referencia) y calcula la correlación usando solamente esos marcadores 
* 🤓 Hace algún tipo de tuneado fino, repitiendo las correlaciones solamente con los genes marcadores de las etiquetas con el mejor score, ayudando a resolver cualquier ambigüedad entre esas etiquetas al eliminar el ruido a partir de marcadores irrelevantes para otras etiquetas



### SingleR incluye varias referencias

[Ver referencias](https://bioconductor.org/packages/3.11/bioc/vignettes/SingleR/inst/doc/SingleR.html#5_available_references)

```{r, warning=FALSE, message=FALSE, eval = FALSE}
# Human
celldex::BlueprintEncodeData()
celldex::DatabaseImmuneCellExpressionData()
celldex::HumanPrimaryCellAtlasData()
celldex::MonacoImmuneData()
celldex::NovershternHematopoieticData()

# Mice
celldex::ImmGenData()
celldex::MouseRNASeqData()
```



### Usando las referencias

```{r, warning=FALSE, message=FALSE}
# if needed install celldex
# create directory? y
library(celldex)
ref <- celldex::BlueprintEncodeData()
```

❓ ¿Qué tipos celulares están disponibles en este dataset de referencia? 



### Usando las referencias integradas

```{r, warning=FALSE, message=FALSE}
library(SingleR)
pred <- SingleR(
    test = sce.pbmc, ref = ref,
    labels = ref$label.main
)
```

* ❓ ¿Qué etiquetas han sido asignadas a los datos single-cell?
* ❓ ¿Cómo usaríamos las etiquetas "finas" con SingleR?


```{r, warning=FALSE, message=FALSE}
plotScoreHeatmap(pred)
```

* 👉 Inspeccionamos los resultados usando un heatmap de los scores por célula y por etiqueta
* 👉 Idealmente, cada célula debería exhibir un score alto en una etiqueta relativa a todas las otras
* 👉 Los scores se muestran antes de cualquier tuneado fino y son normalizadas a [0, 1] dentro de cada célula




### Podado de etiquetas (Label pruning)

```{r, warning=FALSE, message=FALSE}
total_pruned <- sum(is.na(pred$pruned.labels))
plotScoreHeatmap(pred, show.pruned = TRUE)
```

* 👉 SingleR intentará podar aquellas asignaciones de baja calidad marcándolas como NA
* 🤓 El podado se hace calculando la diferencia del score de la etiqueta asignada a partir del score de la mediana dentro de cada célula y entonces podando las células con un valor pequeño de esta diferencia

```{r, warning=FALSE, message=FALSE}
plotScoreDistribution(pred)
```

👉 Distribución de las diferencias del score de la etiqueta asignada a partir del score de la mediana dentro de cada célula



### Identificando los genes con anotación dirigida

* 🤔 ¿Por qué las células en este clúster se etiquetan como el tipo celular X?
* 👉 Examina la expresión de los genes marcadores para cada etiqueta en el dataset de prueba
* 👉 Si una célula en el dataset de prueba está asignado con confianza a una etiqueta en particular, uno esperaría que tenga una fuerte expresión de los marcadores de esa etiqueta (al menos sobreexpresión con respecto a las células asignadas a otras etiquetas)


```{r, warning=FALSE, message=FALSE}
# install gmp, ClusterR, mbkmeans dependencies if needed
sce.pbmc$labels <- pred$labels
all.markers <- metadata(pred)$de.genes
lab <- "B-cells"
# Get top-10 marker genes for B-cells compared to each other cell
# type
top.markers <- Reduce(union, sapply(all.markers[[lab]], head, 10))

plotHeatmap(sce.pbmc,
    order_columns_by = "labels",
    features = top.markers, center = TRUE, zlim = c(-3, 3), main = lab
)
```

❓ Toma otro tipo celular e identifica los genes que dirigen la anotación


### Comparando las etiquetas con los clústeres

```{r, warning=FALSE, message=FALSE}
tab <- table(Assigned = pred$pruned.labels, Cluster = sce.pbmc$cluster)

library(pheatmap)
# Proportion of cells in each cluster assigned to each label
pheatmap(prop.table(tab, margin = 2),
    color = colorRampPalette(c("white", "blue"))(101)
)
# (log-)number of cells in each cluster assigned to each label
# Adding a pseudo-count of 10 to avoid strong color jumps with just
# 1 cell.
pheatmap(log2(tab + 10),
    color = colorRampPalette(c("white", "blue"))(101)
)
```

### Voilà

```{r, warning=FALSE, message=FALSE}
plotTSNE(sce.pbmc, colour_by = "labels", text_by = "labels")

plotTSNE(sce.pbmc, colour_by = "cluster", text_by = "labels")
```

#### Aventura en el tiempo

<blockquote class="twitter-tweet"><p lang="en" dir="ltr">Hmm. We can&#39;t figure out why these two plots are different with the same code, though different R versions.<br><br>BioC 3.11? <a href="https://twitter.com/PeteHaitch?ref_src=twsrc%5Etfw">@PeteHaitch</a> <a href="https://t.co/8DJM6pagdj">https://t.co/8DJM6pagdj</a><br><br>BioC 3.13<a href="https://t.co/t8KiPTyU9D">https://t.co/t8KiPTyU9D</a><a href="https://twitter.com/Bioconductor?ref_src=twsrc%5Etfw">@Bioconductor</a> <a href="https://twitter.com/hashtag/rstats?src=hash&amp;ref_src=twsrc%5Etfw">#rstats</a> <a href="https://twitter.com/yalbi_ibm?ref_src=twsrc%5Etfw">@yalbi_ibm</a> <a href="https://twitter.com/AnaBetty2304?ref_src=twsrc%5Etfw">@AnaBetty2304</a> <a href="https://t.co/P65u7dnVCo">pic.twitter.com/P65u7dnVCo</a></p>&mdash; 🇲🇽 Leonardo Collado-Torres (@lcolladotor) <a href="https://twitter.com/lcolladotor/status/1425242252872409092?ref_src=twsrc%5Etfw">August 10, 2021</a></blockquote> <script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>

* https://github.com/MarioniLab/DropletUtils/issues/67

<blockquote class="twitter-tweet"><p lang="en" dir="ltr">Ok, not identical, but not bad with BioC 3.13 only<br><br>Left: BioC 3.13 with e.out from BioC 3.111<br>Middle: Pete&#39;s slides with BioC 3.11<br>Right: BioC 3.13 only<a href="https://twitter.com/hashtag/rstats?src=hash&amp;ref_src=twsrc%5Etfw">#rstats</a> <a href="https://twitter.com/hashtag/OSCA?src=hash&amp;ref_src=twsrc%5Etfw">#OSCA</a> <a href="https://twitter.com/Bioconductor?ref_src=twsrc%5Etfw">@Bioconductor</a> <a href="https://twitter.com/hashtag/DropletUtils?src=hash&amp;ref_src=twsrc%5Etfw">#DropletUtils</a> <a href="https://t.co/d41LiuXKWn">pic.twitter.com/d41LiuXKWn</a></p>&mdash; 🇲🇽 Leonardo Collado-Torres (@lcolladotor) <a href="https://twitter.com/lcolladotor/status/1425296392222822405?ref_src=twsrc%5Etfw">August 11, 2021</a></blockquote> <script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>

## Resumen de la anotación basada en una referencia (e.g., SingleR)

* ➕ Se centra en aspectos de los datos que se sabe son interesantes, simplifica el proceso de la interpretación biológica 
* ➖ Está restringido por la diversidad y la resolución de las etiquetas disponibles en el dataset de referencia
* 👉 Se pueden suplir referencias personalizadas a SingleR


## Asignando las etiquetas de tipos celulares a partir de marcadores

* 🤔 ¿Cómo podemos hacer uso de nuestros genes marcadores agrupados?
* 🥉 Revisarlos en hojas de cálculo 
* 🥈 Observar heatmaps
* 🥇 Realizar un gene set enrichment analysis



### Gene set enrichment analysis

* 👉 Identifica las rutas y procesos que están (relativamente) activos en cada clúster basado en la sobreexpresión de los genes asociados en comparación con otros clústeres
* ➕ Un método confiable para determinar si las rutas están sobre- o sub- expresadas entre clúesteres 
* ➕ Existen un montón de herramientas para gene set enrichment analysis
* ➖ Todas las conclusiones son relativas a otros clústeres, haciéndolo más difícil para determinar la identidad celular si alguno no está presente en el mismo estudio [más info](https://osca.bioconductor.org/cell-type-annotation.html#assigning-cluster-labels-from-markers)

### Calculando las actividades de los conjuntos de genes

* 👉 Calcular el promedio de la expresión en log en todos los genes, en un conjunto de genes para cada célula y examinar los clústeres con valores altos (_gene set activities_)
* 👉 Se necesita proveer de conjuntos de genes 
* ➖ No todos los genes en el conjunto pueden exhibir el mismo patrón de diferencia y los genes no-DE añadirán ruido, "diluyendo" la fuerza de cualquiera de las diferencias comparadas a un análisis que se centra directamente en genes DE 
* 👉 Es más una visualización útil que la base para cualquier análisis estadístico real [más info](https://osca.bioconductor.org/cell-type-annotation.html#computing-gene-set-activities)


## Resumen y recomendaciones

* 👉 La anotación de tipos celulares "automática", como SingleR, es mejor cuando funciona (i.e. cuando hay un dataset de referencia apropiado)
* 👉 Usualmente necesitaremos usar un método manual, como aquellos basados en agrupar los genes marcadores  (e.g., gene set enrichment analysis)
* 👉 La anotación del tipo celular ofrecerá una reconsideración inmediata de los parámetros del agrupamiento y/o algunos retoques manuales a los clústeres 

## Detalles de la sesión de R

```{r}
## Información de la sesión de R
Sys.time()
proc.time()
options(width = 120)
sessioninfo::session_info()
```

## Patrocinadores {-}

Agradecemos a nuestros patrocinadores:

<a href="https://comunidadbioinfo.github.io/es/post/cs_and_s_event_fund_award/#.YJH-wbVKj8A"><img src="https://comunidadbioinfo.github.io/post/2021-01-27-cs_and_s_event_fund_award/spanish_cs_and_s_award.jpeg" width="400px" align="center"/></a>

<a href="https://www.r-consortium.org/"><img src="https://www.r-consortium.org/wp-content/uploads/sites/13/2016/09/RConsortium_Horizontal_Pantone.png" width="400px" align="center"/></a>