Workshop: Análisis de datos de RNA-Seq

class: title-slide, middle, center
background-image: url(data:image/png;base64,#figures/Clases_RNASeq_Slide1.png) 
background-position: 90% 75%, 75% 75%, center
background-size: 1210px,210px, cover

####Dra. Evelia Coss, Sofia Salazar y Diego Ramirez 
#### 28/02/2024
]

.left[.footnote[R-Ladies Theme[R-Ladies Theme](https://www.apreshill.com/project/rladies-xaringan/)]]

---
background-image: url(data:image/png;base64,#figures/liigh_unam_logo.png) 
background-position: 10% 10%
background-size: 10%
class: middle, center

# Sobre mi equipo
----
.left-col[ 
### <svg aria-hidden="true" role="img" viewBox="0 0 448 512" style="height:1em;width:0.88em;vertical-align:-0.125em;margin-left:auto;margin-right:auto;font-size:inherit;fill:#181818;overflow:visible;position:relative;"><path d="M416 0c17.7 0 32 14.3 32 32c0 59.8-30.3 107.5-69.4 146.6c-28 28-62.5 53.5-97.3 77.4l-2.5 1.7c-11.9 8.1-23.8 16.1-35.5 23.9l0 0 0 0 0 0-1.6 1c-6 4-11.9 7.9-17.8 11.9c-20.9 14-40.8 27.7-59.3 41.5H283.3c-9.8-7.4-20.1-14.7-30.7-22.1l7-4.7 3-2c15.1-10.1 30.9-20.6 46.7-31.6c25 18.1 48.9 37.3 69.4 57.7C417.7 372.5 448 420.2 448 480c0 17.7-14.3 32-32 32s-32-14.3-32-32H64c0 17.7-14.3 32-32 32s-32-14.3-32-32c0-59.8 30.3-107.5 69.4-146.6c28-28 62.5-53.5 97.3-77.4c-34.8-23.9-69.3-49.3-97.3-77.4C30.3 139.5 0 91.8 0 32C0 14.3 14.3 0 32 0S64 14.3 64 32H384c0-17.7 14.3-32 32-32zM338.6 384H109.4c-10.1 10.6-18.6 21.3-25.5 32H364.1c-6.8-10.7-15.3-21.4-25.5-32zM109.4 128H338.6c10.1-10.7 18.6-21.3 25.5-32H83.9c6.8 10.7 15.3 21.3 25.5 32zm55.4 48c18.4 13.8 38.4 27.5 59.3 41.5c20.9-14 40.8-27.7 59.3-41.5H164.7z"/></svg> Sofía Salazar Magaña

Estudiante de la [Licenciatura en Ciencias Genómicas](https://www.enesjuriquilla.unam.mx/?page_id=2136) @ ENES-J UNAM

Asistente de investigación @ [RegGenoLab](https://liigh.unam.mx/), LIIGH UNAM

Asistente de investigación @ [Im Lab](https://hakyimlab.org/), The University of Chicago

[<svg aria-hidden="true" role="img" viewBox="0 0 496 512" style="height:1em;width:0.97em;vertical-align:-0.125em;margin-left:auto;margin-right:auto;font-size:inherit;fill:currentColor;overflow:visible;position:relative;"><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"/></svg> GitHub: SophSM](https://github.com/SophSM)
]

.center-col[
### <svg aria-hidden="true" role="img" viewBox="0 0 448 512" style="height:1em;width:0.88em;vertical-align:-0.125em;margin-left:auto;margin-right:auto;font-size:inherit;fill:#181818;overflow:visible;position:relative;"><path d="M416 0c17.7 0 32 14.3 32 32c0 59.8-30.3 107.5-69.4 146.6c-28 28-62.5 53.5-97.3 77.4l-2.5 1.7c-11.9 8.1-23.8 16.1-35.5 23.9l0 0 0 0 0 0-1.6 1c-6 4-11.9 7.9-17.8 11.9c-20.9 14-40.8 27.7-59.3 41.5H283.3c-9.8-7.4-20.1-14.7-30.7-22.1l7-4.7 3-2c15.1-10.1 30.9-20.6 46.7-31.6c25 18.1 48.9 37.3 69.4 57.7C417.7 372.5 448 420.2 448 480c0 17.7-14.3 32-32 32s-32-14.3-32-32H64c0 17.7-14.3 32-32 32s-32-14.3-32-32c0-59.8 30.3-107.5 69.4-146.6c28-28 62.5-53.5 97.3-77.4c-34.8-23.9-69.3-49.3-97.3-77.4C30.3 139.5 0 91.8 0 32C0 14.3 14.3 0 32 0S64 14.3 64 32H384c0-17.7 14.3-32 32-32zM338.6 384H109.4c-10.1 10.6-18.6 21.3-25.5 32H364.1c-6.8-10.7-15.3-21.4-25.5-32zM109.4 128H338.6c10.1-10.7 18.6-21.3 25.5-32H83.9c6.8 10.7 15.3 21.3 25.5 32zm55.4 48c18.4 13.8 38.4 27.5 59.3 41.5c20.9-14 40.8-27.7 59.3-41.5H164.7z"/></svg> Diego Ramirez

Estudiante de la [Licenciatura en Ciencias Genómicas](https://www.enesjuriquilla.unam.mx/?page_id=2136) @ ENES-J UNAM

]

.right-col[
<img src="figures/LAB-logo.png", height = "150">
<img src = "figures/lupus-logo.png">
<img src = "figures/jaguar-logo.png">
]
---

---

## Pipeline bioinformática para RNA-seq

.pull-right[
<img src="data:image/png;base64,#figures/pipeline1.png" width="80%" style="display: block; margin: auto;" />
]

.left[.footnote[.black[
Imagen proveniente de [mRNA-Seq data analysis workflow](https://biocorecrg.github.io/RNAseq_course_2019/workflow.html)
]]]

---

---

### 1. Descarga de datos públicos de RNA-seq con `wget`

### 2. Análisis de control de calidad de lecturas

### 3. Trimming

### 4. Alineamiento con el genoma de referencia mediante STAR

---

---
### ¿Cómo podemos conseguir datos públicos de RNA-seq?

La forma más simple es ir a repositorios de datos públicos, como [GEO (Gene Expression Omnibus)](https://www.ncbi.nlm.nih.gov/geo/), en donde encontraremos los archivos de datos **crudos** y a veces también las **matrices de cuentas ya procesadas** ó a [Recount3](https://rna.recount.bio/) (aquí podemos encontrar datos ya procesados).

**Para esta clase, usaremos las muestras del estudio [GSE155516](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE155516)** provenientes de [Kodali, *et al*. 2021. *Mol Psychiatry*](https://www.nature.com/articles/s41380-020-00955-5).

### Descarga de los datos con `wget`

Una vez que seleccionamos el estudio que vamos a usar, entramos a la página de [ebi](https://www.ebi.ac.uk/). En el buscador, tecleamos el ID del estudio.

Seleccionamos el estudio donde dice **Nucleotide sequences**, y en la tabla inferior, en la columna **Generated FASTQ files:FTP**, seleccionamos las muestras que queremos y damos click izquierdo en donde dice **Get download script**.

---
### Descarga de los datos con `wget`

Esto descargará un script BASH `.sh` como [este](https://github.com/EveliaCoss/RNAseq_classFEB2024/blob/main/Practica_Dia2/scripts/download.sh), el cual utilizar para descargar las muestras.

Para correr este script y descargar las muestras, debemos ir a la carpeta donde las queremos guardar y ahí guardamos el script. Supongamos que yo renombré mi script a [`download.sh`](https://github.com/EveliaCoss/RNAseq_classFEB2024/blob/main/Practica_Dia2/scripts/download.sh). Y lo tengo en una carpeta llamada `scripts`.

```bash
cd /mnt/Guanina/bioinfo24/data/Clase_RNASeq2024/scripts/
chmod +x downloads.sh
./download.sh
```

Para correr el script, primero me doy permisos de ejecución con `chmod +x` y después ejecuto el script con `./`

> NOTA: Recuerda darle permisos a tu companero:

```bash
setfacl -R -m u:usuario:rwx CARPETA
```

---
### Descarga de los datos con `wget`

Ya que las descargas duran mucho es conveniente mandar el script a la cola de ejecución usando un Sun Grid Engine (SGE) `.sge`. Así el script seguirá ejecutándose aunque tu no estes.

Para lo cual hay 2 formas:
A)	Tener un solo script donde todos los datos se descargue. 
  - Ejemplo: [download_all_rawData.sge](https://github.com/EveliaCoss/RNAseq_classFEB2024/blob/main/Practica_Dia2/scripts/download_all_rawData.sge)
B)	Tener varios script y agilizar la descarga múltiple. 
  - Ejemplos: [download_control_rawData.sge](https://github.com/EveliaCoss/RNAseq_classFEB2024/blob/main/Practica_Dia2/scripts/download_control_rawData.sge),  [download_15min_rawData.sge](https://github.com/EveliaCoss/RNAseq_classFEB2024/blob/main/Practica_Dia2/scripts/download_15min_rawData.sge), [download_30min_rawData.sge](https://github.com/EveliaCoss/RNAseq_classFEB2024/blob/main/Practica_Dia2/scripts/download_30min_rawData.sge) y [download_4h_rawData.sge](https://github.com/EveliaCoss/RNAseq_classFEB2024/blob/main/Practica_Dia2/scripts/download_4h_rawData.sge).

---

## Symlink

Los datos contenido en `data/` se encuentran presentes dentro de la carpeta con ayuda de un symlink:

```bash
ln -s /mnt/Citosina/amedina/dramirez/clase/data/fastq/* .
```

Es una manera de trabajar con los datos sin tener que hacer copias de los mismos. No pesan, pero funcionan como los datos originales.

---

---
# 1. Análisis de control de calidad de lecturas

Para hacer el análisis de control de calidad **QC**. Utilizaremos los programas `fastqc` y `multiqc`

### 1. Fastqc

Este programa va a realizar un análisis de control de calidad en cada una de los archivos `.fastq.gz` y nos va a dar un reporte en forma de un archivo tipo `.html`.

Para más información visita la pagina https://www.bioinformatics.babraham.ac.uk/projects/fastqc/.

### 2. Multiqc

Este programa toma todos los archivos `.html` que arrojó `fastqc` y nos dará un reporte combinado de todas las muestras.

Para más información visita la pagina https://multiqc.info/docs/getting_started/installation/

---

### Correr `fastqc`

En mi caso, debo de cargar el módulo de `fastqc` primero para poder utilizar los comandos. Si estás trabajando en tu computadora local y ya está instalado `fastqc`, tal vez esto no sea necesario.

También creamos una carpeta para los outputs de `fastqc`

```bash
qlogin # pedir un nodo de prueba, default = 2GB
module load fastqc/0.11.3
module load multiqc/1.5
cd /mnt/Guanina/bioinfo24/data/Clase_RNASeq2024/
mkdir quality1
```

Después, usaremos un `for loop` para hacer `fastqc` a todos los archivos que tengan terminación `.fastq.gz`

```bash
for file in data/*.fastq.gz; do fastqc $file -o quality1; done
```

El comando para correr `fastqc` en un solo archivo es:

```bash
fastqc nombre.fastq.gz -o /directorio/de/salida
```

---

## Analicemos el output de `fastQC`

Veamos la informacion contenida en el [`SRR12363092_1_fastqc.html`](https://eveliacoss.github.io/RNAseq_classFEB2024/Practica_Dia2/FastQC_Reports/SRR12363092_1_fastqc.html).

---
## Pausa: Mi carpeta se ve (más o menos) así:

```bash
Clase_RNASeq2024/
├── data/
│   ├── SRR12363092_1.fastq.gz
│   ├── SRR12363092_2.fastq.gz
│   ├── SRR12363093_1.fastq.gz
│   ├── SRR12363093_2.fastq.gz
...
├── quality1/
│   ├── SRR12363092_1_fastqc.html
│   ├── SRR12363092_1_fastqc.zip
│   ├── SRR12363092_2_fastqc.html
├── scripts/
...
```

Mira [aquí](https://github.com/EveliaCoss/RNASeq_Workshop_Nov2023/blob/main/Dia2_Alineamiento/out_logs/qc1.o369176) cómo se ve la ejecución de este comando.

---
### Correr `multiqc`

Multiqc reconoce los outputs de `fastqc` por lo que el comando para utilizarlo es muy sencillo

```bash
multiqc quality1
```

#### **NOTA: siempre es mejor utilizar direcciones absolutas a relativas, para evitar que tus outputs se guarden en un directorio no deseado: **

```bash
multiqc /mnt/Guanina/bioinfo24/data/Clase_RNASeq2024/quality1/
```

---

## Script y ejecución en el cluster

Si desean ver como se analizaron los datos empleando los programas `FastQC` y `multiQC` dentro del cluster DNA, les dejo los siguientes scripts:

- [`qc1.sh`](https://github.com/EveliaCoss/RNAseq_classFEB2024/edit/main/Practica_Dia2/scripts/qc1.sh) 
  - [`qc1.sge`](https://github.com/EveliaCoss/RNAseq_classFEB2024/blob/main/Practica_Dia2/scripts/qc1.sge)

Mira [aquí](https://github.com/EveliaCoss/RNAseq_classFEB2024/blob/main/Practica_Dia2/out_logs/qc1.o369176) la salida del programa.

---

## Analicemos el output de `multiqc`

Veamos la informacion contenida en el [`multiqc_report.html`](https://eveliacoss.github.io/RNAseq_classFEB2024/Practica_Dia2/FastQC_Reports/multiqc_report.html).

---

---

## Trimming

Para hacer **trimming** de las lecturas que no tuvieron una buena calidad, utilizaremos la herramienta `trimmomatic`. Este programa tiene muchas opciones que nos permiten hacer trimming de formas distintas, aquí muestro el comando que utilizaremos para nuestras necesidades. Pero asegúrate de leer el [manual](http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/TrimmomaticManual_V0.32.pdf) para tus análisis personales.

**Para paired-end necesitamos:**

- Las dos lecturas paired end por muestra: `SRRxxxxx_1.fastq.gz` y `SRRxxxx_2.fastq.gz`

- Un archivo con los adaptadores que vamos a cortar: `TruSeq3-PE-2.fa`

Descarga los adaptadores de [aquí](https://github.com/timflutre/trimmomatic/blob/master/adapters/TruSeq3-PE-2.fa)

```bash
wget https://raw.githubusercontent.com/timflutre/trimmomatic/master/adapters/TruSeq3-PE-2.fa
```

Trimmomatic nos dará 4 outputs: Las secuencias que quedaron sin par que eran originalmente del archivo "1": `_1_unpaired.fastq.gz`, las secuencias sin par que eran del archivo "2": `_2_unpaired.fastq.gz` y las secuencias que aun están pareadas: `_1_trimmed.fastq.gz` y `_2_trimmed.fastq.gz`.

---

## Correr Trimmomatic

Creamos una carpeta para los resultados

```bash
mkdir TRIM_results
```

Usamos un `for loop` para hacer trimmomatic a cada par de lecturas `SRRxxxxx_1.fastq.gz` y `SRRxxxx_2.fastq.gz`

```bash
module load trimmomatic/0.33
cd /mnt/Guanina/bioinfo24/data/Clase_RNASeq2024/data
for i in *_1.fastq.gz;
do echo
trimmomatic PE -threads 8 -phred33 $i "${i%_1.fastq.gz}_2.fastq.gz" \
/mnt/Guanina/bioinfo24/data/Clase_RNASeq2024/TRIM_results/"${i%_1.fastq.gz}_1_trimmed.fq.gz" \
/mnt/Guanina/bioinfo24/data/Clase_RNASeq2024/TRIM_results/"${i%_1.fastq.gz}_1_unpaired.fq.gz" \
/mnt/Guanina/bioinfo24/data/Clase_RNASeq2024/TRIM_results/"${i%_1.fastq.gz}_2_trimmed.fq.gz" \
/mnt/Guanina/bioinfo24/data/Clase_RNASeq2024/TRIM_results/"${i%_1.fastq.gz}_2_unpaired.fq.gz" \
ILLUMINACLIP:/mnt/Guanina/bioinfo24/data/Clase_RNASeq2024/TruSeq3-PE-2.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:20 MINLEN:80 \

done
```

---

# Documentación de las opciones

`ILLUMINACLIP:<fastaWithAdaptersEtc>:<seed mismatches>:<palindrome clip threshold>:<simple clip threshold>`

- seedMismatches: specifies the maximum mismatch count which will still allow a full match to be performed

- palindromeClipThreshold: specifies how accurate the match between the two 'adapter ligated' reads must be for PE palindrome read alignment.

- simpleClipThreshold: specifies how accurate the match between any adapter sequence must be against a read.

`SLIDINGWINDOW:<windowSize>:<requiredQuality>`

- windowSize: specifies the number of bases to average across

- requiredQuality: specifies the average quality required.

```bash
ILLUMINACLIP:/mnt/Guanina/bioinfo24/data/Clase_RNASeq2024/TruSeq3-PE-2.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:20 MINLEN:80
```

---

# Documentación de las opciones

`LEADING:<quality>`

Remove low quality bases from the beginning

quality: Specifies the minimum quality required to keep a base.

`TRAILING:<quality>`

Remove low quality bases from the end.

quality: Specifies the minimum quality required to keep a base

```bash
ILLUMINACLIP:/mnt/Guanina/bioinfo24/data/Clase_RNASeq2024/TruSeq3-PE-2.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:20 MINLEN:80
```

Mira [aquí](https://github.com/EveliaCoss/RNAseq_classFEB2024/blob/main/Practica_Dia2/out_logs/trim.o369193) cómo se ve la ejecución de este comando.

---

## QC del Trimming

### ¿Qué tan buena se encuentra la calidad posterior a la limpieza de adaptadores y secuencias de baja calidad?

Corramos `fastqc` y `multiqc` de nuevo

```bash
mkdir quality2
for file in /mnt/Guanina/bioinfo24/data/Clase_RNASeq2024/TRIM_results/*.fq.gz; do fastqc $file -o /mnt/Guanina/bioinfo24/data/Clase_RNASeq2024/quality2; done
```

```bash
multiqc /mnt/Guanina/bioinfo24/data/Clase_RNASeq2024/quality2
```

---

## Script y ejecución en el cluster

Si desean ver como se analizaron los datos empleando los programas `FastQC` y `multiQC` dentro del cluster DNA, les dejo los siguientes scripts:

- [`qc2.sh`](https://github.com/EveliaCoss/RNAseq_classFEB2024/blob/main/Practica_Dia2/scripts/qc2.sh) 
  - [`qc2.sge`](https://github.com/EveliaCoss/RNAseq_classFEB2024/blob/main/Practica_Dia2/scripts/qc1.sge)

---

### Pausa: Mi carpeta se ve (más o menos) así:

```bash
claseSTAR
├── data
│   ├── SRR12363092_1.fastq.gz
...
├── quality1
│   ├── multiqc_data
│   │   ├── multiqc_data.json
│   │   ...
│   ├── multiqc_report.html
...
├── quality2
│   ├── multiqc_data
│   │   ...
│   ├── SRR12363092_1_trimmed.fq_fastqc.html
│   ├── SRR12363092_1_unpaired.fq_fastqc.html
│   ├── SRR12363092_2_trimmed.fq_fastqc.html
...
├── TRIM_results
│   ├── SRR12363092_1_trimmed.fq.gz
│   ├── SRR12363092_1_unpaired.fq.gz
├── scripts/
```
---

### Analicemos el output del `FastQC` y `multiqc`

Veamos la información contenida en los reportes:

- `FastQC` - [SRR12363092_1_trimmed.fq_fastqc.html](https://eveliacoss.github.io/RNAseq_classFEB2024/Practica_Dia2/FastQC_Reports/SRR12363092_1_trimmed.fq_fastqc.html)

- `multiqc` - [`multiqc_report2.html`](https://eveliacoss.github.io/RNAseq_classFEB2024/Practica_Dia2/FastQC_Reports/multiqc_report2.html).

---

---
# Pero antes: ¿Qué es el alineamiento?

La **alineación del genoma** es un proceso bioinformático que consiste alinear las secuencias de ADN o ARN de uno o más genomas. El objetivo principal de la alineación del genoma es identificar *regiones de similitud u homología* entre las secuencias, lo que puede proporcionar información valiosa sobre diversos procesos biológicos, como la identificación de genes, el análisis evolutivo y la anotación funcional.

---

## Existen Diversas formas de alinear en RNA-seq

.left[.footnote[.black[
[Conesa, *et al*, 2016. *Genome Biology*](https://genomebiology.biomedcentral.com/articles/10.1186/s13059-016-0881-8)
]]]

---

## 1. Alineamiento y ensamblaje de lecturas guiado por el **genoma de referencia**

.pull-left[
<img src="data:image/png;base64,#figures/alignment_genoma.png" width="60%" style="display: block; margin: auto;" />
]

.pull-right[
El alineación del genoma de referencia implica mapear las lecturas de RNA-Seq en un **genoma de referencia** conocido.

- Nos permite identificar y cuantificar la expresión de **genes conocidos** y sus isoformas. Además, podemos anotar nuevos transcritos y genes.

- De preferencia contar con un archivo de anotación (GFF).

- La especie debe tener un genoma de buena calidad.

- Empleado normalmente en un organismo modelo (humano, ratón, *Arabidopsis*, etc.).
]

---

## 2. Ensamblaje de **transcriptoma guiado**

.pull-left[
<img src="data:image/png;base64,#figures/alignment_transcriptoma.png" width="60%" style="display: block; margin: auto;" />
]

.pull-right[
Las lecturas de RNA-Seq se asignan a un transcriptoma de referencia, que es una colección de transcritos.

- Veremos expresión de **genes**, pero no isoformas.

- NO hay anotación de nuevos transcritos.

- Si no está en el archivo de anotación (tz2gene/kallisto) no lo veremos.

- Es necesario un archivo de anotación con buena calidad.
]

---

## 2. Ensamblaje de **transcriptoma guiado**

### [**Kallisto**](https://pachterlab.github.io/kallisto/manual)

- **Pseudoalineamiento**.

- Es rápido.

- Se puede ejecutar el programa desde tu computadora.

- Se basa en los grafos de Brujin Graph (T-DBG) .

- Los Nodos (v1,v2,v3) son *k-mers*.

- Omite pasos redundantes en el T-DBG.

]

.pull-right[
<img src="data:image/png;base64,#figures/alignment_kallisto.png" width="80%" style="display: block; margin: auto;" />
]

Si quieren intentarlo les dejo el curso que di en 2023 - [RNAseq_classFEB2023](https://github.com/EveliaCoss/RNAseq_classFEB2023/tree/main/RNA_seq).

---

## 3. Ensamblaje ***de novo***

.pull-left[
<img src="data:image/png;base64,#figures/alignment_denovo.png" width="50%" style="display: block; margin: auto;" />
]

.pull-right[
Ideal para una especie con **genoma de mala calidad o sin referencia**, como **organismos NO modelos**, además de si no contamos con un archivo de anotación bueno.

- Es recomendado utilizar lecturas *Paired-end*.
]

---

# En esta clase utilizaremos el software de alineamiento **STAR**, de tipo alineamiento a **genoma de referencia**

#### Seguiremos 2 sencillos pasos

1. Indexar el genoma de referencia creando un índice de STAR

2. Alinear y contar con STAR

[Manual de STAR](https://github.com/alexdobin/STAR/blob/master/doc/STARmanual.pdf)

**STAR** nos permite, además de alinear las lecturas, hacer un conteo en el mismo paso. El paso de conteo, puede ser separado y con otras herramientas, pero en esta clase, te enseñaré la forma sencilla en la que STAR también puede contar.

> **NOTA: Para obtener la matriz de cuentas, NECESITAMOS UN ARCHIVO DE ANOTACIÓN.**

---

## ¿Qué es un archivo de anotación?

- Un **archivo de anotación GFF** (General Feature File) es un formato de archivo estándar utilizado en bioinformática para almacenar y representar información genómica y de anotación para diversas características dentro de un *genoma*, como **genes, transcritos, exones y otros elementos genómicos**.

**Utilizamos estos archivos junto al genoma de referencia para poder encontrar ubicaciones de genes, transcritos, etc.**

- Los archivos GFF pueden venir en diferentes versiones, siendo GFF3 uno de los formatos más utilizados.

Podemos encontrar estos archivos en bases de datos genómicas como GENCODE para humanos y ratones [aquí](www.gencodegenes.org)

---

## ¿Qué es indexar un genoma de referencia?

Es una forma computacional de crear una "estructura de datos" para el genoma de referencia, mediante **índices**, de tal forma que podramos **accesar a las partes del mismo** de una forma más eficiente al alinear.

El genoma de referencia sirve como plantilla contra la cual se realizan diversos análisis genómicos, como mapeo de lecturas, llamado de variantes y cuantificación de la expresión génica.

**La indexación mejora la velocidad y la eficiencia de estos análisis** al permitir que el software **ubique y acceda rápidamente** a partes relevantes del genoma.

---

## 1. Crear un índice de STAR

Carguemos el módulo de STAR

```bash
module load star/2.7.9a
```

Creamos un directorio para guardar el indice

```bash
mkdir STAR_index
```

El script para crear el indice es el siguiente:

```bash
STAR --runThreadN 12 \
--runMode genomeGenerate \
--genomeDir /mnt/Guanina/bioinfo24/data/Clase_RNASeq2024/STAR_index \
--genomeFastaFiles /mnt/Archives/genome/mouse/mm10/UCSC/chromosomes-1-file/mm10.fa \
--sjdbGTFfile /mnt/Guanina/bioinfo24/data/Clase_RNASeq2024/annotation/mm10.knownGene.gtf.gz \
--sjdbOverhang 149
```

Entremos al 
[Manual de STAR](https://github.com/alexdobin/STAR/blob/master/doc/STARmanual.pdf) para entender las opciones

---

## Script y ejecución en el cluster

Si desean ver como se analizaron los datos empleando el programa `STAR` para realizar el *index del genoma* dentro del cluster DNA, les dejo los siguientes scripts:

- [`STAR_index.sh`](https://github.com/EveliaCoss/RNAseq_classFEB2024/blob/main/Practica_Dia2/scripts/STAR_index.sh) 
  - [`index.sge`](https://github.com/EveliaCoss/RNAseq_classFEB2024/blob/main/Practica_Dia2/scripts/index.sge)

Mira [aquí](https://github.com/EveliaCoss/RNAseq_classFEB2024/blob/main/Practica_Dia2/out_logs/index.o369188) cómo se ve la ejecución de este comando.

---

## 2. Alinear y **CONTAR** con STAR

Si revisamos el [Manual de STAR](https://github.com/alexdobin/STAR/blob/master/doc/STARmanual.pdf) notarás que hay una opción para **Contar lecturas por genes** (Sección 8). Con la opción de `--quantMode`. De hecho, estas cuentas coinciden con las cuentas que nos daría `htseq-count`

```bash
index=/mnt/Guanina/bioinfo24/data/Clase_RNASeq2024/STAR_index
FILES=/mnt/Guanina/bioinfo24/data/Clase_RNASeq2024/TRIM_results/*_1_trimmed.fq.gz
for f in $FILES
do
    base=$(basename $f _1_trimmed.fq.gz)
    echo $base
    STAR --runThreadN 12 --genomeDir $index --readFilesIn $f /mnt/Guanina/bioinfo24/data/Clase_RNASeq2024/TRIM_results/$base"_2_trimmed.fq.gz" \
    --outSAMtype BAM SortedByCoordinate \
    --quantMode GeneCounts \
    --readFilesCommand zcat \
    --outFileNamePrefix /mnt/Guanina/bioinfo24/data/Clase_RNASeq2024/STAR_output/$base
done
```

Mira [aquí](https://github.com/EveliaCoss/RNAseq_classFEB2024/blob/main/Practica_Dia2/out_logs/align.o369662) cómo se ve la ejecución de este comando.

---

### SCRIPTS empleados en la clase

Todos los scripts usados en esta clase están en el [GitHub](https://github.com/EveliaCoss/RNAseq_classFEB2024/tree/main/Practica_Dia2/scripts)

---

### Las carpetas contenidas por equipo deberan ser:

```bash
Equipo1
├── annotation    # Archivo de anotacion del organismo (GFF)
├── data          # raw Data
├── genome        # Genoma de Referencia o symlink del genoma de referencia (opcional)
├── metadata.csv  # Metadata
├── quality1      # FastQC y multiQC de raw Data
├── quality2      # FastQC y multiQC de los datos despues del Trimming
├── results       # Resultados obtenidos de DEG
├── scripts       # Todos los scripts
├── STAR_index    # Index del genoma de referencia
├── STAR_output   # Salida de STAR, cuentas y BAM
├── TRIM_results  # Salida del Trimming
└── TruSeq3-PE-2.fa #Adaptadores PE de Illumina

```

---

Deben entregarme un documento con la **Descripción de los datos** por equipo al correo ecoss@liigh.unam.mx.

---

Gracias por tu atención, respira y coméntame tus dudas.