ゲノムアセンブリにおいて最適なk-merを推定するKmerGenieを試してみた

KmerGenieという，シーケンスデータのみを用いてアセンブラに依存しない形でアセンブルに指定する最適なk-merの値を求めることができるソフトウェアが最近出てきたので，ちょっと使ってみた．

KmerGenie

そういえば以前VelvetKという，これも同じようにk-merパラメータをデータセットから推定するツールについて書いたが，これは名前の通りVelvetに特化したソフトウェアだった．2007年くらいに公開され使われ始めたVelvetは，NGSの普及と相まってかなり長い間ゲノムアセンブラとして広く使われてきた印象がある．その後ABySSやMIRA，ALLPATHS-LG，Rayなどの様々なゲノムアセンブラが出てきたが，それらに対してVelvetKのようなk-mer最適化のようなツールは無かったようで，そのあたりアセンブラに非依存というのがKmerGenieの売りのようだ．

KmerGenieの仕組み

KmerGenieがどういう仕組みで最適なk-merを推定しているかは，以下のスライドやarXiv.orgで公開されている論文に詳しく書かれている．

これによると，単純に言ってしまえば最適なk-merの出現頻度というものはきれいな正規分布に従うと仮定した上で，k-merの値を色々変えて出現頻度のヒストグラムを作成し，生成モデルを立ててフィッティングをして最適なk-merの値を推定しているようだ．実際には，正規分布を仮定したゲノム配列の分布とパレート分布を仮定したエラーの配列の分布が混ざった混合分布を考えて，その混合の重みなども考慮したパラメータ最適化をしているらしい．

KmerGenieのインストール

まずはKmerGenieをインストールする．ダウンロードしたファイルを展開すると既にディレクトリにkmergenieという実行ファイルがあるが，makeをしないと利用することができないので注意．

1
2
3

$ wget http://kmergenie.bx.psu.edu/kmergenie-1.5397.tar.gz
$ cd kmergenie-1.5397
$ make

KmerGenieの実行

それでは実際にfastqファイルに対してKmerGenieを実行してみる．

1. ゲノムアセンブリのチュートリアルで使用されたE. coliのゲノムシーケンス

まずは2013年の6月に行われたMSU NGS Summer Course 2013のチュートリアルで使用されたサンプルデータを使ってみる．これはE. coliのゲノムシーケンスで，元はChitsaz et al.が出したデータの一部のようだ．中身は約4.7Mの70bpのシングルリード．

Assembling E. coli sequences with Velvet — ANGUS 2.0 documentation

1
2

$ curl -O https://s3.amazonaws.com/public.ged.msu.edu/ecoli_ref-5m-trim.fastq.gz
$ gunzip ecoli_ref-5m-trim.fastq.gz

それではecoli_ref-5m-trim.fastqに対してKmerGenieを実行する．実行後に作られる各k-merでのカウントデータや分布をプロットしたpdfはカレントディレクトリにそのまま出力されるので，もし必要ならディレクトリを別途作ってそこで実行したほうが良いだろう．

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
32
33
34
35

$ kmergenie ecoli_ref-5m-trim.fastq
running histogram estimation
Linear estimation: ~130 M distinct 41-mers are in the reads
K-mer sampling: 1/100
| processing                                                                                         |
[going to estimate histograms for values of k: 61 51 41 31 21 
-----------------------------------------------------------------------------------------------------------------------------Total time Wallclock  77.7066 s
fitting model to histograms to estimate best k
fitting histogram for k = 21
fitting histogram for k = 31
fitting histogram for k = 41
fitting histogram for k = 51
fitting histogram for k = 61
estimation of the best k so far: 51
refining estimation around [45; 57], with a step of 2
running histogram estimation
Linear estimation: ~139 M distinct 39-mers are in the reads
K-mer sampling: 1/100
| processing                                                                                         |
[going to estimate histograms for values of k: 57 55 53 51 49 47 45 
-----------------------------------------------------------------------------------------------------------------------------Total time Wallclock  56.4315 s
fitting model to histograms to estimate best k
fitting histogram for k = 21
fitting histogram for k = 31
fitting histogram for k = 41
fitting histogram for k = 45
fitting histogram for k = 47
fitting histogram for k = 49
fitting histogram for k = 51
fitting histogram for k = 53
fitting histogram for k = 55
fitting histogram for k = 57
fitting histogram for k = 61
table of predicted num. of genomic k-mers: histograms.dat
best k: 55

今回の場合，KmerGenieによる最適なk-merの値は55となった．出力結果で示されるk-mer頻度のプロットを並べて見ると，確かにkが小さいときは少しいびつになっており，k=55あたりで一番正規分布っぽくなっている…??ことがわかる．

ただ，このサンプルーデータを解析しているチュートリアルではk-merを31〜35にしてvelvetでアセンブルしていることから，今回の55という推定値は少し大きめのような気がする．

2. Assemblathon2で使用されたセキセイインコのゲノムシーケンス

次はもう少し大きいデータで試してみよう．Assemblathon2というアセンブリツールの性能を競うコンペで使用されたセキセイインコのゲノムシーケンスの中からDuke Illumina GAIIx runsというデータを使ってみる．中身は76bpのペアエンドで合計で126M readsある．

Index of /Data/hcbxz0i7kg/Parrot/illumina_duke_runs

1
2

$ lftp -e "mirror" http://bioshare.bioinformatics.ucdavis.edu/Data/hcbxz0i7kg/Parrot/illumina_duke_runs/ 
$ cat *.fastq > budgie.fastq

Assemblathonのサイトによるとセキセイインコのゲノムは二倍体らしいので，今回はKmerGenieの「–diploid」のパラメータを付けて実行する．

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
32
33
34
35
36
37

$ kmergenie --diploid budgie.fastq
running histogram estimation
Linear estimation: ~3721 M distinct 46-mers are in the reads
K-mer sampling: 1/1000
| processing                                                                                         |
[going to estimate histograms for values of k: 71 61 51 41 31 21
---------------------------------------------------------------------------------------------------------------------------Total time Wallclock  3035.49 s
fitting model to histograms to estimate best k
fitting histogram for k = 21
fitting histogram for k = 31
fitting histogram for k = 41
fitting histogram for k = 51
fitting histogram for k = 61
fitting histogram for k = 71
estimation of the best k so far: 21
refining estimation around [15; 27], with a step of 2
running histogram estimation
Linear estimation: ~6335 M distinct 24-mers are in the reads
K-mer sampling: 1/1000
| processing                                                                                         |
[going to estimate histograms for values of k: 27 25 23 21 19 17 15
---------------------------------------------------------------------------------------------------------------------------Total time Wallclock  2305.61 s
fitting model to histograms to estimate best k
fitting histogram for k = 15
fitting histogram for k = 17
fitting histogram for k = 19
fitting histogram for k = 21
fitting histogram for k = 23
fitting histogram for k = 25
fitting histogram for k = 27
fitting histogram for k = 31
fitting histogram for k = 41
fitting histogram for k = 51
fitting histogram for k = 61
fitting histogram for k = 71
table of predicted num. of genomic k-mers: histograms.dat
best k: 21

データ量が多いと計算時間もかなり大きくなるようで，出力にも少し書かれている通り，今回は1時間半くらいかかった．

さて，今回の場合はKmerGenieによる最適なk-merの値は21となった．出力結果で示されるk-mer頻度のプロットを並べて見てみると，kが大きい値の時には確かに正規分布として近似できないくらいなだらかな曲線になっているものの，K=21の最適値付近でもその様子はあまり変わっていないような感じがする．そもそも可視化の都合もあって横幅が揃えられていなかったりと色々と問題はあるが，とりあえずは今回のデータから推定された最適値は21だった．

KmerGenieで使用できるパラメータ

以上のように，KmerGenieは幾つかのk-merの値で実行して最適と思われる値を自動で出力するが，他にも手動でパラメータを指定することでk-merの範囲や刻みを変更することができる．

1
2
3
4
5
6
7
8
9
10
11
12
13
14

$ kmergenie
KmerGenie

Usage:
    kmergenie <read_file> [options]

Options:
    --diploid    use the diploid model
    --one-pass   skip the second pass
    -k <value>   largest k-mer size to consider (default: 121)
    -l <value>   smallest k-mer size to consider (default: 15)
    -s <value>   interval between consecurive kmer sizes (default: 10)
    -e <value>   k-mer sampling (default: auto-detected power of 10)"
    -o <prefix>  prefix of the output files (default: histograms)"

まとめ

このように，KmerGenieはシーケンスデータのみを用いてアセンブラに依存しない形で，アセンブルに指定する最適なk-merの値を求めることができる．推定された値が最良のアセンブルに繋がるかは実際に確認してみないと判断できないが，大まかな目安であったりk-mer決定プロセスの理由付けなどに力を発揮するだろう．今回は結果に対してそれほど詳細な評価を行わなかったが，arXiv.orgで公開されている論文の方ではGAGEのデータセットに対してアセンブル結果を含めた性能比較をしているので，そちらも参考にしていただきたい．

参考