R言語で核酸配列データを扱ってみた件【その2: FASTAフォーマット・ファイルの読み込み】

はじめに
seqinrパッケージのインストール
Biostringsパッケージのインストール
seqinr::read.fastaを使って、HIV-1 全ゲノムのFASTAファイルの読み込み場合
Biostrings::readDNAStringSetを使って、HIV-1 全ゲノムのFASTAファイルの読み込み場合
- HIV-1 (ヒト免疫不全ウイルス-1) の全ゲノム RefSeq

はじめに

Rを用いた、DNA配列データの基本的な処理の方法を概説します。

今回は、seqinrあるいはBiostringsパッケージを用いて FASTA形式のファイルを読み込んで、DNA配列データとして扱ってみます。

配列データを元にして、相補配列を作成します。

seqinrパッケージのインストール

まずは、Rを起動して、seqinrパッケージをセットアップします。

#インストール
install.packages("seqinr")

#ロード
library(seqinr)

Biostringsパッケージのインストール

次に、Biostringsパッケージをセットアップします。

#インストール
install.packages("BiocManager")
BiocManager::install("Biostrings")

#ロード
library(Biostrings)

seqinr::read.fastaを使って、HIV-1 全ゲノムのFASTAファイルの読み込み場合

Human immunodeficiency virus 1 の全ゲノム(NC_001802.1)の FASTAファイルを用意します。そのFASTAファイル(DNA配列)を読み込みます。

FASTA形式とは、核酸配列やアミノ酸配列を記述するためのデータ形式のひとつです。

この形式では、「>」で始まる行に配列の名前を記載して、次の行以降に配列情報を示した形式です。「>」で始まる行を複数設定することで、1つのFASTAファイル内に、複数の配列情報を列挙できます。

seqinr::read.fastaを使って、FASTAを「ベクトル」として読み込む

#「ベクトル」として読み込む
HIV <- seqinr::read.fasta("https://gist.githubusercontent.com/kumeS/aead5195f323560ab906aff290f51b3f/raw/dedf6599d4c57abc9924c8fa19a2a06c3295cb54/HIV1.fasta")

str(HIV)
#List of 1
# $ NC_001802.1: 'SeqFastadna' chr [1:9181] "g" "g" "t" "c" ...
#  ..- attr(*, "name")= chr "NC_001802.1"
#  ..- attr(*, "Annot")= chr ">NC_001802.1 Human immunodeficiency virus 1, complete genome"
  
#文字列に変換
HIVs <- seqinr::c2s(HIV$NC_001802.1)

str(HIVs)
# chr "ggtctctctggttagaccagatctgagcctgggagctctctggctaactagggaacccactgcttaagcctcaataaagcttgccttgagtgctt#caagtagtgtgtgccc"| __truncated__

seqinr::read.fastaを使って、FASTAを「文字列」として読み込む

#「文字列」として読み込む
HIV_str <- seqinr::read.fasta("https://gist.githubusercontent.com/kumeS/aead5195f323560ab906aff290f51b3f/raw/dedf6599d4c57abc9924c8fa19a2a06c3295cb54/HIV1.fasta", as.string=TRUE)

str(HIV_str)
#List of 1
# $ NC_001802.1: 'SeqFastadna' chr #"ggtctctctggttagaccagatctgagcctgggagctctctggctaactagggaacccactgcttaagcctcaataaagcttgccttgagtgcttcaagt#agtgtgtgccc"| __truncated__
#  ..- attr(*, "name")= chr "NC_001802.1"
#  ..- attr(*, "Annot")= chr ">NC_001802.1 Human immunodeficiency virus 1, complete genome"

#ベクトルに変換
HIVc <- seqinr::s2c(HIV_str$NC_001802.1)

str(HIVc)
#chr [1:9181] "g" "g" "t" "c" "t" "c" "t" "c" "t" "g" "g" ...

読み込んだ配列データから相補配列を求める配列データ処理の基本

#ベクトル部分の取り出し
HIVseq <- HIV$NC_001802.1

#相補配列の結果を返す
HIVcomp <- seqinr::comp(HIVseq, forceToLower = FALSE)

#headで表示
head(HIVseq, n=10)
#[1] "g" "g" "t" "c" "t" "c" "t" "c" "t" "g"

head(HIVcomp, n=10)
#[1] "c" "c" "a" "g" "a" "g" "a" "g" "a" "c"

Biostrings::readDNAStringSetを使って、HIV-1 全ゲノムのFASTAファイルの読み込み場合

Biostringsパッケージでは、FASTAファイル(DNA配列)の読み込みには、readDNAStringSet関数を用います。

#FASTAファイルの読み込み
HIV_bio <- Biostrings::readDNAStringSet(filepath="https://gist.githubusercontent.com/kumeS/aead5195f323560ab906aff290f51b3f/raw/dedf6599d4c57abc9924c8fa19a2a06c3295cb54/HIV1.fasta",
                 format="fasta")

HIV_bio
#DNAStringSet object of length 1:
#    width seq                            names               
#[1]  9181 GGTCTCTCTGGTTA...CCTTGAGTGCTTC NC_001802.1 Human...

#相補配列を求める
Biostrings::complement(HIV_bio)
#DNAStringSet object of length 1:
#    width seq                            names               
#[1]  9181 CCAGAGAGACCAAT...GGAACTCACGAAG NC_001802.1 Human...