Qualité des séquences

Introduction à Bioconductor en R

Paula Andrea Martinez, PhD.

Data Scientist

Scores de qualité – table Phred

Valeur de qualité	Probabilité d'erreur	Exactitude (%)
10	1 sur 10	90
20	1 sur 100	99
30	1 sur 1000	99,9
40	1 sur 10000	99,99
50	1 sur 100000	99,999

Encodage – Phred +33

# quality encoding
encoding(quality(fqsample))

Caractères d'encodage et leurs scores

 !  "  #  $  %  &  '  (  )  *  +  ,  -  .   #  encoding
 0  1  2  3  4  5  6  7  8  9 10 11 12 13   #  score

 /  0  1  2  3  4  5  6  7  8  9  :  ;  <   #  encoding 
14 15 16 17 18 19 20 21 22 23 24 25 26 27   #  score 

 =  >  ?  @  A  B  C  D  E  F  G  H  I      #  encoding
28 29 30 31 32 33 34 35 36 37 38 39 40      #  score

Qualité fastq

library(ShortRead)
quality(fqsample)

class: FastqQuality
A BStringSet instance 

# La qualité est représentée par des caractères ASCII 
[1]    40 ?@@DDDDDHDFDHE>AHFEGFIIEBGDBHH<3FEBEEEEG
[2]    40 BCCDFFFFHHHHHJJJJJJJJJJEHHGHIJJJJJJJJJJJ
[3]    40 BCCFFFFFHFHHHJJJJJJIIJJIIIIIGIIJJIJGIJII
[4]    40 CCCFFFFFHHHHHJJJJJJJJJJIJJJJJJJJJJJJJJJJ

library(ShortRead)

sread(fqsample)[1]

# La qualité est représentée par des caractères ASCII 
quality(fqsample)[1]

50 GTCCCATTTACCTCTGACTCTTTTGATGCTGCAATTGCTGCTCATATACT

50 ?@@DDDDDHDFDHE>AHFEGFIIEBGDBHH<3FEBEEEEGGIGIIGHGHC

## PhredQuality instance
pq <- PhredQuality(quality(fqsample))

# transformer l'encodage en scores 
qs <- as(pq, "IntegerList")
qs # afficher les scores

30 31 31 35 35 35 35 35 39 35 37 35 39 36 29 32 39 37 36 38 37 40 40 36 33 38 35 33 39 39 27 18 37 36 33 36
36 36 36 38 38 40 38 40 40 38 39 38 39 34

Évaluation de la qualité

library(ShortRead)
# évaluation de la qualité
qaSummary <- qa(fqsample, lane = 1)    # voie facultative

# class: ShortReadQQA(10)
# Noms accessibles dans le sommaire d'évaluation de la qualité
names(qaSummary)

 [1] "readCounts"           "baseCalls"            "readQualityScore"          "baseQuality"
 [5] "alignQuality"         "frequentSequences"    "sequenceDistribution"      "perCycle"  
 [9] "perTile"              "adapterContamination"
# Les éléments QA sont accessibles avec qa[["name"]]

# Obtenir un rapport HTML
browseURL(report(qaSummary))

library(ShortRead)

# alphabet des séquences
alphabet(sread(fullSample))

A,C,G,T,M,R,W,S,Y,K,V,H,D,B,N,-,+,.

abc <- alphabetByCycle(sread(fullSample))

# Chaque observation est une lettre et chaque variable, un cycle. D'abord, sélectionner les 4 premières lignes : nucléotides A, C, G, T 
# Puis transposer
nucByCycle <- t(abc[1:4,])

nucByCycle <- nucByCycle %>% as_tibble() %>% # convertir en tibble
                      mutate(cycle = 1:50) # ajouter les numéros de cycle
nucByCycle

    A     C     G     T cycle
16839 16335 16740 10878     1     
13056 13327 12064 22389     2     
13666 15617 13198 18355     3     
14723 15439 14239 16435     4

Êtes-vous emballé(e) ?

Introduction à Bioconductor en R