Mise en pratique des méthodes d'analyse de l'expression génique grâce aux données utilisées par le package airway: création des index nécessaires aux outils d'alignement

• Ici, on va avoir notre première application de l'utilisation de SLURM La construction des index demande à la base un ensemble de fichiers contenant les informations décrivant notre génome d'intérêt. Il y a plusieurs sources potentiels de tels fichiers mais un bon point de départ lorsqu'on travaille avec H. sapiens est le site du projet Gencode. Gencode est un effort international de construction des annotations provenant de diverses sources publiques en y mettant une couche de curation humaine et pas simplement construit programmatiquement. Comme ce processus d'annotation est en mutation constante, Gencode crée les fichiers sous la forme de versions (releases) qui sont publiés sur une base semi-régulière.
• Deuxième application de SLURM sur grappe de calcul: construire les index. Nous allons utiliser trois approches d'alignement: via HISAT2, via STAR et via Salmon. Dans tous les cas, créer les index pour chaque outil profite de l'usage d'une grappe de calcul: la création est accélérée par une exécution multi-fils (multi-threaded) et dans tous les cas, demande une très grande quantité de mémoire pour s'exécuter. Comme il faut créer les index à chaque nouvelle version de Gencode (ou de tout autre fichier d'annotation), mettre ce processus dans un script facile à éditer facilite le travail.

Une fois les index créés, on va vers l'obtention des fichiers de séquence bruts du projet.

• Pour commencer, il faut aller chercher les fichiers d'annotation sur le site du projet Gencode. En partant de notre espace personnel ($HOME), on se dirige là où l'on désire mettre les index pour la suite:

# On assume que l'on se dirige vers /shares/data pour créer
# les répertoires nécessaires
# Ajuster en fonction de votre environnement, évidemment!
% mkdir /shares/data/indexes
# Remarquez que les données d'annotation sont mises dans
# une section distincte de celui des index
% mkdir /shares/data/annotations
% cd /shares/data/annotations
#
# Ici, on utilise une logique sous la forme source/version
# sous /shares/data/annotations
#
% mkdir gencode
# Dernière version en date
# Avril 2026
% mkdir gencode/r49
% cd gencode/r49
# Fichier avec les annotations complètes en format GTF
% curl -L -O https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_49/gencode.v49.annotation.gtf.gz
# Fichier avec les annotations complètes en format GFF3
% curl -L -O https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_49/gencode.v49.annotation.gff3.gz
# Fichier FASTA des transcrits pour protéines, utilisé par Salmon
% curl -L -O https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_49/gencode.v49.pc_transcripts.fa.gz
# Fichier FASTA du genome humain, GRCh38
% curl -L -O https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_49/GRCh38.primary_assembly.genome.fa.gz
# Decompression des fichiers
% for i in `ls | grep gz`; do gunzip $i; done

• Si une nouvelle version (disons r50) est publiée, simplement créez un répertoire sous /shares/data/annotations/gencode et téléchargez les nouvelles versions.

Créer les index pour les outils d'alignement est un processus exigeant et le premier exemple où l'utilisation d'une grappe de calcul et d'un système de gestion des tâches qui y réside est un atout. Dans le cas de la création des index, on est confronté à quelques problèmes:

• C'est un processus qui prend du temps… Les outils sont capables d'opérations multi-fils pour accélérer la tâche alors plus on a de coeurs de calcul, plus ça se fera rapidement.
• C'est un processus qui construit des structures de données qui habitent en premier en mémoire vive et ça demande beaucoup, beaucoup de mémoire vive.

Si on prend pour exemple la grappe de calcul Rorqual, on y trouve des serveurs qui ont tout ce que nous aurons de besoin, il faut simplement leur dire comment faire. En passant, faire cette construction sur une grappe SuperClafoutis n'est pas une option: certaines opérations d'indexation demandent des centaine de Gb de mémoire vive Cependant, une fois les index construits, on peut les télécharger pour utilisation locale car ces fichiers sont indépendants de la plateforme où ils ont été calculés.

Comme ce sont des tâches qui demandent du temps (plusieurs minutes à plusieurs heures), c'est le moment idéal pour apprendre à utiliser SLURM. À la base, SLURM est un système où un usager soumet des actions (interactives ou scriptées) via une application (srun ou sbatch) à une queue de traitement qui gèrera la suite des choses. Une fois la tâche assignée à un noeud de calcul, elle s'exécute indépendamment de votre attention (dans le cas de sbatch).

• HISAT2 construit une série de fichiers qui constituera l'index servant à chercher le meilleur alignement pour chaque séquence sur un génome de référence. Il profite des informations sur la position des exons et des sites d'épissage contenus dans un fichier d'annotations en format GTF et, dans le cas de H. sapiens, il peut aussi utiliser un fichier contenant les infos des variations génétiques, soit en format VCF ou bien en provenance de UCSC avec son format particulier. Si ces fichiers sont disponibles (et ils le sont dans les archives publiques), on utilise des scripts écrit en Python pour créer des fichiers intermédiaires qui seront utiliser par l'utilitaire hisat2-build pour construire les fichiers nécessaires. Quant aux ressources nécessaires, la création des index avec les infos concernant les transcrits (sites d'épissage avec bordures intron-exon et liste des exons) est possible un ordinateur avec au moins 32 Gb de mémoire vive. Si vous voulez ajouter les infos sur les variations génétiques, il vous faudra plusieurs centaines de Gb de mémoire vive…

• Pour commencer, on passe par les scripts Python fournis avec HISAT2 pour la création des fichiers intermédiaires:

# On crée l'arborescence nécessaire sous /shares/data
% cd /shares/data/indexes
#
# Remarquez que nous reprenons la structure logique
# de l'organisation des fichiers retrouvés sous
# /shares/data/annotations
#
# Construction des index utilisés spécifiquement par HISAT2
% mkdir hisat2
% mkdir hisat2/gencode
% mkdir hisat2/gencode/r49 
% mkdir hisat2/gencode/r49/interim_files
% cd hisat2/gencode/r49/interim_files
#
# Ici, on assume que HISAT2 est sur le $PATH; les scripts devraient se trouver 
# au même niveau que l'application.
# Il faut rediriger la sortie vers un fichier sinon ça sort sur STDOUT...
% hisat2_extract_exons.py /shares/data/annotations/gencode/r49/gencode.v49.annotation.gtf > ./gencode_r49_exons.txt
% hisat2_extract_splice_sites.py /shares/data/annotations/gencode/r49/gencode.v49.annotation.gtf > ./gencode_r49_ss.txt
# Etape facultative: extraction des infos pour les variations génétiques
# Dernière version disponible des données UCSC
% curl -L -O http://hgdownload.soe.ucsc.edu/goldenPath/hg38/database/snp151Common.txt.gz
% gunzip snp151Common.txt.gz
# Ici, le script est un peu plus exigeant... Ça prend plus de temps que les étapes précédentes mais ça se fait sur un RPi.
# A la fin, vous obtiendrez deux fichiers: un avec .haplotype et un autre .snp en suffixe avec gencode_r49_snp_haplo comme préfixe.
% hisat2_extract_snps_haplotypes_UCSC.py /shares/data/annotations/gencode/r49/GRCh38.primary_assembly.genome.fa ./snp151Common.txt gencode_r49_snp_haplo

• Ok, pour la suite, ça nous prend un serveur bien doté en mémoire vive et SLURM pour lancer cette tâche. Il nous faut écrire un script bash qui sera donner en entrée à l'application sbatch pour soumettre à l'ordonnanceur de tâches de SLURM. On doit écrire dans l'entête du script des instructions qui seront lues par l'ordonnanceur et qui permettront d'utiliser les bonnes ressources pour l'exécution. Le script qui suit est simplement un exemple ; votre propre grappe de calcul nécessitera fort probablement des instructions différentes…

• Au minimum, on veut avoir dans l'index les informations pour les transcrits via les fichiers d'exons et de bordures exon-intron que nous avons construit ci-dessus. À l'aide de nano, écrire le texte suivant:

#!/bin/bash -l
#
##
# Nous sommes sous /shares/data/indexes/hisat2/gencode/r49
# (ou à l'endroit approprié de votre système de fichiers dans 
# une grappe de calcul)
## 
#SBATCH --account=nom_compte
#SBATCH --job-name=hisat2-build_%j
#SBATCH --output=hisat2-build_%j.out
#SBATCH --error=hisat2-build_%j.err
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=4
#SBATCH --mem=256G
#SBATCH --time=02:00:00
#
# On assume que l'application hisat2-build est sur le $PATH
#
# Sur Rorqual, il faut charger l'application; retirez le
# dièse de la prochaine ligne pour que le script fonctionne
#module load hisat2

hisat2-build -p 4 \
--exon ./interim_files/gencode_r49_exons.txt \
--ss ./interim_files/gencode_r49_ss.txt \
/shares/data/annotations/gencode/r49/GRCh38.primary_assembly.genome.fa \
./h_sapiens_gencode_r49_exons_ss
# Sauvegarder votre fichier avec un nom comme hisat2-builder.sh

• Ok, on a plusieurs lignes qui à première vue, semblent bizarres et pas des commandes… En fait, ça fait partie de l'exécution du script car ce sont des paramètres d'exécution.
  • #!/bin/bash -l: cette ligne dicte que le script sera exécuté via l'interpréteur bash pour comprendre le reste du script.
  • #SBATCH –account=nom_compte: paramètre nécessaire (par exemple, pour le cluster Rorqual) pour la comptabilité des tâches selon les groupes d'usagers.
  • #SBATCH –job-name=hisat2-build_%j: chaque ligne commençant par #SBATCH est en fait une instruction pour sbatch, dans ce cas, la tâche portera le nom hisat2builder_%j où %j est une variable spécifique à sbatch attribuée au moment de la soumission et désignant le numéro de la tâche dans la queue de traitement.
  • #SBATCH –output=hisat2-build_%j.out: nom du fichier de sortie sur STDOUT. Dans notre cas, c'est pratique d'inclure le numéro de la tâche pour faire le suivi.
  • #SBATCH –error=hisat2-build_%j.err: nom du fichier de sortie des erreurs sur STDERR.
  • #SBATCH –ntasks=1: nous avons une seule tâche à accomplir.
  • #SBATCH –cpus-per-task=4: cette tâche demande 4 CPU pour en assurer l'exécution. Assurez-vous que cette valeur soit la même que celle du paramètre -p utilisé par HISAT2
  • #SBATCH –mem=256G: la tâche devrait être attribuée à un serveur avec au moins 256Gb de mémoire vive disponible (eh oui, c'est beaucoup mais nécessaire pour construire les fichiers d'index en utilisant les infos sur les exons et les bordures exon-intron).
    • Si vous voulez inclure les infos sur les variations génétiques, il vous en faudra encore plus, autour du triple de cette valeur.
    • Vu la quantité de mémoire vive nécessaire, une grappe de calcul équipé de machines musclées, comme Rorqual, est absolument nécessaire.
  • #SBATCH –time=02:00:00: la tâche sera abandonnée si elle dure plus de 2 heures. En théorie, ça devrait en prendre moins avec 4 coeurs de calcul pour l'exécution.

• Deux autres paramètres peuvent être utiles ou nécessaire:
  • #SBATCH –mail-user=courriel@institution.org: une adresse de courriel valide pour envoyer les états de la tâche sur le cluster.
  • #SBATCH –mail-type=BEGIN,END,FAIL: le type d'information à communiquer à l'usager. Ici, on demande d'envoyer un courriel au début de l'exécution, à la fin de l'exécution et en cas d'erreur.

• En restant sous /shares/data/indexes/hisat2/gencode/r49, on soumet le script hisat2-builder.sh à la grappe de calcul avec la commande suivante:

% sbatch hisat2-builder.sh
# Vous recevrez sur le terminal le numéro de votre tâche.

• Pour suivre le déroulement de l'exécution, on utilise la commande squeue pour voir:

# Une option:
% squeue | grep nom_usager
# Une autre option avec le numéro de la tâche:
% squeue | grep job_id

• La colonne de la sortie de squeue à surveiller est celle appelée ST pour Status. Si elle indique que la tâche est PD (pour Pending), ça veut dire que la tâche est en attente d'exécution; si elle indique R (pour Running, duh!), ça veut dire qu'elle est en exécution.

• À partir du moment où la tâche s'exécute, vous verrez apparaitre dans le répertoire où se trouve le script différents fichiers, se terminant par l'extension .ht2 qui à la fin de l'exécution constitueront votre index. Lors de l'appel de HISAT2 pour les tâches d'alignement, on utilisera comme source de l'index /shares/data/indexes/hisat2/gencode/r49/h_sapiens_gencode_r49_exons_ss sans l'extension .ht2.

• Vous pouvez télécharger tous ces fichiers sur votre serveur local; exécuter HISAT2 sur cet index demande moins de 16Gb de mémoire vive, ce qui ne devrait pas être un problème à avoir de nos jours HISAT2 s'exécute sans problème sur un RPi 5 doté de 16Gb de mémoire vive. Pour la suite des tâches à faire, on assumera que les fichiers de l'index sont sous /shares/data/indexes/hisat2/gencode/r49.

• Dans le cas de STAR, on utilise une fonction du programme lui-même sans le besoin d'appeler une application auxiliaire. Cependant, il demande à la fois beaucoup de ressources et de temps pour construire l'index (curieuse impression de déjà vu…) et donc. nous aurons besoin d'une grappe de calcul et d'un script à exécuter via SLURM.

• On assume encore que la position des index en général est sous /shares/data/indexes (adaptez selon votre situation) et on doit créer l'arborescence comme pour HISAT2. Un point de différence: on doit aussi créer un répertoire dans lequel seront déposés les fichiers constituant l'index.

% cd /shares/data/indexes
% mkdir star
% mkdir star/gencode
% mkdir star/gencode/r49
% mkdir star/gencode/r49/h_sapiens_r49_index

• Avec nano, on écrit le script nécessaire; disons que nous l'appellerons star_index-builder.sh:

% cd /shares/data/indexes/star/gencode/r49
% nano star_index-builder.sh

• Écrire dans ce fichier le texte qui suit, en adaptant évidemment selon vos circonstances:

#!/bin/bash -l
# 
#SBATCH --account=votre_compte
#SBATCH --job-name=star-build_%j
#SBATCH --output=star-build_%j.out
#SBATCH --error=star-build_%j.err
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=4
#SBATCH --mem=64G
#SBATCH --time=04:00:00
#SBATCH --mail-user=votre@adresse.courriel
#SBATCH --mail-type=ALL
#
# On assume que l'application STAR  est sur le $PATH
#
# Sur Rorqual, il faut charger l'application; retirez le
# dièse de la prochaine ligne pour que le script fonctionne :-)
#module load star

STAR \
--runThreadN 4 \
--runMode genomeGenerate \
--genomeDir ./h_sapiens_r49_index \
--genomeFastaFiles /shares/data/annotations_genome/gencode/r49/GRCh38.primary_assembly.genome.fa \
--sjdbGTFfile /shares/data/annotations_genome/gencode/r49/gencode.v49.annotation.gtf \
--sjdbOverhang 99

• Comment interpréter les paramètres demandés à STAR?
  • –runThreadN 4: spécifie le nombre de coeurs de calcul à utiliser pour exécuter la commande. La valeur doit évidemment être la même que pour le paramètre –mem inscrit dans l'entête du script.
  • –runMode genomeGenerate: la fonction de création d'index de STAR
  • –genomeDir ./h_sapiens_r49_index: le répertoire où seront écrit les divers fichiers constituant l'index
  • –genomeFastaFiles /shares/data/annotations_genome/gencode/r49/GRCh38.primary_assembly.genome.fa: le fichier, en format FASTA, du génome contre lequel faire les alignements
  • –sjdbGTFfile /shares/data/annotations_genome/gencode/r49/gencode.v49.annotation.gtf: le fichier, en format GTF, contenant toutes les annotations disponibles sur les gènes du même génome. Les noms des chromosomes doivent être identiques entre les deux fichiers.
  • –sjdbOverhang 99: correspond à la longueur des lectures contenues dans les fichiers FASTQ à aligner moins 1. Cette valeur correspond à la longueur de génome sur lequel une jonction exon-intron contenue dans le fichier d'annotation doit s'aligner.

• Soumettez votre script à sbatch:

% sbatch star_index-builder.sh

• Après complétion de la construction, vous devriez voir une dernière ligne dans le fichier se terminant en .out contenant le texte suivant: ….. finished successfully. Vous être fin prêt pour la suite.

• Salmon est un logiciel capable soit de faire des alignements directement à partir des fichiers FASTQ ou bien des fichiers BAM alignés via HISAT2 ou STAR et d'en faire l'énumération des transcrits présents. Il est vraiment très rapide car il ne cherche pas à trouver de nouveaux transcrits ou bien de nouveaux sites dépistage. Pour faire son travail avec les fichiers FASTQ, il a besoin de construire un index utilisant en combinaison le fichier du génome, le fichier contenant toutes les séquences des transcrits ainsi que le fichier GTF contenant les annotations.

• C'est un travail en deux temps… En premier, il faut créer des infos afin de diminuer notre risque de faux positifs en créant un ensemble de séquences decoys, permettant de limiter les alignements sur des séquences répétitives présentes dans le génome mais non-représentatives de gènes réels. Pour réussir ça, il nous faut utiliser un script, generateDecoyTranscriptome.sh qui appelle une application appelée MashMap pour réussir ça

# On crée l'arborescence nécessaire sous /shares/data
% cd /shares/data/indexes
#
# Remarquez que nous reprenons la structure logique
# de l'organisation des fichiers retrouvés sous
# /shares/data/annotations
#
% mkdir salmon
% mkdir salmon/gencode
% mkdir salmon/gencode/r49 
% mkdir salmon/gencode/r49/interim_files
% cd salmon/gencode/r49/interim_files

• Il nous faut obtenir le script generateDecoyTranscriptome.sh en téléchargeant le package SalmonTools qui se trouve sur Github:

% git clone https://github.com/COMBINE-lab/SalmonTools.git
# Le script se trouve dans le répertoire ./SalmonTools/scripts; il faut 
# le rendre executable
% chmod u+x ./SalmonTools/scripts/generateDecoyTranscriptome.sh

• L'exécution de ce script demandera à la fois beaucoup de mémoire et beaucoup de temps… L'usage d'une grappe de calcul et de SLURM est donc nécessaire Il naus faut donc écrire un nouveau script pour le soumettre via sbatch:

# On retourne au niveau supérieur
% cd ..
# On utilise nano pour la composition du script
% nano salmon_decoy-builder.sh

• Mettre le texte suivant dans le script; rappelez-vous que vous devrez l'adapter à votre environnement de travail:

#!/bin/bash -l
#
#SBATCH --account=votre_compte
#SBATCH --job-name=salmon_decoy-builder_%j
#SBATCH --output=salmon_decoy-builder_%j.out
#SBATCH --error=salmon_decoy-builder_%j.err
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=1
#SBATCH --mem=384G
#SBATCH --time=08:00:00
#SBATCH --mail-user=votre@adresse.courriel
#SBATCH --mail-type=ALL
#
# On assume que bedtools et mashmap sont sur le $PATH
# Sinon, dans une. grappe de calcul, retirer 
# le # devant les 2 prochaines lignes
#module load bedtools
#module load mashmap

./SalmonTools/scripts/generateDecoyTranscriptome.sh \
-a /shares/data/annotations_genome/gencode/r49/gencode.v49.annotation.gtf \
-g /shares/data/annotations_genome/gencode/r49/GRCh38.primary_assembly.genome.fa \
-t /shares/data/annotations_genome/gencode/r49/gencode.v49.transcripts.fa \
-o ./interim_files

• Il faut maintenant soumettre le script à la queue de traitement comme ci-dessus:

% sbatch salmon_decoy-builder.sh

• Une fois le script exécuté avec succès, vous devriez avoir dans le répertoire interim_files les fichiers suivants: decoys.txt,exons.bed,gentrome.fa et reference.masked.genome.fa, qui seront utilisés à la prochaine étape.

• Ce qui est fascinant, c'est que la création de l'index pour salmon demande pas mal moins de ressources que la création des fichiers decoys: ça se fait sur un Raspberry Pi 5 avec 8Gb de mémoire vive . On peut donc le faire localement en commençant par téléchargeant les fichiers decoys sur notre serveur local.

#
# On assume que nous sommes localement dans le répertoire
# /shares/data/indexes/salmon/gencode/r49
#
# Ajuster le lien distant en consequence de votre propre 
# environnement 
#
% scp -r mon_compte@grappe_calcul:salmon/gencode/r49/interim_files .

• En utilisant les fichiers décrits à la section précédente et en restant sous /shares/data/indexes/gencode/r49, on exécute la fonction de création des index de salmon via un script pour soumission via sbatch. On crée donc un nouveau script (appelons le salmon_index-builder.sh) via nano:

#!/bin/bash -l
#
# Parametres sbatch pour une exécution locale de sbatch sur
# un cluster SuperClafoutis ou une instance locale similaire
#
#SBATCH --job-name=salmon_index-builder_%j
#SBATCH --output=salmon_index-builder_%j.out
#SBATCH --error=salmon_index-builder_%j.err
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=4
#SBATCH --mem=7G
#
# On assume que salmon est sur le $PATH
# Sinon, dans une. grappe de calcul, retirer 
# le # devant la prochaine ligne
#module load salmon

salmon index \
-p 4 \
-t ./interim_files/gentrome.fa \
-i h_sapiens_r49_index \
-k 31
--decoys ./interim_files/decoys.txt

• Vous devriez obtenir un répertoire appelé h_sapiens_r49_index contenant une série de fichiers constituant l'index en question. Lors de l'exécution à venir, on spécifie le répertoire comme site de l'index à utiliser.