Différences

Ci-dessous, les différences entre deux révisions de la page.

--- fr:impilopedia:genex:rnaseq:align:bam_qual_analysis [2016/12/07 09:55] – créée foisys
+++ fr:impilopedia:genex:rnaseq:align:bam_qual_analysis [2021/05/29 15:35] (Version actuelle) – modification externe 127.0.0.1
@@ Ligne 1: / Ligne 1: @@
+====== Comment trouver les valeurs de qualité et en exploiter l'information ======
+===== Introduction =====
+Comme mentionné ailleurs, les aligneurs de séquence pour RNASeq ont une [[:fr:impilopedia:genex:rnaseq:align:filter_bad_align|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:
+^  Type d'aligneur  ^  Valeur correspondant aux "bons" alignements  |
+^  HISAT2  |  60 (depuis la version 2.0.5; 255 avant)  |
+^  Tophat2  |  50  |
+^  STAR  |  255  |
+===== Protocole =====
+  * En utilisant un fichier SAM ou BAM lu par ''samtools'', capturez les valeurs de la 5ème colonne:
+<code bash>
+% samtools view un_fichier_sequences.bam | cut -f5 > un_fichier_sequences_qual.txt
+</code>
+  * 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:
+<code rsplus>
+> data.r<-read.csv(file="un_fichier_sequences_qual.txt",header=FALSE)
+</code>
+  * On doit transformer cet objet en liste:
+<code rsplus>
+> data.l<-unlist(data.r)
+</code>
+  * On applique la magie de la méthode ''hist()'':
+<code rsplus>
+> hist(data.l,col="blue")
+</code>
+  * On remarque que la très grande majorité des lectures se retrouve dans la barre des 50 (exemple avec Tophat2):
+{{ :fr:impilopedia:genex:rnaseq:align:bam_qual:rplot_qual_rnaseq_bam.png?800 |}}
+  * On peut extraire des données descriptives de cet histogramme en interceptant les infos silencieusement créées par la méthode:
+<code rsplus>
+> 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"
+</code>
+  * 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!