#過去記事
前アカウントで書いた記事です。
Bioconductor 解説 Biostrings 編 その１
Bioconductor 解説 Biostrings 編 その２
もう前記事書いたの２年前ですね。。。責任をもって終わらせます。
#はじめに
今回は、前回紹介したBString、DNAString、RNAString、AAStringオブジェクト（XStringオブジェクト）をリスト化したオブジェクトである BStringSet, DNAStringSet, RNAStringSet, AAStringSetオブジェクト（XStringSetオブジェクト）を紹介したいと思います。
そして次の回で前回予告したこれらオブジェクトを引数として要求する諸関数についてみていきたいと思います。
#BStringSet, DNAStringSet, RNAStringSet, AAStringSet
一つの配列情報を記録するXString オブジェクトですが、時として一つでなく複数の配列を一度に処理する場面もありえます。そういった場合に一括処理するオブジェクトがこのXStringSetオブジェクトです。
オブジェクト（インスタンス）を作成するコンストラクタは以下のとおりです。
BStringSet(x=c( "ABCDEFGHIJK", "LMNOPQRSTUVWXYZ" ) )
DNAStringSet( x=c( "ATCGGCAAT", "TTGTACCC" ) )
RNAStringSet( x=c( "AUCGGCAAU", "UUGUACCC" ) )
AAStringSet( x=c( "MARKSLEMSIR", "SRRRIKSMLE" ) )
出力されるインスタンスは以下の通りです。
A BStringSet instance of length 2
width seq
[1] 11 ABCDEFGHIJK
[2] 15 LMNOPQRSTUVWXYZ
A DNAStringSet instance of length 2
width seq
[1] 9 ATCGGCAAT
[2] 8 TTGTACCC
A RNAStringSet instance of length 2
width seq
[1] 9 AUCGGCAAU
[2] 8 UUGUACCC
A AAStringSet instance of length 2
width seq
[1] 11 MARKSLEMSIR
[2] 10 SRRRIKSMLE
引数は以下のようになっています。
#####x：{文字列型を要素としてもつベクトル}
：配列データのベクトル
#####start：{整数}
：読み込み開始地点（負の数で末尾からの参照）
#####end：{整数}
：読み込み終了地点
#####width：{整数}
：読み込む文字の数
#####use.names：{論理値}
：TRUEの場合、ベクトルの要素名をそのまま使用する
ほとんどXStringオブジェクトの場合と同じです。
ただxがベクトル型であることと、use.namesで配列名を使用できる点が異なっています。
###引数について
まず引数xに関して詳しくみていきましょう。
さっきも言いましたが、XStringオブジェクトでは引数xは文字列型（例："xxx"）で、XStringSetオブジェクトでは引数xは文字列ベクトル（例：c( "xxx", "yyy" )）となっています。
配列の数に関してはメモリがある限りはどれだけでも格納できそうです。
試しに１千万個の配列を格納してみます。
BStringSet( x=rep("A", 1E+7) )
A BStringSet instance of length 10000000
width seq
[1] 1 A
[2] 1 A
[3] 1 A
[4] 1 A
[5] 1 A
... ... ...
[9999996] 1 A
[9999997] 1 A
[9999998] 1 A
[9999999] 1 A
[10000000] 1 A
１千億個の配列だと私のPCのメモリ限界でした。
BStringSet( x=rep("A", 1E+10) )
Error: vector memory exhausted (limit reached?)
ちなみにベクトル内の要素が複数である必要はありません。
以下は単一の配列を入力した例です。
DNAStringSet( x=c( "ATCGGCAAT") )
A DNAStringSet instance of length 1
width seq
[1] 9 ATCGGCAAT
このような性質からある程度はXStringオブジェクトの代替として使うことができます。たとえばXStringオブジェクトでは配列に名前をつけて管理できるようにはなっていないので、名前をつけた状態で管理するためにXStringSetオブジェクトを使う例が考えられます。
DNAStringSet( x=c( my.seq="ATCGGCAAT") )
A DNAStringSet instance of length 1
width seq names
[1] 9 ATCGGCAAT my.seq
逆に名前を消して登録する場合は、use.names=FALSEにします。
DNAStringSet( x=c( my.seq="ATCGGCAAT"), use.names=FALSE )
A DNAStringSet instance of length 1
width seq
[1] 9 ATCGGCAAT
ちなみに後付けで命名したい場合はnames()を使えばできます。
my.seq <- DNAStringSet( x=c( "ATCGGCAAT", "CTT" ) )
message("配列の名前を表示します。")
names(my.seq) #名前がない状態
message("命名します。")
names(my.seq) <- c( "new.seq1", "new.seq2" ) # 名前を文字列ベクトルで登録
message("配列の名前を表示します。")
names(my.seq) #名前がある状態
配列の名前を表示します。
NULL
命名します
配列の名前を表示します。
[1] "new.seq1" "new.seq2"
この名前の登録の仕方は、ベクトル型やリスト型で命名する時と同じですね。
Rになれない人は、関数(<オブジェクト>) <- "<名前>"という形に違和感を覚えるかもしれません。Rで要素名や列名を記載する際はこのような記法でSetter的な動きをします。
###XStringSetのメソッド・その他操作に関して
####length()に気を付ける
一部のXStringSetオブジェクトのメソッドは、XStringオブジェクトのメソッドとよく似ています。ただ微妙に仕様が違う点に注意してください。その一例としてlength()があります。
XStringオブジェクトにlength()を適用した場合は、配列長を返します。
length( DNAString("ATCGGCAAT") )
[1] 9
しかし、XStringSetオブジェクトにlength()を適用した場合は、配列の数を返します。
length( DNAStringSet( x=c( A="ATCGGCAAT", B="CTT" ) ) )
[1] 2
XStringSetオブジェクトの各配列長を参照する場合は、nchar()もしくはwidth()を利用します。
nchar( DNAStringSet( x=c( A="ATCGGCAAT", B="CTT" ) ) )
[1] 9 3
ちなみにXStringオブジェクトの場合も、nchar()で配列長を参照できます。ややこしいので「配列長の参照はnchar()だけを使う」というように決めておくほうが混乱しなくてすむかもしれませんね。
nchar( DNAString("ATCGGCAAT") )
[1] 9
####subseq()で一括切り出し
XStringSetを使うことの一番わかりやすい利点は一括で配列を処理できる点にあります。
次の例ではsubseq()を使って１番目から３番目までの配列を切り出しています。
ori.seq <- DNAStringSet( x=c( A="ATCGGCAAT", B="CTT" ) )
subseq( x=ori.seq, start=1, end=3 )
A DNAStringSet instance of length 2
width seq names
[1] 3 ATC A
[2] 3 CTT B
この処理は高速ですので、たとえ一千万個程度の配列でも一秒程度で処理できます。
numerous.seq <- BStringSet( x=rep("ATTTCGCGTTTGCGATGA", 1E+7) )
subseq( x=numerous.seq, start=1, end=15 )
A BStringSet instance of length 10000000
width seq
[1] 15 ATTTCGCGTTTGCGA
[2] 15 ATTTCGCGTTTGCGA
[3] 15 ATTTCGCGTTTGCGA
[4] 15 ATTTCGCGTTTGCGA
[5] 15 ATTTCGCGTTTGCGA
... ... ...
[9999996] 15 ATTTCGCGTTTGCGA
[9999997] 15 ATTTCGCGTTTGCGA
[9999998] 15 ATTTCGCGTTTGCGA
[9999999] 15 ATTTCGCGTTTGCGA
[10000000] 15 ATTTCGCGTTTGCGA
####配列の抽出・変換
XStringSetオブジェクトにはいった配列の一部を取り出したい場合は、[<インデックスのベクトル>]で取り出し可能です。取り出された配列は同じくXStringSetオブジェクトとして生成されます。
ori.seq <- DNAStringSet( x=c( A="ATCGGCAAT", B="CTT", C="AATCGGGTGT" ) )
sub.seq <- ori.seq[ c(1, 3) ] # 1番目と3番目の配列を取り出す。
sub.seq
A DNAStringSet instance of length 2
width seq names
[1] 9 ATCGGCAAT A
[2] 10 AATCGGGTGT C
配列をXStringオブジェクトとして取り出したい場合は、[[<インデックス>]]とすればできます。
ちなみにこの場合で取り出せる配列は一つだけという点には注意してください。
ori.seq <- DNAStringSet( x=c( A="ATCGGCAAT", B="CTT", C="AATCGGGTGT" ) )
one.seq <- ori.seq[[ 1 ]] # 1番目の配列を取り出す。
one.seq
9-letter "DNAString" instance
seq: ATCGGCAAT
文字列ベクトルとして取り出したい場合は、as.character()が使えます。
ori.seq <- DNAStringSet( x=c( A="ATCGGCAAT", B="CTT", C="AATCGGGTGT" ) )
ori.vec <- as.character( ori.seq )
ori.vec
A B C
"ATCGGCAAT" "CTT" "AATCGGGTGT"
配列のマトリックスとして表示したい場合は、as.matrix()が使えます。しかしこの場合はXStringSetオブジェクト内の全ての配列長が同じである必要があります。
ori2.seq <- DNAStringSet( x=c( A="ATC", B="CTT", C="AAT" ) )
ori.mat <- as.matrix( ori2.seq )
ori.mat
[,1] [,2] [,3]
A "A" "T" "C"
B "C" "T" "T"
C "A" "A" "T"
このように手軽に変換できるので、処理方法に策定に行き詰まったら一度変換して処理することも考えてみてください。ただ文字列ベクトル等の基本型で処理した場合はやや処理速度が遅くなると思うので、大規模な配列を処理する場合はまず専用のパッケージがBioconductor内に存在しないか探してみるのが良いかと思います。
#さいごに
XStringSetオブジェクトを使うことで複数の配列を能率的に処理することができます。ぜひ利用してみてください。
次回はXStringオブジェクトとXStringSetオブジェクトを引数とする諸関数を紹介していきたいと思います。