Taller de análisis de datos de expresión diferencial

# Taller de análisis de datos de expresión diferencial
### Rodolfo LCD
### 12 January 2022

---

.remark-slide-content {
  font-size: 28px;
  padding: 20px 80px 20px 80px;
}
.remark-code, .remark-inline-code {
  background: #f0f0f0;
}
.remark-code {
  font-size: 24px;
}
.huge .remark-code { /*Change made here*/
  font-size: 200% !important;
}
.tiny .remark-code { /*Change made here*/
  font-size: 80% !important;
}

.center2 {
  margin: 0;
  position: absolute;
  top: 50%;
  left: 50%;
  -ms-transform: translate(-50%, -50%);
  transform: translate(-50%, -50%);
}

/* From https://github.com/garthtarr/sydney_xaringan */
blockquote, .blockquote {
  display: block;
  margin-top: 0.1em;
  margin-bottom: 0.2em;
  margin-left: 5px;
  margin-right: 5px;
  border-left: solid 10px #0148A4;
  border-top: solid 2px #0148A4;
  border-bottom: solid 2px #0148A4;
  border-right: solid 2px #0148A4;
  box-shadow: 0 0 6px rgba(0,0,0,0.5);
  /* background-color: #e64626; */
  color: #e64626;
  padding: 0.5em;
  -moz-border-radius: 5px;
  -webkit-border-radius: 5px;
}
/* From https://github.com/garthtarr/sydney_xaringan */
.content-box-blue,
.content-box-gray,
.content-box-grey,
.content-box-army,
.content-box-green,
.content-box-purple,
.content-box-red,
.content-box-yellow {
  box-sizing: border-box;
  border-radius: 15px;
  margin: 0 0 15px;
  overflow: hidden;
  padding: 0px 20px 0px 20px;
  width: 100%;
}
</style>

## Material

Este material está basado en el pipeline de [Michael Love](https://www.bioconductor.org/packages/devel/workflows/vignettes/rnaseqGene/inst/doc/rnaseqGene.html) de Bioconductor.

---
class: inverse, center, middle

# 1: Introducción al RNA-seq

---
## Transcriptómica

Estudio de la expresión genética

- ¿Cuánto RNA hay?

---
## RNA-seq

---

## Aplicaciones

---

## Bases de datos

---

# 2: Organización de un proyecto bioinformático

---

## Organización del proyecto

Almacenar el proyecto en .blue[GitHub] <svg viewBox="0 0 496 512" style="height:1em;position:relative;display:inline-block;top:.1em;fill:blue;" xmlns="http://www.w3.org/2000/svg">  <path d="M165.9 397.4c0 2-2.3 3.6-5.2 3.6-3.3.3-5.6-1.3-5.6-3.6 0-2 2.3-3.6 5.2-3.6 3-.3 5.6 1.3 5.6 3.6zm-31.1-4.5c-.7 2 1.3 4.3 4.3 4.9 2.6 1 5.6 0 6.2-2s-1.3-4.3-4.3-5.2c-2.6-.7-5.5.3-6.2 2.3zm44.2-1.7c-2.9.7-4.9 2.6-4.6 4.9.3 2 2.9 3.3 5.9 2.6 2.9-.7 4.9-2.6 4.6-4.6-.3-1.9-3-3.2-5.9-2.9zM244.8 8C106.1 8 0 113.3 0 252c0 110.9 69.8 205.8 169.5 239.2 12.8 2.3 17.3-5.6 17.3-12.1 0-6.2-.3-40.4-.3-61.4 0 0-70 15-84.7-29.8 0 0-11.4-29.1-27.8-36.6 0 0-22.9-15.7 1.6-15.4 0 0 24.9 2 38.6 25.8 21.9 38.6 58.6 27.5 72.9 20.9 2.3-16 8.8-27.1 16-33.7-55.9-6.2-112.3-14.3-112.3-110.5 0-27.5 7.6-41.3 23.6-58.9-2.6-6.5-11.1-33.3 2.6-67.9 20.9-6.5 69 27 69 27 20-5.6 41.5-8.5 62.8-8.5s42.8 2.9 62.8 8.5c0 0 48.1-33.6 69-27 13.7 34.7 5.2 61.4 2.6 67.9 16 17.7 25.8 31.5 25.8 58.9 0 96.5-58.9 104.2-114.8 110.5 9.2 7.9 17 22.9 17 46.4 0 33.7-.3 75.4-.3 83.6 0 6.5 4.6 14.4 17.3 12.1C428.2 457.8 496 362.9 496 252 496 113.3 383.5 8 244.8 8zM97.2 352.9c-1.3 1-1 3.3.7 5.2 1.6 1.6 3.9 2.3 5.2 1 1.3-1 1-3.3-.7-5.2-1.6-1.6-3.9-2.3-5.2-1zm-10.8-8.1c-.7 1.3.3 2.9 2.3 3.9 1.6 1 3.6.7 4.3-.7.7-1.3-.3-2.9-2.3-3.9-2-.6-3.6-.3-4.3.7zm32.4 35.6c-1.6 1.3-1 4.3 1.3 6.2 2.3 2.3 5.2 2.6 6.5 1 1.3-1.3.7-4.3-1.3-6.2-2.2-2.3-5.2-2.6-6.5-1zm-11.4-14.7c-1.6 1-1.6 3.6 0 5.9 1.6 2.3 4.3 3.3 5.6 2.3 1.6-1.3 1.6-3.9 0-6.2-1.4-2.3-4-3.3-5.6-2z"></path></svg> para el control de versiones

---

## Organización del proyecto

- La función del **README** es documentar el proyecto y que otros usuarios entiendan su contenido
--

- Los datos deben de almacenarse en un directorio especial (considerar espacio en disco):
--

- Considerar almacenar los datos en un repositorio como respaldo [OSF](https://osf.io/) <svg viewBox="0 0 512 512" style="height:1em;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <g label="icon" id="layer6" groupmode="layer">    <path id="path2" d="m 255.9997,7.9999987 c -34.36057,0 -62.21509,27.8545563 -62.21509,62.2151643 0,20.303056 9.87066,38.160947 24.91769,49.517247 0.18814,-20.457899 16.79601,-36.993393 37.29685,-36.993393 20.50082,0 37.11091,16.535494 37.29909,36.993393 15.04533,-11.3563 24.9177,-29.212506 24.9177,-49.517247 C 318.21272,35.854555 290.35915,7.9999987 255.99915,7.9999987 Z M 293.29654,392.2676 c -0.18814,20.4601 -16.79601,36.99338 -37.29684,36.99338 -20.50082,0 -37.10922,-16.53551 -37.29684,-36.99338 -15.04759,11.35627 -24.91769,29.21246 -24.91769,49.51722 0,34.36059 27.85453,62.21518 62.2151,62.21518 34.36056,0 62.21508,-27.85459 62.21508,-62.21518 0,-20.30306 -9.87066,-38.16095 -24.91767,-49.51722 z M 441.78489,193.78484 c -20.30301,0 -38.16309,9.87068 -49.51717,24.91769 20.45786,0.18819 36.99333,16.79605 36.99333,37.29689 0,20.50085 -16.53547,37.11096 -36.9911,37.29916 11.35634,15.04533 29.21249,24.91769 49.51721,24.91769 C 476.14549,318.21327 504,290.35948 504,255.99942 504,221.6394 476.14549,193.78425 441.78489,193.78425 Z M 82.738898,255.99997 c 0,-20.50139 16.535509,-37.11096 36.993392,-37.29689 -11.35632,-15.04756 -29.214164,-24.91769 -49.517197,-24.91769 -34.36057,0 -62.2150945,27.85455 -62.2150945,62.21517 0,34.3606 27.8545245,62.21516 62.2150945,62.21516 20.303033,0 38.160877,-9.87068 49.517197,-24.91773 -20.457883,-0.18818 -36.993391,-16.796 -36.993391,-37.29688 z M 431.3627,80.636814 c -24.29549,-24.295544 -63.68834,-24.295544 -87.9844,0 -14.35704,14.357057 -20.00298,33.963346 -17.39331,52.633806 -0.0824,0.0809 -0.18198,0.13437 -0.26434,0.21491 -14.578,14.57799 -14.578,38.21689 0,52.79488 14.57797,14.57799 38.21681,14.57799 52.79484,0 0.0824,-0.0824 0.13455,-0.18198 0.21732,-0.26434 18.66819,2.60796 38.27445,-3.03799 52.63151,-17.39336 24.29378,-24.29778 24.29378,-63.68837 -0.003,-87.986153 z M 186.2806,378.51178 c 14.57798,-14.57799 14.57798,-38.21461 0,-52.79319 -14.57798,-14.57853 -38.21683,-14.57798 -52.79481,0 -0.0825,0.0824 -0.13448,0.18199 -0.21476,0.26215 -18.67046,-2.60795 -38.276723,3.03572 -52.63376,17.39505 -24.297753,24.29552 -24.297753,63.6884 0,87.98449 24.29551,24.29552 63.68833,24.29552 87.98439,0 14.35702,-14.35703 20.00297,-33.96333 17.39333,-52.63386 0.0848,-0.0786 0.18364,-0.13228 0.26672,-0.21505 z m 0,-245.02583 c -0.0826,-0.0824 -0.18198,-0.13436 -0.26445,-0.21494 2.60795,-18.66823 -3.038,-38.27452 -17.39332,-52.633811 -24.29777,-24.295544 -63.68832,-24.295544 -87.984405,0 -24.297747,24.297781 -24.297747,63.688381 0,87.986151 14.357042,14.35706 33.963315,20.00301 52.631515,17.39336 0.0808,0.0824 0.13447,0.18199 0.21475,0.26434 14.57799,14.57799 38.21684,14.57799 52.79482,0 14.57797,-14.57802 14.57797,-38.21689 0,-52.79488 z m 245.0821,209.89048 c -14.35703,-14.35703 -33.96329,-20.00301 -52.63378,-17.39505 -0.0809,-0.0824 -0.13228,-0.18199 -0.21506,-0.26215 -14.57797,-14.57799 -38.21685,-14.57799 -52.79482,0 -14.57797,14.57799 -14.57797,38.21461 0,52.79316 0.0827,0.0828 0.18198,0.13455 0.26434,0.21505 -2.60796,18.67053 3.03802,38.27683 17.39334,52.63386 24.29552,24.29552 63.68834,24.29552 87.98439,0 24.29775,-24.29552 24.29775,-63.68841 0.003,-87.98451 z" style="stroke-width:0.07717"></path>  </g></svg>
--

- **El código es una de las partes más importantes del proyecto**

---
.scroll-output[
## Organización del proyecto
## Código

- Enumerar los scripts de acuerdo al orden en que serán ejecutados

```
+-- 1.QualityControl.sh
+-- 2.Trimming.sh
+-- 3.ReadAlignment.sh
+-- 4.DifferentialExpression.R
```

- Utilizar rutas relativas hacia los archivos *input*

```
fastqc ../Data/*.fastqc -o ../Results/
```

- Comentar el script <svg viewBox="0 0 320 512" style="height:1em;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M285.476 272.971L91.132 467.314c-9.373 9.373-24.569 9.373-33.941 0l-22.667-22.667c-9.357-9.357-9.375-24.522-.04-33.901L188.505 256 34.484 101.255c-9.335-9.379-9.317-24.544.04-33.901l22.667-22.667c9.373-9.373 24.569-9.373 33.941 0L285.475 239.03c9.373 9.372 9.373 24.568.001 33.941z"></path></svg> facilita que sea entendido por humanos

```
#Change the location of bam files
mv ../Data/*.bam ../Results/
```

Consulta el paquete [`here`](https://github.com/jennybc/here_here) desarrollado por Jenny B.

- ¡Mejorar el código! <svg viewBox="0 0 320 512" style="height:1em;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M285.476 272.971L91.132 467.314c-9.373 9.373-24.569 9.373-33.941 0l-22.667-22.667c-9.357-9.357-9.375-24.522-.04-33.901L188.505 256 34.484 101.255c-9.335-9.379-9.317-24.544.04-33.901l22.667-22.667c9.373-9.373 24.569-9.373 33.941 0L285.475 239.03c9.373 9.372 9.373 24.568.001 33.941z"></path></svg> *Less is more*
]

---

---
## Organización del proyecto
<img src="./Images/Proyecto.png" width="450px" style="display: block; margin: auto;" />

---
## Experimento

Número de acceso GEO: GSE136011

---

# 3: Flujo de trabajo para el análisis de RNA-seq
---
## Flujo general de trabajo

---

## Control de calidad

---

## Alineamiento y cuantificación

---

## Expresión diferencial y análisis funcionales

---

## Flujo general de trabajo

---

# 4: Tipos de archivos

---

## FastQ

- Derivado del formato FASTA

- Id secuencia + **Calidad**

- Generalmente no está comprimido (.fastq) <svg viewBox="0 0 320 512" style="height:1em;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M285.476 272.971L91.132 467.314c-9.373 9.373-24.569 9.373-33.941 0l-22.667-22.667c-9.357-9.357-9.375-24.522-.04-33.901L188.505 256 34.484 101.255c-9.335-9.379-9.317-24.544.04-33.901l22.667-22.667c9.373-9.373 24.569-9.373 33.941 0L285.475 239.03c9.373 9.372 9.373 24.568.001 33.941z"></path></svg> comprimir (.fastq.gz)

- Comprende 4 líneas por secuencia:
  * @ ID del read + información de la corrida
  * Secuencia 
  * Símbolo "+"
  * Calidad (Escala **Phred** y código **ASCII**)
  
---

## FastQ  
## Identificador

- El identificador de la secuencia es la primera línea del archivo .fastq

- Plataformas de Illumina superiores a 1.8 contienen la siguiente información:

**@machneID:run number:flowCell ID:lane:tile:read:control number:index**

---

## FastQ
## Calidad

La calidad de cada base secuenciada, en la escala **Phred**, se calcula de acuerdo la siguiente ecuación:

Q = Calidad

p = Probabilidad de que el *base call* sea incorrecto

Valores de **p** cercanos a cero <svg viewBox="0 0 320 512" style="height:1em;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M285.476 272.971L91.132 467.314c-9.373 9.373-24.569 9.373-33.941 0l-22.667-22.667c-9.357-9.357-9.375-24.522-.04-33.901L188.505 256 34.484 101.255c-9.335-9.379-9.317-24.544.04-33.901l22.667-22.667c9.373-9.373 24.569-9.373 33.941 0L285.475 239.03c9.373 9.372 9.373 24.568.001 33.941z"></path></svg> Alta calidad

---
## FastQ
## Calidad

---
## FastQ
## Calidad

- La calidad se codifica en formato **[ASCII](https://www.ascii-code.com/)**

- Phred + 33 <svg viewBox="0 0 320 512" style="height:1em;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M285.476 272.971L91.132 467.314c-9.373 9.373-24.569 9.373-33.941 0l-22.667-22.667c-9.357-9.357-9.375-24.522-.04-33.901L188.505 256 34.484 101.255c-9.335-9.379-9.317-24.544.04-33.901l22.667-22.667c9.373-9.373 24.569-9.373 33.941 0L285.475 239.03c9.373 9.372 9.373 24.568.001 33.941z"></path></svg>  Plataformas de Illumina a partir de 1.8 (símbolos especiales + alfanuméricos)

- Phred + 64 <svg viewBox="0 0 320 512" style="height:1em;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M285.476 272.971L91.132 467.314c-9.373 9.373-24.569 9.373-33.941 0l-22.667-22.667c-9.357-9.357-9.375-24.522-.04-33.901L188.505 256 34.484 101.255c-9.335-9.379-9.317-24.544.04-33.901l22.667-22.667c9.373-9.373 24.569-9.373 33.941 0L285.475 239.03c9.373 9.372 9.373 24.568.001 33.941z"></path></svg> Plataformas de Illumina anteriores a 1.8

---
## FastQ

---

## FastQ

---
class: inverse, center, middle

# 5: Control de calidad de los datos
---
## ¿Porqué es importante revisar los datos crudos?
--

- Los datos influencian los análisis posteriores

.pull-left[
<img src="https://miro.medium.com/max/470/1*DVLpt3zCNaqeZOGNpgME7Q.png" width="150px" style="display: block; margin: auto;" /><img src="https://i.pinimg.com/originals/30/66/ac/3066ac69ae68ac200ace0ca8fe3882c3.jpg" width="150px" style="display: block; margin: auto;" />
]

]

---

## FastQC

Usaremos el programa de [FastQC](https://www.bioinformatics.babraham.ac.uk/projects/download.html#fastqc) para analizar la calidad de los datos.

Podemos usar la interfaz gráfica o la línea de comandos, en este caso usaremos la línea de comandos.

---
class: center, middle
background-image: url("./Images/hands_on.png")
background-size: 90% 80%
---

### Pasos para correr FastQC:

1. Corre `FastQC`

```bash
##La opción -o es para elegir el directorio de salida
fastqc ../data/fastq_files/*.fastq.gz -o ../results/fastqc_output/
```

2. Revisa los archivos generados

```bash
#Busca los archivos generados por fastqc en el directorio de resultados
ls ../results/fastqc_output/
```

---

## Buena calidad

El eje x = longitud de la read, y =  calidad (phred score)

Los datos a continuación muestran lecturas de buena calidad (Verde = buena, amarilla = aceptable, roja = baja).

---

## Mala calidad

¿Hasta donde los datos son de buena calidad? ¿Estos datos los usarías?

---

## Ejercicio

¿Cómo se ven nuestros datos?

---

## MultiQC

[MultiQC](https://multiqc.info/) te permite juntar los archivos de salida para visualizarlos de forma óptima.

1. Corre MultiQC

```bash
##Nuevamente, la opción -o es para especificar la ruta donde se almacenarán los archivos de salida
multiqc ../results/fastqc_output/*.zip -o ../results/multiqc_output/
```

2. Revisar los datos de salida creados por MultiQC

---

# 6: Limpieza de las cuentas crudas

---
## Limpieza

En ocasiones las cuentas crudas presentan problemas de calidad en distintos niveles. Para RNAseq (expresión diferencial) prestar atención en:

---

## Herramientas

[Trimmomatic](http://www.usadellab.org/cms/?page=trimmomatic)

[Trimgalore](https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/)

[Cutadapt](https://cutadapt.readthedocs.io/en/stable/)

--- 
class: center, middle
background-image: url("./Images/hands_on.png")
background-size: 90% 80%
---

## Limpieza

1. Ejecuta `cutadapt` para una sola muestra:

```bash
cutadapt -a AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT ##Para el read forward
-A AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT ##Para el read reverse
-o ../Output/Cutadapt/samp1_clean_R1.fastq.gz ##Archivo de salida para las lecturas forward
-p ../Output/Cutadapt/samp1_clean_R2.fastq.gz ##Archivo de salida para las lecturas reverse
../data/samp1_R1.fastq.gz ##Ruta hacia las lecturas crudas forward
../data/samp1_R2.fastq.gz ##Ruta hacia las lecturas crudas reverse

```

]

---

## Limpieza

De manera posterior a la limpieza de las lecturas, es muy importante analizar la calidad de las lecturas procesadas:

---
class: inverse, center, middle

# 7: Alineamiento
---

---

## Alineamiento
<img src="./Images/ReadAlignment.png" width="600px" style="display: block; margin: auto;" />
Buscar la región del genoma/transcriptoma a partir de la cual se originaron las lecturas

---

## Alineamiento
## Spliced reads

---

## Splice aware

- TopHat 2
<img src="./Images/SpliceAware.png" width="600px" style="display: block; margin: auto;" />

--
- HISAT
<img src="./Images/HISAT.png" width="600px" style="display: block; margin: auto;" />

--
- STAR
<img src="./Images/STAR.png" width="400px" style="display: block; margin: auto;" />

---

## Alineamiento
## Pseudo-alineadores

Pseudo-alineadores (quasi-alineadores):

- Kallisto

- Sailfish

- **Salmon**

![](https://combine-lab.github.io/salmon/images/SalmonLogo.png)

---

## Alineamiento
## Pseudo-alineadores

---

## Alineamiento
## Pseudo-alineadores

--
<img src="./Images/PseudoAlignment2.png" width="500px" style="display: block; margin: auto;" />

---

## Alineamiento
## Pseudo-alineadores

---
class: center, middle
background-image: url("./Images/hands_on.png")
background-size: 90% 80%
---
## Generar índice del transcriptoma
.scroll-output[
Para ello vamos a requerir:

- Archivo **Fasta** del transcriptoma de referencia del humano.

Descargado del sitio de [GeneCode](https://www.gencodegenes.org/human/)

- Código para genear el índice

```bash
salmon index -t ../reference/Homo_sapiens.GRCh38.cdna.all.fa ##Ruta hacia el archivo fasta del transcriptoma de referencia
-i ../reference/hsa_v38_gencode ##Ruta para almacenar los archivos generados del índice
```

`-t`: Ubicación al archivo Fasta del transcriptoma

`-i`: Ubicación para salvar el índice

`--genecode`: El Fasta del transcriptoma de referencia está en formato de GENECODE

]
---
## Cuantificar la abundancia de los transcritos
.scroll-output[
Requerimientos:

- Archivos **Fastq** de las lecturas -> Paired end R1 y R2

Ubicados en el folder de `Output/Cutadapt`

- Índice generado en el paso anterior

- Código para producir cuantificar los transcritos

```bash
salmon quant -i ../reference/hsa_v38_gencode
-l A ##Opción para permitir a salmon deducir de forma automática las características de la libreria
-1 ../Output/Cutadapt/samp1_clean_R1.fastq.gz
-2 ../Output/Cutadapt/samp1_clean_R2.fastq.gz
-p 8 ##Número de cores para paralelizar el proceso
--validateMappings 
-o ../Output/salmon_quants/samp1_quant
```

`-i`: Ubicación del índice

`-l`: Tipo de librería

`-1 y -2`: Ubicación a las lecturas R1 y R2

`-o`: Ubicación para almacenar los resultados

]
---

# 8: Tximeta y Tximport
---
## Recordatorio

¿Dónde estamos?

---

## Recordatorio

¿Dónde estamos?

<img src="Images/worflowTximeta.png" width="450px" style="display: block; margin: auto;" />
 
---

## Tximeta

Paquetería que conserva

---

## Importar datos

Utilizamos `tximeta` para importar los datos.

Para ello requerímos de una tabla de metadatos con información de las muestras

```bash

#La tabla se ubica en el directorio de salmon_quants con el nombre de metadata.txt
ls ../data/salmon_quants
```

```
## metadata.txt
## samp10_quant
## samp11_quant
## samp12_quant
## samp1_quant
## samp2_quant
## samp3_quant
## samp4_quant
## samp5_quant
## samp6_quant
## samp7_quant
## samp8_quant
## samp9_quant
```

---
class: center, middle
background-image: url("./Images/hands_on.png")
background-size: 90% 80%
---

## Importar datos

```r
##Lee la tabla de metadatos en R
coldata <- read.table(here("results/salmon_quants/metadata.txt"), he = T, stringsAsFactors = T, sep = "\t")
##Nombra las filas de la tabla empleando la columna de los id únicos
rownames(coldata) <- coldata$unique_id
##Crea una nueva columna que contenga la ruta hacia los archivos quant.sf
coldata <- dplyr::mutate(coldata, files = file.path(here("results/salmon_quants"), paste0(coldata$key, "_quant"), "quant.sf"))
##Genera una columna para verificar que las rutas hacia los archivos quant.sf fueron creadas correctamente
coldata <- mutate(coldata, file_exist = file.exists(coldata$files))
##Asigna una columna llamada names que contenga el nombre de los directorios de cada muestra
coldata <- mutate(coldata, names = paste0(coldata$key, "_quant"))
```

Importar los datos usando `tximeta`

]

---
class: inverse, center, middle

# 9: Objeto Summarized Experiment
---

## Summarized Experiment

---
## Summarized Experiment

---

## Summarized Experiment
## ColData

.scroll-output[
- El cajón o *slot* correspondiente a `ColData` contiene la tabla de metadatos `coldata` empleada para importar las cuentas con `tximeta`

- Para acceder a **.red[ColData]** usar el siguiente comando:

```r
colData(se)
```

```
## DataFrame with 12 rows and 8 columns
##                                                                                   file_name
##                                                                                    <factor>
## siRNA_NRF2_1    SRR10004274_GSM4039376_NRF2_siRNA_rep_1_Homo_sapiens_RNA-Seq_1.fastq.gz    
## siRNA_control_1 SRR10004275_GSM4039377_NRF2_siRNA_cont_rep1_Homo_sapiens_RNA-Seq_1.fastq.gz
## siRNA_NRF2_2    SRR10004276_GSM4039378_NRF2_siRNA_rep_2_Homo_sapiens_RNA-Seq_1.fastq.gz    
## siRNA_control_2 SRR10004277_GSM4039379_NRF2_siRNA_cont_rep2_Homo_sapiens_RNA-Seq_1.fastq.gz
## siRNA_NRF2_3    SRR10004278_GSM4039380_NRF2_siRNA_rep_3_Homo_sapiens_RNA-Seq_1.fastq.gz    
## ...                                                                                     ...
## mut_NRF2_1      SRR10004281_GSM4039383_NRF2_mutant_rep_1_Homo_sapiens_RNA-Seq_1.fastq.gz   
## CT_2            SRR10004282_GSM4039384_NRF2_wild_type_rep_2_Homo_sapiens_RNA-Seq_1.fastq.gz
## mut_NRF2_2      SRR10004283_GSM4039385_NRF2_mutant_rep_2_Homo_sapiens_RNA-Seq_1.fastq.gz   
## CT_3            SRR10004284_GSM4039386_NRF2_wild_type_rep_3_Homo_sapiens_RNA-Seq_1.fastq.gz
## mut_NRF2_3      SRR10004285_GSM4039387_NRF2_mutant_rep_3_Homo_sapiens_RNA-Seq_1.fastq.gz   
##                     Treatment Replicate       unique_id      key     Geno
##                      <factor> <integer>        <factor> <factor> <factor>
## siRNA_NRF2_1    siRNA_NRF2            1 siRNA_NRF2_1       samp1       WT
## siRNA_control_1 siRNA_control         1 siRNA_control_1    samp2       WT
## siRNA_NRF2_2    siRNA_NRF2            2 siRNA_NRF2_2       samp3       WT
## siRNA_control_2 siRNA_control         2 siRNA_control_2    samp4       WT
## siRNA_NRF2_3    siRNA_NRF2            3 siRNA_NRF2_3       samp5       WT
## ...                       ...       ...             ...      ...      ...
## mut_NRF2_1           mut_NRF2         1      mut_NRF2_1   samp8       Mut
## CT_2                 CT               2      CT_2         samp9       WT 
## mut_NRF2_2           mut_NRF2         2      mut_NRF2_2   samp10      Mut
## CT_3                 CT               3      CT_3         samp11      WT 
## mut_NRF2_3           mut_NRF2         3      mut_NRF2_3   samp12      Mut
##                 file_exist        names
##                  <logical>  <character>
## siRNA_NRF2_1          TRUE  samp1_quant
## siRNA_control_1       TRUE  samp2_quant
## siRNA_NRF2_2          TRUE  samp3_quant
## siRNA_control_2       TRUE  samp4_quant
## siRNA_NRF2_3          TRUE  samp5_quant
## ...                    ...          ...
## mut_NRF2_1            TRUE  samp8_quant
## CT_2                  TRUE  samp9_quant
## mut_NRF2_2            TRUE samp10_quant
## CT_3                  TRUE samp11_quant
## mut_NRF2_3            TRUE samp12_quant
```
]

---

## Summarized Experiment
## ColData

- *Rownames* de **ColData** corresponden a los *Colnames* en el slot **Assay**- **.red[Importante para análisis con DESeq2]**

]
---

## Summarized experiment

---
## Summarized experiment
## rowRanges

.scroll-output[
- El cajón de `rowRanges` es un objeto de tipo *GRanges* que hace referencia a las coordenadas de cada transcrito y su anotación correspondiente

- Para acceder al **.red[rowRanges]** usar:

```r
rowRanges(se)
```

```
## GRanges object with 184844 ranges and 9 metadata columns:
##                             seqnames              ranges strand |
##                                <Rle>           <IRanges>  <Rle> |
##   ENST00000631435  CHR_HSCHR7_2_CTG6 142847306-142847317      + |
##   ENST00000415118                 14   22438547-22438554      + |
##   ENST00000434970                 14   22439007-22439015      + |
##   ENST00000448914                 14   22449113-22449125      + |
##   ENST00000632524 CHR_HSCHR14_3_CTG1 105866322-105866332      - |
##               ...                ...                 ...    ... .
##   ENST00000444456                  1 247974741-247975653      + |
##   ENST00000416377                  1 248003129-248004068      + |
##   ENST00000318104                  1 248121932-248122882      + |
##   ENST00000606902                  1 248498308-248498649      + |
##   ENST00000427827                  1 248549444-248549762      - |
##                             tx_id             tx_biotype tx_cds_seq_start
##                       <character>            <character>        <integer>
##   ENST00000631435 ENST00000631435              TR_D_gene        142847306
##   ENST00000415118 ENST00000415118              TR_D_gene         22438547
##   ENST00000434970 ENST00000434970              TR_D_gene         22439007
##   ENST00000448914 ENST00000448914              TR_D_gene         22449113
##   ENST00000632524 ENST00000632524              IG_D_gene        105866322
##               ...             ...                    ...              ...
##   ENST00000444456 ENST00000444456 unprocessed_pseudogene             <NA>
##   ENST00000416377 ENST00000416377 unprocessed_pseudogene             <NA>
##   ENST00000318104 ENST00000318104 unprocessed_pseudogene             <NA>
##   ENST00000606902 ENST00000606902 unprocessed_pseudogene             <NA>
##   ENST00000427827 ENST00000427827 unprocessed_pseudogene             <NA>
##                   tx_cds_seq_end         gene_id tx_support_level
##                        <integer>     <character>        <integer>
##   ENST00000631435      142847317 ENSG00000282253             <NA>
##   ENST00000415118       22438554 ENSG00000223997             <NA>
##   ENST00000434970       22439015 ENSG00000237235             <NA>
##   ENST00000448914       22449125 ENSG00000228985             <NA>
##   ENST00000632524      105866332 ENSG00000282455             <NA>
##               ...            ...             ...              ...
##   ENST00000444456           <NA> ENSG00000237492             <NA>
##   ENST00000416377           <NA> ENSG00000232215             <NA>
##   ENST00000318104           <NA> ENSG00000177233             <NA>
##   ENST00000606902           <NA> ENSG00000271934             <NA>
##   ENST00000427827           <NA> ENSG00000227102             <NA>
##                       tx_id_version gc_content         tx_name
##                         <character>  <numeric>     <character>
##   ENST00000631435 ENST00000631435.1    83.3333 ENST00000631435
##   ENST00000415118 ENST00000415118.1    25.0000 ENST00000415118
##   ENST00000434970 ENST00000434970.2    55.5556 ENST00000434970
##   ENST00000448914 ENST00000448914.1    61.5385 ENST00000448914
##   ENST00000632524 ENST00000632524.1    54.5455 ENST00000632524
##               ...               ...        ...             ...
##   ENST00000444456 ENST00000444456.1    44.0307 ENST00000444456
##   ENST00000416377 ENST00000416377.2    45.3191 ENST00000416377
##   ENST00000318104 ENST00000318104.2    47.7392 ENST00000318104
##   ENST00000606902 ENST00000606902.2    45.9064 ENST00000606902
##   ENST00000427827 ENST00000427827.2    45.4545 ENST00000427827
##   -------
##   seqinfo: 455 sequences from GRCh38 genome
```

]
---

## Summarized experiment

---
## Summarized experiment
## Assay

.scroll-output[
- El *slot* `assay` almacena la información de las cuentas para cada transcrito dividida en tres niveles:

```r
assayNames(se)
```

```
## [1] "counts"    "abundance" "length"
```

- Para acceder a la matriz de cuentas estimadas por *Salmon*, correr:

```r
head(assay(se), 5)
```

```
##                 siRNA_NRF2_1 siRNA_control_1 siRNA_NRF2_2 siRNA_control_2
## ENST00000631435            0               0            0               0
## ENST00000415118            0               0            0               0
## ENST00000434970            0               0            0               0
## ENST00000448914            0               0            0               0
## ENST00000632524            0               0            0               0
##                 siRNA_NRF2_3 siRNA_control_3 CT_1 mut_NRF2_1 CT_2 mut_NRF2_2
## ENST00000631435            0               0    0          0    0          0
## ENST00000415118            0               0    0          0    0          0
## ENST00000434970            0               0    0          0    0          0
## ENST00000448914            0               0    0          0    0          0
## ENST00000632524            0               0    0          0    0          0
##                 CT_3 mut_NRF2_3
## ENST00000631435    0          0
## ENST00000415118    0          0
## ENST00000434970    0          0
## ENST00000448914    0          0
## ENST00000632524    0          0
```

]
---
## Summarized experiment
## Assay
- Las matriz de abundancia *(TPM)* puede obtenerse:

```r
## Obtener matriz de TPM
head(se@assays@data$abundance, 5)
```

```
##                 samp1_quant samp2_quant samp3_quant samp4_quant samp5_quant
## ENST00000631435           0           0           0           0           0
## ENST00000415118           0           0           0           0           0
## ENST00000434970           0           0           0           0           0
## ENST00000448914           0           0           0           0           0
## ENST00000632524           0           0           0           0           0
##                 samp6_quant samp7_quant samp8_quant samp9_quant samp10_quant
## ENST00000631435           0           0           0           0            0
## ENST00000415118           0           0           0           0            0
## ENST00000434970           0           0           0           0            0
## ENST00000448914           0           0           0           0            0
## ENST00000632524           0           0           0           0            0
##                 samp11_quant samp12_quant
## ENST00000631435            0            0
## ENST00000415118            0            0
## ENST00000434970            0            0
## ENST00000448914            0            0
## ENST00000632524            0            0
```
]

---
## Summarized experiment

```r
rownames(colData(se))
```

```
##  [1] "siRNA_NRF2_1"    "siRNA_control_1" "siRNA_NRF2_2"    "siRNA_control_2"
##  [5] "siRNA_NRF2_3"    "siRNA_control_3" "CT_1"            "mut_NRF2_1"     
##  [9] "CT_2"            "mut_NRF2_2"      "CT_3"            "mut_NRF2_3"
```

```r
colnames(assay(se))
```

```r
## Comprobar que rownames de colData es igual a colnames de assay
row.names(colData(se)) == colnames(assay(se))
```

```
##  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
```

]
]

---

# 10: Exploración de los datos
---

---
## Exploración de datos
## ¿Por qué es importante explorar los datos?

- Paso previo al análisis de expresión diferencial

- Análisis de calidad de los datos

- Permite conocer la congurencian entre individuos o réplicas

- Mediante gráficos, visualizar comportamiento de los datos <svg viewBox="0 0 320 512" style="height:1em;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M285.476 272.971L91.132 467.314c-9.373 9.373-24.569 9.373-33.941 0l-22.667-22.667c-9.357-9.357-9.375-24.522-.04-33.901L188.505 256 34.484 101.255c-9.335-9.379-9.317-24.544.04-33.901l22.667-22.667c9.373-9.373 24.569-9.373 33.941 0L285.475 239.03c9.373 9.372 9.373 24.568.001 33.941z"></path></svg> Outliers

---
## Exploración de datos

]
--
.pull-right[
<img src="./Images/Heatmap.png" width="400px" />

]

---
## Exploración de datos

---
## PCA

- Análisis de componentes principales

- Método algebraico para reducir la dimensionalidad de juegos de datos complejos (múltiples variables)

- Reducción de dimensionalidad o variables permite un análisis exploratorio más intuitivo

.blockquote[Reducción de variables implica  preservar y captar la mayor información sobre los datos
]

---
## PCA
## Normalización de los datos
--

Escalación de los datos <svg viewBox="0 0 320 512" style="height:1em;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M285.476 272.971L91.132 467.314c-9.373 9.373-24.569 9.373-33.941 0l-22.667-22.667c-9.357-9.357-9.375-24.522-.04-33.901L188.505 256 34.484 101.255c-9.335-9.379-9.317-24.544.04-33.901l22.667-22.667c9.373-9.373 24.569-9.373 33.941 0L285.475 239.03c9.373 9.372 9.373 24.568.001 33.941z"></path></svg> Normalización para hacerlos más comparables

Busqueda de correlación entre las variables <svg viewBox="0 0 320 512" style="height:1em;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M285.476 272.971L91.132 467.314c-9.373 9.373-24.569 9.373-33.941 0l-22.667-22.667c-9.357-9.357-9.375-24.522-.04-33.901L188.505 256 34.484 101.255c-9.335-9.379-9.317-24.544.04-33.901l22.667-22.667c9.373-9.373 24.569-9.373 33.941 0L285.475 239.03c9.373 9.372 9.373 24.568.001 33.941z"></path></svg> correlación positiva o negativa

---
## PCA

--
.blockquote[Los componentes principales son los nuevos ejes que maximizan la distancia de los datos al origen
]

---
## PCA
## Componentes principales

- El .green[primer componente] es aquel en el que los datos presentan la mayor separación (variación)

- El .orange[segundo componente] es el que tiene la segunda mayor separación entre los datos y es perpendicular al primer componente

- ¿Cuantos componentes existen? <svg viewBox="0 0 320 512" style="height:1em;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M285.476 272.971L91.132 467.314c-9.373 9.373-24.569 9.373-33.941 0l-22.667-22.667c-9.357-9.357-9.375-24.522-.04-33.901L188.505 256 34.484 101.255c-9.335-9.379-9.317-24.544.04-33.901l22.667-22.667c9.373-9.373 24.569-9.373 33.941 0L285.475 239.03c9.373 9.372 9.373 24.568.001 33.941z"></path></svg> tantas variables en el juego de datos

---
## PCA
## Componentes principales

- Cada componente principal tiene asociado un eigenvector (vector unitario) y un eigenvalue (cantidad escalar)

- Los eigenvalues son la suma del cuadrado de las distancias de los puntos proyectados sobre dicho componente principal

- Los eigenvalues permiten seleccionar cuál es el componente principal que explica la mayor variación en los datos

---
## PCA

<img src="./Images/Screeplot.png" width="450px" style="display: block; margin: auto;" />
Screeplot

---
## PCA

- En sets de datos de RNAseq, el PCA permite visualizar la distancia o congruencia de los datos

- Generalmente son gráficas de puntos (muestras) en dos dimensiones (dos componentes) que resumen las principales fuentes de varianza

--
.center[
<img src="./Images/PCARNAseq.png" width="400px" />
]

---
class: center, middle
background-image: url("./Images/hands_on.png")
background-size: 90% 80%
---
## Exploración de los datos
### Creación del objeto DDS

1. Crear un objeto `DESeq` (*SummarizedExperiment*) a partir del objeto `gse`

```r
dds <- DESeqDataSet(se = gse, ##Objeto SummarizedExperiment que contiene la matriz de cuentas crudas
                    design = ~ Treatment) ##Variable con los grupos experimentales a contrastar. Forma parte de las columnas de coldata
```

```
## using counts and average transcript lengths from tximeta
```

2. Visualizar el grupo de referencia o **control** de la variable seleccionada

```r
##El  primer nivel o elemento es el grupo control  ("CT")
levels(dds$Treatment)
```

```
## [1] "CT"            "mut_NRF2"      "siRNA_control" "siRNA_NRF2"
```

3. Dependiendo de la pregunta a responder, recodifica el  grupo de referencia:

```r
dds$Treatment <- relevel(dds$Treatment, "siRNA_control")

levels(dds$Treatment)
```

```
## [1] "siRNA_control" "CT"            "mut_NRF2"      "siRNA_NRF2"
```

4. Ya que los genes que muestran bajo número de cuentas (abundancia) no nos interesan en el análisis, selecciona aquellos que muestran mayor abundancia usando un punto de corte

```r
##Recupera los genes que muestran más de 1 cuenta en al menos 3 muestras
keep <- rowSums(counts(dds) >= 1) >=3

##Visualiza cuántos genes pasaron el filtro de abundancia
table(keep)
```

```
## keep
## FALSE  TRUE 
## 16846 21299
```

```r
##Corta los genes que pasaron el filtro de abundancia del objeto dds
dds <- dds[keep, ]
```

]

---
## Exploración de los datos
### PCA

5. Para explorar los datos por medio del PCA necesitamos normalizar las cuentas crudas para que  los datos sean más comparables entre sí. Empleamos la transformación *regularized-logarithm transformation* recomendada para juegos de datos con una n < 30:

```r
rld <-  rlog(dds, blind = F)
```

```
## using 'avgTxLength' from assays(dds), correcting for library size
```

6. Realiza el PCA empleando el paquete `PCAtools`:

```r
##Indicarle a  la función pca que los datos deben de ser escalados para que sean más comparables entre sí
pca <- pca(assay(rld), metadata = colData(rld), scale = T)

##Graficar el resultado del PCA mostrando las etiquetas de acuerdo al nombr + tratamiento + réplica y coloreando los puntos de acuerdo al tipo de tratamiento
biplot(pca, lab = rownames(colData(rld)), colby = "Treatment")
```

![](Slides_files/figure-html/unnamed-chunk-80-1.png)

]

---
class: inverse, center, middle

# 11: Análisis de expresión diferencial
---
## Análisis de expresión diferencial (DE)

Prueba estadística para obtener los genes que varían entre las condiciones experimentales

---

## Pasos en el DE

---
## Normalización

Es el primer paso del análisis de expresión diferencial y es necesario para realizar comparaciones acertadas entre muestras.

Las cuentas crudas están conformadas por un componente **"interesante"** (la expresión de RNA) y componentes "no interesantes" (como los efectos de lote, ruido de la plataforma, etc.).

La normalzación escala las cuentas para tratar de reducir los componentes "no interesantes" y poder comparar las muestras entre si.

---

## Criteríos para normalizar

Se puede normalizar considerando:

- La profundidad (tamaño de librería)

---

## Criteríos para normalizar

Se puede normalizar considerando:

- El tamaño del gen

---

## Criteríos para normalizar

Se puede normalizar considerando:

- Composición de RNA

---
## Métodos comunes

| Método de normalización | Descripción | Factores de evaluación | Recomendaciones de uso |
|:-----------------------------------------------------------------------------------:|:----------------------------------------------------------------------------------------------------------------------------:|:--------------------------------------------------:|:-------------------------------------------------------------------------------------------------------------------:|
| CPM (cuentas por millón): cuentas escalados por el número total de lecturas | profundidad de secuenciación | comparaciones de cuentas de genes entre réplicas del mismo grupo de muestras| NO para comparaciones dentro de la muestra o análisis de DE.| 
| TPM (transcritos por kilobase por millón de lecturas): cuentas por longitud de transcripción (kb) por millón de lecturas mapeadas | profundidad de secuenciación y longitud de genes| comparaciones de cuentas de genes dentro de una muestra o entre muestras del mismo grupo de muestras| NO para análisis de DE. |
| RPKM/FPKM (lecturas/fragmentos por kilobase de exón por millón de lecturas/fragmentos mapeados| similar a TPM, profundidad de secuenciación y longitud del gen | comparaciones de cuentas entre genes dentro de una muestra | NO para comparaciones entre muestras o análisis de DE.|
| Mediana de ratios de DESeq2 | cuentas divididas por factores de tamaño específicos de la muestra determinados por la mediana del ratio de cuentas de genes en relación con la media geométrica por gen | profundidad de secuenciación y composición del RNA | comparaciones de cuentas de genes entre muestras y para el análisis de DE; NO para comparaciones dentro de la muestra |
| La media cortada de los valores M de EdgeR (TMM) | utiliza una media recortada ponderada de los ratios de expresión logarítmica entre las muestras | profundidad de secuenciación | composición de RNA y longitud de los genes. |

]

---
## DESeq2

- Normalización para análisis de expresión diferencial:

- Factor de normalización de `DESeq2`

- Normalización para visualización u otras aplicaciones:

- Variance stabilizing transformation (VST)
  
  - Regularized-logarithm transformation (rlog)

---
## Normalización por ratios de medias geométricas (DESeq2)

`DESeq2`ajusta a un modelo lineal generalizado (GLM) de la familia binomial negativa (NB).

```
## # A tibble: 2 × 3
##   gene  muestraA muestraB
##   <chr>    <dbl>    <dbl>
## 1 gene1     1749      943
## 2 gene2       35       29
```

1. Crea una pseudo-referencia por muestra (promedio geometrico por fila) `sqrt(muestraA * muestra B)`

```
## # A tibble: 2 × 4
## # Rowwise: 
##   gene  muestraA muestraB prom_geom
##   <chr>    <dbl>    <dbl>     <dbl>
## 1 gene1     1749      943    1284. 
## 2 gene2       35       29      31.9
```

2. Se calcula la fración `muestra/pseudo-referencia`

```
## # A tibble: 2 × 6
## # Rowwise: 
##   gene  muestraA muestraB prom_geom muestraA_pseudo_ref muestraB_pseudo_ref
##   <chr>    <dbl>    <dbl>     <dbl>               <dbl>               <dbl>
## 1 gene1     1749      943    1284.                 1.36               0.734
## 2 gene2       35       29      31.9                1.10               0.910
```

3. Se calcula un factor de normalización (size factor) utilizando la `mediana` por columnas.

```
## [1] 1.230235
```

```
## [1] 0.8222689
```

4. Se dividen las `cuentas crudas/size factor` para calcular las cuentas normalizadas.

```
## # A tibble: 2 × 3
## # Rowwise: 
##   gene  muestraA muestraB
##   <chr>    <dbl>    <dbl>
## 1 gene1     1749      943
## 2 gene2       35       29
```

```
## # A tibble: 2 × 2
## # Rowwise: 
##   norm_muestraA norm_muestraB
##           <dbl>         <dbl>
## 1        1422.         1147. 
## 2          28.4          35.3
```

]
---

## Otras transformaciones

Puedes realizar otras transformaciones en `DESeq2` para estabilizar la varianza a través de los differentes valores promedio de expresión.

![](Slides_files/figure-html/unnamed-chunk-91-1.png)![](Slides_files/figure-html/unnamed-chunk-91-2.png)

```r
# variance stabilizing transformation (VST), (Anders and Huber 2010)
vsd <- vst(dds, blind = FALSE)
```

```
## using 'avgTxLength' from assays(dds), correcting for library size
```

```r
head(assay(vsd), 3)
```

```
##                 siRNA_NRF2_1 siRNA_control_1 siRNA_NRF2_2 siRNA_control_2
## ENSG00000000003     11.67435        11.63822     11.72695        11.73399
## ENSG00000000419     11.71082        11.78356     11.74230        11.78642
## ENSG00000000457     10.52426        10.50532     10.46261        10.51811
##                 siRNA_NRF2_3 siRNA_control_3     CT_1 mut_NRF2_1     CT_2
## ENSG00000000003     11.84279        11.68878 11.14964   11.17933 11.09670
## ENSG00000000419     11.80449        11.70709 11.48465   11.46978 11.40275
## ENSG00000000457     10.54239        10.52716 10.46414   10.47452 10.46257
##                 mut_NRF2_2     CT_3 mut_NRF2_3
## ENSG00000000003   11.24472 11.17600   11.12306
## ENSG00000000419   11.46835 11.50914   11.50147
## ENSG00000000457   10.50666 10.46139   10.45044
```

```r
# regularized-logarithm transformation (rlog), (Love, Huber, and Anders 2014)
rld <- rlog(dds, blind = FALSE)
```

```
## using 'avgTxLength' from assays(dds), correcting for library size
```

```r
head(assay(rld), 3)
```

```
##                 siRNA_NRF2_1 siRNA_control_1 siRNA_NRF2_2 siRNA_control_2
## ENSG00000000003    10.876558       10.819047    10.958594       10.969500
## ENSG00000000419    10.944160       11.056010    10.992956       11.060526
## ENSG00000000457     8.345928        8.289919     8.160224        8.327933
##                 siRNA_NRF2_3 siRNA_control_3      CT_1 mut_NRF2_1      CT_2
## ENSG00000000003    11.134274       10.899100  9.957696  10.016126  9.851707
## ENSG00000000419    11.087811       10.938303 10.577168  10.551995 10.436115
## ENSG00000000457     8.398004        8.354244  8.165047   8.197047  8.160009
##                 mut_NRF2_2      CT_3 mut_NRF2_3
## ENSG00000000003  10.140861 10.009398   9.905331
## ENSG00000000419  10.549557 10.618447  10.605587
## ENSG00000000457   8.293887  8.156555   8.122530
```

```r
# Normalización con el factor de normalizacion 
dds <- estimateSizeFactors(dds)
```

```
## using 'avgTxLength' from assays(dds), correcting for library size
```

```r
# Juntar los datos de las tres normalizaciones
df <- bind_rows(
  as_data_frame(log2(counts(dds, normalized=TRUE)[, 1:2]+1)) %>%
         mutate(transformation = "log2(x + 1)"),
  as_data_frame(assay(vsd)[, 1:2]) %>% mutate(transformation = "vst"),
  as_data_frame(assay(rld)[, 1:2]) %>% mutate(transformation = "rlog"))
```

```
## Warning: `as_data_frame()` was deprecated in tibble 2.0.0.
## Please use `as_tibble()` instead.
## The signature and semantics have changed, see `?as_tibble`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
```

```r
# Renombrar columnas
colnames(df)[1:2] <- c("x", "y")

# Nombre de las graficas
lvls <- c("log2(x + 1)", "vst", "rlog")

# Agrupar los tres tipos de normalizacion en grupos como factores
df$transformation <- factor(df$transformation, levels=lvls)

# Plotear los datos
ggplot(df, aes(x = x, y = y)) + geom_hex(bins = 80) +
  coord_fixed() + facet_grid( . ~ transformation)  
```

![](Slides_files/figure-html/unnamed-chunk-92-1.png)

]

---
## Dispersión

La dispersión es una medida de la variabilidad de los datos como (varianza, sd, etc.). En `DESeq2` se utiliza `α`:

$$ \alpha \propto 1/mean $$ 
$$ \alpha \propto variance $$

Así que **La dispersión es mayor para genes con poca abundancia y menor para genes abundantes**.

Y la **dispersión refleja la varianza**.

---

## Ajustando la dispersión

---

## Reduciendo la dispersión

---
class: center, middle
background-image: url("./Images/hands_on.png")
background-size: 90% 80%
---

## Ejercicio

La librería de `DESeq2` normaliza y realiza el análisis de expresión diferencial en una sola función. Así que ya realizamos el análisis, vamos a verlo:

1. Corre la función `DESeq()` para normalizar, estimar la dispersión y ajustar los datos al modelo binomial negativo.

```r
dds  <- DESeq(dds)
```

```
## using pre-existing normalization factors
```

```
## estimating dispersions
```

```
## gene-wise dispersion estimates
```

```
## mean-dispersion relationship
```

```
## final dispersion estimates
```

```
## fitting model and testing
```

2. Para obtener los resultados del análisis de expresión diferencial usa la función `results`

```r
##Usemos la función results con los parámetros por default
res <- results(dds)

##Visualiza el objeto que contiene los resultados
res
```

```
## log2 fold change (MLE): Treatment siRNA NRF2 vs siRNA control 
## Wald test p-value: Treatment siRNA NRF2 vs siRNA control 
## DataFrame with 21299 rows and 6 columns
##                    baseMean log2FoldChange     lfcSE       stat      pvalue
##                   <numeric>      <numeric> <numeric>  <numeric>   <numeric>
## ENSG00000000003 1497.757504      0.0982849 0.0902643   1.088857 2.76217e-01
## ENSG00000000419 1789.280526     -0.0107935 0.0702472  -0.153650 8.77885e-01
## ENSG00000000457  304.687311     -0.0243559 0.1151762  -0.211466 8.32523e-01
## ENSG00000000460  467.602306      0.5142010 0.1159843   4.433368 9.27721e-06
## ENSG00000000938    0.618937      1.3390032 2.8076874   0.476906 6.33429e-01
## ...                     ...            ...       ...        ...         ...
## ENSG00000288690    8.913454      0.0326943  0.663455  0.0492789    0.960697
## ENSG00000288694    0.258397      0.8179394  4.365460  0.1873661    0.851374
## ENSG00000288695   11.159112     -0.3769746  0.701013 -0.5377567    0.590745
## ENSG00000288698   23.403296      0.5109625  1.060520  0.4818036    0.629946
## ENSG00000288699    0.869679      1.5093784  3.007207  0.5019203    0.615724
##                        padj
##                   <numeric>
## ENSG00000000003 0.476825323
## ENSG00000000419 0.948127083
## ENSG00000000457 0.924197635
## ENSG00000000460 0.000121719
## ENSG00000000938          NA
## ...                     ...
## ENSG00000288690    0.980632
## ENSG00000288694          NA
## ENSG00000288695    0.768387
## ENSG00000288698    0.795601
## ENSG00000288699          NA
```

```r
##Obten el número de genes sobre- y sub-expresados de acuerdo a la prueba estadística
summary(res)
```

```
## 
## out of 21299 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up)       : 2271, 11%
## LFC < 0 (down)     : 2660, 12%
## outliers [1]       : 26, 0.12%
## low counts [2]     : 5778, 27%
## (mean count < 8)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
```

3. Explora el objeto para tener información del contenido de cada columna:

```r
mcols(res, use.names = TRUE)
```

```
## DataFrame with 6 rows and 2 columns
##                        type            description
##                 <character>            <character>
## baseMean       intermediate mean of normalized c..
## log2FoldChange      results log2 fold change (ML..
## lfcSE               results standard error: Trea..
## stat                results Wald statistic: Trea..
## pvalue              results Wald test p-value: T..
## padj                results   BH adjusted p-values
```
]

---

## Análisis de expresión diferencial

**_Warning_**

.content-box-red[Los parámetros por default de la función `res()` hacen referencia a la siguiente hipótesis nula:]

Si deseamos modificar el *cutoff* o valor de corte del log2FoldChange:

---

```r
##Emplea los argumentos lfcthreshold y alpha para especificar ambos puntos de corte de forma simultánea
res_LFC <- results(dds, lfcThreshold = 1, alpha = 0.05)

res_LFC
```

```
## log2 fold change (MLE): Treatment siRNA NRF2 vs siRNA control 
## Wald test p-value: Treatment siRNA NRF2 vs siRNA control 
## DataFrame with 21299 rows and 6 columns
##                    baseMean log2FoldChange     lfcSE      stat    pvalue
##                   <numeric>      <numeric> <numeric> <numeric> <numeric>
## ENSG00000000003 1497.757504      0.0982849 0.0902643  0.000000  1.000000
## ENSG00000000419 1789.280526     -0.0107935 0.0702472  0.000000  1.000000
## ENSG00000000457  304.687311     -0.0243559 0.1151762  0.000000  1.000000
## ENSG00000000460  467.602306      0.5142010 0.1159843  0.000000  1.000000
## ENSG00000000938    0.618937      1.3390032 2.8076874  0.120741  0.903896
## ...                     ...            ...       ...       ...       ...
## ENSG00000288690    8.913454      0.0326943  0.663455  0.000000  1.000000
## ENSG00000288694    0.258397      0.8179394  4.365460  0.000000  1.000000
## ENSG00000288695   11.159112     -0.3769746  0.701013  0.000000  1.000000
## ENSG00000288698   23.403296      0.5109625  1.060520  0.000000  1.000000
## ENSG00000288699    0.869679      1.5093784  3.007207  0.169386  0.865493
##                      padj
##                 <numeric>
## ENSG00000000003         1
## ENSG00000000419         1
## ENSG00000000457         1
## ENSG00000000460         1
## ENSG00000000938         1
## ...                   ...
## ENSG00000288690         1
## ENSG00000288694         1
## ENSG00000288695         1
## ENSG00000288698         1
## ENSG00000288699         1
```

```r
summary(res_LFC)
```

```
## 
## out of 21299 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 1.00 (up)    : 7, 0.033%
## LFC < -1.00 (down) : 8, 0.038%
## outliers [1]       : 26, 0.12%
## low counts [2]     : 0, 0%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
```

**NUNCA** realices la prueba de expresión diferencial con un valor de corte X (log2FoldChange) y posteriormente filtres los genes diferencialmente expresados empleando otro valor de log2FoldChange.

]
---
class: inverse, center, middle

# 12: Visualización de los resultados
---

---
.center[
### ¿Cuántos tipos de gráficos conoces para visualizar resultados del análisis DE?
]

---
# MAplot
.scroll-output[

Gráfico que representa la distribución de los genes o transcritos en las comparaciones de interés

**M** eje y de *minus*:

**A** eje x de *average* -> Promedio de las cuentas normalizadas para cada gen en las condiciones experimentales de la o las variables de interés

Para generar el **MAplot** usa:

```r
##Es importante que recuerden que la hipótesis nula que se probó fue
##"El lfc del gen n es igual a 0" por lo tanto los genes coloreados son...
plotMA(res)
```

![](Slides_files/figure-html/unnamed-chunk-103-1.png)
]
]

---
# MAplot

- Son poco abundantes

- Tienen alta dispersión

Usemos la función `lfcShrink` de la librería *.red[apeglm]*

```r
res_shrink <- lfcShrink(dds = dds, res = res, coef = "Treatment_siRNA_NRF2_vs_siRNA_control"
, type = "apeglm")
```

```
## using 'apeglm' for LFC shrinkage. If used in published research, please cite:
##     Zhu, A., Ibrahim, J.G., Love, M.I. (2018) Heavy-tailed prior distributions for
##     sequence count data: removing the noise and preserving large differences.
##     Bioinformatics. https://doi.org/10.1093/bioinformatics/bty895
```

```r
summary(res_shrink)
```

```r
plotMA(res_shrink,  ylim = c(-3, 3))
```

![](Slides_files/figure-html/unnamed-chunk-104-1.png)
]

---
# Volcano plot

.scroll-output[
De manera similar al MAplot con el volcano plot visualizamos los genes que muestran expresión diferencial en nuestra condición de interés

- En el eje y se grafica el -log10 de padj

- En el eje x se grafica el lfc o *log2foldchange*

Usa la función personalizada volcanoplotR del script de funciones:

```r
##En los argumentos logfc y padj utiliza los mismos valores empleados en la función results
res_shrink <- as.data.frame(res_shrink)
volcanoplotR(res_shrink, logfc = 0, p.adj = 0.1)
```

![](Slides_files/figure-html/unnamed-chunk-105-1.png)
]

---

## Heatmap

El *h   eatmap* nos permite visualizar la expresión de los genes diferencialmente expresados en terminos de las cuentas normalizadas

Consideraciones:

- Usar los valores de las cuentas normalizadas para una mejor comparación entre muestras

- Escalar los valores de las cuentas (renglones) para visualizar las diferencias en la expresión

Usaremos la librería de *.red[pheatmap]*

---
# Heatmap

1. Obten la matriz de cuentas normalizadas:

```r
norm_counts <- counts(dds, normalized = T)
```

2. Guarda un data.frame con los resultados de los genes mostrando expresión diferencial significativa:

```r
deg <- res_shrink %>% dplyr::filter(log2FoldChange > 0 & padj < 0.1 |
                                                     log2FoldChange < 0 & padj < 0.1)

head(deg)
```

```
##                  baseMean log2FoldChange      lfcSE       pvalue         padj
## ENSG00000000460  467.6023      0.4721205 0.11836789 9.277210e-06 0.0001217192
## ENSG00000001084 1889.3028     -0.2991140 0.09464064 4.368589e-04 0.0031869719
## ENSG00000001167 1068.3113      0.2413441 0.08978195 2.877301e-03 0.0149359373
## ENSG00000001461  564.0667     -0.2490246 0.08409560 1.186278e-03 0.0073086983
## ENSG00000001561  612.2146     -0.2536187 0.10372375 5.197697e-03 0.0240905329
## ENSG00000001617  996.1836      0.2917830 0.09054617 3.973259e-04 0.0029527889
```

3. Usando la librería de `pheatmap()` obten la gráfica

```r
##Genera una paleta de colores 
RedBu <- brewer.pal(n = 10, name = "RdBu")
RedBu <- rev(RedBu)

pheatmap(norm_counts[rownames(deg), 1:6], ##Corta los genes con expresión diferencial de la matriz normalizada
         scale = "row", ##Calcula los valores z por fila (gen)
         color = RedBu, ##Nombre de la paleta de colores
         show_rownames = F, ##No mostrar el nombre de los genes en la gráfica
         show_colnames = F, ##No mostrar el nombre de las muestras en la gráfica
         clustering_distance_rows = "euclidean", ##Especifíca la medida de la distancia entre filas
         clustering_distance_cols = "euclidean", ##O columnas
         clustering_method = "single", ##Método de agrupamiento jerárquico
         border_color = NA, 
         annotation_col = coldata[1:6, c(2, 6)]) ##Añadir información relacionada a las columnas)
```

<img src="Slides_files/figure-html/unnamed-chunk-108-1.png" style="display: block; margin: auto;" />
]

---

---

## Anotación funcional

¿Cuáles son los métodos más comunes para anotar funcionalmente los genes diferencialmente expresados?

---

## ORA

- Primer método para realizar análisis de vías

- Requiere de una lista de genes diferencialmente expresados

- Archivo con la colección de vías o sets de genes curados manualmente

---
class: center, middle
background-image: url("./Images/hands_on.png")
background-size: 90% 90%
---

## ORA

**Pasos para realizar ORA**

1. Obten de  manera independiente la lista de los genes sobre- y sub-expresados

```r
up <- deg %>% dplyr::filter(log2FoldChange > 0)
down <- deg %>% dplyr::filter(log2FoldChange < 0)
```

2. Convierte los ensembl id a entrez_gene id y hgnc_symbol

```r
##Descarga los datos desde el sitio de ensembl usando el paquete de biomaRt
ensembl <- useMart("ensembl")
ensembl <- useDataset("hsapiens_gene_ensembl", mart = ensembl)

##Con la función id_converter() obten los entrezgene_id y hgnc_symbol
##Asigna los nombres de las filas (rownames) a una columna
up <- rownames_to_column(up, var = "ensembl_gene_id")
up <- id_converter(mart = ensembl, ##Objeto con la información de la base de datos  de ensembl
            input = up, ##Lista de genes de entrada
            attributes = c("ensembl_gene_id", "entrezgene_id", "hgnc_symbol"), ##Nombre de los atributos que se desean obtener
            filter = "ensembl_gene_id")
```

3. Carga la colección de vías o set de genes:

```r
go <- read.gmt(here("data/c5.go.bp.v7.4.entrez.gmt"))
```

4. Ejecuta el ORA:

```r
ego_up <- enricher(gene = up$entrezgene_id, ##Vector con los nombres o id de los genes
                   TERM2GENE = go, ##Colección de vías o sets de genes
                   pvalueCutoff = 0.05, ##Punto de corte para obtener vías representadas significativamente
                   pAdjustMethod = "BH") ##Método para ajustar los valores de p
```

]

---

## ORA

Visualización de los resultados:

Las librerías **clusterProfiler** y **enrichplot** cuentan con una variedad de funciones para visualizar los resultados:

No aglomerativos

**Dotplot**:

```r
dotplot(ego_up, showCategory = 10)
```

![](Slides_files/figure-html/unnamed-chunk-113-1.png)

**Barplot**:

```r
barplot(ego_up, showCategory = 10)
```

![](Slides_files/figure-html/unnamed-chunk-114-1.png)

]

Aglomerativos

**emapplot**

```r
emapplot(ego_up, showCategory = 10)
```

![](Slides_files/figure-html/unnamed-chunk-116-1.png)

**Treeplot**

```r
treeplot(ego_up)
```

```
## Warning: The "label" has(have) been found in tree data. You might need to rename the
## variable(s) in the data of "geom_cladelab" to avoid this warning!
```

```
## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.

## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.
```

![](Slides_files/figure-html/unnamed-chunk-117-1.png)

]

---

## ORA

Ventajas:

]

]

---

## GSEA

GSEA surge como una alternativa al ORA ya que resuelve alguna de sus limitaciones. Para ello requiere de:

- La lista completa de genes que se sometieron a un análisis de expresión diferencial. En este caso el objeto `res_shrink`
- El cálculo de una métrica que considere o agregue los valores de la expresión y significancia estadística de cada gen.
- La lista ordenada o *rankeada* los genes de acuerdo a la métrica calculada
- Nuevamente, la colección de vías o set de genes

---

## GSEA

Pasos para correr el GSEA:

1. Convierte los id de los genes, empleando el objeto `res_shrink`

```r
res_shrink <- as.data.frame(res_shrink)
res_shrink <- rownames_to_column(res_shrink, var = "ensembl_gene_id")

res_shrink <- id_converter(mart = ensembl, 
                           input = res_shrink, 
                           attributes = c("ensembl_gene_id", "entrezgene_id", "hgnc_symbol"),  
                           filter = "ensembl_gene_id")
```

2. Calcula la métrica que agrega el log2FoldChange y padj. En este caso usa el producto de -log10(padj)*log2FoldChange

```r
##Agrega una nueva columna al objeto res_shrink que contenga los valores de la métrica
res_shrink <- dplyr::mutate(res_shrink, stat = -log10(res_shrink$padj)*res_shrink$log2FoldChange)
```

3. Ordena de forma descendente, respecto a la métrica, los genes:

```r
res_shrink <- dplyr::arrange(res_shrink, desc(stat))
```

4. Obten un vector nombrado que contenga el valor de la métrica y el id de los genes:

```r
genes <- res_shrink$stat
names(genes) <- res_shrink$entrezgene_id
```

5. Carga la colección de vías:

```r
go <- gmtPathways(here("data/c5.go.bp.v7.4.entrez.gmt"))
```

6. Ejecuta el análisis:

```r
GSEA_res <- fgseaMultilevel(go, ##Objeto con la colección de vías o sets de genes
                            genes, ##Vector nombrado de los genes rankeados
                            maxSize = 500, minSize = 15)
```

```
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (27.25% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
```

```
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize,
## gseaParam, : There are duplicate gene names, fgsea may produce unexpected
## results.
```

]

---

## GSEA

Emplea un barplot para revisar los resultados del GSEA

```r
##Convierte los resultados de GSEA a un objeto de clase tidy y añade una nueva columna que nos indique si la vía está enriquecida de forma significativa
GSEA_res <- as.tibble(GSEA_res) %>%
  arrange(desc(NES)) %>%
  mutate(Significance = ifelse(padj < 0.05, "Significant", "NS"))
```

```
## Warning: `as.tibble()` was deprecated in tibble 2.0.0.
## Please use `as_tibble()` instead.
## The signature and semantics have changed, see `?as_tibble`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
```

```r
##Genera un barplot con los resultados
barplot_GSEA(GSEA_res,  Head = 20, Tail = 20)
```

![](Slides_files/figure-html/unnamed-chunk-124-1.png)
]