01_sesion2.Rmd

# Control de calidad

Dra. Evelia Coss

07 de agosto de 2023

## Material y Diapositivas

[
```{r,echo=FALSE}
knitr::include_url("https://comunidadbioinfo.github.io/cdsb2023/dia1_sesion2_Diapositivas", height = "380px")
```
](https://comunidadbioinfo.github.io/cdsb2023/dia1_sesion2_Diapositivas)

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

Curso 2021 impartido por Leonardo Collado Torres: [Ver información aquí](https://comunidadbioinfo.github.io/cdsb2021_scRNAseq/control-de-calidad.html) 


## Enfoques o Ideas principales

Detectar **células verdaderas y de alta calidad** ✅🍪, de modo que cuando agrupemos nuestras células sea más fácil **identificar poblaciones de tipos celulares distintos**.

Identificar las **muestras fallidas** e intentar **salvar los datos o eliminarlas** del análisis, además de intentar comprender por qué falló la muestra.

## ¿Por qué hay problemas con los datos?

- Problemas con el rompimiento de las células.
- Fallos en la preparación de las bibliotecas (ineficiencias en la transcripción inversa, amplificacion por PCR, etc).
- Más de una célula durante la secuenciación (douplets o multiplets).
- Problemas en el alineamiento.

## Preguntas Básicas en Control de Calidad

Se responden durante el análisis de SingleCell RNA-Seq.

- ¿Cuántos droplets traen más de una célula? (douplets 🍪🍪 o multiplets 🍪🍪🍪). 
- ¿Cuántas células murieron durante el proceso de secuenciacion?

Cada **droplets** debe contener una **sola célula** 🍪.


```{r ImagenControl, echo=FALSE, out.width='80%', fig.align='center'}
knitr::include_graphics("img/quality_control.jpeg")
```
Figura obtenida de [Single-Cell best practices](https://www.sc-best-practices.org/preprocessing_visualization/quality_control.html).


## La mala calidad en los datos puede ser debida a varios factores

Una o varias células podemos encontrar:

- ⚠️ Baja cantidad de cuentas totales
- ⚠️ Baja cantidad de genes expresados
- ⚠️ Alta proporción de cuentas (reads) provenientes de la mitocondria
- ⚠️ Alta proporción de cuentas provenientes de las secuencias control (ERCC spike-in control transcripts)

Mas informacion sobre [ERCC spike-in control transcripts](https://www.agilent.com/cs/library/applications/ERCC%20Spike-in%20App%20Note%205994-1767EN_1.7.pdf).

## Recomendaciones

* Tener un correcto diseño experimental (evitar efecto *bash*). 👾
* Correcta preparación de las muestras. 🎮
* Analizar la calidad de los datos. 🍪


## Tipos de filtros

- **Fixed thresholds:**  Valores de corte fijos y estrictos (FDR < 0.5) ♠️
- **Adaptative thresholds:** Valores de cortes adaptados al comportamiento de nuestros datos.♦️


## Parámetros empleados para evaluar la calidad con la función `addPerCellQC`

La función `addPerCellQC` provenie del paquete `scater`. Agrega las siguientes métricas en cada célula y por cada gen dentro del mismo archivo. 


- `sum`: Número de cuentas (lecturas) totales de cada célula.
- `detected`: Genes expresados con al menos una cuenta.
- `altexps_ERCC_percent`: Porcentaje de cuentas mapeadas de las secuencias control (ERCC spike-in control transcripts).
- `subsets_Mito_percent`: Porcentaje de cuentas mapeadas provenientes de la mitocondria.

## Ejemplo: linea celular 416 en ratón

Realizaremos las siguientes actividades con este ejemplo:

- Eliminar células de baja calidad.
- Comparar entre los tipos de filtros (Fixed thresholds y Adaptative thresholds).
- Filtrar células baja calidad.
- Visualización de datos crudos y filtrados.


```{r ImagenDataset, echo=FALSE, out.width='80%', fig.align='center'}
knitr::include_graphics("img/416Dataset.png")
```
Línea celular de células mieloides progenitoras inmortalizadas de ratón usando [SmartSeq2](https://osca.bioconductor.org/lun-416b-cell-line-smart-seq2.html).


### Paquetes e Importar los datos en R

Los paquetes que vamos a emplear para esta sección son:

```{r Paquetes, message=FALSE, warning = FALSE}
library(scRNAseq) ## para descargar datos de ejemplo
library(DropletUtils) ## para detectar droplets 
library(Matrix) ## para leer datos en formatos comprimidos
library(AnnotationHub) ## para obtener información de genes
library(scater) ## para gráficas y control de calidad
library(BiocFileCache) ## para descargar datos
library(EnsDb.Hsapiens.v86) ## Archivo de anotacion en humanos en Ensembl
library(dplyr) ## Modificacion de archivos dataframe
```

Cargar los datos empleando el paquete `scRNAseq`

```{r Dataset416b}
sce.416b <- LunSpikeInData(which = "416b")
# Conversion a factor, evitemos mensajes de error
sce.416b$block <- factor(sce.416b$block)
```

El paquete `scRNAseq` contiene múltiples data sets compilados en funciones, para saber más da click [aquí](https://www.bioconductor.org/packages/release/data/experiment/vignettes/scRNAseq/inst/doc/scRNAseq.html#ref-kotliarov2020citeseq).

### Anotación de genes

```{r Anotacion}
ah <- AnnotationHub()
query(ah, c("Mus musculus", "Ensembl", "v97"))

# Anotacion de genes con la localizacion de cada cromosoma
ens.mm.v97 <- ah[["AH73905"]] # solo un cromosoma
location <- mapIds(ens.mm.v97, keys=rownames(sce.416b),
    keytype="GENEID", column="SEQNAME")
# deteccion de genes mitocondriales
is.mito <- which(location=="MT")
```

### Análisis de calidad con `addPerCellQC`

Al imprimir la variable `sce.416b` podemos observar que es de tipo *SingleCellExperiment*.


```{r ControlDeCalidad}
# Agregar la informacion de la calidad por cada gen y celula en un mismo archivo
# Identificar genes mitocondriales
sce.416b <- addPerCellQC(sce.416b,
    subsets = list(Mito = is.mito)
)
sce.416b
```

Podemos observar la informacion de este data frame usando `colData`:

```{r InformacionDataset}
# observar la informacion contenida en el dataframe
colData(sce.416b)

# Ver el nombre de las columnas
colnames(colData(sce.416b))
```

### Preguntas sobre los datos

- ¿Cuántas células detectadas se encuentran en el dataframe `sce.416b`? ^[El dataset cuenta con 192 células, `sce.416b`]
- ¿Cuántos genes fueron analizados con al menos una cuenta? ^[El dataset cuenta con 46,604 genes detectados, `sce.416b`]

Podemos conocer estas respuestas observando la informacion contenida en `sce.416b`, en dimensiones. 

### Visualización de los datos crudos

Genes detectados en cada uno de los bloques de secuenciación.

```{r visualizar_qc}
plotColData(sce.416b, x = "block", y = "detected")
```

Genes detectados por cada tratamiento en cada bloque de secuenciación.

```{r visualizar_qc_2}
plotColData(sce.416b,
    x = "block",
    y = "detected",
    other_fields = "phenotype"
) +
    scale_y_log10() +  # Cambiar la escala en Y a logaritmo
    facet_wrap(~phenotype) + # Dividir tratamiento
    labs(x = "Bloques de secuenciación", y="Genes detectados \n(log)", title = "Genes detectados") # Etiquetas en X, Y y titulo
``` 

### Filtro A: Fixed thresholds

Para el filtrado y eliminación de las células de baja calidad emplearemos las variables `sum`, `detected`, `altexps_ERCC_percent` y 
`subsets_Mito_percent`. 

```{r FiltroEstricto}
# --- Fixed thresholds ---
# Valores de corte fijos y estrictos
qc.lib <- sce.416b$sum < 100000 # menos de 100 mil cuentas (lecturas)
qc.nexprs <- sce.416b$detected < 5000 # menos de 5 mil genes
qc.spike <- sce.416b$altexps_ERCC_percent > 10 # 10 % de las cuentas alineando a ERCC
qc.mito <- sce.416b$subsets_Mito_percent > 10 # 10 % alineando al genoma mitocondrial
discard <- qc.lib | qc.nexprs | qc.spike | qc.mito # si falla alguno de estos eliminarlo

# secuencias control = (ERCC spike-in control transcripts)

# Resumen del número de células
# Número de células eliminadas por cada filtro
DataFrame(
    LibSize = sum(qc.lib),
    NExprs = sum(qc.nexprs),
    SpikeProp = sum(qc.spike),
    MitoProp = sum(qc.mito),
    Total = sum(discard)
)
```

### Filtro B: Adaptative thresholds

Valores de cortes adaptados al comportamiento de nuestros datos.

```{r FiltroAdaptativo}
# --- Adaptative thresholds ---
## Usando isOutlier() para determinar los valores de corte
qc.lib2 <- isOutlier(sce.416b$sum, log = TRUE, type = "lower")
qc.nexprs2 <- isOutlier(sce.416b$detected,
    log = TRUE,
    type = "lower"
)
qc.spike2 <- isOutlier(sce.416b$altexps_ERCC_percent,
    type = "higher"
)
qc.mito2 <- isOutlier(sce.416b$subsets_Mito_percent,
    type = "higher"
)
discard2 <- qc.lib2 | qc.nexprs2 | qc.spike2 | qc.mito2

# Extraemos los límites de valores (thresholds)
attr(qc.lib2, "thresholds")

attr(qc.nexprs2, "thresholds")

# Obtenemos un resumen del número de células
# eliminadas por cada filtro
DataFrame(
    LibSize = sum(qc.lib2),
    NExprs = sum(qc.nexprs2),
    SpikeProp = sum(qc.spike2),
    MitoProp = sum(qc.mito2),
    Total = sum(discard2)
)
```

### Preguntas por resolver

- ¿Cuántos células se descartan con cada filtro? 
- ¿Cuántas células se comparten entre ambos filtros? 
- ¿Cúal filtro fue más estricto?
- ¿Existen células que el filtro adaptativo excluyera y que no lo hiciera el filtro estricto?

Para contestar estas preguntas podemos emplear:

```{r Respuesta}
addmargins(table(discard, discard2))
```

## ¿Cómo funciona `isOutlier()`?

Supongamos que la mayor parte del conjunto de datos está formado por células de alta calidad 1

* (Opcional: log-transformar la métrica QC) <sup>1</sup>.
* Calcular la mediana de la métrica QC
* Calcular la desviación absoluta de la mediana (MAD<sup>2</sup>) de la QC
* Identifique los valores atípicos como aquellas celdas con una métrica QC a más de 3 MAD<sup>3</sup> de la mediana en la dirección "problemática".
  - Puede controlar cuántas MAD son aceptables
  - Puede decidir qué dirección es problemática

QC = Quality Control
1 But we'll relax that assumption in a few moments
2 MAD is similar to standard deviation
3 Loosely, 3 MADs will retain 99% of non-outliers values that follow a Normal distribution. 

Figura explicativa, [Diapositiva 24](https://docs.google.com/presentation/d/1pIiA7fZd1GBxaKQpzT2sPn7C6fIZxEJ8NI4p3UtoIOo/edit#slide=id.g7e8a22e342_0_12).


## Consideran el `Batch`


```{r Batch}
## Determino el bloque (batch) de muestras
batch <- paste0(sce.416b$phenotype, "-", sce.416b$block)

## Versión de isOutlier() que toma en cuenta los bloques de muestras
qc.lib3 <- isOutlier(sce.416b$sum,
    log = TRUE,
    type = "lower",
    batch = batch
)
qc.nexprs3 <- isOutlier(sce.416b$detected,
    log = TRUE,
    type = "lower",
    batch = batch
)

qc.spike3 <- isOutlier(sce.416b$altexps_ERCC_percent,
    type = "higher",
    batch = batch
)
qc.mito3 <- isOutlier(sce.416b$subsets_Mito_percent,
    type = "higher",
    batch = batch
)
discard3 <- qc.lib3 | qc.nexprs3 | qc.spike3 | qc.mito3

# Extraemos los límites de valores (thresholds)
attr(qc.lib3, "thresholds")

# Obtenemos un resumen del número de células
# eliminadas por cada filtro
DataFrame(
    LibSize = sum(qc.lib3),
    NExprs = sum(qc.nexprs3),
    SpikeProp = sum(qc.spike3),
    MitoProp = sum(qc.mito3),
    Total = sum(discard3)
)
```

## Visualización gráfica de las células de buena calidad

Usaremos el filtro estricto y visualizaremos la distribución de los datos. 

```{r Visualizar_filtro}
# Reducir el nombre en los fenotipos /  tratamientos
sce.416b_edited <- sce.416b
  
sce.416b_edited$phenotype <- ifelse(grepl("induced", sce.416b_edited$phenotype),
    "induced", "wild type")
# Verificar
unique(sce.416b_edited$phenotype)

# Agregar columna de genes descartados
sce.416b_edited$discard <- discard

# Visualizacion grafica
plotColData(sce.416b_edited, x="block", y="sum", colour_by="discard",
        other_fields="phenotype") + facet_wrap(~phenotype) +
      labs(x = "Bloques de secuenciacion", y="Cuentas totales", title = "Numero de cuentas") # Etiquetas en X, Y y titulo

```

Visualizacion de todas las variables en una sola gráfica.

```{r Visualizar_todo}
# Visualizacion grafica
gridExtra::grid.arrange(
    plotColData(sce.416b_edited, x="block", y="sum", colour_by="discard",
        other_fields="phenotype") + facet_wrap(~phenotype) + 
        scale_y_log10() + labs(x = "Bloques de secuenciacion", y="Cuentas totales", title = "Numero de cuentas"), # Etiquetas en X, Y y titulo
    plotColData(sce.416b_edited, x="block", y="detected", colour_by="discard", 
        other_fields="phenotype") + facet_wrap(~phenotype) + 
        scale_y_log10() + labs(x = "Bloques de secuenciacion", y="Genes expresados", title = "Numero de genes /n(log)"), # Etiquetas en X, Y y titulo
    plotColData(sce.416b_edited, x="block", y="subsets_Mito_percent", 
        colour_by="discard", other_fields="phenotype") + 
        facet_wrap(~phenotype) + labs(x = "Bloques de secuenciacion", y="Porcentaje de genes mitocondriales /n(%)", title = "Contenido mitocondrial"), # Etiquetas en X, Y y titulo
    plotColData(sce.416b_edited, x="block", y="altexps_ERCC_percent", 
        colour_by="discard", other_fields="phenotype") + 
        facet_wrap(~phenotype) + labs(x = "Bloques de secuenciacion", y="Porcentaje de ERCC /n(%)", title = "Contenido de ERCC"), # Etiquetas en X, Y y titulo,
    ncol=1
)
```

Otra visualización gráfica es:

```{r Vizualizacion_agrupada}
plotColData(
    sce.416b_edited,
    x = "sum",
    y = "subsets_Mito_percent",
    colour_by = "discard",
    other_fields = c("block", "phenotype")
) +
    facet_grid(block ~ phenotype)
```

## Identificando droplets vacíos con datos de PBMC

### Conceptos básicos

```{r Reactivo_SingleCell, echo=FALSE, out.width='80%', fig.align='center'}
knitr::include_graphics("img/Reagent_delivery_system.png")
```

Descripción gráfica de la tecnología Next GEM de 10x Genomics. Fuente: [10x Genomics](https://www.10xgenomics.com/technology).

```{r EmptyDrops, echo=FALSE, out.width='80%', fig.align='center'}
knitr::include_graphics("img/emptyDrops_Fig2.png")
```

Opciones algorítmicas para detectar los droplets vacíos. Fuente: [Lun et al, _Genome Biology_, 2019](https://doi.org/10.1186/s13059-019-1662-y).


### Información de los datos

```{r Datos_PMBC, echo=FALSE, out.width='80%', fig.align='center'}
knitr::include_graphics("img/4k_PBMC.png")
```

El dataset proviene de células humanas, con un total de 4,340 células con buena calidad. Para más información consulta la página web en [10Xgenomics](https://www.10xgenomics.com/resources/datasets/4-k-pbm-cs-from-a-healthy-donor-2-standard-2-0-1).

A continuación te presento el reporte de este dataset.

```{r Datos_PMBC_p2, echo=FALSE, out.width='80%', fig.align='center'}
knitr::include_graphics("img/4k_PBMC_summary.png")
```

Captura del reporte de secuenciación del [dataset](https://cf.10xgenomics.com/samples/cell-exp/2.0.1/pbmc4k/pbmc4k_web_summary.html).


10Xgenomics realiza un alineamiento genómico con STAR y emplea CellRanger para realizar la limpieza y cuantificación de los datos. Encontrado la información dividia en dos carpetas importantes:

1) **Gene/cell matrix (raw):** Contiene toda la información posterior al alineamiento, sin limpiar. Contiene toda la información, contemplando la presencia de doplets vacíos.
2) **Gene/cell matrix (filtered):** Contiene los datos filtrados por CellRanger.

Dejo a tu consideración la posibilidad de emplear alguno de estos archivos. Ambas carpetas contienen 3 archivos importantes para el análisis de los datos: 

- `barcode.tsv` = Secuencia que define a cada célula por identificadores de cada célula.
- `feature.tsv` = Genes identificadores, se les midió expresión. 
- `matrix.mtx` = Unión entre las secuencias y genes.


```{r Datos_PMBC_abajo, echo=FALSE, out.width='80%', fig.align='center'}
knitr::include_graphics("img/4k_PBMC_abajo.png")
```

### Importar datos en R

Almacenaremos la información de 10Xgenomics en la memoria Cache de la computadora y posteriormente lo almacenaremos en una variable.

```{r Importar_PBMC}
# Datos crudos, sin procesar
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 Cargar_PBMC}
# Archivos contenidos en Temporales
test_dir <- tempdir()
test_dir <- paste0(test_dir, "/pbmc4k/raw_gene_bc_matrices/GRCh38")
list.files(test_dir) # Encontraremos 3 tipos de archivos: barcodes.tsv, genes.tsv y matrix.mtx  

# Cargar datos
fname <- file.path(tempdir(), "pbmc4k/raw_gene_bc_matrices/GRCh38")
sce.pbmc <- read10xCounts(fname, col.names = TRUE) # Cargar datos de 10Xgenomics de CellRanger
sce.pbmc
```

Contiene 737280 droplets, aún sin filtrar y evaluando un total de 33,694 genes en este dataset.


Si quisieras descargar estos mismos datos pero filtrados, te dejo la ruta:

```{r Cargar_PBMC_v2, eval=FALSE}
# Datos filtrados, sin celulas muertas, sin un alto contenido mitocondrial, sin droplets vacios o con multiples celulas
bfc <- BiocFileCache()
raw.path <- bfcrpath(bfc, file.path(
    "http://cf.10xgenomics.com/samples",
    "cell-exp/2.1.0/pbmc4k/pbmc4k_filtered_gene_bc_matrices.tar.gz"
    
))
untar(raw.path, exdir = file.path(tempdir(), "pbmc4k"))
```

### Visualización de droplets vacíos

```{r Barcode_PBMCs}
bcrank <- barcodeRanks(counts(sce.pbmc))
bcrank
```

```{r Plot_PBMCs}
# Mostremos solo los puntos únicos para acelerar
# el proceso de hacer esta gráfica
uniq <- !duplicated(bcrank$rank)
plot(
    bcrank$rank[uniq],
    bcrank$total[uniq],
    log = "xy",
    xlab = "Rank",
    ylab = "Total UMI count",
    cex.lab = 1.2)


# Agregarle los puntos de corte
abline(
    h = metadata(bcrank)$inflection,
    col = "darkgreen",
    lty = 2)

abline(
    h = metadata(bcrank)$knee,
    col = "dodgerblue",
    lty = 2)

legend(
    "bottomleft",
    legend = c("Inflection", "Knee"),
    col = c("darkgreen", "dodgerblue"),
    lty = 2,
    cex = 1.2)
```

*UMI* =  Secuencia única de cada unión en la bead (unique molecular identifiers). 

Encontremos los droplets vacíos usando `emptyDrops()`. Los siguientes pasos demoran un tiempo.

```{r Eliminar_droplets}
## Usemos DropletUtils para encontrar los droplets
set.seed(100)
e.out <- emptyDrops(counts(sce.pbmc))

# Revisa ?emptyDrops para una explicación de por qué hay valores NA
summary(e.out$FDR <= 0.001)

set.seed(100)

limit <- 100

all.out <- emptyDrops(counts(sce.pbmc), lower = limit, test.ambient = TRUE)

# Idealmente, este histograma debería verse uniforme.
# Picos grandes cerca de cero indican que los _barcodes_
# con un número total de cuentas menor a "lower" no son
# de origen ambiental.
hist(all.out$PValue[all.out$Total <= limit &
                      all.out$Total > 0],
     xlab = "P-value",
     main = "",
     col = "grey80")
```

### Anotación de genes y eliminación de doples vacíos

```{r Anotacion_genes}
# Anotación de los genes
rownames(sce.pbmc) <- uniquifyFeatureNames(
    rowData(sce.pbmc)$ID, rowData(sce.pbmc)$Symbol
)

location <- mapIds(EnsDb.Hsapiens.v86,
    keys = rowData(sce.pbmc)$ID,
    column = "SEQNAME", keytype = "GENEID"
)

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

### Control de calidad

```{r QC_PBMC}
# Obtener las estadisticas en un archivo aparte
stats <- perCellQCMetrics(sce.pbmc,
    subsets = list(Mito = which(location == "MT"))
)
high.mito <- isOutlier(stats$subsets_Mito_percent,
    type = "higher"
)

# Eliminar secuencias con alto contenido de genes mitocondriales
sce.pbmc <- sce.pbmc[, !high.mito]
```

## Visualización de los datos con ISEE

- Github de [iSEE](https://github.com/iSEE/iSEE).
- Workshop de [iSEE](https://isee.github.io/iSEEWorkshop2020/)


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