02_sesion5.Rmd

# Anotación de tipos celulares

Erick Cuevas

08 de agosto de 2023

## Diapositivas

[
```{r,echo=FALSE}
knitr::include_url("https://erickcufe.github.io/anotacion_scRNA/anotacion_celulas.html", height = "380px")
```
](<https://erickcufe.github.io/anotacion_scRNA/anotacion_celulas.html>)

<style>

div.color { 
border-radius: 5px; 
padding: 20px;
margin: 30px 0px 30px;}

div.red { background-color:#f67155; }

div.orange{ background-color:#f0BB51;}

div.pair { 	
display: flex;
flex-direction: row;
justify-content: center;
text-align:center;
padding:0px}

div.inside { width: 49%; padding: 0px}

div.scroll { 
max-height: 400px; 
overflow-y: auto; 
background: #111111;
border-radius:5px;
padding: 10px;
margin: 30px 0px 30px;
color: #999999;}

div.alert{color:#bd475d; background-color:transparent}
</style>



## Introducción

El análisis de secuenciación de ARN a nivel de célula única (scRNAseq) ha revolucionado nuestra capacidad para estudiar la heterogeneidad celular en diversos tejidos y condiciones. A diferencia de la secuenciación de ARN tradicional, que proporciona información promedio de todas las células en una muestra, el scRNAseq permite el estudio detallado de la expresión génica en células individuales. Esto ha abierto puertas a descubrimientos sin precedentes en biología celular, permitiendo identificar nuevos tipos celulares, estados transicionales y vías de señalización específicas de células.

Dentro de este contexto, la anotación de clusters en scRNAseq es esencial. Una vez que se han identificado grupos o "clusters" de células con perfiles de expresión similares, es crucial determinar qué representan estas agrupaciones en términos biológicos. Es aquí donde entra en juego la anotación.

## Motivación

-   **Identificación de Tipos y Subtipos Celulares:** La principal motivación para anotar clusters es identificar y etiquetar tipos y subtipos celulares específicos dentro de una muestra. Esto es fundamental para entender la composición celular de un tejido y cómo esta composición puede cambiar en diferentes condiciones, como en enfermedades.
-   **Entender la Función Celular**: Al anotar clusters, no solo identificamos qué tipo de célula es, sino que también podemos inferir su función actual basándonos en su perfil de expresión génica.
-   **Descubrimiento de Nuevos Tipos Celulares:** En algunos casos, la anotación puede revelar grupos de células que no se ajustan a tipos celulares conocidos, lo que indica la posible existencia de un tipo celular previamente no descrito.
-   **Base para Análisis Posteriores:** Una vez que se han anotado los clusters, se pueden realizar análisis más detallados, como estudios de vías de señalización, análisis de reguladores maestros o estudios de interacción celular.
-   **Comparación Entre Muestras y Condiciones:** La anotación permite comparar la composición celular entre diferentes muestras, tejidos o condiciones, lo que es esencial para estudios comparativos y de enfermedades.

## Aproximaciones para anotar

1. Modo artístico. Usando conocimiento previo de genes marcadores ya publicados.
2. Usando referencias de conjuntos de datos ya anotados.
3. Combinando el modo artístico con referencias. 
4. Anclas con Seurat.

## Paqueterías de R más "famosas" para anotar

- [SingleR](https://bioconductor.org/packages/release/bioc/vignettes/SingleR/inst/doc/SingleR.html)
- [Seurat](https://satijalab.org/seurat/)
- [scCATCH](https://github.com/ZJUFanLab/scCATCH)
- [cellassign](https://github.com/Irrationone/cellassign)
- [SCINA](https://github.com/jcao89757/SCINA)
- [Garnett](https://cole-trapnell-lab.github.io/garnett/)

## Paqueterías de R con conjuntos de datos de referencia

- [SeuratData](https://satijalab.org/seurat/)
- [celldex](http://bioconductor.org/packages/release/data/experiment/vignettes/celldex/inst/doc/userguide.html)
- [scRNAseq](https://bioconductor.org/packages/3.17/data/experiment/vignettes/scRNAseq/inst/doc/scRNAseq.html)

Otros sitios interesantes para obtener conjuntos de datos de referencia: 

- [Azimuth](https://azimuth.hubmapconsortium.org/)
- [Single Cell Expression Atlas](https://www.ebi.ac.uk/gxa/sc/home)
- [GLIASEQ](https://www.liddelowlab.com/gliaseq)


## Preparación del dataset

El dataset tiene el ID en GEO [GSE159677](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE159677)


<div class="alert">

**OJO** PARA FINES DEL TIEMPO DEL CURSO POR AHORA NO GENERAREMOS este dataset desde cero, este proceso realízalo cuando tengas memoria ram libre y no utilices tu equipo para otra cosa

</div> 

```{r, eval=FALSE, echo=TRUE}
library(cowplot)
library(ggplot2)
library(scater)
library(scds)
library(SingleCellExperiment)

#Carga de datos y re formato
#Carga de las cuentas crudas
#Recuerda cambiar el PATH o directorio a la ubicación de los datos
fastq_dirs <- list.dirs("../datos_atherosclerosis/tom_alsaigh_2022/", recursive = FALSE, full.names = TRUE)
names(fastq_dirs) <- basename(fastq_dirs)
sce <- DropletUtils::read10xCounts(fastq_dirs)


#Renombrar row/colData nombre de columnas y SCE dimnames
names(rowData(sce)) <- c("ENSEMBL", "SYMBOL")
names(colData(sce)) <- c("sample_id", "barcode")
sce$sample_id <- factor(basename(sce$sample_id))
dimnames(sce) <- list(
  with(rowData(sce), paste(SYMBOL, sep = ".")),
  with(colData(sce), paste(barcode, sample_id, sep = ".")))


# Agregamos la metadata

md <- data.frame(
  Sample_ID = c("Patient1_AC", "Patient1_PA",
                "Patient2_AC", "Patient2_PA",
                "Patient3_AC", "Patient3_PA"),
  Age = c(82, 82, 87, 87, 65, 65),
  AHA_Classification = c("Type_VII_Calcified", "Type_VII_Calcified",
                         "Type_VII_Calcified", "Type_VII_Calcified",
                         "Type_VII_Calcified", "Type_VII_Calcified"),
  Region = c("Atherosclerotic_Core", "Proximal_Adjacent",
             "Atherosclerotic_Core", "Proximal_Adjacent",
             "Atherosclerotic_Core", "Proximal_Adjacent")
)


m <- match(sce$sample_id, md$Sample_ID)
sce$region <- md$Region[m]
sce$AHA_Classification <- md$AHA_Classification[m]
sce$sample_id <- md$Sample_ID[m]
sce$Age <- md$Age[m]

```

Hacemos ahora un pequeño proceso de limpieza y normalización

```{r, eval=FALSE, echo=TRUE}
#Remover genes no detectados
sce <- sce[Matrix::rowSums(counts(sce) > 0) > 0, ]
dim(sce)

#Remover doublets
#dividir SCE por muestra OJO esto podria ser por lotes
cs_by_s <- split(colnames(sce), sce$sample_id)
sce_by_s_2 <- lapply(cs_by_s, function(cs) sce[, cs])

#correr 'scds'
sce_by_s_2 <- lapply(sce_by_s_2, function(u)
  scds::cxds_bcds_hybrid(scds::bcds(scds::cxds(u))))

#Eliminando doublets
sce_by_s <- lapply(sce_by_s_2, function(u) {
  #Calcula numero de doublets (10x)
  n_dbl <- ceiling(0.01 * ncol(u)^2 / 1e3)
  #Elimina 'n_dbl' celulas con mayor scpre de doublet
  o <- order(u$hybrid_score, decreasing = TRUE)
  u[, -o[seq_len(n_dbl)]]
})

#Integrar nuevamente al objeto SCE
sce <- do.call(cbind, sce_by_s)

#Calcular metricas QC 
(mito <- grep("MT-", rownames(sce), value = TRUE))
# sce <- perCellQCMetrics(sce, subsets = list(Mito=mito))

sce <- addPerCellQCMetrics(sce, subsets = list(Mito=mito))
sce <- addPerFeatureQCMetrics(sce)

```

Filtración de células

```{r, eval=FALSE, echo=TRUE}
#Obtener outliers
cols <- c("total", "detected", "subsets_Mito_percent")
log <- c(TRUE, TRUE, FALSE)
type <- c("both", "both", "higher")

# Con esto decidimos que barcode desechar
drop_cols <- paste0(cols, "_drop")
for (i in seq_along(cols)){
  colData(sce)[[drop_cols[i]]] <- isOutlier(sce[[cols[i]]],
                                            nmads = 2.5, type = type[i],
                                            log = log[i], batch = sce$sample_id)
}

# Muestra un resumen de los barcodes que se eliminaran
sapply(drop_cols, function(i)
  sapply(drop_cols, function(j)
    sum(sce[[i]] & sce[[j]])))

cd <- data.frame(colData(sce))
ps <- lapply(seq_along(cols), function (i) {
  p <- ggplot(cd, aes_string(x = cols[i], alpha = drop_cols[i])) +
    geom_histogram(bins = 100, show.legend = FALSE) +
    scale_alpha_manual(values = c("FALSE" = 1, "TRUE" = 0.4)) +
    facet_wrap(~sample_id, ncol = 1, scales = "free") +
    theme_classic() + theme(strip.background = element_blank())
  if (log[i])
    p <- p + scale_x_log10()
  return(p)
})
plot_grid(plotlist = ps, ncol = 3)

layout(matrix(1:2, nrow = 1))
ol <- Matrix::rowSums(as.matrix(colData(sce)[drop_cols])) != 0
x <- sce$total
y <- sce$detected
LSD::heatscatter(x, y, log="xy", main = "unfiltered",
                 xlab = "Total counts", ylab = "Non-zero features")
LSD::heatscatter(x[!ol], y[!ol], log="xy", main = "filtered",
                 xlab = "Total counts", ylab = "Non-zero features")

#Generar resumen de celulas a mantener
ns <- table(sce$sample_id)
ns_fil <- table(sce$sample_id[!ol])
print(rbind(
  unfiltered = ns, filtered = ns_fil,
  "%" = ns_fil / ns * 100))

#Eliminar celulas outlier
sce <- sce[, !ol]
dim(sce)

#count > 1
sce <- sce[Matrix::rowSums(counts(sce) > 1) >= 20, ]
dim(sce)

```

Agrupamiento, toma en cuenta que esta parte podría demorar bastante.

```{r, eval=FALSE, echo=TRUE}
#Load packages
library(cowplot)
library(Seurat)
library(SingleCellExperiment)
library(ggplot2)

#INTEGRATE
#Crear SeuratObject
so <- CreateSeuratObject(
  counts = counts(sce),
  meta.data = data.frame(colData(sce)),
  project = "Alsaigh_10x_data")

#Dividir por muestra
cells_by_sample <- split(colnames(sce), sce$sample_id)
so <- lapply(cells_by_sample, function(i)
  subset(so, cells = i))

#Normalizar, encontrar genes variables, escalar
so <- lapply(so, NormalizeData, verbose = FALSE)
so <- lapply(so, FindVariableFeatures, nfeatures = 2e3,
             selection.method = "vst", do.plot = FALSE, verbose = FALSE)
so <- lapply(so, ScaleData, verbose = FALSE)

#Encontrar anclas
# Estas anclas nos ayudaran despues para integrar todo con la funcion IntegrateData.
as <- FindIntegrationAnchors(so, verbose = TRUE)
so <- IntegrateData(anchorset = as, dims = seq_len(30), verbose = TRUE)

#Escalar datos integrados
DefaultAssay(so) <- "integrated"
so <- ScaleData(so, display.progress = FALSE)

#Reducción de dimension
so <- RunPCA(so, npcs = 100, verbose = FALSE)

#Cambiar numbero de PCs usedos
so <- RunTSNE(so, reduction = "pca", dims = seq_len(20),
              seed.use = 1, do.fast = TRUE, verbose = FALSE)
so <- RunUMAP(so, reduction = "pca", dims = seq_len(20),
              seed.use = 1, verbose = FALSE)

#CLUSTERING
so <- FindNeighbors(so, reduction = "pca", dims = seq_len(20), verbose = FALSE)
for (res in c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1)){
  so <- FindClusters(so, resolution = res, random.seed = 1, verbose = FALSE)
}

```

<div class = "red color">

Para los ejercicios de la presentación descarga los siguientes archivos

[https://drive.google.com/file/d/1HbJ0syxTcxkI1aG6dIQ6hkJ7SrLHCLfl/view?usp=drive_link](https://drive.google.com/file/d/1HbJ0syxTcxkI1aG6dIQ6hkJ7SrLHCLfl/view?usp=drive_link)

[https://drive.google.com/file/d/1gJ2raTr8BhLsjeW3n1sdKko0qLYE2g5n/view?usp=drive_link](https://drive.google.com/file/d/1gJ2raTr8BhLsjeW3n1sdKko0qLYE2g5n/view?usp=drive_link)

</div>