Excercice 3 : manipulation du format SAM/BAM - samtools

1. Attention, se connecter sur un noeud du cluster.

   srun --pty bash
   

2. Afficher l'aide de la commande samtools et parcourir les commandes disponibles.

   Solution

   samtools --help
   

3. Quels sont les différentes valeurs possibles du champ « FLAG » (et leurs effectifs) pour le run SRR7062654 ?
   Aller sur le site picard/explain-flag et vérifier à quoi correspond le champ « FLAG » rencontré le plus grand nombre de fois.

   Solution

   samtools view SRR7062654.bam | cut -f2 | sort -n | uniq -c | sort -rn
   

   492019 99
   492019 147
   491138 83
   491138 163
     4444 97
     4444 145
     4402 81
     4402 161
      433 73
      433 133
      422 181
      422 121
      396 69
      396 137
      392 185
      392 117
      173 177
      173 113
      135 65
      135 129
      104 2195
       82 2131
       76 2147
       72 2209
       71 2211
       69 2145
       63 2193
       59 2129
       23 2179
       23 2115
       20 2227
       20 2163
        5 2161
        3 2169
        3 2113
        2 2185
        2 2121
        1 2177
   

   99 correspond à :

   • read paired (0x1)
   • read mapped in proper pair (0x2)
   • mate reverse strand (0x20)
   • first in pair (0x40)

   147 correspond à :

   • read paired (0x1)
   • read mapped in proper pair (0x2)
   • read reverse strand (0x10)
   • second in pair (0x80)

4. Combien de lectures ont un champ « CIGAR » égal à 90M ?

   Solution

   samtools view SRR7062654.bam | cut -f6 | grep -c '90M'
   

   1861309
   

5. Combien de lectures sont parfaitement alignées ?

   Solution

   samtools view SRR7062654.bam | grep -c 'NM:i:0'
   

   1080285
   

6. Comment expliquer cette différence ?

   Solution

   M dans le champ « CIGAR » = Match ou Mismatch

7. Utiliser la commande calmd pour mettre en évidence les mismatches.

   Solution

   samtools calmd -e SRR7062654.bam Gallus_gallus.Gallus_gallus-5.0.dna.toplevel_chr25-26.fa | more
   

8. Quelle est la valeur de qualité de mapping maximum et combien de lectures ont cette valeur ?

   Solution

   samtools view SRR7062654.bam | cut -f5 | sort -n | uniq -c | tail -1
   

   1915573 60
   

9. Convertir le fichier BAM du run SRR7062654 en SAM puis en BAM (une seule ligne de commande).

   Solution

   samtools view SRR7062654.bam | samtools view -b - > tmp_output.bam
   

10. En cas d'erreur, attention au header !

    Solution

    samtools view -h SRR7062654.bam | samtools view -b - > tmp_output.bam
    

11. Supprimer le fichier tmp_output.bam.

    Solution

    rm -f tmp_output.bam
    

12. Trier le bam du run SRR7062654 par nom de reads et non par position (tri par défaut).

    Solution

    samtools sort -n SRR7062654.bam -o SRR7062654_byName.bam
    

13. Supprimer le fichier SRR7062654_byName.bam.

    Solution

    rm -f SRR7062654_byName.bam
    

14. Indexer les fichiers bam de votre répertoire de travail.

    Solution

    Bash

    for bam in *.bam; do echo "samtools index $bam"; done > 2_index.jobs
    

    Perl

    \ls -1 *.bam | perl -lne 'print "samtools index $_"' > 2_index.jobs
    

    sarray -J 2_index -e LOGS/%x_%j.err -o LOGS/%x_%j.out 2_index.jobs
    

15. Utiliser la commande faidx de samtools pour extraire les 1 000 premières bases du chr25.

    Solution

    samtools faidx Gallus_gallus.Gallus_gallus-5.0.dna.toplevel_chr25-26.fa 25:1-1000
    

16. Combien de lectures du run SRR7062654 s'alignent entre les coordonnées 1 et 1000 du chr25 ?

    Solution

    samtools view SRR7062654.bam 25:1-1000 | wc -l
    

    74
    

17. Calculer les statistiques d'alignement de chaque fichier bam (commande flagstat de samtools, exécution sur le cluster).

    Solution

    Bash

    for bam in *.bam; do echo "samtools flagstat $bam > ${bam}.flagstat"; done > 3_flagstat.jobs
    

    Perl

    \ls -1 *.bam | perl -lne 'print "samtools flagstat $_ > $_.flagstat"' > 3_flagstat.jobs
    

    sarray -J 3_flagstat -e LOGS/%x_%j.err -o LOGS/%x_%j.out 3_flagstat.jobs
    

Questions avancées :

1. En utilisant le champ « CIGAR », lister les tailles de délétions et leurs effectifs dans le fichier bam du run SRR7062654 ?

   Solution

   Awk

   samtools view SRR7062654.bam | awk -F"\t" 'match($6,/[0-9]+D/){print substr($6,RSTART,RLENGTH-1)}' | sort -n | uniq -c
   

   Perl

   samtools view SRR7062654.bam | perl -lane 'if($F[5]=~/(\d+)D/) {print $1}' | sort -n | uniq -c
   

2. En utilisant le champ « MD », afficher le nombre de total de mismatches dans le fichier bam du run SRR7062654 ?

   Solution

   Awk

   samtools view SRR7062654.bam | awk 'match($0,/MD:Z:([0-9]+(([A-Z]|\^[A-Z]+)[0-9]+)*)/){md=substr($0,RSTART+5,RLENGTH-5); gsub(/\^[ATGC]+/,"",md); gsub(/\^[ATGC]+/,"",md); gsub(/[0-9]+/,"",md); t+=length(md)}END{print t}'
   

   Perl

   samtools view SRR7062654.bam | perl -lne 'if($_=~/MD:Z:([0-9]+(([A-Z]|\^[A-Z]+)[0-9]+)*)/) {$md=$1; $md=~s/\^[ATGC]+//g; $md=~s/\d+//g; $t+=length($md);} END{print $t}'
   

   1589125