01_sesion1.Rmd

# Estructura e importe de datos

Joselyn Cristina Chávez Fuentes

07 de agosto de 2023

## Diapositivas

[
```{r,echo=FALSE}
knitr::include_url("https://comunidadbioinfo.github.io/cdsb2023/dia1_sesion1_slides.html", height = "380px")
```
](https://comunidadbioinfo.github.io/cdsb2023/dia1_sesion1_slides.html)

## Bulk RNAseq vs single-cell RNAseq

```{r, echo=FALSE, out.width='55%', fig.align='center'}
knitr::include_graphics("img/experiment.jpg")
```
Tomado de [Yu X, et al.](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7771369/#:~:text=The%20difference%20between%20single%20cell,within%20the%20same%20tissue%20type.)

- **Bulk RNAseq**

```{r, echo=FALSE, out.width='90%', fig.align='center'}
knitr::include_graphics("img/bulk_matrix.png")
```
Tomado de [Hyeongseon Jeon, et al.](https://www.mdpi.com/2218-273X/13/2/221)

- **scRNAseq**

```{r, echo=FALSE, out.width='100%', fig.align='center'}
knitr::include_graphics("img/single_matrix.png")
```


## Consideraciones experimentales

### Tecnologías de scRNAseq

- **Droplet-based:**

Son las tecnologías más usadas debido a su buen rendimiento con un relativo bajo costo.

- 10X Genomics 
- inDrop seq
- Drop-seq 


- **Plate-based**

Tienen la capacidad de capturar información adicional, como la morfología. Tienen más opciones para personalizar el diseño de experimento.

- Con UMIs:
    - CEL-seq 
    - MARS-seq 

- Con reads: 
    - Smart-seq2


- **Otros**
    - sciRNA-seq

### ¿UMIs o Reads?

Los métodos basados en reads proveen una cobertura de transcriptoma completo, lo que puede ser útil para ciertas aplicaciones, como el análisis de splicing o mutaciones en exones. 

Los métodos basados en UMIs suelen ser más populares, ya que eliminan el ruido causado por la amplificación durante el PCR.

### ¿Cuántas células y profundidad de secuenciación necesito?

Depende...

Si tu objetivo es el estudio de subgrupos raros o poco abundantes de células, necesitarás un mayor número de células por muestra. 

Si tu objetivo es estudiar diferencias sutiles en la expresión de genes, necesitarás una mayor profundidad de secuenciación.

Hasta el momento, las tecnologías basadas en droplets capturan entre 10,000 y 100,000 células, con un aproximado de 1,000 a 10,000 UMIs por célula. 

Existen algunas variaciones entre el rendimiento y la tasa de doublets que pueden afectar la eficiencia real de secuenciación. 

Encontrar un balance entre el número de células por muestra, la profundidad de secuenciación, el número de condiciones y réplicas a secuenciar depende mucho de la aplicación y el presupuesto.


## Generación de la matriz de cuentas

El método para generar la matriz depende de la tecnología utilizada.

- 10X Genomics:

  - [Cellranger](https://pubmed.ncbi.nlm.nih.gov/28091601/) es el programa más popular. Utiliza STAR para alinear los reads con el genoma de referencia y después cuenta el número de UMIs únicos mapeados con cada gene.
  
  - Como alternativa, métodos de pseudo-alineamiento como [alevin](https://genomebiology.biomedcentral.com/articles/10.1186/s13059-019-1670-y) pueden ser usados. No requiere un alineamiento con genoma de referencia, lo que reduce el tiempo de cómputo y los requerimientos de memoria.
  
- El paquete scPipe provee un análisis general. Utiliza Rsubread para alinear los reads y después cuenta reads o UMIs por gene.

- CEL-seq o CEL-seq2 

  - El paquete scruff provee un pipeline para la cuantificación.
  
De manera general, los protocolos basados en reads pueden reusar los métodos desarrollados para bulk RNA-seq.

Si el grupo de datos involucra transcritos spike-in, las secuencias spike-in deben proveerse junto al genoma de referencia durante el alineamiento y cuantificación.


## Nombres de los genes

En todos los casos, los identificadores de genes deben definirse considerando los nombres de Ensembl o Entrez. Esto mantiene un mapeo sin lugar a errores para la identificación de genes en la matriz.

Estos identificadores pueden ser reemplazados por el nombre común del gene durante el análisis, pero considera que estos nombres pueden cambiar con el tiempo. 

Consideraciones adicionales:

- Algunas herramientas de conteo, como HTSeq, incluyen un reporte dentro de la matriz de cuentas con el número lecturas sin alinear. Estos valores pueden ser útiles para el control de calidad, pero deben ser removidas de la matriz de cuentas y guardadas en otro lugar antes de continuar con el análisis para que no sean confundidas con los valores de expresión de genes.

- Generalmente, las secuencias spike-ins son desarrolladas por el External RNA Controls Consortium (ERCC) y suelen tener nombres como ERCC-00002. En las muestra de humano, debemos evitar confundir estos identificadores con los genes de la familia ERCC, que suelen tener nombres similares ERCC1. Podemos evitar estos problemas si usamos los identificadores de Ensembl.

## Importando los datos

### Datos tabulares

Descarguemos los datos de cáncer de páncreas de Muraro et al (2016) [GSE85241](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE85241).

Solamente necesitamos el archivo "GSE85241_cellsystems_dataset_4donors_updated.csv"

- Una forma sencilla

```{r, eval=FALSE}
mat <- read.delim(
  "GSE85241_cellsystems_dataset_4donors_updated.csv")

mat <- as.matrix(mat)

dim(mat)

mat[1:3,1:3]
```

- Una opción más eficiente. Solamente se leen los datos diferentes a cero.

```{r, eval=FALSE}
sparse.mat <- scuttle::readSparseCounts(
  "GSE85241_cellsystems_dataset_4donors_updated.csv")

dim(sparse.mat)

sparse.mat[1:3,1:3]
```


## Actividad

<div class = "orange color">

- Comparemos la clase y la cantidad de memoria que utiliza cada opción.

- Utiliza las funciones class() y lobstr::obj_size() para comparar mat y sparse.mat

</div>

### Datos de cellRanger

Cuando usamos datos de 10X Genomics, Cellranger genera un directorio que contiene 3 archivos: cuentas, anotaciones de features y barcodes. Podemos usar la ruta a este directorio y la función read10xCounts() del paquete DropletUtils.

Utilicemos como ejemplo un dataset de [células sanguíneas periféricas](https://support.10xgenomics.com/single-cell-gene-expression/datasets/2.1.0/pbmc4k?).

```{r, eval=FALSE}
sce <- DropletUtils::read10xCounts(
  "raw_gene_bc_matrices/GRCh38/")
```


## Actividad

<div class = "orange color">

- Evalúa la clase de sce. ¿Habías escuchado antes sobre este tipo de objeto?

- Imprime el objeto sce en la consola. ¿Cómo se ve el resultado?

</div>


### Datos con formato HDF5

El formato Hierarchical Data Format version 5 (HDF5) permite guardar tanto los valores de expresión asociados a los genes, como la anotación de los tipos celulares dentro de un mismo archivo.

Podemos leer este archivo con un objeto SingleCellExperiment utilizando el paquete zellkonverter.

```{r, eval=FALSE}
library(zellkonverter)
demo <- system.file("extdata", 
                    "krumsiek11.h5ad", 
                    package = "zellkonverter")
sce <- readH5AD(demo)
```


### Datos con formato loom

Los archivos con formato Loom son una variante de HDF5, pueden ser leídos como un objeto SingleCellLoomExperiment usando el paquete LoomExperiment. Este objeto es una extensión del SingleCellExperiment.

```{r, eval=FALSE}
library(LoomExperiment)
demo <- system.file("extdata", 
                    "L1_DRG_20_example.loom", 
                    package = "LoomExperiment")
scle <- import(demo, type="SingleCellLoomExperiment")
```


## El objeto SingleCellExperiment

```{r, echo=FALSE, out.width='80%', fig.align='center'}
knitr::include_graphics("img/SingleCellExperiment.png")
```

## Construyendo un objeto SingleCellExperiment

Para crear un SingleCellExperiment rudimentario, solamente necesitamos el slot assays.

Este slot necesita la matriz de cuentas, donde los renglones corresponden a los genes (features) y las columnas corresponden a las células.

```{r, eval=FALSE}
mat <- read.delim(
  "GSE85241_cellsystems_dataset_4donors_updated.csv")
mat <- as.matrix(mat)
sce <- SingleCellExperiment::SingleCellExperiment(
  assays = list(counts = mat))
```

## Accediendo a los elementos del objeto 

Para poder ver el contenido del slot assays tenemos dos opciones:

- Forma general, usando el nombre del assay que queremos ver:

```{r, eval=FALSE}
m <- assay(sce, "counts")
```

- Una opción más corta, pero sólo funciona cuando el assay se llama "counts"

```{r, eval=FALSE}
m <- counts(sce)
```

## Agregando más assays

Se pueden guardar varias versiones de los datos de cuentas.

```{r, eval=FALSE}
sce <- scuttle::logNormCounts(sce)

sce
```

Observa el slot assays, notas algo diferente?

También se pueden conocer los assays disponibles:

```{r, eval=FALSE}
assayNames(sce)
```

## Un vistazo al flujo de trabajo

```{r, echo=FALSE, out.width='40%', fig.align='center'}
knitr::include_graphics("img/workflow.png")
```