# ** U5 Análisis de expr dif con DESeq2**
####
---
## Acceso a [OMICA](http://omica.cicese.mx/docs/index.html)
Acceder al servicor con tu cuenta de omica.
Utilizar `-YC` cuando se requiera despliegue gráfico. Los símbolos ## indican el número de curso asignado al usuario.
```javascript=
ssh Usuario@omica
```
En la carpeta de `Transcriptomica`, crear la carpeta `DIFF_EXP` y acceder.
```javascript=
mkdir DIFF_EXP
cd DIFF_EXP
```
Crear un `link simbólico` del ensamble de buena calidad obtenido con Transrate `good.Trinity.fasta`.
```javascript=
ln -s ../TRANSRATE/gcontigs/Trinity/good.Trinity.fasta ./
ls
```
Crear un `link simbólico` de las lecturas pareadas y limpias obtenidas con trimmomatic.
```javascript=
ln -s ../TRIMMOMATIC/TRIM_P/*.fq.gz ./
ls
```
Crear un archivo `samples_file.txt``
```javascript=
nano samples_file.txt
```
La estructura del archivo se muestra a continuación:
```javascript=
IFM IFM_1 INT5_1P.fastq INT5_2P.fastq
IFM IFM_2 INT6_1P.fastq INT6_2P.fastq
IFM IFM_3 INT15_1P.fastq INT15_2P.fastq
I30 I30_1 INT8_1P.fastq INT8_2P.fastq
I30 I30_2 INT9_1P.fastq INT9_2P.fastq
I30 I30_3 INT23_1P.fastq INT23_2P.fastq
I60 I60_1 INT26_1P.fastq INT26_2P.fastq
I60 I60_2 INT27_1P.fastq INT27_2P.fastq
```
Verificar la estructura tabular del archivo `samples_file.txt`. La tabulación se muestra con el simbolo `^I` y el final de la línea con simbolo `$`.
```javascript=
cat -TE samples_file.txt
```
Generar el script al cual llamaremos `mapping_rsem.slrm` para someter al manejador de tareas del Cluster omica.
```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_ensamble2.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
```
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).
Con los datos de los conteos se procedera a generar las matrices de abundancia. Para ello, generar el script al cual llamaremos `abundance_matrix.slrm` para realizar la matriz de conteos.
```javascript=
nano abundance_matrix.slrm
```
El slrm debe contener lo siguiente:
```javascript=
#!/bin/bash
#SBATCH -p cicese
#SBATCH --job-name=abundance_est
#SBATCH --output=abundance_est-%j.log
#SBATCH --error=abundance_est-%j.err
#SBATCH -N 1
#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/
export PATH=$PATH:/LUSTRE/apps/bioinformatica/trinityrnaseq-v2.15.1/util/support_scripts/
module load conda-2023
module load trinityrnaseq-v2.15.1
UTIL=/LUSTRE/apps/bioinformatica/trinityrnaseq-v2.15.1
cd $SLURM_SUBMIT_DIR
$UTIL/abundance_estimates_to_matrix.pl \
--est_method RSEM \
--name_sample_by_basedir \
./IFM_1/RSEM.isoforms.results \
./IFM_2/RSEM.isoforms.results \
./IFM_3/RSEM.isoforms.results \
./I30_1/RSEM.isoforms.results \
./I30_2/RSEM.isoforms.results \
./I30_3/RSEM.isoforms.results \
./I60_1/RSEM.isoforms.results \
./I60_2/RSEM.isoforms.results \
--gene_trans_map good.Trinity_ensamble2.fasta.gene_trans_map \
exit 0
```
En caso de que el archivo `.err` marque un error de permiso denegado para acceder a las carpetas, se deberan cambiar los permisos de las carpetas generadas con RSEM de la siguiente manera:
```javascript=
chmod 777 FM* SPC30* SPC60* GLiv*
```
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 FM_rep1
:::
Cargar modulo de R de la siguiente manera y someter slrm a manejador de tareas:
```javascript=
module load R-3.5.0
sbatch abundance_matrix.slrm
```
Con esto, se habran generado las matrices de abundancia.
Los siguentes pasos enel procesamiento de datos se realizan en R studio, para ello se recominda utilizar la version 4.0.4 o 4.2.1 de R.
En la
Posteriormente, ger=nerar un archivo de texto con la leyenda de las muestras al cual denominaremos `samples_described.txt` y se generara con la siguiente estructura.
```javascript=
IFM IFM_1
IFM IFM_2
IFM IFM_3
I30 I30_1
I30 I30_2
I30 I30_3
I60 I60_1
I60 I60_2
```
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
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.
:::
Realizar el análisis de expresión per se. Para esto, nos movemos al directorio `DIFF_EXP`,
Crear un script al que denominaremos `diff_express.slrm`.
```javascript=
#!/bin/sh
#SBATCH -p cicese
#SBATCH --job-name=DESeq2
#SBATCH --output=DE-%j.log
#SBATCH --error=DE-%j.err
#SBATCH -N 1
#SBATCH --mem=100GB
#SBATCH --ntasks-per-node=24
#SBATCH -t 6-00:00:00
#SBATCH --exclusive
DE=/LUSTRE/apps/bioinformatica/trinityrnaseq-v2.15.1/Analysis/DifferentialExpression
cd $SLURM_SUBMIT_DIR
$DE/run_DE_analysis.pl \
--matrix RSEM.isoform.counts.matrix \
--method DESeq2 \
--samples_file samples_described.txt \
exit 0
```
Cargar modulo de R
```javascript=
module load R-3.5.0
```
:::info
En caso de que el script de algùn error con respecto a paqueterias de R, utilizar el siguiente script directo en la terminal.
```javascript=
R
BiocManager::install(c("edgeR", "limma", "DESeq2", "ctc", "Biobase", "gplots", "ape", "argparse"))
```
:::
Esto generará una carpeta 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.
Para extraer los transcritos expresados diferencialmente con valores de P y magnitud de cambio (fold change) específicos, nos movemos al directorio
DESeq2 y ejecutamos el script `significant_DE.slrm`.
```javascript=
cd DESeq2.#####.dir
nano significant_DE.slrm
```
El slrm debera contener lo siguiente:
```javascript=
#!/bin/sh
#SBATCH -p cicese
#SBATCH --job-name=sign_DE
#SBATCH --output=sign_DE-%j.log
#SBATCH --error=sign_DE-%j.err
#SBATCH -N 1
#SBATCH --mem=100GB
#SBATCH --ntasks-per-node=24
#SBATCH -t 6-00:00:00
#SBATCH --exclusive
DE=/LUSTRE/apps/bioinformatica/trinityrnaseq-Trinity-v2.8.5/Analysis/DifferentialExpression
cd $SLURM_SUBMIT_DIR
$DE/analyze_diff_expr.pl \
--matrix ../RSEM.isoform.TMM.EXPR.matrix \
--samples ../samples_described.txt \
-P 0.01 -C 1 \
exit 0
```
Someter a manejador de tareas del cluster omica.
```javascript=
sbatch significant_DE.slrm
```
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). Esto generar el heatmap correspondiente en formato `.pdf`.
De las comparaciones pareadas podemos extraer la lista de los transcritos UP y DOWN para ingresarlos en la página de VENNY y obtener los diagramas de VENN. Por ejemplo, si queremos extraer los transcritos UP y DOWN del tratamiento FM100_PREB con respecto al control que es FM100, ejecutamos lo siguiente:
```javascript
cat RSEM#FM100_vs_FM100_PREB#FM100_PREB-UP.subset | cut –f 1 >
FM100_PREB_UP.txt
cat RSEM#FM100_vs_FM100_PREB#FM100-UP.subset | cut –f 1 >
FM100_PREB_DOWN.txt
```
Y así sucesivamente para los otros tratamientos. Estos archivos los copias a tu computadora personal. El símbolo # en los comandos anteriores sólamente se usaron para no poner todo el nombre completo en el texto.
###### tags: `Clase NGS`
Sign in with Wallet
Connect another wallet