02_sesion4.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

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

```{r, warning=FALSE, message=FALSE}
library(DropletUtils)
library(Matrix)
fname <- file.path(tempdir(), "pbmc4k/raw_gene_bc_matrices/GRCh38")
sce.pbmc <- read10xCounts(fname, col.names = TRUE) # 737280 cols
```

Dataset "Células mononucleares humanas de sangre periférica" de 10X Genomics (Zheng, G.X.Y. et al, 2017) ^[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)
rowData(sce.pbmc)$location <- mapIds(EnsDb.Hsapiens.v86,
    keys = rowData(sce.pbmc)$ID,
    column = "SEQNAME", keytype = "GENEID"
)

# Information on the chromosome location of the genes

rowData(sce.pbmc)$location[is.na(rowData(sce.pbmc)$location)] <- "NA"
rowData(sce.pbmc)
rowData(sce.pbmc)[rowData(sce.pbmc)$location == "1",]

# 33694 rows y 737280  cols

# Detección de células
set.seed(100)
e.out <- emptyDrops(counts(sce.pbmc))
summary(e.out$FDR <= 0.001) # La mayoría se fue porque las que tienen menos de 100 se consideraron vacías
sce.pbmc <- sce.pbmc[, which(e.out$FDR <= 0.001)] # 4300 cols!!!!!
```

### Control de calidad

```{r, warning=FALSE, message=FALSE}
# Control de calidad
stats <- perCellQCMetrics(sce.pbmc, 
  subsets=list(Mito=which(rowData(sce.pbmc)$location=="MT")))
high.mito <- isOutlier(stats$subsets_Mito_percent, 
                       type="higher")
sce.pbmc <- sce.pbmc[,!high.mito] # 3896 cols 
```

### Normalización de los datos

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

colLabels(sce.pbmc) <- factor(clust)

```

## Motivación - continuación

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

```{r, warning=FALSE, message=FALSE}
# El gen 1 está asociado con el clustering?
plotExpression(sce.pbmc,
    features = rownames(sce.pbmc)[1],
    x = "cluster", colour_by = "cluster"
)
```

```{r, warning=FALSE, message=FALSE}
# El gen 2 está asociado con el clustering?
plotExpression(sce.pbmc,
    features = rownames(sce.pbmc)[2],
    x = "cluster", colour_by = "cluster"
)
```

```{r, warning=FALSE, message=FALSE}
# El gen 2512 está asociado con el clustering?
plotExpression(sce.pbmc,
    features = rownames(sce.pbmc)[2512],
    x = "cluster", colour_by = "cluster"
)
# PTPRC también es conocido como CD45 que antes se llamaba leukocyte common antigen (LCA). Esto sugiere que casi todos los clusters son leucocitos (excepto el 9)
```

```{r, warning=FALSE, message=FALSE}
# El gen CD3E está asociado con el clustering?
plotExpression(sce.pbmc,
    features = "PF4",
    x = "cluster", colour_by = "cluster"
)
# El gen PF4 (platelet factor 4) solo se expresa en el cluster 9, sugiere que ese cluster son plateletas (trombocitos) células que forman coágulos sólidos y frenan hemorragias y cierran heridas
```

```{r, warning=FALSE, message=FALSE}
# El gen CD3E está asociado con el clustering?
plotExpression(sce.pbmc,
    features = "CD3E",
    x = "cluster", colour_by = "cluster"
)
# El gen PF4 (platelet factor 4) solo se expresa en el cluster 9
```

* 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 parece 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) ^[Soneson C, Robinson MD. Bias, robustness and scalability in single-cell differential expression analysis. Nat Methods. 2018 Apr;15(4):255-261. doi: 10.1038/nmeth.4612. Epub 2018 Feb 26. PMID: 29481549.]

* ➕ Las comparaciones pareadas proveen un _log-fold change_ para indicar cuáles clústeres 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, echo=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", "...", "1 vs 18", "2 vs 1")
logFC <- c("1.50", "-0.08", "...", "1.39", "-1.50")
Pval <- c("1.7e-198", "0.11", "...", "2.2e-9", "1.7e-198")

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 CD3E 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() # Para realizar un prueba t-test Welch (pares)
# scran::combineMarkers() # Para realizar las compaaciones que involucran clusters, aquí está la info de que significa any, all y some
```

- 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, warning=FALSE, message=FALSE}
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 al clúster 9 de cualquier otro clúster.

El conjunto de genes con Top <= X es la unión del top X genes (ordenados por p-value)
de cada comparación de pares que involucran al cluster 9. Por ejemplo, el conjunto de todos los genes con valores Top 1 contienen el gen con el p-value mas pequeño de cada comparación. El conjunto de genes Top valor menor o igual a 10 contiene los 10 genes top de cada comparación. Se consolida en un ranking para cada clúster.

## Usando el _log-fold change_

### Sin especificar 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}
# set direction = "up" to only consider upregulated genes as potential markers
markers.pbmc.up <- findMarkers(sce.pbmc, groups=sce.pbmc$cluster,
  test.type="t", direction="up", pval.type="any")
interesting.up <- markers.pbmc.up[[chosen]]
```

### Usando 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) ^[McCarthy DJ, Smyth GK. Testing significance relative to a fold-change threshold is a TREAT. Bioinformatics. 2009 Mar 15;25(6):765-71. doi: 10.1093/bioinformatics/btp053. Epub 2009 Jan 28. PMID: 19176553; PMCID: 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 

```{r, warning=FALSE, message=FALSE, echo=FALSE}
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))
```

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

Considera todos los genes que están diferencialmente expresados en todas las comparaciones pareadas incluyendo al clúster de interés

```{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. 

Un gen solo logrará un bajo nivel de P-value combinado si está fuertemente DE en todas las comparaciones - cuando funciona es muy efectivo - te da un pequeño conjunto de marcadores candidatos. 

Un gen que esté al mismo nivel que los demás clústeres no será detectado. 
No nos serviría para diferenciar en las poblaciones CD4 CD8.

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

# Ni CD4, o CD8 es un gen marcador específico de un cluster en este escenario
```

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

* 👉 Si `pval.type="all"` es muy estricto, entonces `pval.type="any"` es muy generoso
* 🤓 Podemos entonces usar _some_. Aplica la corrección Holm-Bonferroni a los P-values y toma el mejor valor de en medio como el P-value combinado

Prueba la hipótesis nula de que al menos 50% de las comparaciones pareadas no tienen DE
Se aplica un rango para obtener el conjunto de los mejores marcadores

* ⚠️ 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, warning=FALSE, message=FALSE, fig.show='hide'}
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 binomial

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

```{r, warning=FALSE, message=FALSE, echo=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://bioconductor.org/books/3.17/OSCA.advanced/marker-detection-redux.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()`


### Nota sobre la nueva versión en OSCA

```{r, warning=FALSE, message=FALSE, fig.show='hide'}
library(scran)
marker.info <- scoreMarkers(sce.pbmc, colLabels(sce.pbmc))
marker.info
```

Compara cada par de clusters, y calcula scores cuantificando las diferencias en las distribuciones de la expresión entre clusters

```{r, warning=FALSE, message=FALSE, echo=FALSE}
colnames(marker.info[["1"]]) # statistics for cluster 1.
```

```{r, warning=FALSE, message=FALSE, echo=FALSE}
chosen <- marker.info[["1"]]
ordered <- chosen[order(chosen$mean.AUC, decreasing=TRUE),]
head(ordered[,1:4]) # showing basic stats only, for brevity.
```

```{r, warning=FALSE, message=FALSE, echo=FALSE}
library(scater)
plotExpression(sce.pbmc, features=head(rownames(ordered)), 
    x="label", colour_by="label")
```

## Detalles de la sesión de R

```{r 'reproducibility', cache = TRUE, dependson=knitr::all_labels()}
options(width = 120)
sessioninfo::session_info()
```

## Patrocinadores {-}

La presentación fue hecha con el paquete de R [**xaringan**](https://github.com/yihui/xaringan) y configurada con [**xaringanthemer**](https://github.com/gadenbuie/xaringanthemer).

Este curso está basado en el libro [**Orchestrating Single Cell Analysis with Bioconductor**](https://bioconductor.org/books/release/OSCA/) 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](https://www.sib.swiss/raphael-gottardo-group), 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/).

Puedes encontrar los archivos para este taller en  [comunidadbioinfo/cdsb2023](https://github.com/comunidadbioinfo/cdsb2023). 

Descarga los materiales con `usethis::use_course('comunidadbioinfo/cdsb2023')` o revísalos en línea vía [**comunidadbioinfo.github.io/cdsb2023**](http://comunidadbioinfo.github.io/cdsb2023).