Workshop: Análisis de datos de RNA-Seq

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

####Dra. Evelia Coss, Sofia Salazar y Diego Ramirez 
#### 16/04/2026
]

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

---

# Repaso

### Contraste

<img src="data:image/png;base64,#figures/diseno_experimental.png" alt="" width="80%" style="display: block; margin: auto;" />
]

### Diseño experimental

<img src="data:image/png;base64,#figures/Experimental_design.jpg" alt="" width="50%" style="display: block; margin: auto;" />
]

.left[.footnote[.black[
Figura obtenida del [NBIS SciLifeLab](https://scilifelab.github.io/courses/ngsintro/1910/slides/rnaseq/presentation.html#6)
]]]

---

## Pipeline bioinformática para RNA-seq

.pull-right[
<img src="data:image/png;base64,#figures/pipeline1.png" alt="" 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)
]]]

---

# Contenido de la clase

- 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 [GSE250023](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE250023)** provenientes de [Kavita, *et al*. 2024. *JCIinsight*](https://insight.jci.org/articles/view/176556).

- Bioproject: PRJNA1051620, [NIH](https://www.ncbi.nlm.nih.gov/bioproject/PRJNA1051620) y [EBI](https://www.ebi.ac.uk/ena/browser/view/PRJNA1051620)

---

### Descarga de los datos con `wget`

.pull-left[
1. Selecciona el estudio que vas a utilizar.
2. Ingresa a la página de [EBI](https://www.ebi.ac.uk/).
3. En el buscador, escribe el ID del estudio y presiona Enter.
4. En los resultados, ubica la sección **"Nucleotide sequences"** y selecciona el estudio correspondiente.
5. En la tabla inferior, dirígete a la columna **"Generated FASTQ files: FTP"**.
6. Selecciona las muestras que deseas descargar.
7. Haz clic izquierdo en **"Get download script"** para obtener el script de descarga.
]

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

---

### Bioproject, Biosample y SRA

En pocas palabras, **SRA contiene los datos crudos de secuenciación** obtenidos de **muestras biológicas** (BioSamples) dentro de un **estudio** (BioProject).

.left[.footnote[.black[
Imagen proveniente de [NCBI](https://www.ncbi.nlm.nih.gov/sra/docs/submitbio/)
]]]

---

### Bioproject: PRJNA1051620

En este proyecto se analizaron **18 pacientes con lupus**, contando con 4 tomas en el tiempo posterior a la **vacunación contra COVID-19 (BNT162b2 SARS-CoV-2)**, lo que resultó en un total de 75 biosamples (SRA Experiments).

Para fines de esta práctica, se seleccionaron aquellas correspondientes al **día 7** posterior a la vacunación contra COVID-19 (BNT162b2 SARS-CoV-2) en pacientes con lupus.

<div class="datatables html-widget html-fill-item" id="htmlwidget-926ae72ea4c13870a856" style="width:100%;height:auto;"></div>
<script type="application/json" data-for="htmlwidget-926ae72ea4c13870a856">{"x":{"filter":"none","vertical":false,"caption":"<caption>Información de la metadata<\/caption>","fillContainer":false,"data":[["1","2","3","4","5","6"],[1,8,14,12,17,18],["S003","S008","S018","S013","S023","S028"],["SLE_R","SLE_R","SLE_R","SLE_NR","SLE_NR","SLE_NR"],["R","R","R","NR","NR","NR"],[7,7,7,7,7,7],["d7","d7","d7","d7","d7","d7"],["p1.S003.day7","p8.S008.day7","p14.S018.day7","p12.S013.day7","p17.S023.day7","p18.S028.day7"],["SAMN38793001","SAMN38792997","SAMN38792988","SAMN38792993","SAMN38792983","SAMN38792978"]],"container":"<table class=\"display\">\n  <thead>\n    <tr>\n      <th> <\/th>\n      <th>individual<\/th>\n      <th>sample<\/th>\n      <th>condition<\/th>\n      <th>reponder<\/th>\n      <th>timepoint<\/th>\n      <th>timelabel<\/th>\n      <th>identifiers<\/th>\n      <th>BioSample<\/th>\n    <\/tr>\n  <\/thead>\n<\/table>","options":{"pageLength":8,"columnDefs":[{"className":"dt-right","targets":[1,5]},{"orderable":false,"targets":0},{"name":" ","targets":0},{"name":"individual","targets":1},{"name":"sample","targets":2},{"name":"condition","targets":3},{"name":"reponder","targets":4},{"name":"timepoint","targets":5},{"name":"timelabel","targets":6},{"name":"identifiers","targets":7},{"name":"BioSample","targets":8}],"order":[],"autoWidth":false,"orderClasses":false,"lengthMenu":[8,10,25,50,100]}},"evals":[],"jsHooks":[]}</script>

---

### Descarga de los datos con `wget`

Esto descargará un script BASH `.sh` como [este](https://github.com/EveliaCoss/RNAseq_classFEB2025/blob/main/Practica_Dia2/scripts/PRJNA1051620_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 [`PRJNA1051620_download.sh`](https://github.com/EveliaCoss/RNAseq_classFEB2026/blob/main/Practica_Dia2/scripts/PRJNA1051620_download.sh). Y lo tengo en una carpeta llamada `scripts`.

``` bash
cd /mnt/data/bioinfo-estadistica-2/RNAseq_2026/BioProject_2026/scripts/
```

> NOTA: Recuerda darle permisos a tus compañer@s:

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

---

### Descarga de los datos con `wget`

Ya que las descargas tardan mucho, es conveniente enviar el script a la cola de ejecución usando **Slurm (Simple Linux Utility for Resource Management)** con un archivo `.sh`. De esta manera, el script seguirá ejecutándose aunque tú no estés.

Es necesario tener un solo script donde todos los datos se descargue.

- Ejemplo: [download_all_rawData.sh](https://github.com/EveliaCoss/RNAseq_classFEB2026/blob/main/Practica_Dia2/scripts/download_all_rawData.sh)

---

## NO se trabaja en el **nodo principal**

📢 **Regla importante:** ¡No trabajes directamente en el nodo principal del clúster!

🚫 **No ejecutes scripts pesados ni análisis directamente en el nodo de acceso.**

✔️ **Siempre usa el sistema de colas (slurm)** para enviar trabajos o solicita un nodo de prueba si necesitas ejecutar algo interactivo. Usa `qsub` para enviar trabajos.

Si necesitas hacer pruebas o ejecutar código interactivamente, solicita un nodo de prueba con `srun -p interactive --pty bash`.

``` bash
# pedir un nodo de prueba
# Nos asigna 1 tarea, 1 CPU y memoria de 2 GB
# Ver limite de la particion
sinfo -p interactive
# ver cuanto nos asignaron
scontrol show job $SLURM_JOB_ID
```

Para salir del nodo interactivo, usa:

``` bash
exit
```

---

## Ejemplo de uso interactivo

Si quieres entrar a una sesión interactiva con más control, puedes pedir explícitamente recursos:

``` bash
srun -p interactive --ntasks=1 --cpus-per-task=1 --mem=2G --time=01:00:00 --pty bash
```

Esto nos da: 1 tarea, 1 CPU, memoria 2GB de RAM, 1 hora de tiempo máximo.

---

## Analogía: “Trabajadores en una mesa”

- **Tarea** (`--ntasks`) → es como un trabajador independiente. Cada tarea es un proceso que Slurm lanza.

- **CPU por tarea** (`--cpus-per-task`) → son las manos que tiene cada trabajador. Más manos = puede hacer varias cosas en paralelo.

- **Memoria** (`--mem`) → es el tamaño de la mesa donde ese trabajador pone sus materiales. Si la mesa es pequeña, se le caen los libros; si es grande, puede trabajar cómodo.

Ejemplo:

``` bash
#SBATCH --ntasks=1            # 1 trabajador
#SBATCH --cpus-per-task=1.    # 1 mano
#SBATCH --mem=2G              # en una mesa pequeña
```

---

---

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

Vamos a analizar el contenido del script `qc1.sh`:

``` bash
# Modulos requeridos:
module load fastqc/0.11.3

# FastQC analysis
for file in /mnt/data/bioinfo-estadistica-2/RNAseq_2026/BioProject_2026/data/*.fastq.gz; do fastqc $file -o /mnt/data/bioinfo-estadistica-2/RNAseq_2026/BioProject_2026/quality1; done
```

Para analizar  **todos los archivos** que tengan terminación `.fastq.gz` podemos emplear un `for loop`:

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

---

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

``` bash
/mnt/data/bioinfo-estadistica-2/RNAseq_2026/BioProject_2026/
├── 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/
...
```

---

### Correr `multiqc`

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

``` bash
module load anaconda3/2025.06 # Cargar el módulo de Anaconda
conda activate multiqc-1.5 # Activar el ambiene dentro de Conda
cd /mnt/data/bioinfo-estadistica-2/RNAseq_2026/BioProject_2026/data/raw/
multiqc .
```

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

``` bash
multiqc /mnt/data/bioinfo-estadistica-2/RNAseq_2026/BioProject_2026/data/raw/
```

---

## 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 el siguiente script:

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

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

---

## Analicemos el output del `FastQC` y `multiqc`
### Primer control de calidad

- Veamos la información contenida en el [`SRR12363092_1_fastqc.html`](https://eveliacoss.github.io/RNAseq_classFEB2025/Practica_Dia2/FastQC_Reports/SRR12363092_1_fastqc.html).

- Veamos la informacion contenida en el [`multiqc_report.html`](https://eveliacoss.github.io/RNAseq_classFEB2026/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
cd /mnt/data/bioinfo-estadistica-2/RNAseq_2026/BioProject_2026/adapters
# Paired-end
wget https://raw.githubusercontent.com/timflutre/trimmomatic/master/adapters/TruSeq3-PE-2.fa
# single-end
wget https://raw.githubusercontent.com/timflutre/trimmomatic/master/adapters/TruSeq3-SE.fa
```

---

## Correr Trimmomatic

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/data/bioinfo-estadistica-2/RNAseq_2026/BioProject_2026/data/raw

for i in *_1.fastq.gz;
do echo
trimmomatic PE -threads 8 -phred33 $i "${i%_1.fastq.gz}_2.fastq.gz" \
/mnt/data/bioinfo-estadistica-2/RNAseq_2026/BioProject_2026/data/processed/"${i%_1.fastq.gz}_1_trimmed.fq.gz" \
/mnt/data/bioinfo-estadistica-2/RNAseq_2026/BioProject_2026/data/processed/"${i%_1.fastq.gz}_1_unpaired.fq.gz" \
/mnt/data/bioinfo-estadistica-2/RNAseq_2026/BioProject_2026/data/processed/"${i%_1.fastq.gz}_2_trimmed.fq.gz" \
/mnt/data/bioinfo-estadistica-2/RNAseq_2026/BioProject_2026/data/processed/"${i%_1.fastq.gz}_2_unpaired.fq.gz" \
ILLUMINACLIP:/mnt/data/bioinfo-estadistica-2/RNAseq_2026/BioProject_2026/adapters/TruSeq3-PE-2.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:20 MINLEN:80 \
done
```

---
### Trimmomatic

## Script y ejecución en el cluster

- [`trimm_adapters.sh`](https://github.com/EveliaCoss/RNAseq_classFEB2026/blob/main/Practica_Dia2/scripts/trimm_adapters.sh)

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

## Salida

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

---

### Trimmomatic
## Documentación de las opciones

``` bash
ILLUMINACLIP:/mnt/data/bioinfo-estadistica-2/RNAseq_2026/BioProject_2026/adapters/TruSeq3-PE-2.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:20 MINLEN:80
```

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

- Recorta adaptadores de las lecturas usando el archivo de secuencias de adaptadores especificado.
- `fastaWithAdaptersEtc`: Archivo de adaptadores
- `seed mismatches`: número máximo de desajustes permitidos entre la secuencia de adaptador y la lectura.
- `palindrome clip threshold`: umbral para considerar una coincidencia con el adaptador.
- `simple clip threshold`: tamaño de la ventana para evaluar coincidencias con adaptadores.

LEADING:`<quality>`

- Recorta bases al inicio de la lectura si tienen una calidad menor al valor especificado.
- `quality`: Calidad minima requerida para retener la base.

---

### Trimmomatic
## Documentación de las opciones

``` bash
ILLUMINACLIP:/mnt/data/bioinfo-estadistica-2/RNAseq_2026/BioProject_2026/adapters/TruSeq3-PE-2.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:20 MINLEN:80
```

TRAILING:`<quality>`

- Recorta bases al final de la lectura si tienen una calidad menor al valor especificado.
- `quality`: Calidad minima requerida para retener la base.

SLIDINGWINDOW:`<windowSize>:<requiredQuality>`

- Escanea la lectura con una ventana deslizante de <tamaño_ventana> bases.
- Si el promedio de calidad dentro de la ventana es menor que <calidad>, recorta la lectura desde ese punto.
- `windowSize`: especifica el número de bases que hay que promediar
- `requiredQuality`: especifica la calidad media exigida.

---

### Trimmomatic
## Documentación de las opciones

``` bash
ILLUMINACLIP:/mnt/data/bioinfo-estadistica-2/RNAseq_2026/BioProject_2026/adapters/TruSeq3-PE-2.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:20 MINLEN:80
```

MINLEN: `MINLEN:<longitud>`
- Descarta lecturas con una longitud menor al valor especificado después del recorte.

---

## QC del Trimming
### ¿Qué tan buena es la calidad después de eliminar adaptadores y secuencias de baja calidad?

Corramos `fastqc` y `multiqc` de nuevo

``` bash
mkdir quality2
for file in /mnt/data/bioinfo-estadistica-2/RNAseq_2026/BioProject_2026/data/processed/*.fq.gz; do fastqc $file -o /mnt/data/bioinfo-estadistica-2/RNAseq_2026/BioProject_2026/quality2; done
```

``` bash
multiqc /mnt/atgc-d1/bioinfoII/rnaseq/BioProject_2026/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_classFEB2026/blob/main/Practica_Dia2/scripts/qc2.sh) 
  
- Mira [aquí](https://github.com/EveliaCoss/RNAseq_classFEB2026/blob/main/Practica_Dia2/out_logs/qc.o369176) cómo se ve la ejecución de este comando.

---

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

``` bash
/mnt/data/bioinfo-estadistica-2/RNAseq_2026/BioProject_2026/
├── data/raw
│   ├── SRR12363092_1.fastq.gz
...
├── data/processed
│   ├── SRR12363092_1_trimmed.fq.gz
│   ├── SRR12363092_1_unpaired.fq.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
...
├── scripts/
```

---

## Analicemos el output del `FastQC` y `multiqc`
### Segundo control de calidad

Veamos la información contenida en los reportes:

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

- `multiqc` - [`multiqc_report2.html`](https://eveliacoss.github.io/RNAseq_classFEB2025/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

.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" alt="" 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" alt="" 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 (tx2gene/kallisto) no lo veremos.

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

---

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

.pull-left[
<img src="data:image/png;base64,#figures/alignment_denovo.png" alt="" 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*.
]

---

**Subir en la tarea de Google Classroom**

---

### Tarea
## Cada equipo debe entregarme un documento en formato Rmd y HTML que contenga la siguiente información:

- Descripción de los datos (entregable 1, fecha límite: 16 de marzo).
- Explicación de las gráficas generadas por MultiQC.
- Conclusión sobre los datos: 
    + ¿Son viables para continuar con el análisis? 
    + ¿Qué pasos deben seguirse para mejorar la calidad de los datos?

Ejemplo de entregable: [Reporte](https://eveliacoss.github.io/RNAseq_classFEB2026/Entregable_ejemplo2.html)

**Subir en la tarea de Google Classroom**

---

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