# **Unidad 4. Análisis de expresion diferencial con DESeq2**
#### Dra. Laura Liliana López Galindo (laura.lopez11@uabc.edu.mx)
---
## Paso 1: flujo de trabajo en [OMICA](http://omica.cicese.mx/docs/index.html).
Acceder al servidor con tu cuenta de omica.
Utilizar `-YC` cuando se requiera despliegue gráfico.
```javascript=
ssh lilopez@omica
```
#### 1. Mapeo de lecturas vs el transcriptoma *de novo*
En la carpeta de `Transcriptomica`, crear la carpeta `DIFF_EXP` y acceder.
```javascript=
cd Cg_micropest/Transcriptomica
mkdir DIFF_EXP
cd DIFF_EXP
```
Crear un `link simbólico` del ensamble de buena calidad obtenido con Transrate `good.Trinity_ensamble.Trinity.fasta`.
```javascript=
ln -s ../TRANSRATE/gcontigs/Trinity_ensamble.Trinity/good.Trinity_ensamble.Trinity.fasta ./
ls
```
Crear un `link simbólico` de las lecturas pareadas y limpias obtenidas con trimmomatic.
```javascript=
ln -s ../TRIMMOMATIC/TRIM_P/*.fastq ./
ls
```
Copiar el archivo `samples_file.txt` que se encuentra en la carpeta de `ENSAMBLE`
```javascript=
cp ../ENSAMBLE/De_novo/samples_file.txt ./
ls
```
:::success
Estructura del archivo **samples_file.txt** que debe contener la informacion de cada muestra.
:::
Generar el script al cual llamaremos `mapping_rsem.slrm` para someter al manejador de tareas del Cluster omica. Esto nos permitira generar el mapeo de nuestras lecturas contra el emsamble *de novo* generado previamente.
```javascript=
nano mapping_rsem.slrm
```
El slrm debe contener lo siguiente:
```javascript=
#!/bin/bash
#SBATCH -p cicese
#SBATCH --job-name=rsemCT
#SBATCH --output=rsem-%j.log
#SBATCH --error=rsem-%j.err
#SBATCH -N 2
#SBATCH --mem=100GB
#SBATCH --ntasks-per-node=24
#SBATCH -t 6-00:00:00
#SBATCH --exclusive
#export some soft
export PATH=$PATH:/LUSTRE/apps/bioinformatica/bowtie2/bin/
export PATH=$PATH:/LUSTRE/apps/bioinformatica/samtools-1.17/bin/
export PATH=$PATH:/LUSTRE/apps/bioinformatica/jellyfish-2.3.0/
export PATH=$PATH:/LUSTRE/apps/bioinformatica/salmon/bin/
export PATH=$PATH:/LUSTRE/apps/bioinformatica/RSEM/bin
module load conda-2023
module load trinityrnaseq-v2.15.1
UTIL=/LUSTRE/apps/bioinformatica/trinityrnaseq-v2.15.1/util
cd $SLURM_SUBMIT_DIR
$UTIL/align_and_estimate_abundance.pl \
--seqType fq \
--samples_file samples_file.txt \
--transcripts good.Trinity_ensamble.Trinity.fasta \
--est_method RSEM \
--aln_method bowtie2 \
--trinity_mode \
--prep_reference \
--output_dir RSEM_outdir \
--thread_count=24 \
exit 0
```
Someter script a manejador de tareas de omica
```javascript=
sbatch mapping_rsem.slrm
squeue
```
Este script generará varios archivos y las carpetas con la leyenda de cada condición donde se encuentran los conteos de los transcritos expresados en TPM (Transcripts per million) y FPKM (Fragments per thousand nucleotides per million mapped reads).
Se deberan cambiar los permisos de los directorios generados con RSEM de la siguiente manera:
```javascript=
chmod 777 S0* S8*
```
Las carpetas deberan aparecer de modo subrayado y al revisarlas deberan verse con los siguientes permisos:
```javascript=
ls -l
```
:::info
drwxrwxrwx 2 lilopez clarissa 4096 dic 9 07:58 S0_R1
:::
se renombrara el archivo RSEM.isoforms.results en cada carpeta con el nombre del tratamiento corrwspondiente (es decir el mismo nombre que aparece en la carpeta). Para ello se utilizara el siguiente comando:
```javascript=
cd S0_R1
mv RSEM.isoforms.results S0_R1.isoforms.results
cd ..
```
Este paso se dedera realizar con cada carpeta.
Posteriormente se deberan comprimir los archivos a los cuales les cambiamos el nombre utilizando el siguiente comando:
```javascript=
gzip S*/*.isoforms.results
```
---
### Paso 2: Flujo de trabajo en la terminal de nuestra PC.
:::warning
**NOTA:** Para poder realizar esta seccion de flujo de trabajo deberemos tener instalado [Trinity](https://github.com/trinityrnaseq/trinityrnaseq/wiki/Installing-Trinity) en nuestra PC y **python 2.7**. [Python 2.7](https://netslovers.com/post/how-install-python-27-pip-ubuntu-22-04/) debe ser predeterminado como default.
:::
#### 1. Matrices de abundancia
Abriremos una nueva terminal en nuestra PC y se crearan una serie de directorios anidados denominados `Transcriptomica`, `DIFF_EXP` y `Matrix`. Utilizar los siguientes comandos:
```javascript=
mkdir Transcriptomica
cd Transcriptomica
mkdir DIFF_EXP
cd DIFF_EXP
mkdir Matrix
```
Se procedera a bajar cada archivo `isoforms.results.gz` de nuestros directorios del servidor OMICA a nuestra PC en el directorio que creamos denominado `Matrix`. Para ello se usara el siguiente comando:
```javascript=
scp -r lilopez@omica:/LUSTRE/bioinformatica_data/genomica_funcional/Laura/proyectos/Cg_Micropest/Transcriptomica/DIFF_EXP/S*/*.gz ./
```
Descomprimir los archivos `isoforms.results.gz` con el siguiente comando:
```javascript=
gunzip *.gz
```
Descargar el archivo `good.Trinity_ensamble.Trinity.fasta.gene_trans_map` del servidor omica a nuestro directorio `Matrix`. Utilizar el siguiente comando:
```javascript=
scp lilopez@omica:/LUSTRE/bioinformatica_data/genomica_funcional/Laura/proyectos/Cg_Micropest/Transcriptomica/DIFF_EXP/good.Trinity_ensamble.Trinity.fasta.gene_trans_map ./
```
Con los archivos `isoforms.results` y nuestro archivo `good.Trinity_ensamble.Trinity.fasta.gene_trans_map` se realizara la matriz de conteos. El script a utilizar es el siguiente:
```javascript=
abundance_estimates_to_matrix.pl --est_method RSEM S0_R1.isoforms.results S0_R2.isoforms.results S0_R3.isoforms.results S8_C_T1.isoforms.results S8_C_T5.isoforms.results S8_C_T12.isoforms.results S8_MIX_T4.isoforms.results S8_MIX_T6.isoforms.results S8_MIX_T8.isoforms.results S8_MP_T2.isoforms.results S8_MP_T3.isoforms.results S8_MP_T10.isoforms.results S8_P_T7.isoforms.results S8_P_T9.isoforms.results S8_P_T11.isoforms.results --gene_trans_map ./good.Trinity_ensamble.Trinity.fasta.gene_trans_map
```
Este procedimiento generara las matrices de conteos (abundancia) para genes e isoformas, asi como los archivos TMM y TPM correspondientes. Se generaran un total de 10 archivos que iniciaran con la palabra RSEM.
:::success
Archivos generados:
:::
Posteriormente, generar un archivo de texto, el cual debera contener en la primer columna el tratamiento y en la segunda las replicas de cada tratamiento, y al cual denominaremos `samples_described.txt` y se generara con la siguiente estructura tabular.
```javascript=
S0 S0_R1
S0 S0_R2
S0 S0_R3
S8_C S8_C_T1
S8_C S8_C_T5
S8_C S8_C_T12
S8_MIX S8_MIX_T4
S8_MIX S8_MIX_T6
S8_MIX S8_MIX_T8
S8_MP S8_MP_T2
S8_MP S8_MP_T3
S8_MP S8_MP_T10
S8_P S8_P_T7
S8_P S8_P_T9
S8_P S8_P_T11
```
Verificar la estructura tabular del archivo `samples_described.txt`. La tabulación se muestra con el simbolo `^I` y el final de la línea con simbolo `$`.
```javascript=
cat -TE samples_described.txt
```
:::info
**RECOMENDACION/OPCIONAL:**
En este punto, es decir, antes de empezar un análisis de expresión diferencial per se, se recomienda realizar un análisis de control de calidad asociado con las muestras, las réplicas y su agrupación. Además, con base en las matrices de abundancia, podemos calcular los valores
de ExN50 como estrategia adicional para evaluar la calidad del ensamble. Para ello, creamos los siguientes directorios:
```javascript=
mkdir –p SAMPLES_REPLICATES_QC/{SAMPLES,PCA_CORRELATION,ExN50}
```
El primer paso es verificar el # de fragmentos por réplica biológica para cada una de las condiciones. Para ello, nos movemos a la carpeta SAMPLES y ejecutamos las instrucciones.
```javascript=
cd SAMPLES_REPLICATES_QC/ SAMPLES
PtR –m ../../matrix.counts.matrix –s ../../samples described.txt --log2 --compare_replicates &
```
Explorar los archivos con extensión `.pdf`
Ahora procederemos a calcular los valores de correlación entre las muestras y como están agrupadas vía un análisis de componentes principales. Nos movemos a la carpeta PCA_CORRELATION.
```javascript=
cd ../PCA_CORRELATION
```
Ejecutamos las siguientes instrucciones:
```javascript=
PtR –m ../../matrix.TMM.EXPR.matrix –s ../../samples_described.txt --log2 --sample_cor_matrix &
PtR –m ../../matrix.TMM.EXPR.matrix –s ../../samples_described.txt --log2 --CPM --prin_comp 3 &
```
Revisar los archivos en formato pdf. Al revisar los archivos, hay que recordar que este set de datos es una muestra muy pequeña del total de lecturas y que sólo representa un ejemplo ilustrativo.
Finalmente, procederemos a calcular los valores de ExN50 y obtener su gráfica correspondiente. Para esto, nos movemos a dicho directorio y ejecutamos las directivas.
```javascript=
cd ../ExN50
contig_ExN50_statistic.pl ../../matrix.TMM.EXPR.matrix ../../good.Trinity.fasta > ExN50.stats
plot_ExN50_statistic.Rscript ExN50.stats
```
En la Figura obtenida se muestran los valores de ExN50 e indican los valores del estadístico N50 de los contigs pero restringido a los transcritos altamente expresados.
:::
#### 2. Analisis de expresion diferencial de genes
Realizar el análisis de expresión per se. Para esto, crearemos dentro de `Matrix` un directorio denominado `DE` y accederemos al mismo.
```javascript=
mkdir DE
cd DE
```
Crear un link simbolico al archivo `RSEM.isoform.counts.matrix` y copiar el archivo `samples_described.txt`.
```javascript=
ln -s ../RSEM.isoform.counts.matrix ./
cp ../samples_described.txt ./
```
Posteriormente, se realizara el analisis de expresion diferencial ejecutando la siguiente linea de comando:
```javascript=
run_DE_analysis.pl --matrix RSEM.isoform.counts.matrix --method DESeq2 --samples_file samples_described.txt
```
Esta linea de comando generará un directorio con el nombre `DESeq2.######.dir` que contiene la lista de los genes expresados diferencialmente, los volcano plots y las matrices asociadas a cada una de las comparaciones pareadas.
:::success
Archivos a obtener por cada comparacion pareada.
:::
#### 3. Generacion de heatmap de genes expresados diferencialmente
Para extraer los transcritos expresados diferencialmente con valores de P y magnitud de cambio (fold change) específicos, nos movemos al directorio
`DESeq2.######.dir`.
```javascript=
cd DESeq2.######.dir
```
Posteriormente ejecutaremos el siguiente script:
```javascript=
analyze_diff_expr.pl --matrix ../../RSEM.isoform.TMM.EXPR.matrix --samples ../samples_described.txt -P 0.01 -C 1
```
Con C, representando log2, estamos seleccionando aquellos transcritos sobre o subexpresados con fold change mínimo de 2. El valor *P* 0.01 representa el nivel de significancia para la tasa de falso descubrimiento (FDR), introducida por Benjamini y Hochberg (1995). El valor de P y C puede ser modificado en relacion al numero de genes obtenidos, tomando en cuenta que un valor de *P* menor es mas estricto (*P* 0.05, *P* 0.01, *P* 0.001). El script generara el heatmap correspondiente en formato `.pdf` y la matriz de genes expresados diferencialmente `.txt`.
:::success
Archivos generados:
:::
De las comparaciones pareadas podemos extraer la lista de los transcritos UP y DOWN para ingresarlos en la página de [Venny 2.0](https://bioinfogp.cnb.csic.es/tools/venny/index2.0.2.html) y obtener los diagramas de VENN o tambien se puede ejecutar dentro de la paqueteria de R.
Por ejemplo, si queremos extraer los transcritos UP y DOWN del tratamiento S8_MIX con respecto al control que es S8_C, ejecutamos lo siguiente:
```javascript=
mkdir UP-DOWN_lists
cat RSEM.isoform.counts.matrix.S8_C_vs_S8_MIX.DESeq2.DE_results.P0.01_C1.S8_MIX-UP.subset | cut –f 1 >
./UP-DOWN_lists/S8_MIX_UP.txt
cat RSEM.isoform.counts.matrix.S8_C_vs_S8_MIX.DESeq2.DE_results.P0.01_C1.S8_C-UP.subset | cut –f 1 >
./UP-DOWN_lists/S8_MIX_DOWN.txt
```
Eliminar la primer fila de cada archivo (SAMPLE A) para dejar unicamente los ID de los transcritos.
```javascript=
nano S8_MIX_UP.txt
```
Y así sucesivamente para las otras comparaciones pareadas de nuestros tratamientos vs nuestro control.
###### tags: `Clase NGS`
Sign in with Wallet
Connect another wallet