02_sesion2.Rmd

# Reducción de dimensiones

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

## Diapositivas de Peter Hickey

Contenido adaptado de: [aquí](https://docs.google.com/presentation/d/12CjvQ1beZVeCcMQqD6ptzd9YUHbHOksvy3lt_rvWyRs/edit)

## Motivación

El siguiente paso en el análisis de scRNA-seq usualmente consiste en identificar grupos de células "similares".

Por ejemplo: un análisis de clustering busca identificar células con un perfil transcriptómico similar al calcular distancias entre ellas.

Si tuviéramos un dataset con dos genes podríamos hacer una gráfica de dos dimensiones para identificar clusters de células.

<div>
<p style = 'text-align:center;'>
<img src="img/cluster.gif" width="200px">
</p>
</div>


Pero... tenemos decenas de miles de genes : **Reducción de dimensionalidad**

## Reducción de dimensionalidad

Es posible porque la expresión de diferentes genes estaría correlacionada si estos genes son afectados por el mismo proceso biológico.

Por lo tanto, no necesitamos almacenar información independiente para genes individuales. Podemos comprimir múltiples "features" (genes) en una única dimensión.

**Ventajas:**

- Reduce trabajo computacional en análisis posteriores.
- Reduce el ruido al "promediar" mútiples genes obteniendo una representación más precisa de los patrones en los datos.
- Permite una graficación efectiva en dos dimensiones.

## Dataset ilustrativo: Zeisel


```{r, warning=FALSE, message=FALSE}
library(scRNAseq)
sce.zeisel <- ZeiselBrainData(ensembl = TRUE)

# Estos datos contienen tipos celulares previamente anotados
table(sce.zeisel$level1class)
```

Estudio de tipos celulares en el cerebro de ratón (oligodendrocitos, microglia, neuronas, etc) procesados con el sistema STRT-seq (similar a CEL-Seq)

Descripción [aquí](https://osca.bioconductor.org/zeisel-mouse-brain-strt-seq.html)

*Zeisel, A. et al. Brain structure. Cell types in the mouse cortex and hippocampus revealed by single-cell RNA-seq. Science 347, 1138-1142 (2015)*

```{r , warning=FALSE, message=FALSE}
# Quality control
# Descartar celulas con alto contenido mitocondrial o con alto porcentaje de spike-ins
library(scater)
is.mito <- which(rowData(sce.zeisel)$featureType == "mito")
stats <- perCellQCMetrics(sce.zeisel,
    subsets = list(Mt = is.mito)
)
qc <- quickPerCellQC(stats,
    percent_subsets = c("altexps_ERCC_percent", "subsets_Mt_percent")
)
sce.zeisel <- sce.zeisel[, !qc$discard]
```

```{r , warning=FALSE, message=FALSE}
# normalization
# encontramos unos clusters rápidos para las células y usamos esa información para calcular los factores de tamaño
library(scran)
set.seed(1000)
clusters <- quickCluster(sce.zeisel)
sce.zeisel <- computeSumFactors(sce.zeisel,
    cluster = clusters
)
sce.zeisel <- logNormCounts(sce.zeisel)
```


```{r , warning=FALSE, message=FALSE}
# variance-modelling
dec.zeisel <- modelGeneVarWithSpikes(
    sce.zeisel,
    "ERCC"
)
top.zeisel <- getTopHVGs(dec.zeisel, n = 2000)
```

- ¿Cómo se está modelando la relación media varianza?
- ¿Cómo se están obteniendo los HGVs?


## Dataset ilustrativo: 10x PBMC4k no filtradas


```{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"))
```

```{r , warning=FALSE, message=FALSE}
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)*

```{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 }
# variance modelling
set.seed(1001)
dec.pbmc <- modelGeneVarByPoisson(sce.pbmc)
top.pbmc <- getTopHVGs(dec.pbmc, prop = 0.1)
```

- ¿Cómo se está modelando la relación media varianza?
- ¿Cómo se están obteniendo los HGVs?



## Análisis de Componentes Principales

PCA es el arma principal de la reducción de dimensionalidad.

**PCA descubre las combinaciones (lineales) de "features" que capturan la cantidad más grande de variación**

En un PCA, la primer combinación lineal (componente principal) se elige tal que permite capturar la mayor varianza a través de las células. El siguiente PC se elige tal que es "ortogonal" al primero y captura la cantidad más grande de la variación restante, y así sucesivamente...

### PCA aplicado a datos de scRNA-seq

Podemos realizar reducción de dimensionalidad al aplicar PCA en la matriz de cuentas transformadas (log-counts matrix) y restringiendo los análisis posteriores a los primeros PCs **(top PCs)**.

* Esto puede reducir nuestro dataset de 20,000 dimensiones a, digamos, 10, sin perder demasiada información.
* La técnica de PCA tiene muchas propiedades teóricas bien estudiadas. 
* Hay varias formas rápidas de realizar PCA en datasets grandes.

### Suposiciones de PCA aplicadas a los datos de scRNA-seq

* Los procesos biológicos afectan múltiples genes en una manera coordinada.
* Los primeros PCs probablemente representan la estructura biológica dado que más variación puede ser capturada considerando el comportamiento correlacionado de muchos genes.
* Se espera que el ruido técnico azaroso afecte cada gen independientemente.

**Consideración: Los primeros PCs capturarón "batch effects" (efectos de lote) que afectan muchos genes en una manera coordinada**


```{r , warning=FALSE, message=FALSE}
## Aplicar PCA en la matriz de cuentas transformadas (log-counts matrix)
library(scran)
library(scater)
set.seed(100)
sce.zeisel <- runPCA(sce.zeisel,
    subset_row = top.zeisel
)
```


**¿Estamos corriendo el análisis sobre todos los genes de nuestro dataset?**

Por default, **runPCA()** usa un método rápido aproximado que realiza simulaciones, por lo tanto, es necesario *'configurar la semilla'* para obtener resultados reproducibles.

### Eligiendo el número de PCs

**Esta elección en análoga a la elección del número de highly variable genes (HGV). Elegir más PCs evitará descartar señal biológica a expensas de retener más ruido**

* Es común seleccionar un número de PCs "razonable" pero **arbitrario** (10-50), continuar con el análisis y regresar para checar la robustez de los resultados en cierto rango de valores.

Ahora exploraremos algunas estrategias guiadas por los datos (*data-driven*) para hacer esta selección.

#### Usando el punto del codo

```{r, warning=FALSE, message=FALSE, fig.dim = c(5, 4)}
library(PCAtools)
percent.var <- attr(reducedDim(sce.zeisel), "percentVar")
chosen.elbow <- PCAtools::findElbowPoint(percent.var)
plot(percent.var, xlab = "PC", ylab = "Variance explained (%)")
abline(v = chosen.elbow, col = "red")
```

Una heurística simple es elegir el número de PCs basado en el **porcentaje de varianza explicado** por PCs sucesivos.

#### Basados en la estructura de la población

Esta es una aproximación heurística más sofisticada que usa el número de clusters como un *proxy* del número de subpoblaciones.

```{r }
choices <- getClusteredPCs(reducedDim(sce.zeisel))
chosen.clusters <- metadata(choices)$chosen

plot(choices$n.pcs, choices$n.clusters,
    xlab = "Number of PCs", ylab = "Number of clusters"
)
abline(a = 1, b = 1, col = "red")
abline(v = chosen.clusters, col = "grey80", lty = 2)
```


Supongamos que esperamos *d* subpoblaciones de células, en ese caso, necesitamos *d-1* dimensiones para garantizar la separación de todas las subpoblaciones.

Pero... en un escenario real realmente no sabes cuántas poblaciones hay...

- Intenta con un rango para *d* y únicamente considera valores que produzcan a lo más *d+1* clusters
- Cuando se seleccionan más clusters con menos dimensiones se produce 'overclustering'.
- Elige una *d* que maximice el número de clusters sin caer en 'overclustering'.

**Ventaja**: Es una solución pragmática que soluciona el equilibrio sesgo-varianza en los análisis posteriores (especialmente en el análisis de clustering).

**Desventaja**: Hace suposiciones fuertes sobre la naturaleza de las diferencias biológicas entre los clusters, y de hecho supone la existencia de clusters, los cuales podrían no existir en algunos procesos biológicos como la diferenciación.

### Juntando todo

```{r }
set.seed(100)
# Calcula y almacena todos los posibles PCs
sce.zeisel <- runPCA(sce.zeisel, subset_row = top.zeisel)

# Selecciona y almacena el conjunto reducido de PCs:
# ... por el método del codo
reducedDim(sce.zeisel, "PCA_elbow") <- reducedDim(
    sce.zeisel, "PCA"
)[, 1:chosen.elbow]
# ... basado en la estructura poblacional 
reducedDim(sce.zeisel, "PCA_clusters") <- reducedDim(
    sce.zeisel, "PCA"
)[, 1:chosen.clusters]
```

### EJERCICIO

1. Realiza un PCA para los datos **sce.pbmc**.
2. Elige el número de PCs a conservar utilizando el método del codo.
3. Elige el número de PCs a conservar utilizando la estructura de la población.
4. Agrega esta información al objeto sce.pbmc.




### Usando el ruido técnico

Otra técnica de reducción de dimensiones consiste en **conservar todos los PCs hasta que el % de variación explicado alcance algun límite** (por ejemplo, basado en la estimación de la variación técnica).

**denoisePCA()** automáticamente selecciona el número de PCs.

Por default, denoisePCA() realiza algunas simulaciones, por lo tanto necesitamos *'configurar la semilla'* para obtener resultados reproducibles.

```{r, warning=FALSE, message=FALSE }
library(scran)
set.seed(111001001)
denoised.pbmc <- denoisePCA(sce.pbmc,
    technical = dec.pbmc, subset.row = top.pbmc
)
```


```{r }
dim(reducedDim(denoised.pbmc, "PCA"))
```

La dimensionalidad del output es el límite inferior para el número de PCs requeridos para explicar toda la variación biológica. Lo que significa que cualquier número menor de PCs definitivamente descartará algún aspecto de la señal biológica.

**Esto no grantiza que los PCs retenidos capturen toda la señal biológica**

**Esta técnica usualmente retiene más PCs que el método del punto del codo**

**scran::denoisePCA()** internamente limita el numero de PCs, por default 5-50, para evitar la selección de excesivamente pocos PCs (cuando el ruido técnico es alto relativo al ruido biológico) o excesivamente muchos PCs (cuando el ruido técnico es demasiado bajo).

#### ¿Qué pasa si calculamos la relación media-varianza con la función modelGeneVar para el dataset sce.pbmc (anteriormente usamos la función modelGeneVarByPoisson para este propósito)?

```{r }
dec.pbmc2 <- modelGeneVar(sce.pbmc)
denoised.pbmc2 <- denoisePCA(sce.pbmc,
    technical = dec.pbmc2, subset.row = top.pbmc
)
dim(reducedDim(denoised.pbmc2))
```

**scran::denoisePCA()** tiende a funcionar mejor cuando la relación media-varianza refleja el ruido técnico verdadero, *i.e* estimado por **scran::modelGeneVarByPoisson()** o **scran::modelGeneVarWithSpikes()** en vez de  **scran::modelGeneVar()**.

 *El dataset PBMC está cerca de este límite inferior: el ruido técnico es alto relativo al ruido biológico*.



#### ¿Qué pasa si calculamos el número de PCs usando el ruido técnico para el dataset sce.zeisel?

```{r }
set.seed(001001001)
denoised.zeisel <- denoisePCA(sce.zeisel,
    technical = dec.zeisel, subset.row = top.zeisel
)
dim(reducedDim(denoised.zeisel))
```


*Los datos de cerebro de Zeisel están cerca de este límite superior: el ruido técnico es demasiado bajo*

```{r }
dim(reducedDim(denoised.zeisel, "PCA"))
dim(reducedDim(denoised.zeisel, "PCA_elbow"))
dim(reducedDim(denoised.zeisel, "PCA_clusters"))
```

## Reducción de dimensionalidad para visualización

### Motivación

Los algoritmos de clustering, así como la mayoría de los algoritmos, operan fácilmente sobre 10-50 (a lo más) PCs, pero ese número es aún demasiado para la visualización.

Por lo tanto, necesitamos estrategias adicionales para la reducción de dimensionalidad si queremos visualizar los datos.

### Visualizando con PCA

```{r, fig.dim = c(6, 4) }
plotReducedDim(sce.zeisel, dimred = "PCA")
```

```{r , fig.dim = c(6, 4)}
plotReducedDim(sce.zeisel,
    dimred = "PCA",
    colour_by = "level1class"
)
```

PCA es una técnica lineal, por lo tanto, no es eficiente para comprimir diferencias en más de 2 dimensiones en los primeros 2 PCs.

### Retos y resumen de la visualización con PCA

```{r, fig.dim = c(6, 4) }
plotReducedDim(sce.zeisel,
    dimred = "PCA",
    ncomponents = 4, colour_by = "level1class"
)
```

**Ventajas:**

- PCA es predecible y no introducirá estructura artificial en los datos.
- Es determínistico y robusto a cambios pequeños en los valores de entrada.

**Desventajas:**

- Usualmente la visualización no es suficiente para visualizar la naturaleza compleja de los datos de scRNA-seq.

### Visualización con t-SNE

```{r, fig.dim = c(5, 4) }
set.seed(00101001101)
sce.zeisel <- runTSNE(sce.zeisel, dimred = "PCA")
plotReducedDim(sce.zeisel, dimred = "TSNE", colour_by = "level1class")
```

*t-stochastic neighbour embedding (t-SNE)* es la visualización por excelencia de datos de scRNA-seq. **Intenta encontrar una representación (no-lineal) de los datos usando pocas dimensiones que preserve las distancias entre cada punto y sus vecinos en el espacio multi-dimensional**

#### Retos de la visualización con t-SNE

```{r, fig.dim = c(6, 4) }
set.seed(100)
sce.zeisel <- runTSNE(sce.zeisel,
    dimred = "PCA",
    perplexity = 30
)
plotReducedDim(sce.zeisel,
    dimred = "TSNE",
    colour_by = "level1class"
)
```

#### Preguntas

- ¿Qué pasa si vuelves a correr **runTSNE()** sin especificar la semilla?

- ¿Qué pasa si especificas la semilla pero cambias el valor del parámetro *perplexity*?


#### Continuando


- Baja perplejidad favorece la resolución de la estructura fina, posiblemente al grado de que la visualización parece ruido random.


```{r, fig.width = 21 }
set.seed(100)

sce.zeisel <- runTSNE(sce.zeisel, dimred = "PCA", perplexity = 5)
p1 <- plotReducedDim(sce.zeisel, dimred = "TSNE", colour_by = "level1class")

sce.zeisel <- runTSNE(sce.zeisel, dimred = "PCA", perplexity = 20)
p2 <- plotReducedDim(sce.zeisel, dimred = "TSNE", colour_by = "level1class")

sce.zeisel <- runTSNE(sce.zeisel, dimred = "PCA", perplexity = 80)
p3 <- plotReducedDim(sce.zeisel, dimred = "TSNE", colour_by = "level1class")

library("patchwork")
p1 + p2 + p3
```



El siguiente [foro](http://distill.pub/2016/misread-tsne/) discute la selección de parámetros para t-SNE con cierta profundidad.


- No sobreinterpretes los resultados de t-SNE como un 'mapa' de las identidades de las células individuales.
- Algunos componentes aleatorios y la selección de parámetros cambiarán la visualización.
- La interpretación puede ser engañada por el tamaño y posición de los clusters.
- t-SNE infla clusters densos y comprime clusters escasos.
- t-SNE no está obligado a preservar las localizaciones relativas de clusters no-vecinos (no puedes interpretar distancias no locales).

**Aún así: t-SNE es una herramienta probada para visualización general de datos de scRNA-seq y sigue siendo muy popular**

### Visualización con UMAP

*Uniform manifold approximation and project (UMAP)* es una alternativa a t-SNE.

Así como t-SNE, UMAP **intenta encontrar una representación (no lineal) de pocas dimensiones de los datos que preserve las distancias entre cada punto y sus vecinos en el espacio multi-dimensional**.

t-SNE y UMAP están basados en diferentes teorías matemáticas.

```{r, fig.dim = c(6, 4) }
set.seed(1100101001)
sce.zeisel <- runUMAP(sce.zeisel, dimred = "PCA")
plotReducedDim(sce.zeisel,
    dimred = "UMAP",
    colour_by = "level1class"
)
```

Comparado con t-SNE:

- UMAP tiende a encontrar clusters visualmente más compactos.
- Intenta preservar más de la estructura global que t-SNE.
- Tiende a ser más rápido que t-SNE, lo cual puede ser importante para datasets grandes. La diferencia desaparece cuando se aplican sobre los primeros PCs.


#### Preguntas

- ¿Qué pasa si vuelves a correr **runUMAP()** sin especificar la semilla?

- ¿Qué pasa si especificas la semilla pero cambias el valor del parámetro *n_neighbors*?

#### Continuando

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

- Igual que para t-SNE, es necesario configurar una semilla y diferentes valores para los parámetros cambiaron la visualización.

- Si el valor para los parámetros *n_neighbors* o *min_dist* es demasiado bajo entonces el ruido aleatorio se interpretará como estructura de alta-resolución, si son demasiado altos entonces se perderá la estructura fina.

**TIP: Trata un rango de valores para cada parámetro para asegurarte de que no comprometen ninguna de las conclusiones derivadas de la gráfica UMAP o t-SNE**

### Interpretando las gráficas

**Recuerda:**


- Reducción de dimensionalidad para la visualización de los datos necesariamente involucra descartar información y distorsionar las distancias entre las células.

- No sobre interpretes las gráficas bonitas.

### Resumen y recomendaciones

- Las gráficas de t-SNE y UMAP son herramientas diagnóstico importantes, por ejemplo: para checar si dos clusters son realmente subclusters vecinos o si un cluster puede ser dividido en más de un cluster.

- Es debatible cuál visualización, t-SNE o UMAP, es más útil o estéticamente agradable.

- Está bien elegir aquella que funcione mejor para tu análisis (tomando en cuenta que tratarás la gráfica únicamente como una herramienta de visualización/diagnóstico y que no llegarás a ninguna conclusión fuerte basado únicamente en la gráfica).

## Dónde estamos


<div>
<p style = 'text-align:center;'>
<img src="img/resumen-reduccion.png" width="500px" heigth="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()
```