先前曾介紹 R 讀取 FASTQ 檔的方式,本文則要分別介紹如何在 R 環境,使用 seqnir 套件或 baseR 讀取 FASTA 檔。
FASTA 檔
FASTA 格式用於紀錄核酸或胜肽的序列。以
基於增幅子的資料處理為例,諸如 unique reads、OTU/ASV list、SILVA/RDP reference sequence 皆是以 FASTA 格式儲存。
在 FASTA 格式 ,一條序列的紀錄由兩部分組成:
- 以 “>” 起頭的描述列,用以區隔不同序列的紀錄,並保存序列的名稱、辨識碼或來源等資訊。
- 由鹼基或胺基酸字符組成的序列字串(鹼基和胺基酸字符的 IUPAC 命名可參考維基百科的簡介)。
典型的 FASTA 格式如下,所屬文件通常會以 .fa、.fasta 或 .fna 為副檔名,但實際並無統一規範。
>OTU1
AATTCGCGACTACTATGGC
>OTU2
AATACGCTACTACTCTGGC
與 FASTQ 檔不同,FASTA 檔的序列字串未必在同一列。依慣例,序列字串長度不會超過 80 個字元(這可能跟方便在終端介面閱讀有關),所以較長的序列會被換行符分隔成好幾列。
Seqinr
Seqinr 是專門用於處理核酸和蛋白質序列的套件,其功能枚舉如下:
- 自Genbank等資料庫下載序列
- 讀寫 FASTA 檔
- 計算各式鹼基統計量
- 序列比對
- 轉錄或轉譯等字符抽換
由於 seqinr 存放在 CRAN 庫,所以可以直接輸入install.packages("seqinr")
安裝。
# install and read seqinr package
install.packages("seqinr")
library(seqinr)
先來看看我們的示範資料 dna.fasta
。為了凸顯 FASTA 檔的序列字串不一定在同一行的特質,我故意將示範資料段設定得這麼混亂。然而,正常的資料雖然列數可能隨序列長度而異,但每一列的字串長度是一致的。
>seq1
ATCGATCGATCGATCGATCGAAAAAATCTCAAT
>seq2
GTCGATCTCGCTTGCTCCA
GCGGTAGC
>seq3
TGATCGAGATCGAT
TTAG
CAA
>seq4
TGATCGAGATCGATTTAGCAA
1. 讀取 fasta 檔
以 seqinr 的 read.fasta 能把 FASTA 檔讀為由字串向量構成的 list。
> dna <- read.fasta("dna.fasta")
$seq1
[1] "a" "t" "c" "g" "a" "t" "c" "g" "a" "t" "c" "g" "a" "t" "c" "g" "a"
"t" "c" "g" "a" "a"
[23] "a" "a" "a" "a" "t" "c" "t" "c" "a" "a" "t"
attr(,"name")
[1] "seq1"
attr(,"Annot")
[1] ">seq1"
attr(,"class")
[1] "SeqFastadna"
(下略)
2. 序列簡易操作
統計序列長度分布
# 列出各序列的長度
> sapply(dna, length)
seq1 seq2 seq3 seq4
33 27 21 21
# 統計各長度序列的數量
> table(sapply(dna, length))
21 27 33
2 1 1
計算 GC content
> sapply(dna, GC)
seq1 seq2 seq3 seq4
0.3636364 0.6296296 0.3809524 0.3809524
將 DNA 序列轉換為胺基酸序列
> translate(dna$seq1)
[1] "I" "D" "R" "S" "I" "D" "R" "K" "N" "L" "N"
繪製序列比對的 dot plot
> dotPlot(dna$seq1, dna$seq3)
選取特定序列的方法
dna[sapply(dna, GC) > 0.5] \\ 尋找 GC content > 0.5 者
na[names(dna) %in% c("seq1", "seq2")] \\
3. 儲存 FASTA 檔
write.fasta(dna, names = names(dna), file.out = "dna_rw.fasta", nbchar = 5)
names
:設定序列的名稱,也就是 “>” 後的文字file.out
:儲存的檔名與路徑nbchar
:序列幾個字符一列
Biostrings
除了保存於 CRAN 的 seqinr,存放於 bioconductor 的 Biostrings 也能讀取 FASTA 檔。於 R 終端輸入以下指令即可安裝。
# install and read Biostrings package
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("Biostrings")
library(Biostrings)
此處,我使用 readDNAStringSet 讀取FASTA檔為 DNAStringSet 物件
> dna <- readDNAStringSet("dna.fasta")
> dna
DNAStringSet object of length 4:
width seq names
[1] 33 ATCGATCGATCGATCGATCGAAAAAATCTCAAT seq1
[2] 27 GTCGATCTCGCTTGCTCCAGCGGTAGC seq2
[3] 21 TGATCGAGATCGATTTAGCAA seq3
[4] 21 TGATCGAGATCGATTTAGCAA seq4
除了readDNAStringSet,也可以針對序列的類別分別使用 readRNAStringSet 和 readAAStringSet 讀取 FASTA 檔。讀取的方式會影響能套用的 function 數據的存儲與呈現。例如若以 readRNAStringSet 讀取 DNA 序列會自動掠過不符類別的字符。
> readRNAStringSet("dna.fasta")
RNAStringSet object of length 4:
width seq names
[1] 25 ACGACGACGACGACGAAAAAACCAA seq1
[2] 20 GCGACCGCGCCCAGCGGAGC seq2
[3] 15 GACGAGACGAAGCAA seq3
[4] 15 GACGAGACGAAGCAA seq4
Warning message:
In .Call2("fasta_index", filexp_list, nrec, skip, seek.first.rec, :
reading FASTA file dna.fasta: ignored 27 invalid one-letter sequence
codes
另外,亦可使用 readBStringSet 將 FASTA 檔視為一般字串讀取,只是這樣就無法使用某些 function,例如求序列的 reverse complement:
# read .fasta into DNA string set
> readDNAStringSet("dna.fasta")
> reverseComplement(dna$seq1)
33-letter DNAString object
seq: ATTGAGATTTTTTCGATCGATCGATCGATCGAT
# read .fasta into ordinary string set
> readBStringSet("dna.fasta")
> reverseComplement(dna$seq1)
Error in (function (classes, fdef, mtable) :
unable to find an inherited method for function ‘reverseComplement’ for
signature ‘"BString"’
選取特定序列則需要調出 DNAStringSet 物件的 range 資料,range 資料是一個 data.frame,紀錄了序列的長度和名稱等資訊。
> dna@ranges
group start end width names
1 1 1 33 33 seq1
2 1 34 60 27 seq2
3 1 61 81 21 seq3
4 1 82 102 21 seq4
dna[dna@ranges@width > 25] \\ 保留長度 > 25 bp 的序列
dna[dna@ranges@NAMES == "seq1"] \\ 保留 seq1
而使用 writeXStringSet 儲存 FASTA 檔時,會把示範資料中分為多列的序列貼到同一列。
writeXStringSet(dna, "dna_rw.fasta", format = "fasta")
Biostrings 其他功能可參考 Biostrings Quick Overview
baseR
若不想為了簡單的序列統計和操作而裝套件,也可以使用 baseR 將 FASTA 檔。
1. 讀取 fasta 檔
簡言之,以 readLines 讀取 FASTA 檔所有內容並存為字串向量。接著,使用 grep 搜尋以 “>” 開頭的描述列之索引號,並以 gsub 移除 “>” 取回序列辨識碼。最後利用 paste連接位於描述列的索引號之間的鹼基字串。
> read_fasta <- function (path) {
lines <- readLines(path)
ind <- grep(">", lines)
start <- ind + 1
end <- c((ind - 1)[-1], length(lines))
seqs <- lapply(1:length(ind), function(i) paste(lines[start[i]:end[i]],
collapse = ""))
names(seqs) <- gsub(">", "", lines[ind])
return(seqs)
}
> dna <- read_fasta("dna.fasta")
> dna
$seq1
[1] "ATCGATCGATCGATCGATCGAAAAAATCTCAAT"
$seq2
[1] "GTCGATCTCGCTTGCTCCAGCGGTAGC"
(下略)
2. 序列簡易操作
統計序列長度分布
# 列出各序列的長度
> nchar(dna)
seq1 seq2 seq3 seq4
33 27 21 21
# 統計各長度序列的數量
> table(nchar(dna))
21 27 33
2 1 1
繪製序列長度數量分布直方圖
> hist(nchar(dna))
辨識 unique sequence
> unique(dna)
[[1]]
[1] "ATCGATCGATCGATCGATCGAAAAAATCTCAAT"
[[2]]
[1] "GTCGATCTCGCTTGCTCCAGCGGTAGC"
[[3]]
[1] "TGATCGAGATCGATTTAGCAA"
選取特定序列
dna[nchar(dna) > 25] \\ 保留長度 > 25 bp 的序列
dna["seq1"] \\ 保留 seq1
3. 儲存 FASTA 檔
即用 names 取出 list 的名稱(也是序列的名稱),再貼上 “>” 依序寫到文件之中。
> write_fasta <- function (obj, path) {
lapply(names(obj), function (x) {
write(paste0(">", x), path ,append = TRUE)
write(obj[[x]], path, append = TRUE)
})
}
write_fasta(dna, "dna_rw.fasta")
小結
# 以 seqinr 讀寫 FASTA 檔
library(seqinr)
read.fasta("dna.fasta")
write.fasta(dna, names = names(dna), file.out = "dna_rw.fasta")
# 以 Biostrings 讀寫 FASTA 檔
library(Biostrings)
readDNAStringSet("dna.fasta")
writeXStringSet(dna, "dna_rw.fasta")
# 以 baseR 讀寫 FASTA 檔
(參考前述自訂 function)
沒有留言:
張貼留言