Automatiser l’analyse RNA-seq avec un cluster 🚀

2. Automatiser l’analyse RNA-seq avec un cluster 🚀#

2.1. DĂ©couvrir le cluster#

Un cluster de calcul met Ă  disposition de ses utilisateurs des processeurs et de la mĂ©moire vive pour rĂ©aliser des calculs. Ces ressources sont rĂ©parties sur plusieurs machines, appelĂ©es nƓuds de calcul. Chaque nƓud de calcul contient plusieurs dizaines de processeurs (ou cƓurs) et plusieurs centaines de Go de mĂ©moire vive.

Pour connaütre le nombre de cƓurs disponibles sur le cluster de calcul, ouvrez un terminal dans JupyterLab et entrez la commande suivante :

$ sinfo -O cpusstate -p fast

Rappel

Ne tapez pas le caractÚre $ en début de ligne et faites attention aux majuscules et au minuscules.

Vous obtenez 4 nombres qui correspondent, dans l’ordre, au nombre de cƓurs utilisĂ©s pour le calcul (Allocated), disponibles pour le calcul (Idle), utilisĂ©s pour autre chose (Other) et le nombre total de cƓurs (Total). Normalement, vous devriez avoir suffisamment de cƓurs disponibles (Idle) pour rĂ©aliser vos calculs.

Voici un exemple de sortie :

$ sinfo -O cpusstate -p fast
CPUS(A/I/O/T)       
3018/2687/2523/8228

Sur les 8228 cƓurs du cluster, 3018 sont utilisĂ©s actuellement pour le calcul et 2687 sont disponibles.

Plusieurs utilisateurs utilisent simultanément un cluster de calcul. Pour vous en rendre compte, entrez la commande suivante qui liste tous les calculs (appelés job) en cours sur le cluster :

$ squeue -t RUNNING

Vous voyez que vous n’ĂȘtes pas seul ! Comptez maintenant le nombre de jobs en cours d’exĂ©cution en chaĂźnant la commande prĂ©cĂ©dente avec wc -l :

$ squeue -t RUNNING | wc -l

Vous pouvez aussi compter le nombre de jobs en attente d’exĂ©cution :

$ squeue -t PENDING | wc -l

Et vous alors ? Avez-vous lancé un calcul ? Vérifiez-le en entrant la commande suivante pour lister les calculs associés à votre compte :

$ squeue -u $USER

La colonne ST indique le statut de votre job :

  • CA (cancelled) : le job a Ă©tĂ© annulĂ©

  • F (failed) : le job a plantĂ©

  • PD (pending) : le job est en attente que des ressources soient disponibles

  • R (running) : le job est en cours d’exĂ©cution

Bizarre ! Vous avez un job (appelĂ© sys/dash) avec le statut running en cours d’exĂ©cution alors que vous n’avez, a priori, rien lancĂ© đŸ€”

En fait, le JupyterLab dans lequel vous ĂȘtes est lui-mĂȘme un job lancĂ© sur le cluster. C’est d’ailleurs pour cela qu’avant de lancer JupyterLab, vous avez dĂ» prĂ©ciser le compte Ă  utiliser (2501_duo) et choisir le nombre de processeurs et la quantitĂ© de mĂ©moire vive dont vous aviez besoin. Finalement, vous Ă©tiez dans la matrice sans mĂȘme le savoir đŸ˜±.

Comme plusieurs utilisateurs peuvent lancer des jobs sur un cluster, un gestionnaire de jobs (ou ordonnanceur) s’occupe de rĂ©partir les ressources entre les diffĂ©rents utilisateurs. Sur le cluster de l’IFB, le gestionnaire de jobs est Slurm. DĂ©sormais le lancement de vos analyses se fera via Slurm.

2.2. Analyser un Ă©chantillon#

DĂ©placez-vous dans le rĂ©pertoire /shared/projects/2501_duo/$USER/rnaseq puis vĂ©rifiez que vous ĂȘtes dans le bon rĂ©pertoire avec la commande pwd.

Supprimez les rĂ©pertoires qui contiennent les rĂ©sultats d’une Ă©ventuelle prĂ©cĂ©dente analyse :

$ rm -rf genome_index reads_qc reads_map counts

Téléchargez dans un terminal de JupyterLab un premier script Bash (script_cluster_0.sh) avec la commande wget :

$ wget https://raw.githubusercontent.com/pierrepo/unix-tutorial/master/content/tuto3/script_cluster_0.sh

Ce script correspond au script script_local_2.sh adaptĂ© pour une utilisation sur un cluster. Ouvrez le script script_cluster_0.sh avec l’éditeur de texte de JupyterLab et essayez d’identifier les diffĂ©rences avec script_local_2.sh.

  1. Au début du script, deux lignes sont nouvelles :

    #SBATCH --mem=2G
#SBATCH --cpus-per-task=8

    Ces lignes commencent par le caractĂšre # et indiquent qu’il s’agit de commentaires pour Bash, elles seront donc ignorĂ©es par le shell. Par contre, elles ont un sens trĂšs particulier pour le gestionnaire de jobs du cluster Slurm. Ici, ces lignes indiquent Ă  Slurm que le job a besoin de 2 Go de mĂ©moire vive et de 8 processeurs pour s’exĂ©cuter. Remarquez ensuite qu’on demande ici 8 cƓurs / processeurs alors notre session JupyterLab n’en a que 2. C’est tout Ă  fait normal, car le script va ĂȘtre lancĂ© comme un nouveau job, indĂ©pendant de celui de JupyterLab.

  2. Un peu plus loin, on indique explicitement les modules (les logiciels) Ă  charger avec leurs versions :

    module load fastqc/0.11.9
module load star/2.7.10b
module load samtools/1.15.1
module load htseq/0.13.5
module load cufflinks/2.2.1

  3. L’appel aux diffĂ©rents programmes (STAR, fastqc, samtools
) est prĂ©fixĂ© par l’instruction srun qui va explicitement indiquer Ă  Slurm qu’il s’agit d’un sous-job. Nous en verrons son utilitĂ© plus tard.

  4. Pour STAR, l’indication du nombre de cƓurs Ă  utiliser est dĂ©fini sous la forme :

    srun STAR --runThreadN "${SLURM_CPUS_PER_TASK}" \

    La variable SLURM_CPUS_PER_TASK est utilisĂ©e pour indiquer Ă  STAR le nombre de processeurs Ă  utiliser. Cette variable est automatiquement dĂ©finie par Slurm lors de l’exĂ©cution du script et la lecture de l’instruction #SBATCH --cpus-per-task=8. Dans le cas prĂ©sent, elle vaut 8.

Lancez enfin le script avec la commande suivante :

$ sbatch -A 2501_duo script_cluster_0.sh

Vous devriez obtenir un message du type Submitted batch job 33389786. Ici, 33389786 correspond au numéro du job. Notez bien le numéro de votre job.

  • L’instruction sbatch est la maniĂšre de demander Ă  Slurm de lancer un script et crĂ©er un job indĂ©pendant de celui de JupyterLab.

  • L’option -A 2501_duo spĂ©cifie quel projet utiliser (facturer) pour cette commande. Un mĂȘme utilisateur peut appartenir Ă  plusieurs projets. Le nombre d’heures de calcul attribuĂ©es Ă  un projet Ă©tant limitĂ©, il est important de savoir quel projet imputer pour telle ou telle commande. Pensez-y pour vos futurs projets.

VĂ©rifiez que votre script est en train de tourner avec la commande :

$ squeue -u $USER

Le statut du job devrait ĂȘtre RUNNING (R).

Et pour avoir plus de détails, utilisez la commande :

$ sacct --format=JobID,JobName,State,Start,Elapsed,CPUTime,NodeList -j JOBID

avec JOBID (en fin de ligne) le numéro du job à remplacer par le vÎtre.

Voici un exemple de sortie que vous pourriez obtenir :

$  sacct --format=JobID,JobName,State,Start,Elapsed,CPUTime,NodeList -j 33389786
       JobID    JobName      State               Start    Elapsed    CPUTime        NodeList 
------------ ---------- ---------- ------------------- ---------- ---------- --------------- 
33389786     script_cl+    RUNNING 2023-05-12T10:25:21   00:00:46   00:06:08     cpu-node-11 
33389786.ba+      batch    RUNNING 2023-05-12T10:25:21   00:00:46   00:06:08     cpu-node-11 
33389786.0         STAR  COMPLETED 2023-05-12T10:25:23   00:00:08   00:01:04     cpu-node-11 
33389786.1       fastqc    RUNNING 2023-05-12T10:25:31   00:00:36   00:04:48     cpu-node-11

Quand la colonne State est Ă  COMPLETED, cela signifie que le sous-job est terminĂ©. Un sous-job est crĂ©Ă© Ă  chaque fois que la commande srun est utilisĂ©e dans le script de soumission. Ici, on voit que le sous-job STAR est terminĂ©, mais que le sous-job fastqc est toujours en cours d’exĂ©cution.

Relancez rĂ©guliĂšrement la commande prĂ©cĂ©dente pour suivre l’avancement de votre job.

Note

Pour rappeller une commande précédente, vous pouvez utiliser la flÚche du haut de votre clavier. Cela vous évitera de retaper la commande à chaque fois.

Notez le moment oĂč vous avez terminĂ© l’alignement sur le gĂ©nome de rĂ©fĂ©rence, c’est-Ă -dire que vous avez obtenu un second sous-job STAR avec le statut COMPLETED dans la sortie renvoyĂ©e par sacct. Par exemple :

$  sacct --format=JobID,JobName,State,Start,Elapsed,CPUTime,NodeList -j 33389786
       JobID    JobName      State               Start    Elapsed    CPUTime        NodeList 
------------ ---------- ---------- ------------------- ---------- ---------- --------------- 
33389786     script_cl+    RUNNING 2023-05-12T10:25:21   00:03:25   00:27:20     cpu-node-11 
33389786.ba+      batch    RUNNING 2023-05-12T10:25:21   00:03:25   00:27:20     cpu-node-11 
33389786.0         STAR  COMPLETED 2023-05-12T10:25:23   00:00:08   00:01:04     cpu-node-11 
33389786.1       fastqc  COMPLETED 2023-05-12T10:25:31   00:01:23   00:11:04     cpu-node-11 
33389786.2         STAR  COMPLETED 2023-05-12T10:26:54   00:01:36   00:12:48     cpu-node-11 
33389786.3     samtools    RUNNING 2023-05-12T10:28:30   00:00:16   00:02:08     cpu-node-11

Nous ne souhaitons pas poursuivre plus longtemps cette analyse, car cela prendrait trop de temps.

Stoppez le job en cours avec la commande :

$ scancel JOBID

avec JOBID le numéro de votre job.

VĂ©rifiez avec la commande squeue -u $USER que votre job ne tourne plus.

Note

Le fichier slurm-JOBID.out (avec JOBID le numĂ©ro de votre job) contient le « journal » de votre job. C’est-Ă -dire tout ce qui s’affiche habituellement Ă  l’écran a Ă©tĂ© enregistrĂ© dans ce fichier.

Ouvrez ce fichier avec l’éditeur de texte de JupyterLab pour vous en rendre compte.

2.3. Analyser 50 Ă©chantillons (ou pas loin)#

Maintenant que nous savons comment analyser un Ă©chantillon avec un cluster de calcul, nous allons automatiser l’analyse de 50 Ă©chantillons. Mais pour cela nous allons prendre plusieurs prĂ©cautions :

  • Nous allons sĂ©parer les Ă©tapes suivantes :

    1. Indexer le gĂ©nome de rĂ©fĂ©rence – ne doit ĂȘtre fait qu’une seule fois.

    2. ContrĂŽler la qualitĂ©, aligner et quantifier les reads – doit ĂȘtre fait pour chaque Ă©chantillon, donc 50 fois.

    3. Normaliser les comptages de tous les Ă©chantillons ensemble – ne doit ĂȘtre fait qu’une seule fois.

  • L’étape 2 ne doit pas se faire successivement pour chaque Ă©chantillon, mais en parallĂšle. C’est-Ă -dire qu’on souhaite idĂ©alement lancer l’analyse des 50 Ă©chantillons en mĂȘme temps.

2.3.1. Préparer les données#

Nous avons téléchargé pour vous les 50 échantillons (fichiers .fastq.gz) ainsi que le génome de référence et ses annotations dans le répertoire /shared/projects/2501_duo/data/rnaseq_scere. Vérifiez son contenu avec la commande :

$ tree /shared/projects/2501_duo/data/rnaseq_scere

2.3.2. Lancer le script de l’étape 1 (indexer le gĂ©nome de rĂ©fĂ©rence)#

Supprimez les rĂ©pertoires qui contiennent les rĂ©sultats d’une Ă©ventuelle prĂ©cĂ©dente analyse :

$ rm -rf genome_index reads_qc reads_map counts slurm*.out

Téléchargez dans un terminal de JupyterLab le script Bash (script_cluster_1.sh) avec la commande wget :

$ wget https://raw.githubusercontent.com/pierrepo/unix-tutorial/master/content/tuto3/script_cluster_1.sh

Ouvrez ce script avec l’éditeur de texte de JupyterLab et essayez d’identifier comment sont dĂ©finies les variables base_dir et data_dir.

Puis lancez le script en n’oubliant pas d’indiquer explicitement le compte à utiliser (2501_duo) :

$ sbatch -A 2501_duo script_cluster_1.sh

Vérifiez avec les commandes squeue et sacct que le job est lancé et se termine correctement (normalement rapidement).

2.3.3. Lancer le script de l’étape 2 (traiter les Ă©chantillons)#

Téléchargez dans un terminal de JupyterLab le script Bash (script_cluster_2.sh) avec la commande wget :

$ wget https://raw.githubusercontent.com/pierrepo/unix-tutorial/master/content/tuto3/script_cluster_2.sh

Ouvrez ce script avec l’éditeur de texte de JupyterLab. Notez en dĂ©but de script la ligne :

#SBATCH --array=0-49

qui nous permettra de crĂ©er un job array, c’est-Ă -dire un ensemble de 50 (de 0 Ă  49 inclus) jobs qui vont ĂȘtre lancĂ©s en parallĂšle.

Pour chacun des jobs, on lui attribue un Ă©chantillon avec les commandes suivantes :

# liste de tous les fichiers .fastq.gz dans un tableau
fastq_files=(${fastq_dir}/*fastq.gz)
# extraction de l'identifiant de l'Ă©chantillon
# Ă  partir du nom de fichier : /shared/projects/2501_duo/data/rnaseq_scere/reads/SRR3405783.fastq.gz
# on extrait : SRR3405783
sample=$(basename -s .fastq.gz "${fastq_files[$SLURM_ARRAY_TASK_ID]}")

Cette Ă©tape est importante, car elle permet de savoir quel Ă©chantillon traiter pour chaque job. Par exemple, le job 0 va ĂȘtre associĂ© Ă  l’échantillon SRR3405783 (le premier par ordre alphabĂ©tique).

Pour ne pas emboliser le cluster et pour que tout le monde puisse obtenir des résultats rapidement, modifiez la ligne

#SBATCH --array=0-49

par

#SBATCH --array=0-2

Vous n’analyserez que 3 Ă©chantillons dans 3 jobs en parallĂšle. Pensez Ă  sauvegarder vos modifications.

Lancez enfin le script :

$ sbatch -A 2501_duo script_cluster_2.sh

Suivez la progression de vos jobs avec la commande :

$ sacct --format=JobID,JobName,State,Start,Elapsed,CPUTime,NodeList -j JOBID

avec JOBID le numéro de votre job.

PlutĂŽt que de relancer en permanence la commande sacct vous pouvez demander Ă  la commande watch d’afficher les Ă©ventuels changements :

$ watch sacct --format=JobID,JobName,State,Start,Elapsed,CPUTime,NodeList -j JOBID

Vous devriez obtenir un affichage du type :

Every 2.0s: sacct --format=JobID,JobName,State,Start,Elapsed,CPUTime,NodeList -j ...  Wed May 10 22:22:53 2023

       JobID    JobName      State               Start    Elapsed    CPUTime        NodeList
------------ ---------- ---------- ------------------- ---------- ---------- ---------------
33361021_0   script_cl+    RUNNING 2023-05-10T22:21:54   00:00:59   00:07:52     cpu-node-20
33361021_0.+      batch    RUNNING 2023-05-10T22:21:54   00:00:59   00:07:52     cpu-node-20
33361021_0.0     fastqc    RUNNING 2023-05-10T22:21:56   00:00:57   00:07:36     cpu-node-20
33361021_1   script_cl+    RUNNING 2023-05-10T22:21:54   00:00:59   00:07:52     cpu-node-23
33361021_1.+      batch    RUNNING 2023-05-10T22:21:54   00:00:59   00:07:52     cpu-node-23
33361021_1.0     fastqc    RUNNING 2023-05-10T22:21:56   00:00:57   00:07:36     cpu-node-23
33361021_2   script_cl+    RUNNING 2023-05-10T22:21:54   00:00:59   00:07:52     cpu-node-25
33361021_2.+      batch    RUNNING 2023-05-10T22:21:54   00:00:59   00:07:52     cpu-node-25
33361021_2.0     fastqc    RUNNING 2023-05-10T22:21:56   00:00:57   00:07:36     cpu-node-25

On apprend ici que 3 jobs sont en cours d’exĂ©cution (RUNNING) et qu’ils sont appelĂ©s 33361021_0, 33361021_1 et 33361021_2. Chacun de ces jobs a Ă©tĂ© lancĂ© sur un noeud de calcul diffĂ©rent (cpu-node-20, cpu-node-23 et cpu-node-25) mais cela aurait pu ĂȘtre le mĂȘme. Chaque job est dĂ©composĂ© en sous-jobs, pour le moment Ă  l’étape fastqc.

Le temps que les 3 jobs se terminent, profitez-en pour faire une pause cafĂ© ☕ bien mĂ©ritĂ©e et rĂ©aliser Ă  quel point la bioinformatique et ses possibilitĂ©s d’automatisation sont cools.

Hint

Utilisez la combinaison de touches Ctrl + C pour arrĂȘter la commande watch.

Quand tous les jobs sont Ă  COMPLETED, comparez le temps d’exĂ©cution de chacun (dans la colonne Elapsed au niveau des lignes script_cl+). Ce temps d’exĂ©cution devrait ĂȘtre compris entre 15 et 20 min, ce qui est Ă©quivalent au temps d’exĂ©cution du script script_local_2.sh pour un seul Ă©chantillon. C’est tout l’intĂ©rĂȘt de lancer des jobs en parallĂšle.

Remarquez que le temps CPU (CPUTime) qui correspond au temps de calcul consommĂ© par tous les processeurs est supĂ©rieur au temps « humain » (Elapsed). Ainsi, si vous avez lancĂ© un job avec 8 coeurs qui se termine en 1 heure, la consommation CPU sera de 8 coeurs x 1 heure = 8 heures. Sur un cluster, c’est toujours le temps CPU qui est facturĂ©.

2.3.4. Lancer le script de l’étape 3 (normaliser les comptages)#

Vous allez maintenant normaliser les différents comptages.

Téléchargez dans un terminal de JupyterLab le script Bash (script_cluster_3.sh) avec la commande wget :

$ wget https://raw.githubusercontent.com/pierrepo/unix-tutorial/master/content/tuto3/script_cluster_3.sh

Ouvrez ce script avec l’éditeur de texte de JupyterLab. Retrouvez les lignes spĂ©cifiques Ă  l’utilisation d’un cluster de calcul et qui dĂ©butent par :

  • #SBATCH

  • module load

  • srun

Puis lancez le script :

$ sbatch -A 2501_duo script_cluster_3.sh

Affichez l’avancement de votre job avec la commande sacct :

$ sacct --format=JobID,JobName,State,Start,Elapsed,CPUTime,NodeList -j JOBID

ou

$ watch sacct --format=JobID,JobName,State,Start,Elapsed,CPUTime,NodeList -j JOBID

avec JOBID le numéro de votre job.

L’exĂ©cution de ce script devrait moins d’une minute.

Le fichier qui contient le comptage normalisé des transcrits est counts/genes.count_table.

2.4. Bilan#

Vous avez lancĂ© une analyse RNA-seq complĂšte en utilisant le gestionnaire de ressources (Slurm) d’un cluster de calcul. Bravo ✹

Un peu plus tard, nous vous inviterons à reprendre cette analyse, mais cette fois sur les 50 échantillons. Pour cela, il vous faudra modifier le script script_cluster_2.sh en remplaçant la ligne #SBATCH --array=0-2 (pour 3 échantillons) par #SBATCH --array=0-49 (pour 50 échantillons). Pensez aussi à relancer le script script_cluster_3.sh pour normaliser les résultats de comptage.

2.5. Aller plus loin : connecter les jobs#

Note

Vous n’ĂȘtes pas obligĂ© de lancer les commandes ci-dessous. Elles sont lĂ  pour vous montrer comment automatiser le lancement de plusieurs jobs les uns aprĂšs les autres.

Pour cette analyse, il faut lancer 3 scripts Ă  la suite : script_cluster_1.sh, script_cluster_2.sh et script_cluster_3.sh. À chaque fois, il faut attendre que le prĂ©cĂ©dent soit terminĂ©, ce qui peut ĂȘtre pĂ©nible. Slurm offre la possibilitĂ© de chaĂźner les jobs les uns avec les autres avec l’option --dependency. Voici un exemple d’utilisation.

Tout d’abord on lance le script script_cluster_1.sh :

$ sbatch -A 2501_duo script_cluster_1.sh
Submitted batch job 33390286

On récupÚre le job id (33390286) puis on lance immédiatement le script script_cluster_2.sh :

$ sbatch -A 2501_duo --dependency=afterok:33390286 script_cluster_2.sh
Submitted batch job 33390299

On récupÚre le nouveau job id (33390299) puis on lance immédiatement le dernier script script_cluster_3.sh :

$ sbatch -A 2501_duo --dependency=afterok:33390299 script_cluster_3.sh
Submitted batch job 33390315

Dans l’exemple ci-dessus, le job 33390315 (pour script_cluster_3.sh) ne va s’exĂ©cuter que quand le job 33390299 (pour script_cluster_2.sh) sera terminĂ© avec succĂšs. Et le job 33390299 ne va s’exĂ©cuter que quand le job 33390286 (pour script_cluster_1.sh) sera, lui aussi, terminĂ© avec succĂšs.

Voici le résultat obtenu avec squeue :

$ squeue -u ppoulain
             JOBID PARTITION     NAME     USER ST       TIME  NODES NODELIST(REASON)
          33390315      fast script_c ppoulain PD       0:00      1 (Dependency)
          33389748      fast sys/dash ppoulain  R      28:36      1 cpu-node-15
        33390299_0      fast script_c ppoulain  R       1:44      1 cpu-node-2
        33390299_1      fast script_c ppoulain  R       1:44      1 cpu-node-30
        33390299_2      fast script_c ppoulain  R       1:44      1 cpu-node-35

Dans cet exemple le premier job (33390286) est dĂ©jĂ  terminĂ© (et n’apparait pas). Le job 33390299 est en cours d’exĂ©cution (pour 3 Ă©chantillons seulement) et le job 33390315 est en attente que le job 33390299 se termine.

Cette mĂ©thode Ă©vite d’attendre que le job prĂ©cĂ©dent se termine pour lancer le suivant, mais il faut quand mĂȘme les lancer manuellement pour rĂ©cupĂ©rer les job ids. On peut automatiser cela avec les commandes suivantes :

jobid1=$(sbatch -A 2501_duo script_cluster_1.sh | cut -d' ' -f 4)
echo "Running job 1: ${jobid1}"
jobid2=$(sbatch -A 2501_duo --dependency=afterok:${jobid1} script_cluster_2.sh | cut -d' ' -f 4)
echo "Running job 2: ${jobid2}"
sbatch -A 2501_duo --dependency=afterok:${jobid2} script_cluster_3.sh
squeue -u $USER

La commande squeue Ă  la fin permet simplement de vĂ©rifier que le premier job est en cours d’exĂ©cution et que les deux autres sont en attente.

Contents