05-control_calidad.Rmd

# Control de calidad

Instructor: [Leonardo Collado Torres](https://comunidadbioinfo.github.io/es/authors/lcollado/)

```{r cargar_paquetes, message = FALSE}
## Paquetes de este capítulo
library("scRNAseq") ## para descargar datos de ejemplo
library("AnnotationHub") ## para obtener información de genes
library("scater") ## para gráficas y control de calidad
library("BiocFileCache") ## para descargar datos
library("DropletUtils") ## para detectar droplets
library("Matrix") ## para leer datos en formatos comprimidos
```

## Diapositivas de Peter Hickey

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

## Ejercicio: entendiendo addPerCellQC

```{r datos_ejercicio_1}
## Datos
library("scRNAseq")
sce.416b <- LunSpikeInData(which = "416b")
sce.416b$block <- factor(sce.416b$block)

# Descarga los archivos de anotación de la base de datos de Ensembl
# correspondientes usando los recursos disponibles vía AnnotationHub
library("AnnotationHub")
ah <- AnnotationHub()
query(ah, c("Mus musculus", "Ensembl", "v97"))
# Obtén la posición del cromosoma para cada gen
ens.mm.v97 <- ah[["AH73905"]]
location <- mapIds(
    ens.mm.v97,
    keys = rownames(sce.416b),
    keytype = "GENEID",
    column = "SEQNAME"
)

# Identifica los genes mitocondriales
is.mito <- which(location == "MT")
library("scater")
sce.416b <- addPerCellQC(sce.416b,
    subsets = list(Mito = is.mito)
)
## Si quieres guarda los resultados de addPerCellQC() para responder
## las preguntas del ejercicio. Eventualmente si necesitaremos los
## resultados de addPerCellQC() para las secciones posteriores a este
## ejercicio.
```


* ¿Qué cambió en nuestro objeto `sce` después de `addPerCellQC`? ^[Ahora tenemos más información en `colData(sce.416b)`]
* Haz una gráfica de _boxplots_ del número de genes por bloque (block) de células. ^[`with(colData(sce.416b), boxplot(detected ~ block))`]

## Gráficas sobre medidas de control de calidad (QC)

```{r visualizar_qc}
plotColData(sce.416b, x = "block", y = "detected")
plotColData(sce.416b, x = "block", y = "detected") +
    scale_y_log10()
plotColData(sce.416b,
    x = "block",
    y = "detected",
    other_fields = "phenotype"
) +
    scale_y_log10() +
    facet_wrap(~phenotype)
```

### Ejercicio: gráficas QC ERCC

Adapta el código de las gráficas anteriores para otra variable de control de calidad. Por ejemplo, escribe el código para reproducir las siguientes gráficas.

```{r qc_altexps_ERCC_percent, echo = FALSE}
plotColData(sce.416b, x = "block", y = "altexps_ERCC_percent")
plotColData(sce.416b,
    x = "block",
    y = "altexps_ERCC_percent",
    other_fields = "phenotype"
) +
    facet_wrap(~phenotype)
```

* Basado en las gráficas encuentra la variable de `colData(sce.416b)` que contiene la información que queremos gráficar.
* ¡No hay que reemplazar todo lo que diga `phenotype`
* Tengo cuidado con las transformaciones de valores en el eje Y. No aplican para todo tipo de datos.

## Eliminar células de baja calidad

```{r valores_qc}
# Valores de límite ejemplo
qc.lib <- sce.416b$sum < 100000
qc.nexprs <- sce.416b$detected < 5000
qc.spike <- sce.416b$altexps_ERCC_percent > 10
qc.mito <- sce.416b$subsets_Mito_percent > 10
discard <- qc.lib | qc.nexprs | qc.spike | qc.mito

# Obtenemos un resumen del 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)
)

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

## Más pruebas
plotColData(sce.416b,
    x = "block",
    y = "detected",
    other_fields = "phenotype"
) +
    scale_y_log10() +
    facet_wrap(~phenotype)

## 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")
attr(qc.nexprs3, "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)
)
```

## Ejercicio: filtrado de células

* ¿Fue necesario `qc.lib` para crear `discard`? ^[Sí, usando `table(qc.lib , qc.spike)` y `table(qc.lib , qc.mito)`.]
* ¿Cúal filtro fue más estricto? ¿`discard` o `discard2`? ^[`discard` de `table(discard, discard2)`]
* Al considerar el grupo de cada muestra (batch), ¿descartamos más células usando un valor de límite automático? ^[Sí, usando `table(discard, discard2, discard3)`]

## Datos de Grun et al

¿Qué patrón revela esta gráfica?

```{r grun_problema}
sce.grun <- GrunPancreasData()
sce.grun <- addPerCellQC(sce.grun)

## ¿Qué patrón revela esta gráfica?
plotColData(sce.grun, x = "donor", y = "altexps_ERCC_percent")
```

¿Cúal de las siguientes gráficas identifica mejor las células de baja calidad?

```{r grun_isOutlier}
## isOutlier() puede ayudarnos cuando un grupo de muestras
## tuvo más problemas que el resto
discard.ercc <- isOutlier(sce.grun$altexps_ERCC_percent,
    type = "higher",
    batch = sce.grun$donor
)
discard.ercc2 <- isOutlier(
    sce.grun$altexps_ERCC_percent,
    type = "higher",
    batch = sce.grun$donor,
    subset = sce.grun$donor %in% c("D17", "D2", "D7")
)

## isOutlier() tomando en cuenta el batch
plotColData(
    sce.grun,
    x = "donor",
    y = "altexps_ERCC_percent",
    colour_by = data.frame(discard = discard.ercc)
)

## isOutlier() tomando en cuenta batch y muestras que fallaron
plotColData(
    sce.grun,
    x = "donor",
    y = "altexps_ERCC_percent",
    colour_by = data.frame(discard = discard.ercc2)
)
```

## Gráficas de QC extra

Otras gráficas que podemos hacer.

```{r qc_extra_416b}
# Agregamos información sobre que células
# tienen valores extremos
sce.416b$discard <- discard2

# Haz esta gráfica para cada medida de
# control de calidad (QC)
plotColData(
    sce.416b,
    x = "block",
    y = "sum",
    colour_by = "discard",
    other_fields = "phenotype"
) +
    facet_wrap(~phenotype) +
    scale_y_log10()

# Otra gráfica de diagnóstico útil
plotColData(
    sce.416b,
    x = "sum",
    y = "subsets_Mito_percent",
    colour_by = "discard",
    other_fields = c("block", "phenotype")
) +
    facet_grid(block ~ phenotype)
```

## Ejercicio: ERCC Grun et al

Adapta el código de `sce.416b` para los datos de Grun et al y reproduce la imagen siguiente.

```{r qc_extra_grun, echo = FALSE}
sce.grun$discard <- discard.ercc2
plotColData(
    sce.grun,
    x = "altexps_ERCC_detected",
    y = "altexps_ERCC_sum",
    colour_by = "discard",
    other_fields = c("donor")
) +
    facet_grid(~donor)
```

* Fíjate en que variables de `colData()` estamos graficando.
* ¿Existe la variable `discard` en `colData()`?
* ¿Qué variable tiene valores de D10, D17, D2, D3 y D7?

## Identificando droplets vacíos con datos de PBMC


```{r echo=FALSE, fig.cap="Descripción gráfica la tecnología _Next GEM_ de 10x Genomics. Fuente: [10x Genomics](https://www.10xgenomics.com/technology)."}
knitr::include_graphics("https://cdn.10xgenomics.com/image/upload/dpr_2.0,e_sharpen,f_auto,q_auto/v1607106030/Reagent_delivery_system.png")
```

```{r echo=FALSE, fig.cap="Opciones algorítmicas para detecar los droplets vacíos. Fuente: [Lun et al, _Genome Biology_, 2019](https://doi.org/10.1186/s13059-019-1662-y)."}
knitr::include_graphics("img/emptyDrops_Fig2.png")
```


```{r pbmc_qc}
## Descarguemos los datos
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"))

## Leamos los datos en R
library("DropletUtils")
library("Matrix")
fname <- file.path(tempdir(), "pbmc4k/raw_gene_bc_matrices/GRCh38")
sce.pbmc <- read10xCounts(fname, col.names = TRUE)
bcrank <- barcodeRanks(counts(sce.pbmc))

# 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
)
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
)
```
Encontremos los _droplets_ vacíos usando `emptyDrops()`.

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

# Revisa ?emptyDrops para una explicación de porque 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"
)
```

## Ejercicio: detección de droplets vacíos

* ¿Por qué `emptyDrops()` regresa valores `NA`? ^[Debajo de `lower` son considerados _droplets_ vacíos. Solo se usan para la correción estadística de pruebas múltiples.]
* ¿Los valores p son iguales entre `e.out` y `all.out`? ^[No, debido a los `NA`s.]
* ¿Son iguales si obtienes el subconjunto de valores que no son `NA`? ^[Sí: `identical(e.out$PValue[!is.na(e.out$FDR)], all.out$PValue[!is.na(e.out$FDR)])`.]

## Filtrado de expresión mitocondrial adicional

Después de filtar los droplets, el filtrado por expresión mitocondrial nos va a ayudar a eliminar células de baja calidad.

```{r pbmc_chrMT_ayuda}
sce.pbmc <- sce.pbmc[, which(e.out$FDR <= 0.001)]
is.mito <- grep("^MT-", rowData(sce.pbmc)$Symbol)
sce.pbmc <- addPerCellQC(sce.pbmc, subsets = list(MT = is.mito))
discard.mito <-
    isOutlier(sce.pbmc$subsets_MT_percent, type = "higher")
plot(
    sce.pbmc$sum,
    sce.pbmc$subsets_MT_percent,
    log = "x",
    xlab = "Total count",
    ylab = "Mitochondrial %"
)
abline(h = attr(discard.mito, "thresholds")["higher"], col = "red")
```

## Ejercicio avanzado

Volvamos a crear `sce.pbmc` para poder usar `plotColData()` y visualizar la relación entre `total` y los niveles de expresión mitocondrial (en porcentaje) separando lo que pensamos que son droplets vacíos y las células de acuerdo a los resultados que ya calculamos de `emptyDrops()`. El resultado final se verá como en la siguiente imagen.

```{r pbmc_combined, echo = FALSE}
## Leer los datos crudos de nuevo
sce.pbmc <- read10xCounts(fname, col.names = TRUE)

## Almacenar si es o no es célula
sce.pbmc$is_cell <- e.out$FDR <= 0.001

## Filtar los casos para los cuales no tuvimos pruebas estadísticas
sce.pbmc <- sce.pbmc[, which(!is.na(sce.pbmc$is_cell))]

## Identicar chr mitocondrial
is.mito <- grep("^MT-", rowData(sce.pbmc)$Symbol)

## Calcular medidas de QC
sce.pbmc <- addPerCellQC(sce.pbmc, subsets = list(MT = is.mito))

## Almacenar cuales hay que descartar por MT
sce.pbmc$discard_mito <-
    isOutlier(sce.pbmc$subsets_MT_percent, type = "higher")

## Visualizar resultados
plotColData(
    sce.pbmc,
    x = "total",
    y = "subsets_MT_percent",
    colour_by = "discard_mito",
    other_fields = c("is_cell")
) + facet_grid(~ ifelse(is_cell, "Célula", "Droplet vacío"))
```

* No podemos usar nuestro objeto `sce.pbmc` porque ya eliminamos los droplets vacíos al correr `sce.pbmc <- sce.pbmc[, which(e.out$FDR <= 0.001)]`. Por eso tendremos que volver a usar `sce.pbmc <- read10xCounts(fname, col.names = TRUE)`.
* Una vez que hayamos vuelto a hacer `sce.pbmc`, tenemos que guardar en ese objeto los resultados de `emptyDrops()`. Por ejemplo, con `sce.pbmc$is_cell <- e.out$FDR <= 0.001`.
* Como `e.out$FDR` tiene muchos `NA`, nos conviene filtrar esos datos.
* Tendremos que volver a correr `addPerCellQC()` y guardar los resultados en nuestro objeto `sce.pbmc`.
* Al final usaremos `plotColData()` junto con `facet_grid(~ sce.pbmc$is_cell)`.

## Discusión ¿Conviene eliminar datos?


```{r filtrar_o_no}
# Eliminemos las células de calidad baja
# al quedarnos con las columnas del objeto sce que NO
# queremos descartar (eso hace el !)
filtered <- sce.416b[, !discard2]
# Alternativamente, podemos marcar
# las células de baja calidad
marked <- sce.416b
marked$discard <- discard2
```


* ¿Cúal de estos objetos es más grande? ^[`marked` es más grande que `filtered`]
* ¿Cúal prefieres usar? ^[Yo prefiero usar `marked` si tengo suficiente memoria para usarlo.]

### Un nuevo paquete: ExperimentSubset

En [BioC2021](https://bioc2021.bioconductor.org/) presentaron [`ExperimentSubset`](http://bioconductor.org/packages/ExperimentSubset/) que provee otro camino para resolver este dilema.

```{r echo=FALSE, fig.cap="Descripción gráfica de `ExperimentSubset`. Fuente: [vignette `ExperimentSubset`](http://bioconductor.org/packages/release/bioc/vignettes/ExperimentSubset/inst/doc/ExperimentSubset.html)."}
knitr::include_graphics("img/ExperimentSubset.png")
```

## Explorando datos de forma interactiva con iSEE

<blockquote class="twitter-tweet"><p lang="en" dir="ltr"><a href="https://twitter.com/hashtag/rstats?src=hash&amp;ref_src=twsrc%5Etfw">#rstats</a> / <a href="https://twitter.com/Bioconductor?ref_src=twsrc%5Etfw">@Bioconductor</a> congrats winners of the 1st Shiny Contest: iSEE <a href="https://t.co/oHgGkWqRsJ">https://t.co/oHgGkWqRsJ</a> <a href="https://t.co/vZLFvcMBIS">https://t.co/vZLFvcMBIS</a> !</p>&mdash; Bioconductor (@Bioconductor) <a href="https://twitter.com/Bioconductor/status/1114773873537449984?ref_src=twsrc%5Etfw">April 7, 2019</a></blockquote> <script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>

* http://bioconductor.org/packages/release/bioc/html/iSEE.html
* http://bioconductor.org/packages/release/bioc/vignettes/iSEE/inst/doc/basic.html

```{r isee_basic}
## Hagamos un objeto sencillo de tipo RangedSummarizedExperiment
library("SummarizedExperiment")
## ?SummarizedExperiment

## De los ejemplos en la ayuda oficial
## Creamos los datos para nuestro objeto de tipo SummarizedExperiment
## para 200 genes a lo largo de 6 muestras
nrows <- 200
ncols <- 6

## Números al azar de cuentas
set.seed(20210223)
counts <- matrix(runif(nrows * ncols, 1, 1e4), nrows)

## Información de nuestros genes
rowRanges <- GRanges(
    rep(c("chr1", "chr2"), c(50, 150)),
    IRanges(floor(runif(200, 1e5, 1e6)), width = 100),
    strand = sample(c("+", "-"), 200, TRUE),
    feature_id = sprintf("ID%03d", 1:200)
)
names(rowRanges) <- paste0("gene_", seq_len(length(rowRanges)))

## Información de nuestras muestras
colData <- DataFrame(
    Treatment = rep(c("ChIP", "Input"), 3),
    row.names = LETTERS[1:6]
)

## Juntamos ahora toda la información en un solo objeto de R
rse <- SummarizedExperiment(
    assays = SimpleList(counts = counts),
    rowRanges = rowRanges,
    colData = colData
)

## Exploremos el objeto resultante
rse

## Explora el objeto rse de forma interactiva
library("iSEE")
if (interactive()) {
    iSEE::iSEE(rse)
}
```

### Ejercicio iSEE con sce.416b

Repitamos la imagen que hicimos anteriormente.

<a href="https://comunidadbioinfo.github.io/cdsb2021_scRNAseq/control-de-calidad.html#gr%C3%A1ficas-sobre-medidas-de-control-de-calidad-qc"><img src="https://comunidadbioinfo.github.io/cdsb2021_scRNAseq/05-control_calidad_files/figure-html/visualizar_qc-3.png"/></a>

```{r isee_416b}
## Explora el objeto sce.416b de forma interactiva
if (interactive()) {
    iSEE::iSEE(sce.416b, appTitle = "sce.416b")
}
```

### Datos de LIBD de Tran et al

<a href="https://libd.shinyapps.io/tran2021_AMY/"><img src="https://raw.githubusercontent.com/LieberInstitute/10xPilot_snRNAseq-human/master/screenshot_tran2021_AMY.png"></a>

* Datos de [Tran et al, _bioRxiv_, 2020](https://github.com/LieberInstitute/10xPilot_snRNAseq-human#explore-the-data-interactively).

* Código de R para el sitio web: https://github.com/LieberInstitute/10xPilot_snRNAseq-human/tree/master/shiny_apps/tran2021_AMY. 

### Más detalles de iSEE

<iframe width="560" height="315" src="https://www.youtube.com/embed/bK8D30MqXb8" title="YouTube video player" frameborder="0" allow="accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>

* [Notas en inglés](https://docs.google.com/document/d/19ZMeaFiFqhJUJJHiEj0qsCcH5WMrn32H-PCKzWag1oI/edit).

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