Hatena::Grouplifesciencedb

ゲノム周辺 このページをアンテナに追加 RSSフィード

2009-05-16

tc-subseq、もしくは Tokyo Cabinet でヒトゲノム配列の部分配列取り出しをしてみたら速くて驚いた

|  tc-subseq、もしくは Tokyo Cabinet でヒトゲノム配列の部分配列取り出しをしてみたら速くて驚いた - ゲノム周辺 を含むブックマーク はてなブックマーク -  tc-subseq、もしくは Tokyo Cabinet でヒトゲノム配列の部分配列取り出しをしてみたら速くて驚いた - ゲノム周辺  tc-subseq、もしくは Tokyo Cabinet でヒトゲノム配列の部分配列取り出しをしてみたら速くて驚いた - ゲノム周辺 のブックマークコメント

id:dritoshi 氏がほ乳動物ゲノムの部分配列を頻繁にとりだすソリューションを探していたのと、そのような配列解析の周辺事情に興味があったので、いろいろ調べてみた結果いつのまにか Tokyo Cabinet をつかう高速なツールを作っていました。

ほ乳類ゲノム配列、たとえばヒトゲノム配列の場合では、長さは30億塩基染色体数は23個。染色体では、最長の一番染色体で219,475,005塩基、最小の22番染色体で49,691,432塩基というスケールです。これらから任意の部分配列を取り出すのは、各種ウェブサービスを利用すれば容易いことです。インタラクティブにはすぐ解決する類いの問題です。しかし、3万個くらいの遺伝子の各エクソンごとに周辺の部分配列を10個づつくらい取り出してみたいということを考えると、90万クエリを処理する計算になります。ウェブサービスはつかうのはつらそうです。

となると、ローカルに解析環境をととのえることになります。いろいろ手段やツールがありそうですが、実はスケールしないものが多かったという話です。


dbxfasta と seqret は遅過ぎた

fasta フォーマットの配列データベースを作成するとしたら、とりあえず EMBOSS の dbifasta や dbxfasta が挙げられます。EMBOSS の dbxfasta は B+tree をつかっているので、高速にエントリ取得が可能なインデクスを作成できます。しかし、今回の問題設定では、ヒトゲノムの場合でも23本の染色体分のエントリ数なので、エントリ取得に関してはインデクスはいくら高速でも寄与しません。

dbxfasta で作成したインデクスからのエントリ取り出しは seqret で行う。seqret の USA 引数に部分配列を指定することができる。たとえは、chr1 の 10..100 は dbname:chr1[1:100] という USA で指定できる。ここでの部分配列のとりだしは、計算上では、エントリの全体を取得してから切り出しているので、巨大な chr1 の場合は、エントリ全体をコピーするといった無駄な計算が行われていました。


fastacmd は ID 処理が複雑だった

NCBI Tools の formatdb と fastacmd は NCBIgi のついたファイルのだと手軽に扱えますが、それ以外のものは、思ったように動かすのが難しく、うまくいきません。きっと高速に検索できるのでしょうが、データの準備に手間取ってしまいます。


MySQL に格納するには不利だった

UniProtKB/SwissProt 程度のサイズなら、配列エントリの部分文字列マッチや部分文字列のとりだしは MySQL にテキストとして格納しても実用的に利用できます。text なら SUBSTRING(seq, offset, length) で取り出せます。SwissProt 程度のサイズとは、一レコード配列長がせいぜい数百文字長くとも数万文字という範囲です。ヒトゲノム配列の場合は、一レコードが長いので、そもそもそのままではテーブルに格納できません。


jksrc にあると思っていた

nibFrag が nib フォーマットの配列データから任意の領域をとりだすツールです。あとで調べるてみます。


Tokyo Cabinet のハッシュデータベース(TCH)で十分だった

巨大データ関係で注目していた Tokyo Cabinet のハッシュデータベースをつかって UCSC の hg18 の chromFa.zipダウンロードからデータインポートまでをおこなうスクリプトと部分配列取得スクリプトgithubtc-subseq としておきました。

データのダウンロードUCSC から hg18 chromFa.zipダウンロードし、展開します。 

$ rake data:ucsc:hg18

格納するデータは、塩基配列ファイル(chr1.faなど)が50文字毎に改行がはいっているので、50文字分を値、その行番号からキーを作成し、染色体名_行番号をキーとして配列ファイル(chr1.faなど)から生成しました。

chr1_1	NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN

Tokyo Cabinet のハッシュデータベースインポートに時間がかかったからインデクスファイルを分割した

インデクスファイルを一個にすると、インポートが一晩かかっても終わらなかったので断念。染色体ごとにわけることにしました。この場合はキーから染色体名を取り除いて行番号だけにします。

1	NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN

スクリプトでは、中間tsvファイルは生成せずに、tokyocabinet-ruby 経由で直接インデクスを作成します。

$ time rake tch:create
real 13m34.084s
user 6m16.345s
sys 6m1.639s

rake tch:create 480.11s user 438.57s system 46% cpu 32:35.06 total

インデクス作成は十分くらいで終わりました。

ベンチマークをとってみます。

秒速4827クエリ

$ rake benchmark:tch:chr1
                   user     system      total        real
14928 gets:    1.890000   1.150000   3.040000 (  3.136869)

部分文字列の長さに依存してるようです。

$ rake benchmark:tch:length
                       user     system      total        real
TCH   1 nt/440:    0.030000   0.030000   0.060000 (  0.063573)
TCH  10 nt/440:    0.020000   0.030000   0.050000 (  0.061263)
TCH 100 nt/440:    0.060000   0.040000   0.100000 (  0.094042)
TCH  1k nt/440:    0.090000   0.030000   0.120000 (  0.143346)
TCH 10k nt/440:    0.410000   0.050000   0.460000 (  0.501937)


Tokyo Cabinet の固定長インデクス再び

キーが数字だけなら固定長インデクス(TCF)で良いわけで、TCF をつかってみたら、秒速7000クエリくらいでました。

インデクス作成

$ time rake tcf:create

高速な文字列取得

$ time ruby tcf_subseq chr1:1000,1010

ベンチマークをとってみます。

一番染色体のなかで1000塩基未満の適当な組み合わせの問い合わせ

$ rake benchmark:tcf:chr1
                   user     system      total        real
14928 gets:    1.570000   0.940000   2.510000 (  2.787780)

およそ秒速7000クエリです。


取り出す塩基配列長ごとの比較

$ rake benchmark:tcf:length
                       user     system      total        real
TCF   1 nt/440:    0.000000   0.000000   0.000000 (  0.000108)
TCF  10 nt/440:    0.000000   0.000000   0.000000 (  0.000105)
TCF 100 nt/440:    0.000000   0.000000   0.000000 (  0.000093)
TCF  1k nt/440:    0.000000   0.000000   0.000000 (  0.000156)
TCF 10k nt/440:    0.000000   0.000000   0.000000 (  0.000116)

取り出す長さはあまり影響していないようです。


nibFrag が速かった

UCSC の jksrc のなかの nibFrag が高速でした。nib フォーマットは1バイトで塩基配列二文字表現するバイナリフォーマットです。そこから高速に部分配列をとりだすのがnibFragです。なからずファイルを生成してしまうのですが、秒速300クエリほどです。似たような条件として、Tokyo Cabinet の tcfmgr でファイル生成をしながら実行すると、秒速220クエリでした。

chr1.faファイルからchr1.nibファイル生成に3秒ちょっと。

$ faToNib chr1.fa chr1.nib
Writing 245203898 bases in 122601957 bytes
3.72 real 2.59 user 0.84 sys

nibファイルは元ファイルの半分くらいのサイズ。

$ du -sh chr1.nib 
117M chr1.nib
$ du -sh chr1.fa
239M chr1.fa

部分配列の取り出しは高速

$ time nibFrag chr1.nib 1 2 + out.fa

real 0m0.002s
user 0m0.000s
sys 0m0.002s

nibFrag の引数は offset と offset+length です。offset は塩基配列の座標ではなくて、座標の整数列の間に相当します。たとえば、"ACGT" から G に相当するのは、offset = 2 です。2 なのは 0 オリジンだからです。


まとめ

  1. 簡単そうな問題ですぐに使えるツールがない場合は、簡単に解決できる問題だった。だからツールが無い。グラフと似ている。
  2. tc-subseq は二時間程度で実装が終わった。GitHub - nakao/tc-subseq: human genome on tokyo cabinet for sub-sequence extraction
  3. ヒトゲノム全体を一つのインデクスにするのは、インポート時間がかかりすぎるのでおすすめできない。
  4. nibFrag が高速だった。
  5. TCからnibデータにアクセスできると最強かとおもった。だれかたのむ。
トラックバック - http://lifesciencedb.g.hatena.ne.jp/nakao_mitsuteru/20090516