VelvetKを使ってVelvetを用いたアセンブルに最適なk-merサイズを推定する

(追記, 2013/07/18, Velvetに依らないk-mer推定の記事：ゲノムアセンブリにおいて最適なk-merを推定するKmerGenieを試してみた - Wolfeyes Bioinformatics beta)

VelvetKというスクリプトが公開されたらしい．これはVelvetのアセンブルの際にユーザが指定する必要のあるk-merのkの値を，ゲノムサイズとショートリードのサイズから自動推定するスクリプトのようだ．普通ならばkの値を細かく変えてVelvetを大量に走らせてアセンブル結果を評価するというのが定石だと思うが，もしkの値が自動推定できるなら，そのkの値付近だけを重点的に調べるといったことが出来るので大幅に労力を削減できる．私自身あまりVelvet/Oasesを使った経験が無いのだが，取り敢えず一通り使ってみた．

インストール

インストールは簡単で，Perlスクリプトを落としてくるだけで実行出来る．外部ライブラリの依存が無いのでお手軽．

1

$ wget http://bioinformatics.net.au/velvetk.pl

使い方

VelvetKでkの推定値を算出するには，

  • ゲノム配列 または ゲノムサイズ
  • アセンブルの元になるショートリード

の2つの入力が必要になる．ゲノム配列が既知の場合には--genomeオプションを使って，fasta形式のゲノム配列のファイルを指定する．ショートリードはfastq形式などのファイルを指定することになるが，gzやbz2などの圧縮形式にも対応している．ゲノム配列を指定してVelvetKを実行するには，以下のようにコマンドを実行する．

1

$ perl velvetk.pl --genome chr.fasta R1.fastq R2.fastq

一方でde novoのようなゲノム配列が無い場合には，対象種のおおまかなゲノムサイズを入力することになる．ゲノムサイズの指定の場合には，k/M/Gといった接頭辞を使う事ができる．ゲノムサイズを数値で指定してVelvetKを実行するには，以下のようにコマンドを実行する．

1

$ perl velvetk.pl --size=1M chr.fasta R1.fastq R2.fastq

とする．

擬似データで試してみる

今回はお手製の擬似データで試してみることにする．

擬似データの作り方は簡単で，乱数を使って擬似的に作ったゲノム配列からこれまた乱数を使って擬似的にショートリードを作成するというもの．ゲノムの長さを1Mbとして，100bpのペアエンドの配列を500,000本用意して，カバレッジ100xくらいのゲノムシーケンスを想定した．

その擬似データセットで実行してみると以下のような出力が得られた．実際にはショートリードの長さや本数も不均一になるので，以下のように各値が綺麗な感じにはならないと思う．

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

$ perl velvetk.pl --genome chr.fasta R1.fastq R2.fastq
Estimating target genome size from 'chr.fasta'
Target genome size is 1000000 bp
Desire k-mer coverage of 25
Using cat for R1.fastq
Read 500000 sequences from R1.fastq
Using cat for R2.fastq
Read 500000 sequences from R2.fastq
Considered 1000000 reads with lengths 100..100 bp
K       #Kmers  Kmer-Cov
1       100000000       100.0
3       98000000        98.0
5       96000000        96.0
7       94000000        94.0
9       92000000        92.0
〜〜〜中略〜〜〜
71      30000000        30.0
73      28000000        28.0
75      26000000        26.0
77      24000000        24.0
79      22000000        22.0
81      20000000        20.0
83      18000000        18.0
85      16000000        16.0
87      14000000        14.0
89      12000000        12.0
91      10000000        10.0
93      8000000 8.0
95      6000000 6.0
97      4000000 4.0
99      2000000 2.0

入力ファイルをそれぞれ読み込んだ後，それぞれのkの値に対して#kmersとKmer-Covが表示される．#kmersはショートリードのデータからカウントできるk-merの本数で，例えばK=99の場合では100bpから抜き出せる99merは2つあるので，今回の場合リード数500,000*2の合計2,000,000という感じだろう．次のKmer-Covは，いわゆるK-mer Coverageと呼ばれるもので，実際のカバレッジとk-merで見た時のカバレッジの割合，もしくは変換するときの係数に当たる．

ではどうやってこのリストからVelvetに最適なkの値を推定するかというと，出力の始めの方に書かれている

1

Desire k-mer coverage of 25

のところを見て，Kmer-Covが25を下回った時のkの値が最適なkの値とする．今回の場合だとk=77でKmer-Covが24.0になるので，k=77が最適なk-merのサイズとなる．

それではなぜ25なのかというと，このスクリプトの場合だとデフォルトだからということになってしまうのだが，どうやらk-mer Coverageにはだいたい10~20あれば十分という知見(多分経験則)があるらしい．一般的にはkの値が小さいとコンティグが繋がりにくく，kの値が大きいと細かな違いをコンティグに反映出来無いといったsensitivityとspecificityのトレードオフがある．それらのいい具合の中間点ということで，k-mer Coverageを見て判断しようということらしい．

また，このスクリプトでは閾値となるk-mer Coverageを変更するオプション--covや，リストの出力を抜きにして最適なKの値だけを出力するオプション--bestなども用意されている．

1
2
3
4
5
6
7
8
9

$ perl velvetk.pl --best --cov=20 --size=1M R1.fastq R2.fastq
Target genome size is 1000000 bp
Desire k-mer coverage of 20
Using cat for R1.fastq
Read 500000 sequences from R1.fastq
Using cat for R2.fastq
Read 500000 sequences from R2.fastq
Considered 1000000 reads with lengths 100..100 bp
81


上の場合では，k=81が最適とおもわれるkの推定結果となる．

実はWeb版もある

ちなみに，ここまで書いておいて難だが，実はVelvetKと同作者が作った簡易版「Velvet Advisor」がWebアプリケーションとして既に存在している．ただ，この簡易版は複数の入力ファイルに対応していないので，使えるケースが制限されている．普通は複数のライブラリを読んでアセンブルするとと思うので，まあ今回のVelvetKで使いやすくなったという感じだろうか．

http://dna.med.monash.edu.au/~torsten/velvet_advisor/

まとめ

今回は擬似データを使ってVelvetKを動かしてみたが，これからアセンブルをするという人にとってはパラメータの大雑把な目安となって良いと思う．どれだけの精度で当たるのか不安なところもあるので，実際にVelvetを使ったアセンブルのプロジェクトでVelvetKの推定値がどれだけ当たっているかが知りたいところだけれども，結局目安としてしか使わないのだから，まあ適当でいいんだろう．

参考資料