Filtrage des alignements basé selon leur qualité

Introduction

Il est quasi-impossible de tout aligner: beaucoup de vos séquences proviendront de séquences répétitives, donc difficiles à placer ou bien, si vous utilisez un séquençage paired-end, vos deux lectures associées ne seront pas placé ensembles… Il vous faudra donc retirer tout ça pour éviter de polluer notre ensemble de données pour l'expression différentielle.

Cette recette est inspirée de trois pages très intéressantes:

Procédure

Premièrement, retirer les séquences qui ne sont pas aligner ensembles. Si vous avez fait un séquençage paired-end, certaines paires pourraient ne pas avoir été alignées dans la bonne orientation, par exemple sur deux chromosomes distincts ou bien dans le même sens plutôt que 5'→3' 3'←5' tel qu'attendu. On utilisera une des colonnes du format SAM pour réussir ceci. Les aligneurs tel Tophat 2 ou STAR indiquent dans la 2ème colonne (FLAG) une valeur entière correspondant à la somme des bits d'information applicables à cette séquence. Dans notre cas, nous recherchons une valeur numérique de 2 ou bien sa représentation binaire de 0x2. Selon la définition SAM, ceci correspond à Each segment properly aligned according to the aligner. Non seulement nous voulons cette valeur mais toutes les autres combinaisons incluant celle-ci.

Deuxièmement, on veut retirer les alignements ambiguës, placés à plusieurs endroits possibles sur le génome de référence. À première vue, il ne semble pas aisé de la faire mais les aligneurs spécialisés en RNASeq ont tordu le standard SAM pour y insérer l'information nécessaire… Le format SAM spécifie que la 5ème colonne du fichier contient une valeur numérique, de 0 à 255, contenant la qualité de l'alignement selon la norme MAPQ. Comme Tophat 2 est un wrapper pour Bowtie 2, on pourrait s'attendre à y retrouver les valeurs trouvées par Bowtie 2 mais ce n'est pas le cas… STAR d'ailleurs fait la même chose que Tophat 2, de même que Hisat 2. Quand on regarde la distribution des valeurs de cette colonne, on observe que nous ne voyons pas une distribution régulière des valeurs mais au contraire beaucoup de valeurs vers 0,1,2 et 3 ainsi que vers 50 (Tophat 2) ou 255 (STAR). Pourtant, 255 correspond selon MAPQ à l'absence de valeur…
- Voici une méthode pour trouver les valeurs contenues dans un fichier BAM.

Tophat 2 (comme STAR) a redéfini l'usage de la 5ème colonne: au lieu de contenir une valeur ressemblant à un Q score, les valeurs entières utilisées sont plutôt une façon de faire tenir de l'information pertinente sur l'unicité de l'alignement obtenu (selon cette source, celle-ci et surtout celle-ci). Ainsi:

Pour Tophat:
- 0: lecture alignée à 10 endroits et sur le génome
- 1: lecture alignée à 4-9 endroits sur le génome
- 2: lecture alignée à 3 endroits sur le génome
- 3: lecture alignée à 2 endroits sur le génome
- 50 ou 255: lecture alignée à un seul endroit sur le génome (255 est pour Tophat 1.4 et moins; 50 est pour toutes les versions depuis. Il n'y a pas de raison documentée pour ce changement…)
Pour STAR: même chose que pour Tophat mais la valeur maximale est 255, non 50.

Donc comment faire pour la suite? Heureusement, Samtools à ce qu'il faut! Nous devrons filtrer nos fichiers BAM pour ne garder que les lectures qui comptent. Nous utiliserons deux paramètres de filtrage pour la fonction view de Samtools:
- -f 0x2 nous permet de garder les séquences qui ont le FLAG 0x2 paramètré dans la 2ème colonne (-f 2 fonctionnera aussi!);
- -q int nous permet de garder les séquences qui ont une qualité d'alignement supérieure ou égale à une certaine valeur seuil int.

* Donc pour un fichier mySequences.bam, on utilisera Samtools de 2 manières possibles:

// Première méthode
% samtools view -f 0x2 -q 50 -b -o mySequences_properly_paired_unique_mapping.bam mySequences.bam

// Deuxième méthode
% samtools view -f 0x2 -q 50 -b mySequences.bam > mySequences_properly_paired_unique_mapping.bam