10-genes_marcadores.Rmd

# Identificación de genes marcadores

Instructora: [**Yalbi I. Balderas-Martínez**](http://Yalbibalderas.github.io/).

## Diapositivas de Peter Hickey

Ver las diapositivas originales [aquí](https://docs.google.com/presentation/d/18NArqbfeo3Ujo8bcSJTWXirY8N4lEG-T_tFEqlI-uek/edit)

## Motivación 

* Ahora que hemos obtenido los clústeres, nos preguntamos, pero qué son? (e.g. ¿qué tipo celular es el clúster 1?)

* ¿Cuáles genes están dirigiendo el agrupamiento (e.g., ¿cuáles son los genes diferencialmente expresados entre los clústeres 1 y 2?)

* Idea: Mirar las diferencias en los perfiles de expresión de las células de los diferentes clústeres


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

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

### Genes variables

```{r, warning=FALSE, message=FALSE}
## Identificación de genes altamente variables
set.seed(1001)
dec.pbmc <- modelGeneVarByPoisson(sce.pbmc)
top.pbmc <- getTopHVGs(dec.pbmc, prop = 0.1)
```

### Reducción de dimensiones

```{r, warning=FALSE, message=FALSE}
## Reducción de dimensiones
set.seed(10000)
sce.pbmc <- denoisePCA(sce.pbmc,
    subset.row = top.pbmc,
    technical = dec.pbmc
)
set.seed(100000)
sce.pbmc <- runTSNE(sce.pbmc, dimred = "PCA")
set.seed(1000000)
sce.pbmc <- runUMAP(sce.pbmc, dimred = "PCA")
```

### Clustering

```{r, warning=FALSE, message=FALSE}
# clustering
g <- buildSNNGraph(sce.pbmc, k = 10, use.dimred = "PCA")
clust <- igraph::cluster_walktrap(g)$membership
sce.pbmc$cluster <- factor(clust)
```

## Motivación - continuación

¿Algunos de estos genes están asociados con los resultados de _clustering_?

```{r, warning=FALSE, message=FALSE}
# Is gene 1 associated with the clustering?
plotExpression(sce.pbmc,
    features = rownames(sce.pbmc)[1],
    x = "cluster", colour_by = "cluster"
)
```

```{r, warning=FALSE, message=FALSE}
# Is gene 2 associated with the clustering?
plotExpression(sce.pbmc,
    features = rownames(sce.pbmc)[2],
    x = "cluster", colour_by = "cluster"
)
```

```{r, warning=FALSE, message=FALSE}
# Is gene 2512 associated with the clustering?
plotExpression(sce.pbmc,
    features = rownames(sce.pbmc)[2512],
    x = "cluster", colour_by = "cluster"
)
```


```{r, warning=FALSE, message=FALSE}
# Is gene CD3E associated with the clustering?
plotExpression(sce.pbmc,
    features = "CD3E",
    x = "cluster", colour_by = "cluster"
)
```


* Ver una gráfica como una forma de encontrar los genes marcadores obviamente no nos sirve a gran escala

* Necesitamos un método estadístico para identificar estos genes marcadores 

* 👉 La prueba t de Welch es una opción obvia para probar las diferencias en la expresión entre clústeres



## Prueba t modificada de Welch pareada

* ➕ Rápidas y buenas propiedades estadísticas para un gran número de células [(Soneson and Robinson, 2018)](https://www.ncbi.nlm.nih.gov/pubmed/29481549)

* ➕ Las comparaciones pareadas proveen un _log-fold change_ para indicar cuáles clústerse son distinguidos por cada gen 

* 🤔 ¿Por qué no comparar cada clúster con el promedio de todas las otras células?
  - Sensible a la composición poblacional, una subpoblación dominante sola que dirige la selección de los marcadores _top_ para cualquier otro clúster 


## Ejemplo ilustrativo: CD3E como gen marcador en el dataset PBMC4k 10X

### Pruebas pareadas

```{r, warning=FALSE, message=FALSE}
plotExpression(sce.pbmc,
    features = "CD3E",
    x = "cluster", colour_by = "cluster"
)
```

```{r, warning=FALSE, message=FALSE, echo=FALSE}
library(kableExtra)
comparison <- c("1 vs 2", "1 vs 3", "...", "2 vs 1", "...", "18 vs 17")
logFC <- c("1.50", "-0.08", "...", "1.39", "...", "0.11")
Pval <- c("1.7e-198", "0.11", "...", "1.7e-198", "...", "0.46")
paired_tests <- data.frame(comparison, logFC, Pval)
knitr::kable(paired_tests, format = "html") %>%
    kable_styling(
        bootstrap_options = c("striped", "hover"),
        full_width = F,
        font_size = 12,
        position = "left"
    )
```

* K = 18 clústeres
* K!/(K-2)! = 306 comparaciones
* La mitad de ellas son **redundantes**


### Combinando comparaciones del gen CD3E para el clúster 1

"Me interesa saber si el gen CD3 está diferencialmente expresado entre el clúster 1 y .."

- **cualquier** (**any**) otro clúster = P = 1.3 x 10<sup>-205</sup> (Simes adjusted P-value)
- **todos** (**all**) los otros clústeres = P = 0.11 (Berger's intersection-union test)
- **algunos** (**some**) de los otros clústeres = P = 2.0 x 10<sup>-44</sup> (mediana u otro cuantil, Holm-adjusted P-values)



#### Extendiendo a todos los genes

```{r, warning=FALSE, message=FALSE, eval = FALSE}
scran::pairwiseTTests()
scran::combineMarkers()
```

- M = 33,694 genes
- 🤓 K x M = 10,310,364 pruebas
- Comparaciones involucrando clúster 1...
- Comparaciones involucrando clúster ...
- Comparaciones involucrando clúster 18


## Aplicación estándar

Para cada clúster, usar pruebas **t de Welch** para identificar los genes que están diferencialmente expresados entre éste y **cualquier** (**any**) otro clúster 

```{r, warning=FALSE, message=FALSE}
# scran::findMarkers()
library(scran)
markers.pbmc <- findMarkers(sce.pbmc,
    groups = sce.pbmc$cluster,
    test.type = "t", pval.type = "any"
)
```

### Explorando los resultados

```{r, warning=FALSE, message=FALSE}
chosen <- "9"
interesting <- markers.pbmc[[chosen]]
```

```{r, warning=FALSE, message=FALSE}
plotExpression(sce.pbmc, rownames(interesting)[1:4],
    x = "cluster", colour_by = "cluster"
)
```

### Con un heatmap

```{r, fig.height = 10}
best.set <- interesting[interesting$Top <= 6, ]
logFCs <- as.matrix(best.set[, -(1:3)])
colnames(logFCs) <- sub("logFC.", "", colnames(logFCs))
library(pheatmap)
pheatmap(logFCs, breaks = seq(-5, 5, length.out = 101))
```

👉 Usamos el campo _Top_ para identificar un conjunto de genes que distinguen el clúster 9 de cualquier otro clúster


## Usando el _log-fold change_

### Sin espeficiar el lfc

Para cada clúster, usa pruebas t de Welch para identificar los genes que están sobreexpresados entre éste y **cualquier** otro clúster

```{r, warning=FALSE, message=FALSE}
markers.pbmc.up <- findMarkers(sce.pbmc,
    groups = sce.pbmc$cluster,
    test.type = "t", direction = "up", pval.type = "any"
)
interesting.up <- markers.pbmc.up[[chosen]]
```

### Especificando el lfc

Para cada clúster, usa pruebas t de Welch para identificar los genes que están sobreexpresados con un log-fold change (lfc) o al menos 1 entre éste y **cualquier** otro clúster

```{r, warning=FALSE, message=FALSE}
markers.pbmc.up2 <- findMarkers(sce.pbmc,
    groups = sce.pbmc$cluster,
    test.type = "t", direction = "up", lfc = 1, pval.type = "any"
)
interesting.up2 <- markers.pbmc.up2[[chosen]]
```

* 👉 La prueba t también nos permite especificar un _log-fold change_ diferente de cero como la hipótesis nula 
* 🤓 Es más riguroso que simplemente filtrar por _log-fold change_ [TREAT](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2654802/)



### Heatmap

```{r, fig.height = 8}
best.set <- interesting.up2[interesting.up2$Top <= 5, ]
logFCs <- as.matrix(best.set[, -(1:3)])
colnames(logFCs) <- sub("logFC.", "", colnames(logFCs))
pheatmap(logFCs, breaks = seq(-5, 5, length.out = 101))
```

* 👉 Los promedios están más centrados en un conjunto de genes marcadores candidatos que están sobreexpresados en el clúster 9 
* ⚠️ El incremento del rigor no se da sin costo 
* ⚠️ Si el _lfc_ es muy grande podría descartar genes útiles
	- E.g., un gen sobreexpresado en una proporción pequeña de células en un clúster sigue siendo un marcador efectivo si el foco está en la especificidad más que en la sensibilidad 
	

## Encontrando marcadores específicos de clústeres

* 👉 Por defecto, `scran::findMarkers()` dará un alto rango a genes que están DE en cualquier comparación pareada
* 🤔 Quiero genes que son específicos de cada clúster
* 👉 Tú quieres genes que son DE en todas las comparaciones pareadas


Para cada clúster, usa pruebas t de Welch para identificar genes que están sobreexpresados entre éste y **todos** los otros clústeres 

```{r, warning=FALSE, message=FALSE}
markers.pbmc.up3 <- findMarkers(sce.pbmc,
    groups = sce.pbmc$cluster,
    direction = "up", pval.type = "all"
)
interesting.up3 <- markers.pbmc.up3[[chosen]]
```

🤓 Usa una prueba de unión-intersección para combinar los P-values que es el máximo P-value de todas las comparaciones pareadas 


### Pros/cons de los genes marcadores específicos de los clústeres

```{r, warning=FALSE, message=FALSE, echo=FALSE}
library(kableExtra)
Poblacion <- c("DN(CD4-/CD8-)", "CD4+>", "CD8+>", "DP(CD4+/CD8+)")
Expresion_CD4 <- c("No", "Si", "No", "Si")
Expresion_CD8 <- c("No", "No", "Si", "Si")
poblacion <- data.frame(Poblacion, Expresion_CD4, Expresion_CD8)
knitr::kable(poblacion, format = "html") %>%
    kable_styling(
        bootstrap_options = c("striped", "hover"),
        full_width = F,
        font_size = 12,
        position = "left"
    )
```

### findMarkers con pval.type some

Para cada clúster, usa pruebas t de Welch para identificar los genes que están sobreexpresados entre éste y **algunos** de los otros clústers 

```{r, warning=FALSE, message=FALSE}
markers.pbmc.up4 <- findMarkers(sce.pbmc,
    groups = sce.pbmc$cluster,
    direction = "up", pval.type = "some"
)
interesting.up4 <- markers.pbmc.up4[[chosen]]
```

* 👉 Para cuando `pval.type="all"` es muy estricto todavía `pval.type="any"` es muy generoso
* 🤓 Aplica la corrección Holm-Bonferroni a los P-values y toma el mejor valor de en medio como el P-value combinado
* ⚠️ Perderás algunas garantías ofrecidas por los otros métodos 



## Pruebas alternas

### Motivación

La prueba t no es la única forma de comparar dos grupos de mediciones

* 🤔 Quiero una prueba que pueda ser usada perfectamente para distinguir dos clústeres uno del otro
- 👉 Prueba de rangos Wilcoxon

* 🤔 Quiero identificar genes que son expresados más frecuentemente en un clúster que en otro
- 👉 Prueba Binomial



## Prueba de rangos de Wilcoxon

Evalúa directamente la separación entre la distribución de la expresión de los diferentes clústeres

* 🤓 Es proporcional al área bajo la curva (AUC), que es la probabilidad de que una célula al azar de un clúster tenga mayor que expresión que una célula al azar de otro clúster

* 👉 AUCs de 1 o 0 indican que los dos clústeres tienen distribuciones de expresión separadas 

* 🤓 También se conoce como prueba Wilcoxon-Mann-Whitney (WMW)


### findMarkers para Wilcoxon
Para cada clúster, usa la prueba de rangos de Wilcoxon para identificar genes que están sobreexpresados entre éste y **cualquier** otro clúster 

```{r, warning=FALSE, message=FALSE}
markers.pbmc.wmw <- findMarkers(sce.pbmc,
    groups = sce.pbmc$cluster, test.type = "wilcox",
    direction = "up", pval.type = "any"
)
interesting.wmw <- markers.pbmc.wmw[[chosen]]
```

### Heatmap de genes marcadores con Wilcoxon

```{r, fig.height = 8}
best.set <- interesting.wmw[interesting.wmw$Top <= 5, ]
AUCs <- as.matrix(best.set[, -(1:3)])
colnames(AUCs) <- sub("AUC.", "", colnames(AUCs))
pheatmap(AUCs,
    breaks = seq(0, 1, length.out = 21),
    color = viridis::viridis(21)
)
```


### Resumen de la prueba de rangos de Wilcoxon

* ➕ Ofrece directamente la propiedad deseable de un gen marcador (i.e. que el gen distinga perfectamente entre dos clústeres)
* ➕ Es simétrico con respecto a las diferencias en el tamaño de los grupos comparados
* ➖ Es mucho más lento comparado con la prueba t (aunque esto en general no es un problema en la práctica)



## Prueba binomial

* Es una prueba que identifica los genes que difieren en la proporción de células que se expresan entre clústeres 
* Una definición mucho más estricta de genes marcadores 

* 🤓 Convierte la expresión en una medida binaria de presencia/ausencia, por lo que toda la información cuantitativa es ignorada 
- Desde una perspectiva práctica, puede ser más fácil para validar 



### findMarkers para binomial

Para cada clúster, usa la prueba Binomial para identificar genes que están más frecuentemente expresados (sobreexpresados) en comparación con cualquier otro clúster 

```{r, warning=FALSE, message=FALSE}
markers.pbmc.binom <- findMarkers(sce.pbmc,
    groups = sce.pbmc$cluster, test.type = "binom",
    direction = "up", pval.type = "any"
)
interesting.binom <- markers.pbmc.binom[[chosen]]
```

* 🤓 El efecto en el tamaño se reporta como el _log-fold change_ en la proporción de las células que se expresan entre clústeres 
* 👉 _Log-fold changes_ grandes positivos, indican que el gen está más frecuentemente expresado en un clúster comparado con otro

### Visualizando genes marcadores de la prueba bionomial

```{r, warning=FALSE, message=FALSE}
top.genes <- head(rownames(interesting.binom))
plotExpression(sce.pbmc, x = "cluster", features = top.genes)
```



### Resumen de la prueba binomial

* La prueba Binomial no toma en cuenta la normalización

* ➕ Produce genes marcadores que pueden ser más fáciles de validar

* ➖ Ser más estricto puede llevar a la pérdida de buenos marcadores candidatos  



## Métodos de expresión diferencial personalizados

* 🤔 ¿Por qué no usar edgeR/DESeq2/limma-voom u otros métodos personalizados (e.g., MAST)?
* 👉 Claro que puedes! [Checa OSCA](https://osca.bioconductor.org/marker-detection.html#using-custom-de-methods)

* 👉 Pero éstos son tal vez algo exagerados para identificar genes marcadores
* 🤓 Las células son nuestras "réplicas" para el propósito de identificar genes marcadores
* 🤓 edgeR/DESeq2/limma-voom hacen asunciones más fuertes acerca de los datos que es más probable que no se cumplan para células individuales en scRNA-seq



## Problemas estadísticos


### Invalidez de P-values

Todas las estrategias de DE para detectar genes marcadores entre clústeres son estadísticamente defectuosas de alguna manera

* 🤓 "Dragado de datos": El análisis DE se realiza usando los mismos datos usados para obtener los clústeres 

* 👉 Las pruebas para genes DE entre clústeres producirá inevitablemente algunos resultados significativos y así es como los clústeres serán definidos!

* 👉 Aún cuando los P-values son defectuosos, el efecto no es muy dañino para la detección de genes ya que los P-values solo son usados para los rangos

* 🤓 No se pueden usar P-values para definir "diferencias significativas" entre los clústeres con respecto a un umbral de la tasa de error 



### Naturaleza de la replicación

👉 Idealmente, validar algunos de los marcadores con una población de células independientes (y idealmente usando una técnica diferente, e.g., hibridación fluorescente _in situ_ o qPCR)


### Comentarios adicionales

* 👉 La estrategia de análisis DE es que los marcadores son definidos relativo a subpoblaciones en el mismo dataset 
* 👉 Si un gen se expresa uniformemente a través de la población no servirá como un marcador
- e.g., los marcadores de las células T no serán detectados si solamente hay células T en los datos
- usualmente no es un problema, ya que tenemos idea de las células que se capturaron
* 👉 Existen métodos de machine learning para hacer la identificación de los genes marcadores, pero la humilde prueba t sigue siendo muy buena


## Resumen y recomendaciones

* 👉 Crea múltiples listas de genes marcadores con diferentes niveles de rigor
* 👉 La forma más simple de interpretar los genes marcadores es que son los sobreexpresados de "forma única", o son "genes específicos de clústeres", especialmente si queremos imponer un _log-fold change_ mínimo
* 👉 Puedes requerir hacer una identificación de genes marcadores más enfocada, e.g., subset de los datos de solo 2 clústeres de interés y entonces correr `scran::findMarkers()`


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