01_sesion3.Rmd

# Normalización de datos

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

## Material

- Diapositivas de Peter Hickey: Ve las diapositivas [aquí](https://docs.google.com/presentation/d/1_tCNLiEsQ_TgsqHHf9_1lzXSaM_LunEHxBq3k130dQI/edit#slide=id.g7cc450648d_0_6)

- Capítulo de OSCA: Ve el capítulo del libro OSCA [aquí](https://bioconductor.org/books/release/OSCA/normalization.html)

## Motivación

Al igual que con otras tecnologías, single-cell RNA-seq (scRNA-seq) puede presentar diferencias sistemáticas en la cobertura de las librerías de scRAN-seq.

Típicamente se debe a diferencias técnicas en procesos como la captura de cDNA, la eficiencia de la amplificación por PCR. 

La **normalización** tiene como objetivo remover estas diferencias sistemáticas para que no interfieran cuando comparamos los perfiles de expresión entre células.

Al normalizar los datos, las diferencias observadas entre poblaciones celulares o condiciones son debido a la biología y no por factores técnicos. 


**La normalización es diferente a la corrección por lote **. El proceso de normalización se lleva a cabo sin importar la estructura de los lotes y sólo considera sesgos técnicos, los cuales tienden a afectar los genes en una manera medianamente similar, por ejemplo: la eficiencia de la amplificación por PCR debido al contenido de GC o la longitud del gene.



- ¿Qué es correción por lote (*batch effect correction*)? Da un ejemplo. 

"Large single-cell RNA sequencing (scRNA-seq) projects usually need to generate data across multiple batches due to logistical constraints. However, the processing of different batches is often subject to uncontrollable differences, e.g., changes in operator, differences in reagent quality. This results in systematic differences in the observed expression in cells from different batches, which we refer to as “batch effects”. Batch effects are problematic as they can be major drivers of heterogeneity in the data, masking the relevant biological differences and complicating interpretation of the results" -OSCA


- ¿Cuáles son las diferencias entre correción por lote y normalización?

"Normalization occurs regardless of the batch structure and only considers technical biases, while batch correction - as the name suggests - only occurs across batches and must consider both technical biases and biological differences. Technical biases tend to affect genes in a similar manner, or at least in a manner related to their biophysical properties (e.g., length, GC content), while biological differences between batches can be highly unpredictable" -OSCA



## Datos

Usaremos el dataset de [**Zeisel**](https://bioconductor.org/books/release/OSCA/zeisel-mouse-brain-strt-seq.html).

- Tipos celulares en cerebro de ratón (oligodendrocitos, microglias, neuronas, etc.)

- Procesado con STRT-seq (similar a CEL-seq), un sistema de microfluido. 

- 3005 células y 18441 genes.

- Contiene UMIs.


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

### Primero se debe hacer el control de calidad 

Eliminando células de baja calidad:

- Outliers bajos para log(cuenta total)

- Outliers bajos para log(número de features detectados)

- Outliers altos para el porcentaje de cuentas de conjuntos de genes especificados (mitocondriales, transcritos spike-in)


```{r message=FALSE, warning=FALSE}
# Control de calidad
library("scater")

unfiltered <- sce.zeisel

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

colSums(as.data.frame(qc))
sce.zeisel <- sce.zeisel[, !qc$discard]
```



## Normalización por escalamiento (scaling normalization)

La normalización por escalamiento es la estrategia más simple y usada. 

**Suposición:** Cualquier sesgo específico en cada célula  afecta a **todos los genes de igual manera** a través de escalar por el promedio esperado de cuentas para dicha célula. 

Se realiza de la siguiente manera:

1. Estima un **"size factor"** para cada célula.

2. Este factor representa el sesgo relativo en esa célula.

3. Divide las cuentas de los genes de cada célula entre su correspondiente **"size factor"**.


$$ CuentasNormalizadas = Cuentas / Size factor$$

Los valores de expresión normalizados pueden ser usados por análisis posteriores como *clustering* o *reducción de dimensiones*.


### Normalizacion por tamaño de biblioteca (*Library Size normalization*)

**Tamaño de biblioteca (*Library Size*):** La suma total de las cuentas a tráves de todos los genes en una célula. Asumimos que este valor escala con cualquier sesgo específico en cada célula.


$$Library Size_{cell} = \sum_{n=1}^{j} gene$$
Donde $j$ es el número total de genes y $gene$ es el número de cuentas por gen en la célula $cell$.

```{r echo=FALSE}
knitr::include_graphics(here::here("img/libsize.png"))
```


Para escalar los datos ocuparemos un factor de escalamiento llamado **Library Size factor**. Este valor se calcula a partir de escalar el tamaño de la biblioteca, tal que el promedio de los *Library Size factor* en todas las células sea igual a 1.


$$ mean(Library Size factor) = 1 $$

Lo que nos permite que los valores normalizados están en la misma escala y pueden ser útiles para la interpretación. 

```{r echo=FALSE}
knitr::include_graphics(here::here("img/libfactor.png"))
```


```{r}
# Estimar factores de normalización
lib.sf.zeisel <- librarySizeFactors(sce.zeisel)
# Examina la distribución de los tamaños de librerías
# que acabamos de estimar
summary(lib.sf.zeisel)
hist(log10(lib.sf.zeisel), xlab = "Log10[Library Size factor]", col = "grey80")
```

**¿A qué crees que se deba la variación en los factores de escalamiento?**
  
  
  
  
```{r}
## Calculando el tamaño de las librerias
ls.zeisel <- colSums(counts(sce.zeisel))
plot(ls.zeisel, lib.sf.zeisel, log="xy",
  xlab="Library size", ylab="Size factor")
```

### Puntos finales

- La normalización por tamaño de biblioteca asume que no hay un **sesgo de composición**, es decir, si existe un grupo de genes cuya expresión aumenta entonces debe existir otro grupo de genes para los que la expresión disminuye en la misma magnitud, esto generalmente no es verdad para los experimentos de scRNA-seq.

- Para algunos análisis exploratorios, la precisión de la normalización no es un punto mayor a considerar. La normalizacion por tamaño de biblioteca suele ser suficiente en algunas ocasiones donde se busca identificar clusters y los marcadores de los clusters, dado que el sesgo por composición normalmente no afecta la separación de los clusters, sólo la magnitud. 

- Sin embargo, esta falta de precisión podría afectar cuando se busca obtener estimaciones y estadísticas a nivel de genes individuales.


## Normalización por deconvolución (deconvolution)

El problema causado por el **sesgo de composición** y cómo eliminarlo ha sido estudiado en bulk RNA-seq por métodos como `DESeq2::estimateSizeFactorsFromMatrix()` y `edgeR::calcNormFactors()`. En ambos casos, se asume que la **mayoría de genes no estarán diferencialmente expresados** entre las muestras (en nuestro caso células) y cualquier diferencia entre los genes no diferencialmente expresados representa un sesgo el cual se remueve utilzando el factor de normalización calculado. 

Sin embargo, en single-cell RNA-seq se tienen muchas cuentas bajas y ceros debido a limitaciones en la tecnología, las cuales no necesariamente indican ausencia de expresión. Para este fenómento no están preparados los métodos utilizados en bulk RNA-seq. 

**La normalización por deconvolución (o normalización scran)**:

- Junta las cuentas de varias células (pool) para incrementar el tamaño de las cuentas, reduciendo ceros y cuentas bajas.
- Aplica los métodos de normalización que corrigen el sesgo de composición.
- Repiten esto para conjuntos sobrelapantes de células para obtener **size factors** en pool.
- **Deconvoluciona** estos **size factors** en pool para obtener **size factors** para cada célula.



```{r}
# Normalización por deconvolución (deconvolution)
library("scran")
# Pre-clustering (establece una semilla para obtener resultados reproducibles)
set.seed(100)
clust.zeisel <- quickCluster(sce.zeisel)

# Calcula factores de tamaño para la deconvolución (deconvolution)
deconv.sf.zeisel <-
    calculateSumFactors(sce.zeisel, clusters = clust.zeisel, min.mean = 0.1)
```

El **pre-clustering mejora la estimación de los size factors** al normalizar células similares juntas.

```{r}
# Examina la distribución de los factores de tamaño
summary(deconv.sf.zeisel)
hist(log10(deconv.sf.zeisel),
    xlab = "Log10[Deconvolution size factor]",
    col = "grey80"
)


plot(lib.sf.zeisel,
    deconv.sf.zeisel,
    xlab = "Library size factor",
    ylab = "Deconvolution size factor",
    log = "xy",
    pch = 16, 
    col=as.integer(factor(sce.zeisel$level1class))
)
abline(a = 0, b = 1, col = "red")
```


### Puntos finales

La normalización por deconvolución (deconvolution) mejora los resultados para análisis posteriores de una manera más precisa que los métodos para bulk RNA-seq.

El método `scran::calculateSumFactors` algunas veces calcula factores negativos o ceros lo cual altera la matriz de expresión normalizada. ¡Checa los factores que calculas! Si obtienes factores negativos intenta variar el número de clusters.

```{r}
summary(deconv.sf.zeisel)
```


## Transformación logarítmica

Una vez calculados los factores de normalización con `computeSumFactors()`, podemos calcular las cuentas en escala logarítmica usando `logNormCounts()`. ¿Qué es lo que hace esta función?

- Divide las cuentas de cada gene por el *size factor* correspondiente para esa célula (estos serían los valores normalizados) y entonces le aplica una transformción logarítmica a estos valores normalizados.
- Típicamente agrega una pseudo-cuenta de 1 para evitar calcular log(0).
- Estas cuentas en escala log se conocen como: valores de expresión normalizados en escala logarítmica (*log-transformed normalized expression values*).


```{r}
sce.zeisel

# Normalization
# Esto agrega el "sizeFactor" como parte del objeto
# set.seed(100)
# clust.zeisel <- quickCluster(sce.zeisel)
sce.zeisel <- computeSumFactors(sce.zeisel, cluster=clust.zeisel, min.mean=0.1)

sce.zeisel

# Log transformation
sce.zeisel <- scater::logNormCounts(sce.zeisel)
assayNames(sce.zeisel)
```

### Motivación, ¿Por qué hacemos esta transformación?

¿Qué gen es más interesante?

- *Gen X*: el promedio de expresión en el tipo celular A: 50 y B: 10

- *Gen Y*: el promedio de expresión en el tipo celular A: 1100 y B: 1000

```{r}
50 - 10
1100 - 1000

log(50) - log(10)
log(1100) - log(1000)
```

**Exacto, Gen X**

### Lluvia de ideas

- ¿Qué pasa cuando hacemos la transformación logarítmica?

- ¿Qué es una pseudo-cuenta? ¿Por qué se usa?

- ¿Qué valor de pseudo-cuenta usa `logNormCounts()`?




## Otras normalizaciones

Te invitamos a leer más sobre otras formas de normalizar, para empezar puedes consultar el siguiente curso del [Sanger Institute](https://www.singlecellcourse.org/basic-quality-control-qc-and-exploration-of-scrna-seq-datasets.html#normalization-theory).

Si estas interesad@ en diferencias en el contenido total de RNA en cada célula, checa la normalización por *spike-ins*. La cual asume que los *spike-ins* fueron añadidos en un nivel constante en cada célula.

Si tienes resultados donde el *library size* está asociado a tus datos a pesar de haber normalizado, checa la opción de `downsample=TRUE` dentro de la función de `logNormCounts()`.

Si te gustaría tener factores de normalización específicos para cada gen que tomen en cuenta la varianza, checa las funciones `varianceStabilizingTransformation()` de DESeq2 o `sctransform` de Seurat.

## Dónde estamos


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

## Adicionales

[1] [2018 BioInfoSummer Workshop](https://www.stephaniehicks.com/2018-bioinfosummer-scrnaseq/introduction-to-single-cell-rna-seq.html)

[2] [HBC training](https://hbctraining.github.io/DGE_workshop/lessons/02_DGE_count_normalization.html)

## Agradecimientos 

Este curso está basado en el libro [**Orchestrating Single Cell Analysis with Bioconductor**](https://osca.bioconductor.org/) de [Aaron Lun](https://www.linkedin.com/in/aaron-lun-869b5894/), [Robert Amezquita](https://robertamezquita.github.io/), [Stephanie Hicks](https://www.stephaniehicks.com/) y [Raphael Gottardo](http://rglab.org), además del [**curso de scRNA-seq para WEHI**](https://drive.google.com/drive/folders/1cn5d-Ey7-kkMiex8-74qxvxtCQT6o72h) creado por [Peter Hickey](https://www.peterhickey.org/).

Y en el material de la  [comunidadbioinfo/cdsb2020](https://github.com/comunidadbioinfo/cdsb2020) con el permiso de [**Leonardo Collado-Torres**](http://lcolladotor.github.io/).


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