Les différentes méthodes de séquençage à haut débit ne sont pas totalement infaillibles dans leur capacité à faire le dit séquençage! Elles sont toutes basées sur des processus qui à divers moments de la méthode, peuvent engendrer des problèmes comme des échantillons de qualité discutable, des erreurs de lecture, des séquences partielles dans l'ensemble d'un fichier (rappel: chaque fichier contient des millions de séquences), etc. Il nous faudra dans une étape subséquente retirer toute donnée de séquence qui pourra perturber les reste de notre analyse; pour établir les critères de filtration, il nous faut regarder le profil de nos données brutes. Quoi utiliser? Dans la communauté, dans le moment, les deux applications les plus utilisées sont FastQC et fastP, qui toutes les deux vont créer des rapports sous la forme de pages web qui permettront de vous donner un aperçu de vos données et qui vous permettront de déterminer les meilleurs critères de filtration dans l'étape suivante.

Pour une première fois, on a une situation où l'on doit faire à répétition la même opération sur tous les fichiers FASTQ que nous avons téléchargé. On pourrait faire 32 fois la même opération mais, soyons sérieux, c'est à ça que sert une grappe de calcul et Python Pourquoi ne pas utiliser Python pour écrire tous les fichiers bash à soumettre à SLURM? À la base, il faudrait encoder statiquement tous les paramètres à soumettre dans notre script, ce qui le rend peu flexible… Introduisons un autre outil: un fichier de configuration à lire par notre script. Comment écrire ce fichier de paramètres? En utilisant un fichier texte écrit dans un format structuré qui pourra être lu par le script pour construire une structure de données directement utilisable. Un tel format est le format YAML, très flexible et facilement interprété par Python.

• Dans un premier temps, construisons un fichier texte en format YAML pour y mettre nos paramètres:

#
# qc_airway.yaml - Un fichier de paramètres à soumettre aux scripts
#                  de QC utilisant FASTQC ou FastP
#
# Est-ce qu'on veut créer les scripts sans les soumettre?
# Mettre à false pour lancer les tâches sur SLURM
#
dry_run: true
#
# Ou se trouve les données brutes?
#
path2raw: /shares/data/rnaseq/airway
#
# Ou écrire les données de QC obtenues par
# FASTQC ou FastP?
#
path2qc: /home/username/analysis/rnaseq/airway/1.fastq_raw_qc
#
# Conditions utilisées sur les lignées cellulaires
# Doit être les mêmes que celles du répertoire  
# /shares/data/rnaseq/airway
#
condition:
  - alb
  - alb_dexa
  - dexa
  - dmso
#
# Identification des lignées cellulaires
#
cellline:
    - N052611
    - N061011
    - N080611
    - N61311
#
# Quel outil utilisé? Ce paramètre est obligatoire
# pour la suite
#
tool: fastqc
#tool: fastp
#
# Paramètres spécifiques pour chaque outil à utiliser
# pour chaque soumission. NE PAS METTRE DE PARAMÉTRES 
# SPÉCIFIQUES! Ils seront construit par le script lui-même 
# Assurez-vous d'utiliser les bons!! Vérifier la doc de
# chaque outil au besoin. Si possible, utiliser les mêmes noms
# que ceux attendu par le programme.
#
# Dans l'exemple, FASTQC ne demande que le nombre de CPU pour
# faire l'analyse
#
toolParams:
  #
  # Combien de processeurs utilisés?
  #
  threads: 4

• Ce fichier peut être sauvegarder à l'endroit de votre choix car le script Python qui le lira créera automatiquement les chemins à suivre pour la lecture et l'écriture des fichiers.
• Passons maintenant au script:

#
# fastq_qc_processor.py - Un programme simple pour automatiser l'analyse
#                         QC des séquences téléchargées en utilisant FastQC
#                         sur une grappe de calcul via SLURM.
#
# Ici, on introduit le concept de fichier de paramètres écrit en format
# YAML pour rendre le code indépendant des variables; les variables sont 
# définies dans le fichier YAML. 
#
# Il faut lire chacun des fichiers FASTQ pour créer un rapport FastQC dans
# notre répertoire perso pour l'analyse. Le script est aussi agnostique
# sur l'outil utilisé: FastQC ou fastp. Il faut simplement le spécifier
# dans le fichier YAML
#
# Le code de ce script est hyper-documenté pour la formation ;-)
# 
# 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...
from datetime import datetime
import os
import subprocess
import sys
import yaml
#
# Permet de capturer la position du fichier YAML dans le système
# de fichiers:
#
# % python ./fastq_qc_processor.py qc_airway.yaml 
#
# sys.argv[1] devrait être qc_airway.yaml :-)
# 
with open(sys.argv[1], "r") as flux:
    #
    # Une boucle try...except permet de tenter d'accomplir
    # une tâche en capturant les messages d'erreur en cas
    # d'échec. 
    #
    try:
        #
        # Lire le fichier pour en mettre le contenu dans 
        # un dictionnaire
        #
        data = yaml.safe_load(flux)
    #
    # Une bonne pratique: inclure les messages d'erreur pour
    # info de débuggage
    #
    except yaml.YAMLError as err:
        print(f"Erreur de lecture du fichier YAML: {err}")

path2raw = data["path2raw"]
path2qc =  data["path2qc"]
#
# Liste des répertoires en relation avec les traitements
#
conditions = data["condition"]
#
# Liste des lignées cellulaires du projet
#
cellLines = data["cellline"]
#
# Comme notre script utilise FastQC ou fastP, il faut que la 
# clé "tool" soit dans data.
#
# Ici, on a un exemple de vérification des infos avant de continuer.
# À vous de voir si vous voulez le faire pour les autres variables en
# amont :-)
#
try:
  tool = data["tool"]
except:
  print("L'outil à utiliser n'est pas spécifié; pas de champs tool dans le fichier YAML...") 
  print("Je quitte!!")
  sys.exit()
#
# Si toolParams est présent, les valeurs seront utilisées à la place des valeurs par défaut.
#
params = data["toolParams"]
#
# Quelle logique à employer? On veut reproduire la structure qui se trouve dans
# /shares/data/rnaseq/airway en détaillant étape par étape les procédures
# faites sur les fichiers. Une idée:
# 
# - On a construit /home/user/analyse/rnaseq/airway;
# - On aura comme première étape le contrôle de la qualité des fichiers FASTQ
#   avec FastQC ou fastp; appelons ça 1.fastq_raw_qc donc créons ce répertoire;
# - Une bonne pratique: identifier le plus clairement possible les répertoires
#   alors disons que sous 1.fastq_raw_qc, on créera fastqc;
# - Chaque exécution doit être vue comme une expérience unique; il faut donc
#   trouver une façon de rendre ça unique aussi dans le système de fichiers;
# - On recréera la structure des conditions et sous chaque condition, on crée
#   une structure des lignées cellulaires. 
#
if not os.path.exists(f"{path2qc}"):
  os.mkdir(f"{path2qc}")
  print(f"Le répertoire {path2qc} n'existe pas... On le crée!")
else:
  print(f"Le répertoire {path2qc} existe... On continue!")
#
# Création du répertoire spécifique à l'outil choisi
#
if not os.path.exists(f"{path2qc}/{tool}"):
  os.mkdir(f"{path2qc}/{tool}")
else:
  print(f"Le répertoire {path2qc}/{tool} existe... On continue!")
#
# Création d'un ID unique basé sur le moment présent
#
uniq = datetime.now().strftime("%Y_%m_%d_%H_%M_%S")
#
# Variable pour simplifier notre vie
#
path2rez = f"{path2qc}/{tool}/{uniq}"
#
# Répertoire pour mettre les données QC
#
os.mkdir(path2rez)
#
# Création des répertoires selon traitement et lignée cellulaire
#
# parce que l'on sait que le répertoire path2rez n'existait pas, on
# sait que les répertoires internes n'existeant pas non plus ;-)
# Pas besoin de vérifier leur existence
#
for i in conditions:
  os.mkdir(path2rez+"/"+i)
  for j in cellLines:
    path2rezS = path2rez+"/"+i+"/"+j
    os.mkdir(path2rezS)
    #
    # On a maintenant un répertoire pour y mettre nos analyses QC
    # par lignée cellulaire en y mettant les répertoires pour
    # chaque fichier de cette lignée pour cette condition.
    #
    os.chdir(path2rezS)
    #
    # Prefixe de nos répertoires et de nos fichiers 
    # de données brutes
    #
    nameF = i+"_"+j
    # 
    # 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={tool}_{nameF}_%j \n"
    batch+= f"#SBATCH --output={path2rezS}/{tool}_{nameF}_%j.out \n"
    if "threads" in params:
      batch+= "#SBATCH --ntasks=1 \n"
      batch+= f"#SBATCH --cpus-per-task={params["threads"]}\n"
    #
    # On assume que FastQC et/ou fastp sont sur le $PATH
    #
    # On sait que nos fichiers fastq identifie chaque read par _1 ou _2
    #
    for k in range(1,3):
      #
      # Comment faire si on utilise FASTQC
      # On veut une tâche par fichier fastq.gz à traiter.
      #
      if tool=="fastqc":
        k = str(k)
        #
        # Il faut créer le repertoire de stockage en premier car 
        # FastQC ne le fera pas...
        #
        os.mkdir(nameF+"_"+k)
        #
        # Sur Rorqual, il faut charger l'application; retirez le
        # dièse de la ligne correspondante pour que le script fonctionne
        #batch+="module load fastqc \n"
        batch+=f"fastqc  --noextract --threads {params["threads"]} -o ./{nameF}_{k} {path2raw}/{i}/{nameF}_{k}.fastq.gz \n"
      #
      # Comment faire si on utilise fastp
      #
      # On ne veut pas la fonction de filtration de fastp alors on ne spécifie
      # pas de sortie (-o/-O)
      #
      elif tool=="fastp":
        #
        # Sur Rorqual, il faut charger l'application; retirez le
        # dièse de la ligne correspondante pour que le script fonctionne
        #batch+="module load fastp \n"
        batch+=f"fastp  --thread {params["threads"]} -j ./{nameF}_{k}.json -h ./{nameF}_{k}.html -i {path2raw}/{i}/{nameF}_{k}.fastq.gz \n"
      # 
      # Écriture du script bash à soumettre à SLURM
      #
      scriptB = f"{i}_{j}_{k}.sh"
      with open(scriptB, "w", encoding="utf-8") as scriptL:
        scriptL.write(batch)
      #
      # Soumission du script bash
      # Utilisation du paramètre dry_run: si false, alors soumettre la tache.
      #
      if not data["dry_run"]:
        subprocess.run(["sbatch",scriptB])

• Ici, le nom du fichier est fastq_qc_processor.py mais ça peut-être ce que vous voulez Si tout va bien, vous aurez l'envoi de 32 tâches distinctes à votre grappe de calcul qui s'exécuteront successivement et que vous pouvez suivre grâce à la commande squeue.
• Plus à venir…