Mise en pratique des méthodes d'analyse de l'expression génique grâce aux données utilisées par le package airway: obtenir les fichiers de séquence bruts du projet

Il faut maintenant se concentrer sur la méthodologie à employer pour optimiser l'automatisation des étapes. En construisant une hiérarchie unifiée pour les fichiers source, on pourra réutiliser cette hiérarchie dans les étapes de l'analyse. Dans un premier temps, il nous faut obtenir les fichiers source FASTQ bruts, qu'ils proviennent d'une plateforme de séquençage ou bien d'une base de données. Comment les organisés de manière optimale? Dans le cas précis de notre projet, on a diverses lignées cellulaires humaines qui nous servent à observer les changements de l'expression génique selon divers traitements. Il est par conséquent logique d'assembler les fichiers par traitement plutôt que par lignée.

De plus, plusieurs personnes dans un projet pourraient avoir un besoin de lire les informations contenues dans ces fichiers. Il est alors nécessaire de les mettre dans un endroit du système de fichiers accessible à tous mais en les protégeant en modifiant les permissions afin qu'ils soient accessibles en lecture mais ne peuvent être éditer. Un autre point: on ne veut pas que personne d'autre vienne y mettre d'autres fichiers par accident ou autre

Est-ce que vos fichiers proviennent d'une plateforme de séquençage ou bien d'un entrepôt de données comme Séquence Read Archives? La plupart du temps, ça importe peu car il faudra qu'ils arrivent à votre plateforme via l'Internet. De plus, il faudra télécharger beaucoup de fichiers, jusqu'à 2 par échantillon si il a été séquencé paired-end donc ça fait beaucoup de tâches distinctes de téléchargement… C'est ici que l'automatisation des tâches via un script Python, seul ou en exécution sur une grappe de calcul, sera d'un usage fort utile Nous présenterons deux méthodes: un script pour l'ensemble des tâches de téléchargement ou bien un script qui créera autant de tâches que d'échantillons pour soumettre à une. grappe de calcul.

Pour commencer, il nous faut une liste des fichiers à télécharger à partir d'un serveur distant. Pour notre projet, les auteurs nous disent que les infos se trouvent sur Gene Expression Omnibus sous l'identificateur GSE52778. En inspectant cette page, on voit que les données de séquences brutes se trouvent dans le site Sequence Read Archive avec l'identificateur SRP033351. On y est presque! En allant sur le site de la page de sélection des données du SRA, on fait une recherche avec cet identificateur et on tombe sur la page contenant l'information que nous cherchons!! Dans le panneau Select, on voit deux boutons: Metadata ou bien Accession List; cliquer sur Metadata pour obtenir un fichier CSV plus informatif, SraRunTable.csv, qui nous servira pour la suite. Allez mettre ce fichier dans l'arborescence de votre serveur; par ex., créer un répertoire appelé z.misc.files sous /shares/data/rnaseq/airways.

• La suite: plus à venir…

• Une idée sur la marche à suivre: nous avons construit les fichiers index sous /shares/data, un répertoire accessible pour tous dans le système de fichiers. alors construisons un nouveau répertoire pour y mettre tous les projets de transcriptomique basés sur le séquençage à haut débit:

% mkdir /shares/data/rnaseq
% mkdir /shares/data/rnaseq/airway

• Nous savons (en lisant le papier associé) que les lignées ont subit l'un des traitements suivant: DMSO comme contrôle, albuterol, dexaméthasone ou albuterol+dexaméthasone. On crée alors une hiérarchie avec cette information:

% cd /shares/data/rnaseq/airway
% mkdir dmso
% mkdir alb
% mkdir dexa
% mkdir alb_dexa
% mkdir z.misc.files

• Après les opérations de téléchargement des fichiers, on reviendra retoucher les permissions. En principe, on veut pouvoir voir le contenu des répertoires et lire leur contenu ainsi que lire les fichiers qui s'y trouvent. Comme l'usager est le propriétaire des fichiers, il peut se retirer des permissions d'écriture des répertoires et fichiers pour éviter d'en écraser le contenu et d'y mettre des fichiers pas rapport. Évidemment, il faut voir comment les usagers du groupe d'usagers auquel appartient l'usager créateur et le reste des usagers du système peuvent accéder aux données.

• Dans le cas de l'analyse, le simple fait que chacun peut la faire selon ses propres hypothèses, on va plutôt mettre le tout dans notre répertoire $HOME:

# On reste générique car le repertoire peut contenir
# autre chose que des projets de transcriptomique...
% mkdir $HOME/analysis
% mkdir $HOME/analysis/rnaseq
% cp -r /shares/data/rnaseq/airway $HOME/analysis/rnaseq

• Premier exemple: téléchargement directement dans le serveur stockant les fichiers, sans recours à un système de gestion de tâches tel que SLURM. Noter que le code est fortement documenté à des fins de formation

#
# fetch_fastq_from_csv.py - Un programme simple pour télécharger des
# fichiers FASTQ directement du Sequence Read Archive du NCBI en
# utilisant les outils de la suite SRA Toolkit.
#
# Un point important à considérer quand on fait ce travail, c'est la
# structure des données contenues dans le fichier. Dans notre cas 
# particulier, on va observer que les noms de traitement sont écrit
# au long: Untreated, Albuterol, Dexamethasone, Albuterol_Dexamethasone
# donc différents de notre structure proposé, alors on devra composer 
# avec ça.
#
# Le code de ce script est hyper-documenté pour la formation ;-)
# 
# NOTE: ce script prendra beaucoup de temps pour rouler vu la quantité de
#       données à récolter. Si nécessaire, utiliser la commande nohup devant
#       pour éviter de se faire déconnecter... Ou bien utiliser SLURM si
#       vous avez une grappe Super-Clafoutis ou bien votre grappe permet
#       l'accès à l'Internet pour les noeuds de traitement. Les grappes
#       de Calcul-Québec ne permettent pas l'accès à leurs noeuds de travail...
#
# Inclure les libraries qui seront nécessaires...
import csv
import subprocess
import time
#
# Pointer cette variable vers la source de votre
# fichier SraRunTable.csv 
#
csvFile = "./SraRunTable.csv"
with open(csvFile, mode="r", encoding="utf-8") as file:
    #
    # Lire le fichier pour retourner son contenu sous la forme
    # d'une collection de dictionnaires via la méthode DictReader. 
    # C'est important car avec un dictionnaire, on peut accéder 
    # aux clés, c.-à-d., aux noms de colonnes :-)
    # 
    # Il faut ensuite transformer le contenu en une liste utilisable
    # car la méthode csv.DicReader retourne en fait un objet au contenu
    # inaccessible...  
    #
    sampleList = list(csv.DictReader(file))
#
# Imprimons simplement le premier item de la liste pour voir les noms de 
# colonne pour faire notre filtration et téléchargement.
# 
# Retirez le dièse pour voir le resultat:
#print(sampleList[0])
#
# Vous devriez y trouver des items dans la liste suivante:
#  Run: ID unique dans SRA
#  cell_line: la source du matériel séquencé
#  treatment: quel traitement a été appliqué sur ce matériel
#
# Avec ces champs, on pourra faire le travail :-)
#
# Chemin pour la sauvegarde; adapter selon votre environnement
dlPath = "/shares/data/rnaseq/airway"
#
# Liste des répertoires en relation avec les traitements
#
treatment = [{"Albuterol":"alb"},{"Dexamethasone":"dexa"},{"Untreated":"dmso"},{"Albuterol_Dexamethasone":"alb_dexa"}]
#
# Parcourons la liste des échantillons du projet
#
for i in sampleList:
  for j in treatment: 
    if i["treatment"] in j:
       treat = j.get(i["treatment"])
       saveF = dlPath+"/"+treat
       uid = i["Run"]
       cell = i["cell_line"]
       nameF = treat+"_"+cell+"_"+uid
       #*********************
       # Téléchargement
       #*********************
       # Ça ne fonctionnera que si les applications est sur $PATH
       # L'utilisation de la méthode run assure que les commandes seront
       # effectuée une à la fois.
       # 
       # On utilise la procédure suggérée par le NCBI: passer par la 
       # commande prefetch pour télécharger le fichier binaire compact
       # pour l'échantillon demandé pour ensuite prendre fasterq-dump
       # pour en extraire les fichiers FASTQ.
       # 
       # C'est un méthode plus rapide car prefetch vérifie si le
       # téléchargement à échouer: il le continue s'il a été arrêté et
       # passe par dessus s'il a été complété.
       #  
       # À partir de ce point, l'exécution prendra du temps, possiblement
       # beaucoup de temps... On parle de 16 fichiers à ~1-2Gb chacun; utiliser
       # SLURM si possible ;-)
       #
       subprocess.run(["prefetch","-O",f"{saveF}",uid])
       #
       # Si nécessaire, espaçons nos requetes pour rester dans les bonnes 
       # graces du NCBI ;-)
       #
       # Retirez le dièse pour mettre un délai dans vos requetes de 
       # téléchargements. Pas forcément nécessaire si le script est exécuté
       # sur un seul ordi car les étapes qui suivront vont prendre du temps.
       #
       #time.sleep(120)
       #
       #*******************************
       # Extraction des fichiers FASTQ
       #*******************************
       #
       # Le reste du script ne demande pas d'accès Internet donc pourrait
       # être délégué à une grappe via SLURM; à suivre... 
       #
       # L'exemple fonctionne sur un RPi 4/5 alors on utilise 
       # 4 fils à la fois. Si votre serveur dispose de plus, pourquoi pas 
       # en mettre plus ;-)
       #
       subprocess.run(["fasterq-dump","--split-3","--threads","4","--outdir",saveF,"--outfile",nameF,f"{saveF}/{uid}"])
       #
       # Compressons maintenant les fichiers sur notre stockage
       # Les outils à venir sont tous capables de lire les fichiers 
       # comprimés. 
       #
       # On utilise le progamme bgzip, installé sur Impilo, pour faire
       # une compression multi-fils et accélérer le processus.
       #
       # Comme on fait deux fois la même procédure (sur chaque fichier), 
       # on utilise une boucle. La commande range() part à 0 alors il faut
       # ajouter 1 pour avoir le nom réel du fichier final. 
       #
       for k in range(2):
         k = k+1
         #
         # Pas besoin d'effacer le fichier .fastq; bgzip le fera après la 
         # compression
         #
         subprocess.run(["bgzip","--threads","4",f"{saveF}/{nameF}_{k}.fastq"])

• À la fin de l'exécution, vous devriez avoir 32 fichiers dont l'extension se termine par _x.fastq.gz avec x=1 ou x=2. Les fichiers devraient se trouver dans le bons répertoire avec un nom sous la forme traitement_ligneeCellulaire_SRRxyz, ce qui nous facilitera la vie pour la suite

• De plus, vous devriez avoir pour certains échantillons un autre fichier se terminant simplement par .fastq ou sans extension mais avec le même nom; ce fichier existe car certains échantillons ont généré plusieurs lectures qui n'ont pas d'équivalent entre les deux fichiers pour toute sorte de raison: séquences techniques comme les codes barres ou bien des séquences d'amorces ou autres séquences sans correspondance . On ne veut pas utiliser ces fichiers donc le paramètre –split-3 nous permet de les retirer en amont plutôt que de devoir le faire post-alignement.

• Comme nous avons une liste de tâches répétitives, pourquoi ne pas utiliser une grappe de calcul pour accélérer le processus? On peut mais il faudra avoir des détails sur la grappe à utiliser. Par exemple, sur la grappe rorqual.calculquebec.ca, aucun noeud de traitement n'a accès à l'internet alors il faudra que le script travaille en deux temps: le noeud de connexion devra faire le téléchargement du fichier .sra et le reste (utilisation de fasterq-dump et compression des fichiers FASTQ) pourra se faire sur les noeuds de traitement car aucun recours à une connexion hors du réseau local n'est nécessaire.

• Le script qui suit devrait fonctionner sur une grappe comme Rorqual ou bien sur une grappe Super-Clafoutis.

#
# fetch_fastq_from_csv_slurm_1.py - Un programme simple pour télécharger des
# fichiers FASTQ directement du Sequence Read Archive du NCBI en
# utilisant les outils de la suite SRA Toolkit et une gappe de calcul
# via SLURM.
#
# Un point important à considérer quand on fait ce travail, c'est la
# structure des données contenues dans le fichier. Dans notre cas 
# particulier, on va observer que les noms de traitement sont écrit
# au long: Untreated, Albuterol, Dexamethasone, Albuterol_Dexamethasone
# donc différents de notre structure proposé, alors on devra composer 
# avec ça.
#
# Le code de ce script est hyper-documenté pour la formation ;-)
# 
# NOTE: ce script prendra tout de même un certain temps pour rouler 
#       vu la quantité de données à télécharger. Si nécessaire, utiliser 
#       la commande nohup devant pour éviter de se faire déconnecter...
#
#       La vitesse de réalisation du script dépendra du nombre de noeuds
#       de traitement: plus il y en a, plus ça sera rapide :-)
#
# Inclure les libraries qui seront nécessaires...
import csv
import subprocess
#
# Pointer cette variable vers la source de votre
# fichier SraRunTable.csv 
#
csvFile = "./SraRunTable.csv"
with open(csvFile, mode="r", encoding="utf-8") as file:
    #
    # Lire le fichier pour retourner son contenu sous la forme
    # d'une collection de dictionnaires via la méthode DictReader. 
    # C'est important car avec un dictionnaire, on peut accéder 
    # aux clés, c.-à-d., aux noms de colonnes :-)
    # 
    # Il faut ensuite transformer le contenu en une liste utilisable
    # car la méthode csv.DicReader retourne en fait un objet au contenu
    # inaccessible...  
    #
    sampleList = list(csv.DictReader(file))
#
# Imprimons simplement le premier item de la liste pour voir les noms de 
# colonne pour faire notre filtration et téléchargement.
# 
# Retirez le dièse pour voir le resultat:
#print(sampleList[0])
#
# Vous devriez y trouver des items dans la liste suivante:
#  Run: ID unique dans SRA
#  cell_line: la source du matériel séquencé
#  treatment: quel traitement a été appliqué sur ce matériel
#
# Avec ces champs, on pourra faire le travail :-)
#
# Chemin pour la sauvegarde; adapter selon votre environnement
dlPath = "/shares/data/rnaseq/airway"
#
# Liste des répertoires en relation avec les traitements
#
treatment = [{"Albuterol":"alb"},{"Dexamethasone":"dexa"},{"Untreated":"dmso"},{"Albuterol_Dexamethasone":"alb_dexa"}]
#
# Parcourons la liste des échantillons du projet
#
for i in sampleList:
  for j in treatment: 
    if i["treatment"] in j:
       treat = j.get(i["treatment"])
       saveF = dlPath+"/"+treat
       uid = i["Run"]
       cell = i["cell_line"]
       nameF = treat+"_"+cell+"_"+uid
       #*********************
       # Téléchargement
       #*********************
       # Ça ne fonctionnera que si l'application est sur $PATH
       # L'utilisation de la méthode run assure que les commandes seront
       # effectuée une à la fois.
       # 
       # On utilise la procédure suggérée par le NCBI: passer par la 
       # commande prefetch pour télécharger le fichier binaire compact
       # pour l'échantillon demandé pour ensuite prendre fasterq-dump
       # pour en extraire les fichiers FASTQ.
       # 
       # C'est un méthode plus rapide car prefetch vérifie si le
       # téléchargement à échouer: il le continue s'il a été arrêté et
       # passe par dessus s'il a été complété.
       #  
       # À partir de ce point, l'exécution prendra du temps, possiblement
       # beaucoup de temps... On parle de 16 fichiers à ~1-2Gb chacun; utiliser
       # SLURM si possible ;-)
       #
       subprocess.run(["prefetch","-O",f"{saveF}",uid])
       #
       # Si nécessaire, espaçons nos requetes pour rester dans les bonnes 
       # graces du NCBI ;-)
       #
       # Retirez le dièse pour mettre un délai dans vos requetes de 
       # téléchargements. Pas forcément nécessaire si le script est exécuté
       # sur un seul ordi car les étapes qui suivront vont prendre du temps.
       #
       #time.sleep(120)
       #
       #******************************************************
       # Extraction des fichiers FASTQ - Utilisation de SLURM
       #******************************************************
       #
       # Le reste du script ne demande pas d'accès Internet donc 
       # on passe à SLURM. Il ne faut pas oublier que sbatch demande
       # un script bash et non Python, alors on devra créer un fichier
       # pour ensuite le soumettre via sbatch.
       # 
       # Création de l'entête du fichier bash 
       #
       batch = "#!/bin/bash -l \n"
       batch+= "## \n" 
       batch+= "#SBATCH --account=nom_compte \n"
       batch+= f"#SBATCH --job-name=faster-dump_{uid}_%j \n"
       batch+= f"#SBATCH --output={saveF}/faster-dump_{uid}_%j.out \n"
       batch+= "#SBATCH --ntasks=1 \n"
       batch+= "#SBATCH --cpus-per-task=4\n"
       #
       # On assume que les applications du SRA Toolkit sont sur le $PATH
       #
       # Sur Rorqual, il faut charger l'application; retirez le
       # dièse des deux lignes suivantes pour que le script fonctionne
       #batch+="module load sra-toolkit \n"
       #batch+="module load tabix \n"
       
       # L'exemple fonctionne sur un RPi 4/5 alors on utilise 
       # 4 fils à la fois sur une grappe Super-Clafoutis. Si votre serveur 
       # dispose de plus, pourquoi pas en mettre plus ;-)
       #
       # Remarquer le "\\\n": c'est pour mettre un "\" à la fin de chaque 
       # ligne pour que bash sache que la commande continue sur l'autre
       # ligne :-)
       #
       batch+=f"fasterq-dump --split-3 --threads 4 \\\n--outdir {saveF} \\\n--outfile {nameF} \\\n{saveF}/{uid} &&\n"
       #
       # Compressons maintenant les fichiers sur notre stockage
       # Les outils à venir sont tous capables de lire les fichiers 
       # comprimés.
       #
       # On utilisera bgzip car le programme est multithread par défaut
       # donc une accélération de la compression. 
       #
       # Comme on fait deux fois la même procédure (sur chaque fichier), 
       # on utilise une boucle. La commande range() part à 0 alors il faut
       # ajouter 1 pour avoir le nom réel du fichier final. 
       #
       for k in range(2):
         k = k+1
         if k==2:
           batch+=f"bgzip --threads 4 {saveF}/{nameF}_{k}.fastq \\\n"
         else:
           batch+=f"bgzip --threads 4 {saveF}/{nameF}_{k}.fastq &&\\\n"

       scriptB = f"{saveF}/{nameF}.sh"

       with open(scriptB, "w", encoding="utf-8") as scriptL:
         scriptL.write(batch)
       #
       # Envoyer le script à sbatch
       #
       subprocess.run(["sbatch",scriptB])

• Vous devriez voir apparaître vos fichiers SRA dans votre espace de stockage et un ensemble de fichiers en attente d'exécution via la commande squeue. Comme on utilise la méthode run(), on sait que le fichier SRA sera téléchargé avant de mettre en exécution l'extraction/compression des fichiers FASTQ. Une fois les fichiers SRA téléchargés, il ne reste plus qu'à attendre que les scripts bash créés par le script Python s'exécutent .

• Il faut faire le ménage ensuite: vérifier que tout s'est bien passé en regardant les fichiers .out, effacer les fichiers de séquences non-appariés dans certains cas et qui ne se terminent pas par .fastq.gz et effacer les fichiers *.sh car ils ne sont plus nécessaires. On peut toujours garder les répertoires dont le nom commence par SRR: si nécessaire, on peut y retourner pour récupérer les séquences sous de nouveaux paramètres.