Exercice 2: alignement BWA

1. Rechercher et charger BWA (v0.7.17) et samtools (v1.18) en utilisant la commande module.

   Solution

   search_module bwa
   module load bioinfo/bwa/0.7.17
   search_module samtools
   module load bioinfo/samtools/1.18
   

2. Afficher l'aide des commandes bwa et bwa index.

3. Créer l'index bwa de la référence.

   Solution

   bwa index Gallus_gallus.Gallus_gallus-5.0.dna.toplevel_chr25-26.fa
   

4. Explorer l'aide de la commande bwa mem.

5. Lancer l'alignement en ajoutant les ReadGroup, convertir en BAM et trier.
   Pour l'ajout du readgroup, voir l'option -R.
   Pour convertir en bam, rediriger la sortie vers | samtools - > echantillon.bam.

   Vous pouvez explorer la FAQ > Bioinfo tips > "How to generate an sarray command file with perl ... " ou utiliser la solution ci dessous.

   Solution fastidieuse

   Créer à la main un fichier 1_bwamem.jobs contenant toutes les commandes bwa :

   bwa mem -R "@RG\tID:1\tSM:SRR7062654\tPL:illumina\tLB:SRR7062654\tPU:1" -t4 Gallus_gallus.Gallus_gallus-5.0.dna.toplevel_chr25-26.fa SRR7062654_R1.fastq.gz SRR7062654_R2.fastq.gz | samtools sort - > SRR7062654.bam
   bwa mem -R "@RG\tID:1\tSM:SRR7062655\tPL:illumina\tLB:SRR7062655\tPU:1" -t4 Gallus_gallus.Gallus_gallus-5.0.dna.toplevel_chr25-26.fa SRR7062655_R1.fastq.gz SRR7062655_R2.fastq.gz | samtools sort - > SRR7062655.bam
   

   Lancer l'exécution sur le cluster :

   sarray -J 1_bwamem -e LOGS/%x_%j.err -o LOGS/%x_%j.out -c 5 1_bwamem.jobs
   

   Solution compliquée

   Bash

   for R1 in *_R1.fastq.gz; do echo "bwa mem -R \"@RG\tID:1\tSM:${R1%_*}\tPL:illumina\tLB:${R1%_*}\tPU:1\" -t4 Gallus_gallus.Gallus_gallus-5.0.dna.toplevel_chr25-26.fa $R1 ${R1/_R1.f/_R2.f} | samtools sort - > ${R1%_*}.bam"; done > 1_bwamem.jobs
   

   💡

   • ${R1%_*} : affiche le contenu de la variable R1 en supprimant en fin de chaine le motif défini après le caractère %. Ici _* signifie n'importe quel caractère après _. La chaîne SRR7062654_R1.fastq.gz devient donc SRR7062654.
   • ${R1/_R1.f/_R2.f} : affiche le contenu de la variable R1 en remplaçant le motif défini entre les // par la chaîne de substitution définie après le second /. Ici, nous cherchons à remplaçer le motif _R1.f par la chaîne _R2.f.

   Perl

   \ls -1 *_R1.fastq.gz | perl -lne '$id=$_; $id=~s/\_.*//; $out=$_; $out=~s/\_R1.*//; $r2=$_; $r2=~s/\_R1\.f/\_R2\.f/; print "bwa mem -R \"\@RG\\tID:1\\tSM:$id\\tPL:illumina\\tLB:$id\\tPU:1\" -t4 Gallus_gallus.Gallus_gallus-5.0.dna.toplevel_chr25-26.fa $_ $r2 | samtools sort - > $out.bam"' > 1_bwamem.jobs
   

   💡

   • \ls -1 *_R1.fastq.gz : liste tous les read1 du répertoire (le \ permet d'utiliser la commande ls native)
   • perl -lne : avec un retour à la ligne comme séparateur (-l) et pour chaque ligne en entrée (-n), j'exécute le bout de code perl '...' (-e)
   • $id=$_; : initialisation de la variable id avec le nom du fichier fastq
   • $id=~s/\_.*//;: subtitution de tout ce qui suit le premier '_' par rien dans la variable id
   • $out=$_;: initialisation de la variable out avec le nom du fichier fastq
   • $r2=$_;: initialisation de la variable r2 avec le nom du fichier fastq
   • $r2=~s/\_R1\.f/\_R2\.f/;: subtitution de la chaine _R1 par _R2 dans la variable r2
   • print "ligne de commande en utilisant les variables précédement créées"

   sarray -J 1_bwamem -e LOGS/%x_%j.err -o LOGS/%x_%j.out -c 5 1_bwamem.jobs