02_sesion1.Rmd

# Selección de genes altamente variables

José Antonio Ovando

08 de agosto de 2023

## Diapositivas de Peter Hickey

Ver las diapositivas originales [aquí](https://docs.google.com/presentation/d/19J2FyjKlBQdAkku4Oa6UZ6SA-Y4P7AEKCRIbEQWA9ho/edit#slide=id.ga100bba375887aa_0)


## Motivación

* Usualmente usamos datos scRNA-seq para caracterizar la heterogeneidad entre células

* Para hacer esto, usamos métodos como el clustering y la reducción de dimensionalidad

* Esto involucra resumir las diferencias por gen en una sola medida de (dis)similitud entre un par de células

* **¿Cuáles genes deberíamos usar para calcular esta medida de (dis)similitud?**



## Selección de _features_ (genes)

La elección de los _features_ tiene un mayor impacto en qué tan similares decidimos que son las células 

* ➕ _Features_ que contienen información útil biológica 
* ➖ _Features_ que contienen ruido aleatorio
* 👉 Efectos laterales al reducir la dimensionalidad de los datos 

Deseamos seleccionar los **genes altamente variables** (High Variable Genes **HVGs**). Genes con una variación incrementada en comparación con otros genes que están siendo afectados por ruido técnico u otra variación biológica que no es de nuestro interés.

## Dataset ilustrativo: PBMC4k 10X sin filtrar

### Descargar datos

```{r, warning=FALSE, message=FALSE}
# Usemos datos de pbmc4k
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"))

library(DropletUtils)
library(Matrix)
fname <- file.path(tempdir(), "pbmc4k/raw_gene_bc_matrices/GRCh38")
sce.pbmc <- read10xCounts(fname, col.names = TRUE)
sce.pbmc
```

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).]

### Anotación

```{r, warning=FALSE, message=FALSE}
# Anotación de los genes
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"
)

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

### Control de calidad

```{r, warning=FALSE, message=FALSE}
# Control de calidad
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]

# Normalización de los datos
library(scran)
set.seed(1000)
clusters <- quickCluster(sce.pbmc)
sce.pbmc <- computeSumFactors(sce.pbmc, cluster = clusters)
sce.pbmc <- logNormCounts(sce.pbmc)
```


## Dataset ilustrativo: 416B

```{r, warning=FALSE, message=FALSE}
library(scRNAseq)
sce.416b <- LunSpikeInData(which = "416b")
sce.416b$block <- factor(sce.416b$block)
```

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](https://osca.bioconductor.org/lun-416b-cell-line-smart-seq2.html) ^[Lun, A. T. L., Calero-Nieto, F. J., Haim-Vilmovsky, L., Göttgens, B. & Marioni, J. C. Assessing the reliability of spike-in normalization for analyses of single-cell RNA sequencing data. Genome Res. 27, 1795–1806 (2017)).]

```{r, warning=FALSE, message=FALSE}
# gene-annotation
library(AnnotationHub)
ens.mm.v97 <- AnnotationHub()[["AH73905"]]
rowData(sce.416b)$ENSEMBL <- rownames(sce.416b)
rowData(sce.416b)$SYMBOL <- mapIds(ens.mm.v97,
    keys = rownames(sce.416b),
    keytype = "GENEID", column = "SYMBOL"
)
rowData(sce.416b)$SEQNAME <- mapIds(ens.mm.v97,
    keys = rownames(sce.416b),
    keytype = "GENEID", column = "SEQNAME"
)
library(scater)
rownames(sce.416b) <- uniquifyFeatureNames(
    rowData(sce.416b)$ENSEMBL,
    rowData(sce.416b)$SYMBOL
)
```

```{r, warning=FALSE, message=FALSE}
# quality-control
mito <- which(rowData(sce.416b)$SEQNAME == "MT")
stats <- perCellQCMetrics(sce.416b, subsets = list(Mt = mito))
qc <- quickPerCellQC(stats,
    percent_subsets = c("subsets_Mt_percent", "altexps_ERCC_percent"),
    batch = sce.416b$block
)
sce.416b <- sce.416b[, !qc$discard]

# normalization
library(scran)
sce.416b <- computeSumFactors(sce.416b)
sce.416b <- logNormCounts(sce.416b)
```

## Cuantificando la varianza por gen


### Varianza de los _log-counts_

El enfoque más simple para cuantificar la variación _per-feature_ es simplemente calcular la varianza de los _log-counts_

* ➕ Selección del _feature_  basado en los _log-counts_ (que serán usadas en los análisis más adelante)
* ⚠️ La transformación log no logra la estabilización de la varianza perfecta, así que se requiere modelar la relación de la varianza-media de los _features_.


### Enfoque simple

1. Calcular la varianza de los _log-counts_ para cada gen (ignorando grupos experimentales)
2. Ordenar los genes del más-al-menos variable

### Un enfoque más sofisticado

1. Calcular la varianza de los _log-counts_ para cada gen (ignorando grupos experimentales)
2. Modelar la relación entre la media y la varianza de los _log-counts_ para estimar la variación _técnica_
3. Estimar la varianza _biológica_ sustrayendo la varianza _técnica_ de la varianza total
4. Ordenar los genes de la variable de mayor-a-menor biológicamente


### Supuestos

```{r, warning=FALSE, message=FALSE}
# Varianza de las log-counts
library(scran)
dec.pbmc <- modelGeneVar(sce.pbmc)
```

* 🤓 El supuesto es que a cualquier abundancia dada, la abundancia de los perfiles de expresión de la mayoría de los genes están dominados por el ruido aleatorio _técnico_
* 🤓 Por lo consiguiente, una tendencia representa un estimado del ruido técnico como una función de la abundancia
* 🤓 Podemos entonces descomponer la varianza total de cada gen en un componente _técnico_ y uno _biológico_ 
* 🤓 Genes con una gran varianza _biológica_ son considerados interesantes 

### Visualizando la media y varianza

```{r, warning=FALSE, message=FALSE}
# Visualicemos la relación entre la media y la varianza
fit.pbmc <- metadata(dec.pbmc)
plot(fit.pbmc$mean, fit.pbmc$var,
    xlab = "Mean of log-expression",
    ylab = "Variance of log-expression"
)
curve(fit.pbmc$trend(x), col = "dodgerblue", add = TRUE, lwd = 2)
```


### Ordenando genes interesantes

```{r, warning=FALSE, message=FALSE}
# Ordenemos por los genes más interesantes para checar
# los datos
dec.pbmc[order(dec.pbmc$bio, decreasing = TRUE), ]
```



## Coeficiente de variación de las cuentas

El coeficiente de variación de las cuentas al cuadrado (CV<sup>2</sup>) es una alternativa a la varianza de los _log-counts_

* 👉 Se calcula usando las cuentas en lugar de los _log-counts_

* 🤓 CV es el _ratio_ de la desviación estándar a la media y está muy relacionada con el parámetro de _dispersión_ de la distribución binomial negativa usada en edgeR y DESeq2



### Coeficiente de variación

```{r, warning=FALSE, message=FALSE}
# Coeficiente de variación
dec.cv2.pbmc <- modelGeneCV2(sce.pbmc)
```

* 🤓 Modela la relación de la media de la varianza cuando se considera la relevancia de cada gen
* 🤓 Asume que la mayoría de los genes contienen ruido aleatorio y que la tendencia captura la mayoría de la variación técnica 
* 🤓 Genes con un gran CV<sup>2</sup> que se desvían fuertemente de la tendencia es probable que representen genes afectados por la estructura biológica 
* 🤓 Usa la tasa (en lugar de la diferencia) del CV<sup>2</sup> a la tendencia



### Visualizando el coeficiente de variación

```{r, warning=FALSE, message=FALSE}
# Visualicemos la relación con la media
fit.cv2.pbmc <- metadata(dec.cv2.pbmc)
plot(fit.cv2.pbmc$mean, fit.cv2.pbmc$cv2,
    log = "xy"
)
curve(fit.cv2.pbmc$trend(x),
    col = "dodgerblue",
    add = TRUE, lwd = 2
)
```



### Genes por coeficiente de variación

```{r, warning=FALSE, message=FALSE}
# Ordenemos por los genes más interesantes para checar
# los datos
dec.cv2.pbmc[order(dec.cv2.pbmc$ratio,
    decreasing = TRUE
), ]
```

## Varianza de los _log-counts_ vs coeficiente de variación

**Generalmente se usa la varianza de los _log-counts_** 

* Ambas son medidas efectivas para cuantificar la variación en la expresión génica 
* CV<sup>2</sup> tiende a tener otorgar rangos altos a genes altamente variables (HGVs) con bajos niveles de expresión
  - Éstos son dirigidos por una sobreregulación en subpoblaciones raras 
  - Puede asignar un alto rango a genes que no son de nuestro interés con varianza baja absoluta	
* La variación descrita por el CV<sup>2</sup> de las cuentas es menos relevante para los procedimientos que operan en los _log-counts_

## Cuantificando el ruido técnico

* Previamente, ajustamos una línea de tendencia a todos los genes endógenos y asumimos que la mayoría de los genes no están dominados por ruido técnico

* En la práctica, todos los genes expresados tienen algún nivel de variabilidad biológica diferente de cero (e.g., transcriptional bursting)

* Esto sugiere que nuestros estimados de los componentes técnicos estarán inflados probablemente

* 👉 Es mejor que pensemos estos estimados como una variación _técnica_ más la variación biológica que no es interesante

* 🤔 Pero que tal si un grupo de genes a una abundancia particular es afectado por un proceso biológico?
E.g., fuerte sobre regulación de genes específicos de un tipo celular podrían conllevar a un enriquecimiento de HVGs en abundancias altas. Esto inflaría la tendencia y compromete la detección de los genes relevantes

**¿Cómo podemos evitar este problema?**

Podemos revisar dos enfoques:

1. Cuando tenemos spike-ins
2. Cuando no tenemos spike-ins 

### En la presencia de spike-ins

```{r, warning=FALSE, message=FALSE}
dec.spike.416b <- modelGeneVarWithSpikes(
    sce.416b,
    "ERCC"
)
# Ordering by most interesting genes for
# inspection.
dec.spike.416b[order(dec.spike.416b$bio,
    decreasing = TRUE
), ]
```

* 👉 Ajusta la tendencia solo con los spike-ins (que deberían estar afectados solamente por variación técnica)

```{r, warning=FALSE, message=FALSE}
plot(dec.spike.416b$mean, dec.spike.416b$total,
    xlab = "Mean of log-expression",
    ylab = "Variance of log-expression"
)
fit.spike.416b <- metadata(dec.spike.416b)
points(fit.spike.416b$mean, fit.spike.416b$var,
    col = "red", pch = 16
)
curve(fit.spike.416b$trend(x),
    col = "dodgerblue",
    add = TRUE, lwd = 2
)
```

## En la ausencia de spike-ins

```{r, warning=FALSE, message=FALSE}
set.seed(0010101)
dec.pois.pbmc <- modelGeneVarByPoisson(sce.pbmc)
# Ordering by most interesting genes for inspection.
dec.pois.pbmc[order(dec.pois.pbmc$bio, decreasing = TRUE), ]
```

* 👉 Realiza algunas asunciones estadísticas acerca del ruido
* 🤓 Las cuentas UMI típicamente muestran una variación cercana a Poisson si solo consideramos ruido técnico de la preparación de las librerías y la secuenciación 
* ⚠️ modelGeneVarByPoisson() realiza simulaciones, por lo que necesitamos “ajustar la “semilla” para obtener resultados reproducibles
* 🤓 modelGeneVarByPoisson() pueden también simular una variación binomial negativa (variación de Poisson sobredispersada)

```{r, warning=FALSE, message=FALSE}
plot(dec.pois.pbmc$mean, dec.pois.pbmc$total,
    pch = 16, xlab = "Mean of log-expression",
    ylab = "Variance of log-expression"
)
curve(metadata(dec.pois.pbmc)$trend(x),
    col = "dodgerblue", add = TRUE
)
```

* 🤓 La línea de tendencia basada puramente en ruido técnico tiende a producir componentes “biológicos” más grandes por los genes altamente expresados, que frecuentemente incluyen los genes “house-keeping”
* 🤔 Necesitas considerar si tales genes son “interesantes” o no en tu dataset de interés

## Recordemos propiedades de los datos de sce.416b

Este dataset consiste de células de una línea celular de células inmortalizadas mieloides progenitoras de ratón utilizando SmartSeq2

Una cantidad constante de spike-in ERCC RNA se agregó a cada lisado celular antes de la prepatación de la librería

Descripción [aquí](https://osca.bioconductor.org/lun-416b-cell-line-smart-seq2.html)

_Lun, A. T. L., Calero-Nieto, F. J., Haim-Vilmovsky, L., Göttgens, B. & Marioni, J. C. Assessing the reliability of spike-in normalization for analyses of single-cell RNA sequencing data. Genome Res. 27, 1795–1806 (2017)_

## Considerando factores experimentales 

- Los datos que contienen múltiples *batches* muy seguido presentan **efecto de bloque** que pueden crear HGVs artificiales
- Se debe identificar los HGVs en cada *batch* y combinarlos en una única lista de HGVs


```{r }
# calculando la variacion por bloque
dec.block.416b <- modelGeneVarWithSpikes(sce.416b,
    "ERCC",
    block = sce.416b$block
)
dec.block.416b[order(
    dec.block.416b$bio,
    decreasing = TRUE
), ]
```

Al calcular tendencias específicas por batch se tomarán en cuenta las diferencias en la relación media-varianza entre batches

Se deben obtener estimados de los componentes biológico y técnico para cada gene específicos de cada batch, los cuales se promedian entre los batches para crear una única lista de HVGs


```{r echo=FALSE, eval=FALSE, fig.cap="Factor experimental.", out.width = "100%"}
knitr::include_graphics("img/experimental-factor.png")
```



## Seleccionando genes altamante variables (high-variable genes, HVGs)


Hasta ahora hemos ordenado los genes del más al menos **interesantemente variable**

*¿Qué tanto debemos de bajar en la lista para seleccionar nuestros HVGs?*

Para responder esta pregunta debemos tomar en cuenta lo siguiente: elegir un subset más grande:

- Reduce el riesgo de desechar señal biológica
- Incrementa el ruido por la inclusión de genes irrelevantes

Es difícil determinar el balance óptimo porque el rudio en un contexto podría ser una señal útil en otro contexto

Discutiremos algunas estrategias para seleccionar HVGs


### Seleccionando HVGs sobre la métrica de varianza

La estrategia más simple es seleccionar los **top-X genes** con los valores más grandes para la métrica relevante de varianza, *por ejemplo, la varianza biológica más grande calculada con `scran::modelGeneVar()`*

**Pro**: El usuario puede controlar directamente el número de HVGs

**Contra**: ¿Qué valor de X se debe usar?


```{r }
# Works with modelGeneVar() output
hvg.pbmc.var <- getTopHVGs(dec.pbmc, n = 1000)
str(hvg.pbmc.var)

# Works with modelGeneVarWithSpikes() output
hvg.416b.var <- getTopHVGs(dec.spike.416b, n = 1000)
str(hvg.416b.var)

# Also works with modelGeneCV2() but note `var.field`
hvg.pbmc.cv2 <- getTopHVGs(dec.cv2.pbmc,
    var.field = "ratio", n = 1000
)
str(hvg.pbmc.cv2)
```


#### Estrategias para seleccionar X

Asume que, por ejemplo, no más del 5% de los genes están diferencialmente expresados entre las células de nuestra población:

 - *Establece X como el 5% de los genes*
 
Normalmente no conocemos el número de genes diferencialmente expresados desde antes, por lo tanto, solo hemos cambiado un número arbitrario por otro número arbitrario
 
**RECOMENDACIÓN**: Si decides utilizar los top-X HGVs, elige un valor de X y procede con el resto del análisis con la intención de regresar más adelante y probar otros valores, en vez de dedicarle mucho esfuerzo a encontrar el valor óptimo


### Seleccionando HVGs de acuerdo a su significancia estadística

Establece un límite fijo en alguna métrica de significancia estadística. Por ejemplo: algunos de los métodos reportan un p-valor para cada gene, entonces selecciona todos los genes con un p-valor ajustado menor que 0.05

**Recuerda que las pruebas estadísticas siempre dependen del tamaño de la muestra**

**Ventajas:** 
* Fácil de implementar
* Menos predecible que la estrategia de los top-X

**Desventajas:**
* Podría priorizar genes con significancia estadística fuerte en vez de significancia biológica fuerte


```{r }
# Works with modelGeneVar() output
hvg.pbmc.var.2 <- getTopHVGs(dec.pbmc, fdr.threshold = 0.05)
str(hvg.pbmc.var.2)
# Works with modelGeneVarWithSpikes() output
hvg.416b.var.2 <- getTopHVGs(dec.spike.416b,
    fdr.threshold = 0.05
)
str(hvg.416b.var.2)

# Also works with modelGeneCV2() but note `var.field`
hvg.pbmc.cv2.2 <- getTopHVGs(dec.cv2.pbmc,
    var.field = "ratio", fdr.threshold = 0.05
)
str(hvg.pbmc.cv2.2)
```


### Seleccionando genes por arriba de la tendencia media-varianza

Selecciona todos los genes con **una varianza biológica positiva**

Este es un extremo del equilibrio sesgo-varianza que minimiza el sesgo con el costo de maximizar el ruido

Si seguimos esta aproximación, estamos:

* Dándole a la estructura secundaria de la población una oportunidad de manifestarse
* Capturando más ruido, lo cual puede reducir la resolución de poblaciones bien separadas enmascarando la señal secundaria que intentamos preservar

Funciona mejor si tenemos datasets altamente heterogeneos que contienen muchos tipos celulares diferentes


```{r }
# Works with modelGeneVar() output
hvg.pbmc.var.3 <- getTopHVGs(dec.pbmc, var.threshold = 0)
str(hvg.pbmc.var.3)

# Works with modelGeneVarWithSpikes() output
hvg.416b.var.3 <- getTopHVGs(dec.spike.416b,
    var.threshold = 0
)
str(hvg.416b.var.3)

# Also works with modelGeneCV2() but note `var.field` and
# value of `var.threshold`
hvg.pbmc.cv2.3 <- getTopHVGs(dec.cv2.pbmc,
    var.field = "ratio", var.threshold = 1
)
str(hvg.pbmc.cv2.2)
```

### EJERCICIO: Dibujando los HVGs

Para este ejercicio tendrás que repetir la gráfica que muestra la tendencia de la relación media-varianza (ejeX: media de la expresión, ejeY: varianza de la expresión) incluyendo la línea de tendencia obtenida con alguna de las funciones vistas en la primer parte de la clase (modelGeneVar, modelGeneVarWithSpikes, modelGeneCV2). En esta gráfica, deberás colorear los puntos que corresponden a los HGVs obtenidos con algunos de los enfoques revisados

**RESPUESTA**
```{r }
plot(fit.pbmc$mean, fit.pbmc$var,
    xlab = "Mean of log-expression",
    ylab = "Variance of log-expression"
)
points(fit.pbmc$mean[hvg.pbmc.var], fit.pbmc$var[hvg.pbmc.var], col = "orange")
curve(fit.pbmc$trend(x), col = "dodgerblue", add = TRUE, lwd = 2)
```


### Seleccionando genes de interés *a priori*

Una estrategia contundente es usar sets predefinidos de genes de interés. No hay vergüenza en aprovechar el conocimiento biológivo previo

Sin embargo, limita nuestra capacidad de descubrir aspectos nuevos o inesperados de la variación. Por lo tanto, considera esta como una estrategia complementaria a otros tipo de estrategias de selección de HGVs


También podrías eliminar listas pre-definidas de genes:

* Genes de proteínas ribosomales o genes mitocondriales son conocidos por encontrarse dentro de los genes *más variables* y por interferir con análisis posteriores
 
 Sin embargo, tampoco hay que pecar de prevacido, espera a que estos genes demuestren ser problemáticos para removerlos
 

## Poniendo todo junto


```{r }
# Elegimos el 10% de los genes con con componente biologico de variacion mayor
dec.pbmc <- modelGeneVar(sce.pbmc)
chosen <- getTopHVGs(dec.pbmc, prop = 0.1)
str(chosen)
```

Después de esto tenemos varias opciones para imponer nuestra selección de HGVs durante el resto del análisis:

1. Hacer un subset de SCE para quedarnos únicamente con los HGVs
2. Especificar los HGVs en funciones posteriores
3. Magia (altExps)


### Quedándonos sólo con los HGVs

```{r }
sce.pbmc.hvg <- sce.pbmc[chosen, ]
sce.pbmc.hvg
```

**PRO:** Te aseguras de que los métodos posteriores **sólo** usen estos genes para sus cálculos

**CONTRA:** Los genes no-HGVs son eliminados del nuevo objeto *SingleCellExperiment*, lo cual hace menos conveniente la interrogación del dataset completo sobre genes que no son HGVs


### Especificando los HGVs

```{r }
# Example of specifying HVGs in a downstream function
# Performing PCA only on the chosen HVGs.
library(scater)
sce.pbmc <- runPCA(sce.pbmc, subset_row = chosen)
sce.pbmc
```

Mantiene el objeto *SingleCellExperiment* original y especifica los genes para usar en funciones posteriores mediante un argumento adicional como **subset_row**

**PRO:** Es útil si el análisis usa varios conjuntos de HGVs en diferentes pasos

**CONTRA:** Podría ser inconveniente especificar repetidamente el mismo conjunto de HGVs en diferentes pasos


### Witchcraft (Brujería)

```{r }
# Add the full SCE to the subsetted data SCE
altExp(sce.pbmc.hvg, "original") <- sce.pbmc
sce.pbmc.hvg
altExp(sce.pbmc.hvg, "original")
```

Utilizando el sistema de "experimento alternartivo" en la clase *SingleCellExperiment*

**PRO:** Evita algunos problemas a largo plazo cuando el dataset original no está sincronizado con el conjunto filtrado por HVGs

**CONTRA:** Ralentiza todos los análisis subsecuentes


## Resumen y recomendaciones

* Es fácil atorarse en este paso debido a la variedad de opciones disponibles
* Elige un conjunto de HVGs y continúa con el análisis, **recuerda regresar para probar la robustez de tus resultados usando una forma diferente para seleccionar los HVGs**
* Si tienes spike-ins, trata de usarlos. No obstante, recuerda que los spike-ins no pueden capturar todas las fuentes de variación técnica



## Recomendaciones para empezar

Para CEL-Seq2:

- `scran::modelGeneVarWithSpikes()`

Para 10X:

- `scran::modelGeneVarByPoisson()`


Si quieres irte por el lado de conservar demasiados genes:

- `scran::getTopHVGs(dec, var.threshold=0)`

Y realiza una comparación rápida con, por lo menos, el top-1000 HVGs

Regresa al paso de selección de HVG para **eliminar genes problemáticos** tantas veces como sea necesario


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