06-normalizacion.Rmd

# Normalización de datos

Instructora: [Ana Beatriz Villaseñor Altamirano](https://comunidadbioinfo.github.io/es/authors/anabva/)

## 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 otras tecnologías, single-cell RNA-seq (scRNA-seq) tiene errores sistemáticos y es importante conocerlos. 

Dentro de los más comunes se encuentran la captura de cDNA y la amplificación de PCR que se lleva acabo por célula. Tener poco material de inicio dificulta lograr una preparación de libreria consistente dando lugar a errores sistemáticos. 


```{r echo=FALSE, fig.cap="Figura tomada de [1]"}
knitr::include_graphics("https://scrnaseq-course.cog.sanger.ac.uk/website/figures/RNA-Seq_workflow-5.pdf.jpg")
```

La **normalización** tiene como objetivo remover estás 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 célulares o condiciones son debido a la biología y no por factores técnicos. 

### Ejercicio: Conceptos básicos

- Da ejemplos de sesgos técnicos (TIP: ¿Qué es RPKM?).

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

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



## 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 microfluio. 

- 3005 células y 18441 genes

- Contiene UMIs


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

### Ejercicio: QC

- ¿Cuántos genes son mitocondriales? (TIP: `is.mito`)

- ¿Cuántos genes tienen: bajas cuentas de librería, bajos features, alto porcentaje de expresión de ERCC, alto porcentaje de genes MT? ¿Cuántas células descartamos? (TIP: `perCellQCMetrics` y `quickPerCellQC`)

- Gráfica los resultados


```{r message=FALSE, warning=FALSE, include=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. 

Representa el estimado del sesgo relativo en cada célula. 

Se realiza dividiendo todas las cuentas de cada célula por un factor de escalamiento específico para cada una. Este factor de escalamiento se le conoce como *Library Size factor*.

$$ CuentasNormalizadas = Cuentas / Library Size factor$$

Suposición: Cualquier sesgo específico en cada célula (e.j. eficiencia en la captura o en la amplificación) afecta a **todos los genes de igual manera** a través de escalar por el promedio esperado de cuentas para dicha célula. 

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


### Tamaño de biblioteca (*Library Size*)

**Tamaño de biblioteca (*Library Size*):** La suma total de las cuentas a tráves de todos los genes en una 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 para cada célula.


El valor de *library size* es el que asumimos que escala con cualquier sesgo específico en cada célula.


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



Para escalar los datos ocuparemos un factor de escalamiento llamado *Library Size factor*. 

$$ Library Size \propto Library Size factor $$

Se calcula usando *library size*:

$$ Library Size factor =  {Library Size} / {mean(Library Size)}$$

Y se define de tal manera que el promedio de *Library Size factor*  en todas las células es 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 tamaños de librerías
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")
```



### Ejercicio: *library Size*

* Revisa los detalles (**Details**) en `?scater::librarySizeFactors`

* Calcula *library Size* `ls.zeisel`

* ¿Son idénticos `ls.zeisel` y `lib.sf.zeisel`?

* ¿Son proporcionales?

* Calcula `lib.sf.zeisel` de forma manual. TIP: Checa el [código fuente](https://github.com/LTLA/scuttle/blob/master/R/librarySizeFactors.R)

### Puntos finales

- Normalizar por *Library Size factor* asume que no hay desigualdad en la cantidad de genes differencialmente expresados (DE) entre dos células.  Es decir, que para cada grupo de genes sobre-expresados, debe existir un grupo de genes sub-expresados en la misma magnitud, cuando esto no pasa se le conoce como *sesgo de composición* (Veáse a continuación). 

- Para análisis exploratorios, la precisión de la normalización no es un punto mayor a considerar. El sesgo por composición normalmente no afecta la separación de los clusters, solo la magnitud. 

- La normalización por *Library Size factor* suele ser suficiente en algunas ocasiones donde se busca identificar clusters y los marcadores de los clusters.


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

Un sesgo técnico que es importante considerar es el *sesgo de composición* de RNA (transcriptoma). 

Supongamos que un gen X (o grupo de genes) se expresa en mayor cantidad en la célula A comparado a la célula B. Esto significa que más recursos fueron tomados por el gen X, disminuyendo la covertura de los demás.

¿Qué pasa si escalamos por tamaño de biblioteca?


```{r echo=FALSE, fig.cap="Figura tomada de [2]", fig.height=5}
knitr::include_graphics("https://hbctraining.github.io/DGE_workshop/img/normalization_methods_composition.png")
```

Este problema ha sido estudiado en bulk RNA-seq, `DESeq2::estimateSizeFactorsFromMatrix()` y `edgeR::calcNormFactors()`, contemplan este sesgo. Se assume que la mayoría de genes no estarán DE entre las muestras (en nuestro caso células) y cualquier diferencia entre los genes non-DE representa un sesgo el cual se remueve (calculando un factor de normalización). 

Sin embargo, single-cell RNA-seq tiene muchas cuentas bajas y ceros debido a limitaciones en la tecnología y no necesariamente indica ausencia de expresión. 

El método de `scran` resuelve este problema juntando las cuentas de varias células (pool) para incrementar el tamaño de las cuentas y obtener un factor de estimación que remueva el sesgo de composición de manera más precisa.

Este factor calculado con las cuentas pool se les regresa individualmente a cada célula mediante decircunvolución (deconvolution). Utilizando este factor se normalizan los datos con `scran::calculateSumFactors()`.


```{r}
# Normalización por decircunvolución (deconvolution)
library("scran")
# Pre-clustering
set.seed(100)
clust.zeisel <- quickCluster(sce.zeisel)

# Calcula factores de tamaño para la decircunvolución (deconvolution)
deconv.sf.zeisel <-
    calculateSumFactors(sce.zeisel, clusters = clust.zeisel, min.mean = 0.1)
# 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
)
abline(a = 0, b = 1, col = "red")
```


### Ejercicios: deconvolution


* ¿Cúantos clusters rápidos obtuvimos?

* ¿Cúantas células por cluster obtuvimos?

* ¿Cúantos clusters rápidos obtendríamos si cambiamos el tamaño mínimo a 200? Usa 100 como la semilla (seed).

* ¿Cúantas líneas ves en la gráfica?


### Puntos finales

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

`scran` algunas veces alcula factores negativos o ceros lo cual altera la matrix de expresión normalizada. ¡Checa los factores que calculas!

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

Si obtienes factores negativos intenta variar el número de clusters, checa si incrementar el número de células por cluster te dan factores positivos.

## Transformación logarítmica

### Motivació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)
```


Una vez calculados los factores de normalización con `computeSumFactors()`, podemos calular las cuentas en escala logaritmica usando `logNormCounts()`.

Estos valores resultantes son valores de expresión normalizados transformados en escala logarítmica.

```{r}
# Normalization
# set.seed(100)
# clust.zeisel <- quickCluster(sce.zeisel)
# sce.zeisel <- computeSumFactors(sce.zeisel, cluster=clust.zeisel, min.mean=0.1)

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

### Ejercicio: Transformación logarítmica

- ¿Qué es una pseudo-cuenta?

- ¿Porqué se usa?

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

- ¿Qué es la opción `downsample=TRUE`?



## Otras normalizaciones

Te invitamos a leer más sobre otras formas de normalizar, un lugar para empezar lo puedes encontrar en el 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()`.

### Seurat 

La normalización de Seurat con `NormalizeData()` (tomado de [aquí](https://github.com/satijalab/seurat/issues/3630)):

1. Dividir cada célula por el número total de moléculas medidas en la célula. - ¿Será *library size*?

2. Multiplicar ese número por un *scaling factor* (e.j. 10000)

3. Sumar 1 y tomar el logaritmo natural.

```{r, message=FALSE, warning=FALSE}
library("Seurat")
# Create a Seurat obj
sce <- sce.zeisel
sce <- removeAltExps(sce)
seurat.zeisel <- as.Seurat(sce, counts = "counts", data = NULL)
seurat.zeisel

# Normalize using Seurat function
seurat.zeisel <- NormalizeData(seurat.zeisel, normalization.method = "LogNormalize")
# Compare Total counts per cell after normalization
ls.seurat <- colSums(seurat.zeisel[[SingleCellExperiment::mainExpName(x = sce)]]@data)
## Relacionado a
## https://github.com/satijalab/seurat/blob/9b3892961c9e1bf418af3bbb1bc79950adb481d7/R/objects.R#L1041-L1046
## donde podemos ver como Seurat convierte el objeto de SingleCellExperiment
## a un objeto de Seurat

summary(ls.seurat)
hist(ls.seurat)

# Trying to replicate it
ls.zeisel <- colSums(counts(sce.zeisel))
summary(ls.zeisel)
step1 <- t(counts(sce.zeisel)) / ls.zeisel # matrix(2,2,2,2) /c(1,2)
step2 <- step1 * 10000
step3 <- t(log1p(step2))
ls.steps <- colSums(step3)
summary(ls.steps)
plot(ls.seurat, ls.steps)

# Compare with deconv normalization
ls.log <- colSums(logcounts(sce.zeisel))
plot(ls.seurat, ls.log)
```

Nota: [scanpy](https://scanpy-tutorials.readthedocs.io/en/latest/pbmc3k.html) ocupa un factor de normalización igual que Seurat.

## Notas finales


### Ejercicio: Conceptos básicos

- Da ejemplos de sesgos técnicos (TIP: ¿Qué es RPKM?).

"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)" - [hbctraining](https://hbctraining.github.io/DGE_workshop/lessons/02_DGE_count_normalization.html)

Algunos ejemplos de los sesgos técnicos que se contemplan son **la profundidad de secuenciación** y **la longitud del gen**. 

- Profundidad de secuenciación: Es necesario contemplarlo para comparar entre muestras, en este caso células. 

```{r echo=FALSE, fig.cap="Figura tomada de [2]", fig.height=5}
knitr::include_graphics("https://hbctraining.github.io/DGE_workshop/img/normalization_methods_depth.png")
```

- Longitud del gen: Es necesario contemplarlo para comparar entre genes. 

```{r echo=FALSE, fig.cap="Figura tomada de [2]", fig.height=5}
knitr::include_graphics("https://hbctraining.github.io/DGE_workshop/img/normalization_methods_length.png")
```



- ¿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



### Ejercicio: QC

- ¿Cuántos genes son mitocondriales? (recuerdas `is.mito`)

```{r}
length(is.mito)
```


- ¿Cuántos genes tienen: bajas cuentas de librería, bajos features, alto porcentaje de expresión de ERCC, alto porcentaje de genes MT? ¿Cuántas células descartamos? (TIP: `perCellQCMetrics` y `quickPerCellQC`)

```{r}
colSums(as.data.frame(qc))
```


- Gráfica los resultados

```{r}
# Plots
colData(unfiltered) <- cbind(colData(unfiltered), stats)
unfiltered$discard <- qc$discard

gridExtra::grid.arrange(
    plotColData(unfiltered, y = "sum", colour_by = "discard") +
        scale_y_log10() + ggtitle("Cuentas Totales"),
    plotColData(unfiltered, y = "detected", colour_by = "discard") +
        scale_y_log10() + ggtitle("Features (genes) detectados"),
    plotColData(unfiltered,
        y = "altexps_ERCC_percent",
        colour_by = "discard"
    ) + ggtitle("ERCC %"),
    plotColData(unfiltered,
        y = "subsets_Mt_percent",
        colour_by = "discard"
    ) + ggtitle("Mito %"),
    ncol = 2
)
```


### Ejercicio: *library Size*

* Revisa los detalles (**Details**) en `?scater::librarySizeFactors`

* Calcula *library size* `ls.zeisel`

```{r}
ls.zeisel <- colSums(counts(sce.zeisel))
summary(ls.zeisel)
hist(log10(ls.zeisel), xlab = "Log10[Library size]", col = "grey80")
```


* ¿Son idénticos `ls.zeisel` y `lib.sf.zeisel`? 

```{r}
identical(lib.sf.zeisel, ls.zeisel)
```


* ¿Son proporcionales?

```{r}
# Checamos proporcionalidad
plot(
    ls.zeisel,
    lib.sf.zeisel,
    log = "xy",
    main = "Proporcionalidad",
    xlab = "Library size",
    ylab = " Library size factor"
)
```


* Calcula `lib.sf.zeisel` de forma manual. TIP: Checa el [código fuente](https://github.com/LTLA/scuttle/blob/master/R/librarySizeFactors.R)

```{r}
## Ahora asegurate que su media sea 1 (unity mean)
lib_size_factors <- ls.zeisel / mean(ls.zeisel)
summary(lib_size_factors)
identical(lib_size_factors, lib.sf.zeisel)
```


### Ejercicios: deconvolution


* ¿Cúantos clusters rápidos obtuvimos?

```{r}
levels(clust.zeisel)
```


* ¿Cúantas células por cluster obtuvimos?

```{r}
cells_cluster <- sort(table(clust.zeisel))
cells_cluster
barplot(cells_cluster)
```


* ¿Cúantos clusters rápidos obtendríamos si cambiamos el tamaño mínimo a 200? Usa 100 como la semilla (seed).

```{r}
set.seed(100)
sort(table(quickCluster(sce.zeisel, min.size = 200)))
```


* ¿Cúantas líneas ves en la gráfica?

```{r}
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")
abline(a = -.2, b = 0.95, col = "red")
abline(a = 0.08, b = 1, col = "red")
```

### Ejercicio: Transformación logarítmica

- ¿Qué es una pseudo-cuenta?

Un número que se agrega para poder sacar el logarítmo 

- ¿Porqué se usa?

Por que `log(0) = -Inf` y produce error más adelante.

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

`pseudo.count = 1`

- ¿Qué es la opción `downsample=TRUE`?

[OSCA: Downsampling](https://bioconductor.org/books/release/OSCA/normalization.html#downsampling-and-log-transforming). Para cuando existe un efecto en los valores que se asocia a la *library size* a pesar de haber nomalizado.

Funciones interesantes para después de normalizar

```{r}
# sce.zeisel <- runPCA(sce.zeisel)
# plotPCA(sce.zeisel, colour_by = "level1class")
# plotRLE(sce.zeisel, exprs_values = "logcounts", colour_by = "level1class")
```


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

## 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>