wgsimのコードを眺める #ポエム

はじめに

お正月なので、Heng Liさんのコードでも眺めてみようと思って、wgsim というツールのコードを眺めています。このツールはリファレンス配列を与えると、ショートリードシーケンサーのシミュレーションを行います。13年前のツールで、全部で400行程度。複雑なアルゴリズムが使われているわけでもなさそうなので読みやすいのではないかと思って見始めました。完全に理解できたわけではないですが、いくつか勉強になるところがあったので記事に残します。

この記事は対象読者を想定していません
自分が1年後とか、2年後に読み返すことを想定して書いています。

定番のライブラリ klib

Heng Li さんのツールは、C言語で記述されて、そのライブラリとしてklibが使用されています。

klibはC言語にハッシュテーブル、動的配列、キュー、木構造、文字数付き文字列、ソートアルゴリズム、コマンドラインオプション解析、マルチスレッド機構など、FASTA/FASTQファイル読み込みなど、便利な機能を追加します。このライブラリも多くの部分は Heng Li さんによって作成されています。

世間を見ると、複雑なアルゴリズムを記述した高速なコードの作成ではC++を利用している人が多いように思います。しかしHeng Li製ツールはあくまでC言語とklibの最小構成でツールを作ります。

外部ライブラリである klib の使い方も大変面白いものです。Rubyであれば、Gemspec ファイルを記述することで外部のライブラリを使います。Rustであれば toml ファイルを記述するでしょう。C++であれば git submodule として外部のライブラリをプロジェクトに組み入れるかも知れません。しかし、klibでは、単にヘッダーファイルをプロジェクトの中にコピーして使います。シンプルですね。

wgsimで使われる kseq.h

wgsimでもkseq.hというFASTA/FASTQを読み込むライブラリが使われています。

#include <zlib.h>
#include <stdio.h>
#include "kseq.h"
KSEQ_INIT(gzFile, gzread)

fp_fa = gzopen(fn, "r");           // gzipで圧縮されたfastaを開く
ks = kseq_init(fp_fa);             // 初期化
tot_len = 0;                       // 全塩基の数
n_ref = 0;                         // 染色体の数
while ((l = kseq_read(ks)) >= 0) { // 染色体読み込む(戻り値は塩基数)
	tot_len += l;                  // 塩基数をインクリメント
	++n_ref;                       // 染色体数をインクリメント
    // 
    // ks->name.s
    // ks->seq.s
    // ks->comment.s
    // ks->qual.s
}

ASCII文字列を数値配列に変換　nst_nt4_table[256] の謎

Heng Li作品を見ているとよく登場するのが次の変換テーブルです。最初に見たときには意味がわからなくてぎょっとすると思います。

const uint8_t nst_nt4_table[256] = {
	4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 
	4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 
	4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 5 /*'-'*/, 4, 4,
	4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 
	4, 0, 4, 1,  4, 4, 4, 2,  4, 4, 4, 4,  4, 4, 4, 4, 
	4, 4, 4, 4,  3, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 
	4, 0, 4, 1,  4, 4, 4, 2,  4, 4, 4, 4,  4, 4, 4, 4, 
	4, 4, 4, 4,  3, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 
	4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 
	4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 
	4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 
	4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 
	4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 
	4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 
	4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 
	4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4
};

ほとんどは4ですが、ところどころに、0, 1, 2, 3, 5 があります。
途中に出てくるコメントの /*'-'*/ が顔みたいに見えて、なんとなくかわいいですね。

これは塩基配列を表現する文字を、数値に変換するためのテーブルです。例えば次のようなRubyのコードを実行すると、このテーブルの意味が明らかになります。

arr = [4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
       4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
       4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 4, 4,
       4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
       4, 0, 4, 1, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 4, 4,
       4, 4, 4, 4, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
       4, 0, 4, 1, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 4, 4,
       4, 4, 4, 4, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
       4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
       4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
       4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
       4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
       4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
       4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
       4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
       4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4]

arr.each_with_index do |value, index|
  puts "|#{index}|#{index.chr}|#{value}|" if value != 4
end

インデックス番号	対応する文字	変換後の数値
45	-	5
65	A	0
67	C	1
71	G	2
84	T	3
97	a	0
99	c	1
103	g	2
116	t	3

C言語の世界で、FASTA/FASTQを読み込むと、ASCII文字の配列のポインタとして扱われます。ASCII文字はchar型ですから、内部的には1バイトの整数です。uint8ですね。そこで、上記のテーブルを利用して、数値配列に変換してしまおうというわけです。

同じことを、Rubyでやる場合を考えてみましょう。

seq = "acgtnACGTN-"

seq.chars.map { |c|
  case c
  when 'a', 'A' then 0
  when 'c', 'C' then 1
  when 'g', 'G' then 2
  when 't', 'T' then 3
  when '-' then 5
  else 4
  end
}
# => [0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 5]

とか

seq = "acgtnACGTN-"

table = {'A' => 0, 'C' => 1, 'G' => 2, 'T' => 3, '-' => 5}
table.default = 4

seq.chars.map { |c| table[c.upcase] }
# => [0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 5]

こんな感じです。
しかし、このようなコードですと、大量のデータを処理するバイオインフォマティックスのツールとしては処理が遅くなってしまうのでしょう。

この nst_nt4_table の使い方は非常に簡単で、次のようなコードでmut_t型に変換します。この mut_t は uint16 なのです。

typedef unsigned short mut_t;

c = (mut_t)nst_nt4_table[char_code];

uint16配列の上位12ビットはメタデータが格納できる

さて、この時点で？？って思った人もいるかもしれません。C言語の世界では、ACCII文字列を8ビットの整数として扱います。そうであるならば、わざわざ16ビットの整数に変換する必要などないのでは？

その答えは僕にもはっきり分かりませんが、まず変換先を16ビットにするかどうかはともかく「変換」の必要性については、

a と A など大文字小文字を統一して扱いたいので変換は必要
a c g t - 以外の文字をすべて 4 に統一したいので変換は必要
- （FASTA/FASTQには、actgおよびn以外にも使える文字がいくつかあったはずです。）

ということになると思います。そして、変換先が8ビットではなく16ビットである理由ですが、主に拡張された領域にメタデータを格納しているためだと思います。wgsimの場合は、

下位4ビットに塩基の情報を格納
上位4ビットに deletion, 挿入塩基数などのメタデータ
中位8ビットに insertion 用の塩基データ

としています。挿入塩基とその塩基数のメタデータの埋め込みには、最小でも８ビットが必要であり、Uint8ではまかなうことができません。これが uint16 が利用されている理由の一つだと思います。

さて、上位４ビットには次のようなパターンが与えられています。

enum muttype_t {NOCHANGE = 0, INSERT = 0x1000, SUBSTITUTE = 0xe000, DELETE = 0xf000};

static mut_t mutmsk = (mut_t)0xf000;

名前	２進数での表記	意味
NOCHANGE	0000 0000 0000 0000	変異なし
INSERT	0001 0000 0000 0000	挿入
SUBSTITUTE	1110 0000 0000 0000	一塩基置換
DELETE	1111 0000 0000 0000	欠失
mutmsk	1111 0000 0000 0000	ビットマスク

しかしながら、実際には、INSERT が使われていなかったり、後述するように、この上位4ビットに、insertionの塩基数が書き込まれていたりするので、実際にはゴチャゴチャしている印象ですね。

ビット演算の活用

さて、これらのメタデータがどのように活用されるか見ていきましょう。
ここではビット演算というものが登場します。（僕もビット演算についてはよくわかっておらず、必要になるたびにGoogleに聞いて確認するぐらいのレベルです。）

メタデータの読み取り

読み込みは以下のように行っています。

(c & mutmsk) != NOCHANGE

(c & mutmsk) == SUBSTITUTE)

ビットマスク mutmusk と AND 演算を行いメタデータが保存されている上位4ビットのメタデータのみを取り出します。それからビットパターンと比較を行います。

メタデータの書き込み

書き込みは、OR 演算を行います。たとえば、deletionが伸展するときのコードは次のようになります。

c |= DELETE;

この |= は += と同じで、

c = c | DELETE;

を意味します。DELETEは上述のように

1111 0000 0000 0000

を意味しますので、下位ビットに関わらず、上位4ビットにフラグが立ちます。
ちなみに、deletionの伸展ではなく、最初の欠失の場合は、

c = DELETE;

となっていますので、この場合は元の塩基の情報は消えているようです。

メタデータにinsertion用の塩基配列を書き込む

さて、難しいのが、insertionにおける塩基の扱いです。挿入される塩基配列もメタデータの中に書き込まれます。

int num_ins = 0, ins = 0;
do {
	num_ins++;
	ins = (ins << 2) | (int)(drand48() * 4.0);
} while (num_ins < 4 && drand48() < INDEL_EXTEND);
	if (is_hap || drand48() < 0.333333) { // hom-ins
		ret[0]->s[i] = ret[1]->s[i] = (num_ins << 12) | (ins << 4) | c;
	} else { // het-ins
		ret[drand48()<0.5?0:1]->s[i] = (num_ins << 12) | (ins << 4) | c;
	}
}

num_ins は挿入される塩基数を表します。これは 1, 2, 3 のいずれかです。
num_ins << 12 は挿入数を上位4ビットに配置します。
ins は挿入される塩基を表します。ACTGが2ビットずつ表現されます。挿入数は最大で3であり、このとき6ビットが消費されます。
ins << 4 は挿入された塩基のシーケンスを5-10ビットに配置します。
残りの4ビットには現在の塩基 c が配置されます。

上位4ビットに num_ins が書き込まれています。そして、上で規定した INSERT はコード中で使われていません。このあたりは少し混乱があると思いますがそんな感じになっています。

挿入される塩基の数は、最大で3つです。このことを考えると、上位4ビットはあくまでフラグ専用にして、上位 5-7 に num_ins を格納し、上位8-12（下位 5-11)ビットに、挿入塩基を格納し、下位4ビットに元々の塩基情報を入れてもいいかなと思いますが、これだとカツカツそうですね。

コード中をよくみるとスペースとタブが混在しています。これを見ると、Heng Li さんも締め切りなどに追われて焦っていたのかなという気になります。

Heng Liさんは現役では世界最高のバイオインフォマティシャンの１人だと思いますが、やはりそんな彼も人間で、コードの中には人間らしい部分が残っているようだ、などと思いました（すみません）。

正規分布の乱数

C言語の組み込みの関数では正規分布の乱数を生成しないらしく下のようなコードが用意されています。（実はRubyやCrystalでも標準の関数では正規分布の乱数は用意されていないので外部のライブラリが必要になります）

/* Simple normal random number generator, copied from genran.c */

double ran_normal()
{ 
	static int iset = 0; 
	static double gset; 
	double fac, rsq, v1, v2; 
	if (iset == 0) {
		do { 
			v1 = 2.0 * drand48() - 1.0;
			v2 = 2.0 * drand48() - 1.0; 
			rsq = v1 * v1 + v2 * v2;
		} while (rsq >= 1.0 || rsq == 0.0);
		fac = sqrt(-2.0 * log(rsq) / rsq); 
		gset = v1 * fac; 
		iset = 1;
		return v2 * fac;
	} else {
		iset = 0;
		return gset;
	}
}

まとめ

ここまで見てきた範囲でも、

kseq.h でFASTA/FASTQファイルを読み込む
16ビット整数配列に変換する
上位12ビットにメタデータを配置する
- 読み込みにはビットマスクを用いる
- 書き込みには OR 演算を用いる
insertionのように長さが決まっていない場合は、
- 2ビットで１塩基を表現（a: 0, c: 1, g: 2, t: 3, n: 4)
- 適宜左シフト >> 2 して確保できたスペースに | で書き込み
- さらに length も書き込んでおく

といった技法が明らかになりました。

16に拡張したことでできたスペースにビット演算を利用してメタデータを書き込むというのが、一つの高速化技法ということになるのでしょう。これは、特別な専用の構造体や　ビットフィールドを用意せずに済ましているという点では、Rubyの配列を用いた雑なスクリプトの技法に似ているところもあります。（よくわかってない個人の感想です）

コードを眺めていて簡単に思いつく変更のアイディアは

insertionのサイズが3以下に制約される
insertionの途中からのリードは生成されにくい？？
copy number variationをはじめとした構造変異が考慮されていない
Quality値が一定である
シークエンスのシミュレーションと、変異発生という異なるプロセスが一つのツールにまとめられているので、リファクタリングの余地がある
マルチスレッドに対応していない

などですね。なお、これは13年前のツールであり、現在では多数の素晴らしいシーケンシングシミュレータが存在することに留意ください。

芸術作品としてのコード(ポエム)

逆の見方をすると、wgsimは決して全ての変異を再現しようとはしていません。生物で起きていることを完全に再現しようとはせず、１塩基変異を中心に、あくまで注目している側面についてのみを単純なモデルと小さなデータ構造で端的に表現しています。

私たちはともすれば完璧さを求めて、ツールの機能拡張を行いますが、科学という窓を通して現象を理解しようとする人間の営みという、歴史的な観点から作品を評価するならば、むしろこのような「切り取り方のセンス」の方に、大きな意味があるのかもしれません。

Heng Liさんのツールはいつもミニマムなデータ構造に対する強いこだわりを感じます。これは、バイオインフォマティクス用ソフトウェアの特徴というよりも、Liさんという一個人の個性、作品の特徴ではないかなと思います。私もいつの日か自分の作風を確立したいなーとか思うなどしました。

wgsimはアルゴリズム的な難しさがほとんどないプロジェクトでした。

アルゴリズムは私にとってはとても難しいのでおいおい勉強するとして、ビット演算をもっと直感的に理解できれば、C言語のコードを今よりも読み解けるようになるのかな、と思いました。

実は、もともとwgsimをCrystal言語に移植しようとして、wgsimを調べ始めました。(それは仕事で必要性に迫られたり、実益を目的としたものではなく、単純な興味からです)冒頭にも書いた通り、Heng Li作品の中ではコードの量も少なく、アルゴリズムもほとんどありませんから、ChatGPTを使えば簡単に移植できそうな気がして、入門編にちょうどいい気がしたのです。しかし、wgsimに対する理解が深まるとともに、この作品をどのような方針でCrystalに移植するべきなのか、よくわからなくなってしましたした。

確かにCrystalでビット演算を行うことも不可能ではないでしょう。しかし、それでは例えCrystal言語で書いたとしても、Crystal言語で書いたことにはならない気がします。かといって、冗長なデータ構造で処理を進めるのであれば、元の作品の精神を十分に理解したとは言えない気がします。

これまでプリミティブではないCrystalの型がどう表現されているかについて、あまり気にしていませんでしたが、Crystalもオブジェクト指向言語であるから、多少はメタデータなどのデータがメモリに展開されているはずで、それがどうなっているのか少し気になってしまいました。

コードを観察しながらそんな事を考えた正月でした。

この記事は以上です。

wgsimのコードを眺める