Exercice 2: alignement BWA
Rechercher et charger BWA et samtools (commande module).
Solutionsearch_module bwa module load bioinfo/bwa/0.7.17 search_module samtools module load bioinfo/samtools/1.18
Afficher l'aide des commandes
bwa
etbwa index
.Créer l'index bwa de la référence.
Solutionbwa index Gallus_gallus.Gallus_gallus-5.0.dna.toplevel_chr25-26.fa
Explorer l'aide de la commande
bwa mem
.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 fastidieuseCréer à la main un fichier
1_bwamem.jobs
contenant toutes les commandesbwa
: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éeBash
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 variableR1
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îneSRR7062654_R1.fastq.gz
devient doncSRR7062654
.${R1/_R1.f/_R2.f}
: affiche le contenu de la variableR1
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 commandels
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 variableid
avec le nom du fichier fastq$id=~s/\_.*//;
: subtitution de tout ce qui suit le premier '_' par rien dans la variableid
$out=$_;
: initialisation de la variableout
avec le nom du fichier fastq$r2=$_;
: initialisation de la variabler2
avec le nom du fichier fastq$r2=~s/\_R1\.f/\_R2\.f/;
: subtitution de la chaine _R1 par _R2 dans la variabler2
- 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