02_sesion3.Rmd

# Clustering

Instructora: [Laura Gómez-Romero](https://comunidadbioinfo.github.io/es/authors/lgomez/)

Este contenido está basado en las diapositivas de Peter Hickey. Ve las diapositivas [aquí](https://docs.google.com/presentation/d/1YTmP5QLxRg5JIT_bnQtMADg3eszsKEz7l5u6rOMeQ1A/edit). Y en el curso de OSCA, lee el material [aquí](https://bioconductor.org/books/release/OSCA/clustering.html)

## Dataset ilustrativo: 10X PBMC4k no filtrado

```{r, warning=FALSE, message=FALSE}
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 de células mononucleares de sangre periférica humana (PBMC) 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)*

```{r, warning=FALSE, message=FALSE}
# gene-annotation
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"
)

# cell-detection
set.seed(100)
e.out <- emptyDrops(counts(sce.pbmc))
sce.pbmc <- sce.pbmc[, which(e.out$FDR <= 0.001)]
```

```{r, warning=FALSE, message=FALSE}
# quality-control
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]

# normalization
library(scran)
set.seed(1000)
clusters <- quickCluster(sce.pbmc)
sce.pbmc <- computeSumFactors(sce.pbmc, cluster = clusters)
sce.pbmc <- logNormCounts(sce.pbmc)
```

```{r, warning=FALSE, message=FALSE}
# variance modelling
set.seed(1001)
dec.pbmc <- modelGeneVarByPoisson(sce.pbmc)
top.pbmc <- getTopHVGs(dec.pbmc, prop = 0.1)
```

```{r, warning=FALSE, message=FALSE}
# dimensionality-reduction
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")
```

- **¿Para qué sirve el método que estamos usando para reducir la dimensionalidad de los datos?**

- **¿Los HGVs están almacenados en nuestro objeto sce.pbmc?**


## Motivación

**Clustering** es un procedimiento **no supervisado** par definir grupos de células con perfiles de expresión similares

Su propósito principal es resumir los datos en un formato digerido susceptible a interpretación humana

Nos permite asignar **etiquetas** (por ejemplo, tipos celulares) a las células

## ¿Por qué no realizamos el clustering sobre las coordenadas de t-SNE/UMAP?

Las técnicas de t-SNE/UMAP han comprimido datos altamente multi-dimensionales en dos dimensiones

Esta compresión inevitablemente ha provocado la perdida de información

Por lo tanto, agrupamos sobre los PCs y después visualizamos las identidades de los clusters en la gráfica t-SNE/UMAP

## ¿Cuál es el verdadero clustering?


<div>
<p style = 'text-align:center;'>
<img src="img/real-clustering.png" width="800px">
</p>
</div>

**Un cluster no implica un tipo celular**

Nosotros podemos definir tantos clusters como queramos y podemos utilizar el algoritmo que más nos acomode

**El clustering, como un microscopio, simplemente es una herramienta para explorar los datos**

Preguntar por el mejor clustering es similar a preguntar cuál es la mejor magnificación en un microscopio sin contenido

## Clustering basado en grafos


### Antecedentes

El clustering basado en grafos fue popularizado (más no inventado) por su uso en Seurat

**Objetivo: Construir un grafo en el que cada nodo es una célula que está conectada a sus *vecinos más cercanos* (otras células con perfiles de expresión similares) en el espacio multidimensional**


### Gráfica de los *k* vecinos más cercanos (k-nearest neighbour -KNN- graph)

Ilustremos como funciona para 20 células

<div>
<p style = 'text-align:center;'>
<img src="img/knn1.png" width="400px">
</p>
</div>

<div>
<p style = 'text-align:center;'>
<img src="img/knn2.png" width="400px">
</p>
</div>

<div>
<p style = 'text-align:center;'>
<img src="img/knn3.png" width="400px">
</p>
</div>

<div>
<p style = 'text-align:center;'>
<img src="img/knn19.png" width="400px">
</p>
</div>

### Gráfica de los vecinos más próximos compartidos (SNN)

De una gráfica KNN se puede construir una gráfica SNN. En este tipo de grafo, dos células estarán conectadas por una arista si comparten alguno de sus vecinos más próximos.

<div>
<p style = 'text-align:center;'>
<img src="img/knn-snn-1.png" width="800px">
</p>
</div>

<div>
<p style = 'text-align:center;'>
<img src="img/knn-snn-3.png" width="800px">
</p>
</div>
 
### Gráfica SNN con pesos en las aristas

Podemos asignar **pesos** a cada arista del grafo, basándonos en la similaridad de las células involucradas, dándole pesos más altos a células que están más cercanamente relacionadas

Para ver los distintos esquemas de pesado puedes consultar la documentación de la función **makeSNNGraph** del paquete *bluster*. Algunos ejemplos son:

- **Rango**: El peso entre dos nodos está dado por *k*-*r*/2 donde *r* es la suma más pequeña de los rangos (de proximidad, el vecino más cercano tiene el rango 1) para cualquiera de los vecinos compartidos
- **Número**: el peso entre dos nodos es igual al número de vecinos más próximos compartidos
- **Jaccard**: el peso entre dos nodos es igual a la similaridad de Jaccard entre los conjuntos de vecinos de estos nodos

### Obteniendo comunidades a partir de una gráfica SNN pesada mediante un algoritmo de clustering

A partir de una gráfica SNN pesada podemos aplicar algoritmos para identificar **comunidades** de células

**Una comunidad es un grupo de células que están más conectadas a células en el mismo grupo que lo que están a células de un grupo diferente**

Cada comunidad representa un cluster

<div>
<p style = 'text-align:center;'>
<img src="img/knn-snn-3-com.png" width="800px">
</p>
</div>

### Resumen de clustering basado en grafos

* La construcción y búsqueda de una red KNN es rápida, por lo tanto, es escalable para datasets grandes
* Debes evitar obtener conclusiones fuertes acerca de la forma de los clusters o la distribución de células dentro de cada cluster
* El algoritmo conecta cada célula con un número mínimo de células vecinas lo cual reduce, más no elimina, el riesgo de clusters no informativos 

Después de la construcción del grafo, no se almacena información adicional más alla de las células vecinas. Esto puede producir subclusters artificiales en regiones con muchas células

### Detalles a considerar en la implementación

* ¿Cuántas céulas vecinas debo considerar durante la construcción del grafo?
* ¿Cómo debo pesar las aristas?
* ¿Cuál algoritmo de detección de comunidades se debe usar para definir los clusters?


### Implementación

```{r, warning=FALSE, message=FALSE}
library(scran)
# Construir grafo usando k= 10 vecinos más cercanos en el espacio definido por el PCA
g <- buildSNNGraph(sce.pbmc, k = 10, use.dimred = "PCA")
# Identificar comunidades utilizando el método Walktrap
clust <- igraph::cluster_walktrap(g)$membership
```

```{r, warning=FALSE, message=FALSE, fig.dim = c(6, 4)}
# Visualizar clusters en una gráfica t-SNE 
library(scater)
sce.pbmc$cluster <- factor(clust)
plotReducedDim(sce.pbmc, "TSNE", colour_by = "cluster")
```

**¿Qué pasa si utilizas una k más grande o más pequeña?**

```{r, warning=FALSE, message=FALSE}
library(scran)
#Construir grafo usando k= 50 vecinos más cercanos en el espacio definido por el PCA 
g50 <- buildSNNGraph(sce.pbmc, k = 50, use.dimred = "PCA")
# Identificar comunidades utilizando el método Walktrap
clust50 <- igraph::cluster_walktrap(g50)$membership
```

```{r, warning=FALSE, message=FALSE, fig.dim = c(6, 4)}
# Visualizar clusters en una gráfica t-SNE 
library(scater)
sce.pbmc$cluster50 <- factor(clust50)
plotReducedDim(sce.pbmc, "TSNE", colour_by = "cluster50")
```


**En esta implementación:**

* La construcción de la red KNN se baso en la distancia Euclideana entre células
* La construcción de la red KNN implica que las aristas se crean entre todos los pares de células que comparten por lo menos un vecino
* Se utilizó el esquema de peso de: [Xu and Su (2015)](https://pubmed.ncbi.nlm.nih.gov/25805722/)

### Eligiendo un valor de *k*

* El valor de *k* puede ser toscamente interpretado como el tamaño anticipado de la subpoblación más pequeña

* Si una subpoblación tiene menos que *(k+1)* células entonces el método será forzado a construir aristas entre células de esa subpoblación y células de otras subpoblaciones, incrementando el riesgo de que la subpoblación en cuestión no forme su propio cluster

### Una implementación diferente: estilo Seurat

```{r, warning=FALSE, message=FALSE, fig.dim = c(6, 4)}
# Pesos definidos estilo Jaccard seguidos por clustering de Louvain 
# aka 'clustering estilo Seurat' 
g2 <- buildSNNGraph(sce.pbmc, k = 10, use.dimred = "PCA", type = "jaccard")
clust2 <- igraph::cluster_louvain(g2)$membership
sce.pbmc$cluster2 <- factor(clust2)
plotReducedDim(sce.pbmc, "TSNE", colour_by = "cluster2")
```

### Detalles de las implementaciones más comunes

**Pipelines basados en Seurat:**

* Pesos basados en Jacard
* Clustering Louvain

**Pipelines basados en Scran:**

* Pesos basados en Rangos
* Clustering Walktrap

Para detalles sobre la seleccion de parámetros y comparaciones: [visitar esta página](https://bioconductor.org/books/release/OSCA/clustering.html#clustering-graph).

```{r, fig.cap = "Estilo scran vs estilo Seurat.", fig.width=10}
library("patchwork")

## Estilo scran vs estilo Seurat
plotReducedDim(sce.pbmc, "TSNE", colour_by = "cluster") +
plotReducedDim(sce.pbmc, "TSNE", colour_by = "cluster2")
```

### Otras implementaciones

**Distintas métricas de distancia**

```{r, eval = FALSE}
g.num <- buildSNNGraph(sce.pbmc, use.dimred = "PCA", type = "number")
g.jaccard <- buildSNNGraph(sce.pbmc, use.dimred = "PCA", type = "jaccard")
g.none <- buildKNNGraph(sce.pbmc, use.dimred = "PCA")
```

**Distintos métodos de clustering**

```{r, eval = FALSE}
clust.louvain <- igraph::cluster_louvain(g)$membership
clust.infomap <- igraph::cluster_infomap(g)$membership
clust.fast <- igraph::cluster_fast_greedy(g)$membership
clust.labprop <- igraph::cluster_label_prop(g)$membership
clust.eigen <- igraph::cluster_leading_eigen(g)$membership
```



## Evaluando la separación de los clusters

**Modularidad** es una métrica natural para evaluar la separación entre comunidades/clusters

La modularidad se define como la diferencia (escalada) entre el peso total observado de las aristas entre los nodos en el mismo cluster y el peso total esperado si los pesos fueran distribuidos aleatoriamente entre todos los pares de nodos

Nosotros calcularemos un **score de modularidad** para cada cluster usando las tasas en vez de las diferencias, debido a que las tasas no se ven tan fuertemente influenciadas por el tamaño de los clusters


```{r, warning=FALSE, message=FALSE}
library(bluster)

# obteniendo la métrica de modularidad
ratio <- pairwiseModularity(g, clust, as.ratio = TRUE)
dim(ratio)
```


```{r, warning=FALSE, message=FALSE, fig.dim = c(6, 4)}
library(pheatmap)
pheatmap(log2(ratio + 1),
    cluster_rows = FALSE,
    cluster_cols = FALSE,
    color = colorRampPalette(c("white", "blue"))(100)
)
```

Un dataset que contiene clusters bien separados debería contener la mayoría del peso total observado en las entradas diagonales, *i.e* la mayoría de las aristas ocurren entre células del mismo cluster

Para más detalles sobre evaluación de la separación entre clusters [visite esta página](https://osca.bioconductor.org/clustering.html#assessing-cluster-separation)


## Otros métodos de clustering

**Clustering por k-means**

 - **PRO:** Rápido
 - Se debe especificar el número de clusters de antemano
 - Favorece clusters esféricos
 
 
 **Clustering jerárquico**
 
 - **Produce un dendograma (árbol) representando las células y la similaridad entre subpoblaciones a varias resoluciones**
 - Demasiado lento para correrse en algo más grande que los datasets más pequeños de scRNA-seq
 
## Evaluando la estabilidad de los clusters

Una propiedad deseable de cualquier cluster es que éste sea estable a las perturbaciones en los datos de entrada, de esta manera:

* Pequeños cambios en el procesamiento de los datos no cambiarán el resultado
* Se incrementa la probabilidad de que las conclusiones puedan ser replicadas en un estudio independiente

Uno puede hacer un proceso de **bootstrap** para evaluar la estabilidad de un algoritmo de clustering en un dataset dado y calcular la **coasignación**. 

**La coasignación es la probabilidad de que células elegidas al azar del cluster X y Y sean asignadas al mismo cluster en la réplica del proceso de bootstrap**

```{r, warning=FALSE, message=FALSE}
myClusterFUN <- function(x) {
    g <- buildSNNGraph(x, use.dimred = "PCA", type = "jaccard")
    igraph::cluster_louvain(g)$membership
}

originals <- myClusterFUN(sce.pbmc)
set.seed(0010010100)
coassign <- bootstrapStability(sce.pbmc,
    FUN = myClusterFUN,
    clusters = originals
)
```


```{r, warning=FALSE, message=FALSE, fig.dim = c(6, 4)}
pheatmap(coassign,
    cluster_row = FALSE, cluster_col = FALSE,
    color = rev(viridis::magma(100))
)
```

**Probabilidad alta** de coasignación indica que X es **estable** con respecto a su separación de Y.

**Queremos altas probabilidades de coasignación en la diagonal**

Debes considerar que el bootstraping solo considera el efecto del ruido de muestreo e ignora otros factores que pueden afectar la reproducibilidad (como efectos de batch)

Además, **una pobre separación puede ser altamente estable**

## Subclustering


* Mejora la resolución al repetir el proceso de *feature selection* y *clustering* dentro de un único cluster

* Se enfoca en los HGVs y PCs que son los más relevantes para un cluster específico

**Veamos como se comporta la expresión de ciertos genes en nuestros clusters**

*CD3E, CCR7, CD69, y CD44 son marcadores de células T de memoria*. 

```{r, warning=FALSE, message=FALSE, fig.dim = c(6, 4)}
g.full <- buildSNNGraph(sce.pbmc, use.dimred = "PCA")
clust.full <- igraph::cluster_walktrap(g.full)$membership
sce.pbmc$clust.full <- factor(clust.full)
plotExpression(sce.pbmc,
    features = c("CD3E", "CCR7", "CD69", "CD44"),
    x = "clust.full", colour_by = "clust.full"
)
```

De esta gráfica deducimos que las células T de memoria se encuentran en el cluster 10

Dentro de las  células T de memoria, ¿dónde están las subpoblaciones CD4+ y CD8+?


```{r, warning=FALSE, message=FALSE}
# Repetimos TODO el proceso de clustering en el subconjunto de células que hemos identificado como células T de memoria (cluster 10).
memory <- 10
sce.memory <- sce.pbmc[, clust.full == memory]
dec.memory <- modelGeneVar(sce.memory)
sce.memory <- denoisePCA(sce.memory,
    technical = dec.memory,
    subset.row = getTopHVGs(dec.memory, prop = 0.1)
)

g.memory <- buildSNNGraph(sce.memory, use.dimred = "PCA")
clust.memory <- igraph::cluster_walktrap(g.memory)$membership
sce.memory$clust.memory <- factor(clust.memory)
```

```{r, warning=FALSE, message=FALSE, fig.dim = c(6, 4)}
plotExpression(sce.memory,
    features = c("CD8A", "CD4"),
    x = "clust.memory"
)
```

Expresión de CD4 es bajo, por lo tanto, su cambio es modesto, pero la interpretación es clara

Si tipos celulares o estados celulares se extienden sobre las fronteras de los clusters, entonces un subcluster podría representar contaminación de un tipo celular en un cluster separado

## Resumen y recomendaciones

- **Un cluster no implica un tipo celular**

- Nosotros podemos definir tantos clusters como queramos y podemos utilizar el algoritmo que más nos acomode

- **El clustering, como un microscopio, simplemente es una herramienta para explorar los datos**

- Preguntar por el mejor clustering es similar a preguntar cuál es la mejor magnificación en un microscopio sin contenido

- Clustering basado en grafos es rápido y evita tener que hacer suposiciones fuertes sobre la forma de los clusters o la distribución de las células dentro de cada cluster:

  1. **scran::buildSNNGraph() **
  
  2. **igraph::cluster_walktrap()** o **igraph::cluster_louvain()**
 
- **Modularidad y estabilidad** de los clusters son diagnósticos útiles

- El proceso de **subclustering** podría mejorar la resolución dentro de clusters grandes


## Dónde estamos


<div>
<p style = 'text-align:center;'>
<img src="img/resumen-clustering.png" width="500px" height="400px">
</p>
</div>

## 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()
```