Comment trouver les valeurs de qualité et en exploiter l'information

Comme mentionné ailleurs, les aligneurs de séquence pour RNASeq ont une manière très particulière de "calculer" la qualité des alignements obtenus… Il peut devenir nécessaire de voir la distribution des valeurs de qualité non seulement pour juger du résultat de l'alignement mais aussi pour connaitre la valeur extrême utilisée pour identifier les “bons” alignements.

Pour référence:

• En utilisant un fichier SAM ou BAM lu par samtools, capturez les valeurs de la 5ème colonne:

% samtools view un_fichier_sequences.bam | cut -f5 > un_fichier_sequences_qual.txt

• Ce fichier aura autant de lignes que de séquences, avec une seule valeur sur chaque ligne. Le reste de la magie se fait via R. Donc, après avoir démarrer R, on crée un objet contenant nos données:

> data.r<-read.csv(file="un_fichier_sequences_qual.txt",header=FALSE)

• On doit transformer cet objet en liste:

> data.l<-unlist(data.r)

• On applique la magie de la méthode hist():

> hist(data.l,col="blue")

• On remarque que la très grande majorité des lectures se retrouve dans la barre des 50 (exemple avec Tophat2):

• On peut extraire des données descriptives de cet histogramme en interceptant les infos silencieusement créées par la méthode:

> hist.d<-hist(data.l)
> hist.d
$breaks
 [1]  0  2  4  6  8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44
[24] 46 48 50
 
$counts
 [1]  1188218  1077044        0        0        0        0        0
 [8]        0        0        0        0        0        0        0
[15]        0        0        0        0        0        0        0
[22]        0        0        0 34913044
 
$density
 [1] 0.01597999 0.01448485 0.00000000 0.00000000 0.00000000 0.00000000
 [7] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
[13] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
[19] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
[25] 0.46953516
 
$mids
 [1]  1  3  5  7  9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45
[24] 47 49
 
$xname
[1] "data.l"
 
$equidist
[1] TRUE
 
attr(,"class")
[1] "histogram"

• On peut remarquer qu'on observe que trois classes contenant des valeurs autre que 0: 0-1, 2-3 et 49-50. De plus, il est assez évident que le maximum est 50 Donc, c'est cette valeur à prendre pour filtrer!