使用 R 語言的 ShortRead 套件讀取 FASTQ 檔

發表於 更新於

R 的 ShortRead 套件用於處理高通量定序產出的 FASTQ 檔，其功能包含讀寫 FASTQ 檔、計算序列統計量、裁切與篩選低品質序列等。目前，ShortRead 被保留在 Bioconductor ，並由開發團隊持續更新與維護。

在 R 視窗中，輸入以下指令即可安裝此套件。

1
2
3
4
>if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("ShortRead")

FASTQ 格式簡介

FASTQ 格式以四列字串記錄每一條核酸的辨識碼、序列及定序品質。

1
2
3
4
>@A77AV140207:1:1101:10732:16139/1  \\ sequence identifier
AACGTAGGTCACAAGCGTTGTCCGGAATTACTG   \\ raw sequence letters
+ 				    \\ optional description
?ABBBFFBF5B4BBG4AEEGGB5BEE1=1<E@1   \\ quality code

  • sequence identifier：首列以 @ 起頭，藉此區分每條核酸的紀錄。這一列通常會記錄序列的標識碼、定序儀器代碼、核酸在定序晶片的位置、barcode等資訊，不同定序平台有不同的格式。值得留意的是，在雙端定序的結果中，順向序列會在首列的最後標註 1，逆向序列則會標註 2。

  • raw sequence letters：次列則核酸各位置的鹼基種類，即 A、T、C、G、N（N 表示不明），所以此列長度等同定序讀長。

  • optional description：+ 號列，無含義。（雖然維基百科介紹說，這列可添加額外資訊，但目前還沒看過這樣標示的當案）

  • quality code：與第二列一一對應，反映定序品質的ASCII字符。

序列的定序品質可以用 error probability（定序儀器誤判鹼基種類的機率）評估，error probability 愈低則定序的品質越高，即兩者呈負相關。

為了便於理解和檢視，error probability 需經數值轉換和 ASCII 編碼，才會成為 FASTQ 檔內的 quality code，所以 FASTQ 格式的第四列才會存在各種符號。

簡言之，error probability (P) 會先轉變為整數形式且與品質正相關的 quality value (Q)。雖然曾出現許多種轉換方式，但目前最通用方式為 phred quality score，即

$$Q = -10 log_{10}{P}$$
$$P= 10^{-\frac{Q}{10}}$$​

接著，為了方便存儲，quality value 會再依 ASCII 規範編碼成單一字符。編碼方式則依照定序平台和其版本而異，所以平台間的 quality value 範圍也各不相同，可參考：https://en.wikipedia.org/wiki/FASTQ_format。

1
2
3
4
5
6
7
8
9
10
S - Sanger        Phred+33,  raw reads typically (0, 40)
X - Solexa        Solexa+64, raw reads typically (-5, 40)
I - Illumina 1.3+ Phred+64,  raw reads typically (0, 40)
J - Illumina 1.5+ Phred+64,  raw reads typically (3, 41)
    with 0=unused, 1=unused, 2=Read Segment Quality Control Indicator (bold) 
    (Note: See discussion above).
L - Illumina 1.8+ Phred+33,  raw reads typically (0, 41)
N - Nanopore      Phred+33,  Duplex reads typically (0, 50)
E - ElemBio AVITI Phred+33,  raw reads typically (0, 55)
P - PacBio        Phred+33,  HiFi reads typically (0, 93)

綜上所述，error probability 與 quality code 之間的關係為

1
2
3
>  error probability (float) 
→ quality value      (int) 		
→ quality code 	     (character)

讀取 fastq 檔案

使用 readFastq 即可讀取 FASTQ 並存為 ShortRead 物件。

1
2
3
4
5
>> library(ShortRead)
> rfq <- readFastq("./data/G64589_R1_001.fastq")
> rfq
class: ShortReadQ
length: 13510 reads; width: 175 cycles \\ length 是序列數，width 是序列長

也可以借助檔名的共同模式，一次讀取多個 FASTQ 檔。然而，這項操作會合併所有檔案的序列，所以無法追溯各自源於哪個檔案。

1
2
3
4
> rfqs <- readFastq("./data/", pattern = "_001.fastq")
> rfqs
class: ShortReadQ
length: 51508 reads; width: 175 cycles

檢視序列的標識碼、鹼基種類與品質分數

ShortReadQ 物件包含 @id、 @sread 和 @qulity，分別對應 FASTQ 格式的ID、sequence letters 與 quality code。

瀏覽標識碼

1
2
3
4
5
> rfq@id[1:2]	\\ or id(rfq)[1:2]
BStringSet object of length 2:
    width seq
[1]    32 A77AV140207:1:1101:10732:16139/1
[2]    32 A77AV140207:1:1101:10743:16162/1

id(rfq)這寫法有時會和其他套件的 function 衝突，這時要改寫為 ShortRead::id(rfq)才能順利運作。

瀏覽序列

1
2
3
4
5
> rfq@sread[1:2] \\ or sread(rfq)[1:2]
DNAStringSet object of length 2:
    width seq
[1]   175 AACGTAG...GATAC
[2]   175 AACGTAG...ATATG

瀏覽品質分數

1
2
3
4
5
6
7
> rfq@quality[1:2] \\ or quality(rfq)[1:2]
class: FastqQuality
quality:
BStringSet object of length 2:
    width seq
[1]   175 ?ABBBFF...#####
[2]   175 AAAAACF...#####

查閱每個 ASCII 記號對應的品質分數

1
2
3
> encoding(quality(rfq))[1:15]
 !  "  #  $  %  &  '  (  )  *  +  ,  -  .  / 
 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14

取得定序品質指標

將 FastQuality class 的物件轉換為 quality value matrix。其中，row 為序列， column 則是序列上各鹼基的位置，每一元素都記錄了序列各鹼基的 quality value。 

1
2
3
4
5
6
7
8
9
10
11
12
># transform FastqQuality class into matrix
> rfq_to_matrix <- function (rfq) {
	return(as(quality(rfq), "matrix"))
}
> Qval <- rfq_to_matrix(rfq)
> Qval[1:5, 1:5]
     [,1] [,2] [,3] [,4] [,5]
[1,]   30   32   33   33   33
[2,]   32   32   32   32   32
[3,]   32   32   32   32   32
[4,]   32   32   33   33   33
[5,]   34   34   34   34   34

由於 quality value matrix 不含序列 ID，所以要手動添加

1
>row.names(Qval) <- as.character(rfq@id)

延伸計算

Quality value，依照 Phred score 關係式回推 error probability。

1
2
># Transform phred score into error probability
errProb <- 10^-(Qval/10)

接著加總各列的定序錯誤率，可獲得各序列的 expected errors（錯誤鹼基數量期望值）

1
2
># Calculate expected errors for each read
ee <- rowSums(errProb)

最後使用內建的資料選取方式，移除低品質序列。

1
2
># Discard read with expected errors above 5
> ee[ee <= 5]

小結

1
2
3
4
5
6
7
8
># 讀取 FASTQ 為 ShortRead 物件
rfq <- readFastq(path)
# 轉換 ShortRead 物件為 quality value matrix
Qval <- as(quality(rfq), "matrix")
# 轉換 Phred score 為 error probability
errProb <- 10^-(Qval/10)
# 計算 expected errors
ee <- rowSums(errProb)