はじめに
２年ぐらい前からruby-htslibというものを作っています。
これは、Rubyでゲノム関連のファイルフォーマットを扱うライブラリです。RやPythonのライブラリを活用すればいいのですが、右も左もSAM、VCFもよくわからない状態なので、htslibバインディングを自作しながら勉強しようと考えました。
Rubyは日本語コミュニティが強力で、言語仕様から、コア開発まで、日本語で質問でき、的確な回答が得られるのがいいところです。他の言語では日本語ネイティブのコア開発者が少ないのでなかなかそうはいきません。
Ruby-FFI
まずはFFIでC言語のネイティブ関数を呼べるようにしていきました。これは、そこまで難しいことではありませんでした。
しかし、問題はこのあとです。htslibは、一般的なライブラリよりもやたら難しいのです。最初のうちは「なぜ難しいのかわからない難しさ」を感じていました。
この記事では、難しく感じたことをリストして記事にしていきます。（まとまりのない記事で、個人のメモのようなものになっています）
ファイルの仕様を理解するのが難しい
いきなりですが、ファイルの仕様を理解するのは難しかったです。今もよくわかっているとは言えません。英語の仕様書を読む必要があります。日本語の解説はありません。hts-specはTexやPDFの形式で配布されており、改行が機械翻訳の利用を阻みます。「必要な箇所のhts-specを地道に英語で読む」「普段からゲノムのファイル形式によく触れて慣れる」「他言語のhtslibバインディング実装を使って学ぶ」を地道に続けようと思いました。
また英語であれば、クリアカットに物事が説明されているかと言うと、必ずしもそうではなくて、歴史的経緯などにより用語や概念が少々混乱していると思われる時があります。仕様書、実装、各言語のバインディングによって用語にバリエーションがあったりして、その点も悩ましいところです。しかし、これは日進月歩の科学という世界では、ある意味仕方がないことなのかもしれません。つまり右か左かわからない混沌が存在し、未来は見通せず、事後的にだんだんとはっきりとした形をとっていくという営みの本質を反映して、ライブラリやフォーマットもそのようなものになっているのかもしれないということです。
マクロ関数
最初の関門は、マクロ関数です。htslibでは非常に多くのマクロ関数が登場します。しかも、よく使いそうな処理に限って、マクロ関数が利用されています。マクロ関数は共有ライブラリから呼び出すことができません。言い換えると、Ruby-FFIではマクロ関数を呼び出すことができません。だからRubyで処理を再実装する必要があります。（Rust-htslibでも、マクロ関数はRust言語で再実装が行われています。）FFIではなくC拡張を使えばマクロ関数を呼び出すことができます。しかし、C拡張は採用せずFFIを使い続けました。（今でもhtslibでC拡張を採用しなかった方針は正しかったと思っています。）
なぜマクロ関数が多用されているのか？推測ですが、実行速度を少しでも速くするためだと思います。htslibは巨大なゲノムデータを対象にしたライブラリです。C言語やC++言語から呼び出す目的で作成され、PythonやR、あるいはRustなどの他の言語から呼び出すことは最初の目的ではなかったのだと思います。
ビットフィールド
htslibは、ビットフィールドを多用しています。ビットフィールドは整数型をさらに細かく分割して、フラグの格納場所などとして活用するCの手法です。仕方がないので、RubyのレベルでRuby-FFIを拡張し、ビットフィールドに対応させるライブラリを自作しました。ビット演算が全くわからなかったので、自分の頭で考えつつ、ruby-jpのSlackでも質問して何名かの方にスニペットを書いてもらいました。その成果は分離してffi-bitfieldというgemになりました。今考えると、ビットフィールドを避けてhtslibバインディングを作成することも不可能ではなかったと思います。しかし、やはりビットフィールドに対応していた方が便利であるのと、海外でhtslib以外でもffi-bitfiledを利用している方がいらっしゃるようなので、作って良かったなと思っています。
Enumerableの利用
初期の開発目標は、Bamファイルを読み込んで each メソッドを呼び出すことでした。しかし、Enumerableの使い方に関する理解度が低くて、いくつかの困難に突き当たりました。最初はPythonのhtslibバインディングを参照にしていたのですが、Pythonは外部イテレータ、Rubyは内部イテレータを利用するため実装に違いがありました。そのほか、Rubyのイテレータは一回実行したあと、必ずしも最初に戻る必要はない（StringIO#each)ことをmrknさんから教わりました。（現在は不完全ながらtellやseekを実装しており、元の位置に戻る手段は多少確保しています。その際にBGZFについての考え方を学びました。）__enum__ を使ってブロックがないときはEnumeratorを返却するような書き方があることや、Rubyのイテレータではメソッドチェーンを実現できるように最後にselfをReturnするべきだとか、そういったことを学びました。
メモリの確保・解放の問題
なんとかBam#eachを実装して、イテレータを回してRecordを取得することに成功したのですが、速度が遅くてがっかりしてしまいました。このときのショックでruby-htslibの開発は一時期停滞することになりました。あとから見ると、Samtoolsが速すぎるだけで、実際にはそこまで遅かったわけでもないのですが…。Bam#eachがうまくいっていなかった理由はいくつもあります。たとえば、Record構造体のメモリーのリリースが行われていませんでした。GC時にメモリの解放を行う関数を登録していなかったので、メモリーリークが発生しメモリの使用量が爆増していました。まずは、Bam1, Bcf1, といった構造体に関してはFFI::StructのかわりにFFI::ManagedStructを利用して、メモリ解放用の関数を登録するようにしました。（本当は全ての構造体についてそのようにするべきですが、今のところメモリリークが問題になりそうな一部の構造体しか対応できていません。
ちなみにRuby-FFIのManagedStructの実装は、GCが発火したときにコールされるProcを登録するRubyのメソッドObjectSpace#define_finalizerを利用してメモリを開放します。C言語レベルではなく、Ruby言語レベルでの実装です。サーバー上でruby-htslibを動かし続ける用途は想定していないのですが、サーバー上で大量のhtsファイルを開いたり閉じたりして動かし続けると、メモリの問題が発生するかもしれません。もちろんrubyコマンドが終了すれば、確保されたメモリは全て開放されます。）
マルチスレッド
遅かった他の理由の一つはCPUのマルチスレッドに対応していなかったことです。これについては、こちらの記事を参考にし、さらに hts_set_threads という関数を利用すれば良いことがわかったので、今ではスレッド数を指定するとこの関数を呼び出されるようになっています。
しかし、この話には続きがあります。hts-nimで hts_set_threadsが使われていない理由がわからなかったのですが、どうやらhts_set_threadsには問題があるようです。実際にhts_set_threadsを設定するとruby-htslibがうまく動作しないケースもありそうです。十分に調査できていません。
Crystal言語のバインディングの開発
Crystal言語ってみんなご存知でしょうか？ Rubyとよく似た文法でめちゃくちゃ速い言語です。
ruby-htslibは、私が作っていなかったら今でも全く存在しなかったとかなり自信をもって断言することができます。しかし、Crystal言語についてはわかりません。この言語は成長中で、誰かがhtslibのバインディングを先に作ってしまうかもしれません。（それはそれでとても良いことなんですけどね！）なので、自分でCrystal言語のバインディングを作成してみることにしました。Crystal言語では、crsytal_lib というツールがあって、これである程度バインディングを自動生成することができます。しかし、生成されたファイルをそのまま利用しているわけではなく、手直しは必要です。特にビットフィールドには対応していないので手動での修正が必要です。Crystal言語のバインディングを作っていてわかったのは、Crystal言語では、本当に相当のレベルまでRubyとの互換性を確保することができること、驚くほど高速だということです。（ただし、メタプログラミングを使わないような標準的な書き方のRubyの話です。Crystalでも、マクロを書くことにより、Rubyのメタプログラミングに相当することはかなりの程度まで可能です）
構造的に逐次的にデータを展開する関数
この部分はあまり理解できていないのですが、Bcfでは、構造体に必要な分だけデータを解凍して展開する関数 bcf_unpack というものが存在します。フィールドの値が知りたくなったら、単に構造体のフィールドを呼び出すだけでは不十分、LibHTS.bcf_unpack(@bcf1, LibHTS::BCF_UN_INFO) みたいなことをして展開してからフィールドを呼ばなければならないのです。いままでそういうコードを見たことがなかったので理解するのに時間がかかりました。
これもマクロ関数と同じで、必要のないフィールドは展開しない。処理を減らして高速化する、そういう工夫なんだと思います。結局のところ、htslibの難しさの正体は、高速化のために初心者が見慣れない手法をいろいろ行っているからということになりそうです。
vcf-bench に挑戦・Recordのメモリの使いまわし
Brentp氏は、多種多様な言語のhtslibバインディングを作成されている国際的に著名なバイオインフォマティシャンです。VCFの呼び出しのベンチマークを公開しています。実行してみたら、hts-nimよりもかなり速度が遅いことがわかりました。原因をつきとめるために、人生ではじめてプロファイラというものを起動し、KCachegrindなどのGUIツールを使って調べました。その結果、Nim言語ではRecord構造体を使いまわしていたのに対して、Crystalでは毎回アロケートしている。これが良くないということが判明しました。つまり、高速化するためには、繰り返し処理のたびにメモリを確保するのはやめて、毎回同じメモリを使いまわせばいいということになります。これはRubyだと、eachで、同一のオブジェクト（または同じ構造体オブジェクトを参照する別のオブジェクト）をyieldするということを意味しており、意図しない挙動の原因になります。しかし、値のような即値はコピーされるため、そのようにしても map も含め大抵の場合はうまくいくということがわかりました。デメリットよりも高速化のメリットの方が大きいだろうと判断して、eachでは構造体を使いまわし、each_copyというメソッドを別に用意することにしました。ひょっとすると逆の方がいいかもしれません。というのは、each は Enumerable の他のメソッドでも使い回されるため不具合がでるメソッドもあるのではないかと思うからです。（どこで不具合が出そうかはあまり追究できていません）しかし、今のところはeachは使いまわしということになっています。これにより大幅に速度が向上し、手元のPCではhts-nimよりも高速な結果を出すことに成功しました。しかし、その差はわずかであり、私の手元のコンピュータで測定すると常にCrystalの方がほんの少しだけ高速ですが、環境によってはNimの方が早い結果になることもあるようです。
深いコピー
本当はRubyのオブジェクトはcloneやdupに対応しているべきです。その場合、対応するC言語の構造体もコピーする必要があります。
Crystalの型システム、リターンする型が問題
Crystalはほとんどの場面でRubyとの互換性を保つことができます。しかし、一箇所、あまり予想していないところで問題が発生しました。それは戻り値の型です。Ruby言語の場合は、戻り値にいろいろな型を取るのは普通のことです。その方が便利だからです。しかしながら、Crystal言語では戻り値の型が一意に決まっていないと非常に面倒です。これによって、RubyとCrystalでベストなAPIが変わってしまう傾向があります。たとえば Rubyでは tag というメソッドが、Crystalでは戻り値によって tag_int tag_float などと分割しておいたほうが、ユーザーにとっては使いやすくなると感じます。またNilを戻り値の可能性に含めるかということも問題になります。（これらについてはCrystal言語の経験が浅いため何がベストプラクティスなのかよくわかっていません。）
Bam::Tag, VCF::Info, VCF::Format の読み込みと書き込み
Bam、VCFなどのファイルフォーマットは、レコードに任意の情報を書き込める仕組みを用意しています。それが、BamであればTagであり、VCFであれば、InfoやFormatといったものです。Samtoolsのようなコマンドラインツールではなく、わざわざhtslib使うような人は、これらのTag、Info、Formatをフル活用したいという人が多く、これらにある程度対応しておかないと、htslibバインディングの価値が半減してしまいます。ruby-htslib, htslib.cr これらの追加フィールドの読み込みに関してはそこそこ対応しているつもりです。一方で、書き込みに関してはまだ対応できておらず、今後の課題となっています。
- VCF::Infoはすべてのサンプルで同一のもの。
- VCF::Formatはサンプルごとに異なるデータを入れる。
画像出典：https://upload.wikimedia.org/wikipedia/commons/3/39/Binary_BCF_versus_VCF_format.png
これらのデータは普段は構造体に展開されないので、利用するときは、bcf_unpackを呼ぶ必要があります。また、infoやformatを多用するVCFファイルでは、それに応じてヘッダーの部分も長大になる傾向があり、それらを読み書きできる仕組みも必要ですが、これはあまり対応できていません。
CRAM対応
当初はCRAMをBamと別のクラスにしようと考えていましたが、htslibの構造上、Bamに対応していると、Cramも自動的に対応するため結局Bamクラスにまとめる方針となりました。CRAMの問題点はFastaファイルの絶対パスを必要とすることです。CRAMはFASTAの存在を前提にした圧縮形式であるため。常にFASTAを必要とします。そのためCRAMはポータブルな形式になりにくく、ファイルの欠点となっています。このためGithubActionsでのテストを通すのに少し苦労しました。
VCF、Format の GTは String って書いてあるのに本当は Integer
これは現在対処中で全てを把握しきれていないのですがワナですね。
GTはHeaderでStringと宣言されていても、特別扱いでIntegerとして処理されています。bcftoolsやhtslibのコードを見ても、GTは完全に特別扱いで、GTを検出した場合は別の処理を行わせるようになっています。BcfではGTをIntegerで扱うので、それに合わせて実装したものが使われているようです。
わかっていないこと・今後の課題
大量にあります。むしろわかっていないことだらけで、わかっていることの方がはるかに少ないのが実情です。
いろいろ学ばせてもらっています。Ruby-htslibとhtslib.crに関しては、少しずつテストや機能を追加していきたいと考えています。また、htslib側で処理が失敗した場合、関数は -1 などの値を返します。これらの値を監視して、Rubyのエラーを発生させるべきですがサボっています。一方で、OSSでは単にメンテナンスを継続するだけでも一定のエネルギーと時間が必要なので、とりあえず開発が継続できたらいいなとも思っています。
夢
いつかファイルフォーマットの誕生に携わりたい。この段落は単なるポエムです。いつの時代も人々は最高の方法で、複雑なものを構築します。htslibもそのなかの一つです。しかし、時間が経つともっとよりよい方法が発明されて、これまでのやり方が時代遅れになっていきます。そのような場合でも、まずファイルフォーマットがあり、低レベルの言語のライブラリがあり、それが次第に高レベルな言語になっていくという構造は変わらないと思います。いつの時代も、新しいファイル形式とツールが同時に開発されます。はるか遠い夢になりますが、いつか新しいファイルフォーマットの誕生に立ち会い、何らかの貢献をできたらいいなと願っています。
この記事は以上です。
まとまりのない記事ですが、こういった頭の中でぼんやり考えているような情報は書き残しておかないと消えてしまうので無理やり記事にしました。