【R】 DNA配列データの読み込みと変換

FASTA形式のDNA配列を読み込む

library(tidyverse)
library(ape)
library(Biostrings)
naseq1 <- readDNAStringSet("seq.fasta")         # -> DNAStringSet
naseq2 <- read.dna("seq.fasta", format="fasta") # -> DNAbin
naseq3 <- read.FASTA("seq.fasta")               # -> DNAbin(list)

FASTA形式で書かれたDNA配列を読み込む関数はいくつかありますが、それぞれ出力されるクラスが異なります。class()で調べることができます。naseq2とnaseq3は共にDNAbinですが、後者は内部がlist形式になっていてデータ構造が異なります。便宜的にDNAbin(list)と書き区別することにします。str()でどちらかが分かります。

FASTQ形式のDNA配列を読み込む

naseq4 <- readDNAStringSet("seq.fastq", with.qualities=TRUE, format="fastq") # -> DNAStringSet
naseq5 <- read.fastq("seq.fastq")                       # -> DNAbin(list)

Quality Scoreを反映させることができます。

NEXUS形式のDNA配列を読み込む

naseq6 <- read.nexus.data("seq.nex") # -> list

こちらはlistクラスになります。1文字ずつのベクトルのリストになっています。
次に、系統樹を書くために距離行列を計算してみます。

dist.dna(naseq2)
dist.dna(naseq3)

DNAbinとDNAbin(list)クラスを受け付けてくれます。
DNAStringSetとlistはそのままではエラーになります。

データの変換

以下の関数で変換できます。

as.DNAbin(naseq1)    # DNAStringSet -> DNAbin(list)
nexus2DNAbin(naseq6) # list -> DNAbin(list)
as.DNAbin(naseq6)    # list -> DNAbin(list)
as.list(naseqs2)     # DNAbin -> DNAbin(list)
as.matrix(naseq3)    # DNAbin(list) -> DNAbin
as.character(naseq3) # DNAbin(list) -> list

次にDNAStringSetクラスへの変換です。
DNAStringSetには、DNA配列の文字列(1文字ずつのベクトルではない)のベクトルを渡す必要があります。

# DNAbin -> DNAStringSet
DNAStringSet(map_chr(as.character(as.list(naseq2)), ~str_c(.x, collapse="")))

# DNAbin(list) -> DNAStringSet
DNAStringSet(map_chr(as.character(naseq3), ~str_c(.x, collapse="")))

# list -> DNAStringSet
DNAStringSet(map_chr(naseq6, ~str_c(.x, collapse="")))

ベクトル、リストデータとの変換

(2022.8.5 追記)

プログラム中の配列データの場合は、一旦DNAStringSetクラスに変換すると上記の変換ができます。

naseq.vec <- c(seq1="ATGCATAGT", seq2="ATGAAAGGC")
# 名前付きベクトル -> DNAStringSet
naseq7 <- DNAStringSet(naseq.vec)
# DNAStringSet -> 名前付きベクトル
as.character(naseq7)

naseq.list <- list(seq1="ATGCATAGT", seq2="ATGAAAGGC")
# リストデータ -> DNAStringSet 
naseq8 <- DNAStringSet(unlist(naseq.list))
# DNAStringSet -> リストデータ
as.list(as.character(naseq8))

データフレームとの変換

(2022.8.5 追記)

naseq.df <- data.frame(name=c("seq1", "seq2"), naseq=c("ATGCATAGT", "ATGAAAGGC"))

# data.frame -> DNAStringSet
naseq9 <- DNAStringSet(naseq.df$naseq)
names(naseq9) <- naseq.df$name
# DNAStringSet -> data.frame
data.frame(name=names(naseq9), naseq=as.character(naseq9))