定序品質

R 中的 Bioconductor 入門

Paula Andrea Martinez, PhD.

Data Scientist

品質分數－Phred 對照表

品質值	出錯機率	準確率 (%)
10	1/10	90
20	1/100	99
30	1/1000	99.9
40	1/10000	99.99
50	1/100000	99.999

編碼－Phred +33

# quality encoding
encoding(quality(fqsample))

編碼字元與分數

 !  "  #  $  %  &  '  (  )  *  +  ,  -  .   #  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

fastq 品質

library(ShortRead)
quality(fqsample)

class: FastqQuality
A BStringSet instance 

# 品質以 ASCII 字元表示 
[1]    40 ?@@DDDDDHDFDHE>AHFEGFIIEBGDBHH<3FEBEEEEG
[2]    40 BCCDFFFFHHHHHJJJJJJJJJJEHHGHIJJJJJJJJJJJ
[3]    40 BCCFFFFFHFHHHJJJJJJIIJJIIIIIGIIJJIJGIJII
[4]    40 CCCFFFFFHHHHHJJJJJJJJJJIJJJJJJJJJJJJJJJJ

library(ShortRead)

sread(fqsample)[1]

# 品質以 ASCII 字元表示 
quality(fqsample)[1]

50 GTCCCATTTACCTCTGACTCTTTTGATGCTGCAATTGCTGCTCATATACT

50 ?@@DDDDDHDFDHE>AHFEGFIIEBGDBHH<3FEBEEEEGGIGIIGHGHC

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

# 將編碼轉成分數 
qs <- as(pq, "IntegerList")
qs # 列印分數

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

品質評估

library(ShortRead)
# 品質評估
qaSummary <- qa(fqsample, lane = 1)    # 選填的 lane

# class: ShortReadQQA(10)
# 可從品質評估摘要存取的名稱
names(qaSummary)

 [1] "readCounts"           "baseCalls"            "readQualityScore"          "baseQuality"
 [5] "alignQuality"         "frequentSequences"    "sequenceDistribution"      "perCycle"  
 [9] "perTile"              "adapterContamination"
# 以 qa[["name"]] 取用各 QA 元素

# 產生 HTML 報告
browseURL(report(qaSummary))

library(ShortRead)

# 定序字母表
alphabet(sread(fullSample))

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

abc <- alphabetByCycle(sread(fullSample))

# 每筆觀測是一個字母，每個變數是一個循環。先選前 4 列核苷 A、C、G、T，
# 再轉置
nucByCycle <- t(abc[1:4,])

nucByCycle <- nucByCycle %>% as_tibble() %>% # 轉為 tibble
                      mutate(cycle = 1:50) # 加入循環編號
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

準備好了嗎？

R 中的 Bioconductor 入門