bwaのソースコードをちょっとだけ改造してindexコマンドでもverboseオプションを利用できるようにしてみよう

## はじめに

　この記事は、[bwa](https://github.com/lh3/bwa)というバイオインフォマティクスの領域で利用されているアライナーのソースコードのコマンドラインオプションの部分を少しだけいじって、`--verbose` オプションを `bwa index` コマンドでも使えるようにしてみた記録です。

　bwaは広く使われているリードのアライナーです。bwaを利用する時には、まずリファレンスのインデックスを作る必要があります。

```sh
bwa index poo.fa
```

とすると、たとえば、

```
[bwa_index] Pack FASTA... 0.00 sec
[bwa_index] Construct BWT for the packed sequence...
[bwa_index] 0.00 seconds elapse.
[bwa_index] Update BWT... 0.00 sec
[bwa_index] Pack forward-only FASTA... 0.00 sec
[bwa_index] Construct SA from BWT and Occ... 0.00 sec
[main] Version: 0.7.17-r1188
[main] CMD: bwa index poo.fa
[main] Real time: 0.093 sec; CPU: 0.004 sec
```

こんな感じの表示が得られると思います。（実行時間が短いのは、poo.fa が非常に小さなファイルだからですね。）

```:poo.fa
>chr1 1
AAAAAAAAAA
```

~~Aしかないのかよ。~~

さて，bwa index をしたときのlogですが、[bwaのソースコード](https://github.com/lh3/bwa/blob/3ddd7b87d41f89a404d95f083fb37c369f70d783/bwtindex.c#L267-L274)を参照しますと、bwaindex.cというファイルのなかの `bwa_idx_build` という関数で出力されるようになっています。その一部を引用すると、


```c:bwtindex.c
	{ // nucleotide indexing
		gzFile fp = xzopen(fa, "r");
		t = clock();
		if (bwa_verbose >= 3) fprintf(stderr, "[bwa_index] Pack FASTA... ");
		l_pac = bns_fasta2bntseq(fp, prefix, 0);
		if (bwa_verbose >= 3) fprintf(stderr, "%.2f sec\n", (float)(clock() - t) / CLOCKS_PER_SEC);
		err_gzclose(fp);
	}
```

こんな感じで `bwa_verbose` という変数が３以上の場合だけ出力する仕組みになっていることがわかります。

ところが実際の `index` サブコマンドには `--verbose` 引数は存在しません。だから bwa_verbose はずっと3のままで、この変数をいじることができません。

```
Usage:   bwa index [options] <in.fasta>

Options: -a STR    BWT construction algorithm: bwtsw, is or rb2 [auto]
         -p STR    prefix of the index [same as fasta name]
         -b INT    block size for the bwtsw algorithm (effective with -a bwtsw) [10000000]
         -6        index files named as <in.fasta>.64.* instead of <in.fasta>.* 

Warning: `-a bwtsw' does not work for short genomes, while `-a is' and
         `-a div' do not work not for long genomes.
```

ひょっとしたらヘルプに存在しないだけで、実際には `-v` を指定できる可能性もあるのでは？と思って

```
bwa index -v 0 poo.fa
```

などとしても `index: invalid option -- 'v'` として怒られてしまいます。だから、今回はせっかく用意されている `bwa_verbose` 変数を利用するべく、`-v` コマンドラインオプションを `index` に追加して、エラーログの出力も制御できるようにしてみます。

もちろん、実際にはエラーログをターミナルに表示したくない場合は、リダイレクトを利用して、

```
bwa index -v 0 poo.fa 2>/dev/null
```

としてしまえばよいです。だから、実用においては `--verbose` オプションがなくても大丈夫なのです。しかし、せっかく `bwa_verbose` で処理が切り替えられる形でソースコードが準備されているので、今回はC言語の勉強のためにbwaの本体を少しいじってみようと思いました。

## ソースコードをダウンロードしよう

ソースコードはGithubで公開されているので、gitを使ってcloneします。

```
git clone https://github.com/lh3/bwa
```

変更を本家に還元する場合、つまりプルリクエストを考慮している場合は、フォークしてからcloneする方がよいと思いますが、`bwa` には index のほかにも `--verbose` オプションが使えないコマンドがたくさんあります。bwa_verbose が指定できない仕様は意図的なものだと考えられます（後述）。なので、今回はプルリクエストを出す予定はありませんからフォークはしてもしなくても良いと思います。

### 最新の安定バージョンでブランチを切る（オプション）

今回は遊びや検証の要素が強いため、過去の安定バージョンを使うメリットはないので、マスターにそのまま変更を加えます。しかし、実際に何らかのプロジェクトで使うような場合は、安定バージョンを使いたいと思うので、

```
git checkout -b kojix2-v0.7.17 refs/tags/v0.7.17
```

としてブランチを切って、そこから作業をスタートするのが良いかと思います。

## どこでコマンドラインオプションをさばいているのか探す grep コマンド

今回は練習のために `vim` を使っています。なので、ソースコードの検索も `grep` コマンドを使って探しました。`bwa index` コマンドをうったときに表示されるヘルプの一部の `BWT construction algorithm` というフレーズを探すと。

```sh
grep -rn 'BWT construction algorithm' .
```

```sh
./bwtindex.c:236:		fprintf(stderr, "Options: -a STR    BWT construction algorithm: bwtsw, is or rb2 [auto]\n");
```

bwtindex.c の 236行目あたりに該当箇所があります。ここをvimで開くと、

```
vim bwtindex.c +236
```

すると、すぐ上の211行目〜219行目のあたりで、コマンドラインオプションを処理する記述があります。

```C:bwtindex.c
	int c, algo_type = BWTALGO_AUTO, is_64 = 0, block_size = 10000000;
	char *prefix = 0, *str;
	while ((c = getopt(argc, argv, "6a:p:b:")) >= 0) {
		switch (c) {
		case 'a': // if -a is not set, algo_type will be determined later
			if (strcmp(optarg, "rb2") == 0) algo_type = BWTALGO_RB2;
			else if (strcmp(optarg, "bwtsw") == 0) algo_type = BWTALGO_BWTSW;
			else if (strcmp(optarg, "is") == 0) algo_type = BWTALGO_IS;
			else err_fatal(__func__, "unknown algorithm: '%s'.", optarg);
			break;
		case 'p': prefix = strdup(optarg); break;
		case '6': is_64 = 1; break;
		case 'b':
			block_size = strtol(optarg, &str, 10);
			if (*str == 'G' || *str == 'g') block_size *= 1024 * 1024 * 1024;
			else if (*str == 'M' || *str == 'm') block_size *= 1024 * 1024;
			else if (*str == 'K' || *str == 'k') block_size *= 1024;
			break;
		default: return 1;
		}
```

`getopt`関数を使って、オプションを処理しています。3行目 `"6a:p:b:"` の `:` は、getoptではその直前のオプションが引数を持つことを意味しているそうです

switch-case文の最後に、`fastmap.c:181` を参考に下を書き加えます。

```C
		case 'v': bwa_verbose = atoi(optarg); break;
```

さらに `getopt` の `6a:p:b:` も書き換えて、`v` をコマンドラインオプションに取れるようにします。

```C
while ((c = getopt(argc, argv, "6a:p:b:v:")) >= 0) {
```

一応ヘルプの説明文にもverboseを追加しておきましょう。

```C
		fprintf(stderr, "         -v INT    verbosity level. [%d]\n", bwa_verbose);
```

変更点をまとめるとこんな感じになります。今回は git のディレクトリで作業しているので、`git diff` とすれば差分が把握できます。

`git diff`

```diff
diff --git a/bwtindex.c b/bwtindex.c
index 6a27ae1..12839de 100644
--- a/bwtindex.c
+++ b/bwtindex.c
@@ -210,7 +210,7 @@ int bwa_index(int argc, char *argv[]) // the "index" command
 {
 	int c, algo_type = BWTALGO_AUTO, is_64 = 0, block_size = 10000000;
 	char *prefix = 0, *str;
-	while ((c = getopt(argc, argv, "6a:p:b:")) >= 0) {
+	while ((c = getopt(argc, argv, "6a:p:b:v:")) >= 0) {
 		switch (c) {
 		case 'a': // if -a is not set, algo_type will be determined later
 			if (strcmp(optarg, "rb2") == 0) algo_type = BWTALGO_RB2;
@@ -226,6 +226,7 @@ int bwa_index(int argc, char *argv[]) // the "index" command
 			else if (*str == 'M' || *str == 'm') block_size *= 1024 * 1024;
 			else if (*str == 'K' || *str == 'k') block_size *= 1024;
 			break;
+		case 'v': bwa_verbose = atoi(optarg); break;
 		default: return 1;
 		}
 	}
@@ -237,6 +238,7 @@ int bwa_index(int argc, char *argv[]) // the "index" command
 		fprintf(stderr, "         -p STR    prefix of the index [same as fasta name]\n");
 		fprintf(stderr, "         -b INT    block size for the bwtsw algorithm (effective with -a bwtsw) [%d]\n", block_size);
 		fprintf(stderr, "         -6        index files named as <in.fasta>.64.* instead of <in.fasta>.* \n");
+		fprintf(stderr, "         -v INT    verbosity level. [%d]\n", bwa_verbose);
 		fprintf(stderr, "\n");
 		fprintf(stderr,	"Warning: `-a bwtsw' does not work for short genomes, while `-a is' and\n");
 		fprintf(stderr, "         `-a div' do not work not for long genomes.\n\n");
```


## ビルドする

これが、Heng Li氏のバイオインフォマティックスのツールが愛される一つの理由だと思うのですが、ビルドはものすごく簡単です。

```
make
```

楽ですね。


## 実行してみる

それでは変更したソフトを実行してみましょう。
まずはヘルプの表記を確認します。ディレクトリの中にある実行ファイル `bwa` を実行します。ディレクトリ内のファイルを実行する場合は、ドットとスラッシュをつけて、`./bwa` とする必要があります。

```
./bwa index
```

下記のようにヘルプに、`-v` オプションが追加されていれば成功です。

```

Usage:   bwa index [options] <in.fasta>

Options: -a STR    BWT construction algorithm: bwtsw, is or rb2 [auto]
         -p STR    prefix of the index [same as fasta name]
         -b INT    block size for the bwtsw algorithm (effective with -a bwtsw) [10000000]
         -6        index files named as <in.fasta>.64.* instead of <in.fasta>.* 
         -v INT    verbosity level. [3]

Warning: `-a bwtsw' does not work for short genomes, while `-a is' and
         `-a div' do not work not for long genomes.

```

それでは実際にファイルにインデックスをつけてみましょう。ここでは、`poo.fa` があるディレクトリで実行します。（ディレクトリの場所は各自置き換えてください）

```
cd ~/kojix2/Bio/poo
```

`bwa` コマンドを実行しますが、実行ファイルのパスを直接打ち込みます。もしも`bwa`と打ってしまうと、おそらく、すでにシステムにインストールされている `bwa` が起動されてしまいます。心配なひとは `which bwa` で `bwa` コマンドがすでにインストールされているか確認できます。

```C
~/C/bwa/bwa index -v 0 poo.fa
```

うまくいきました。[bwa_index] のログが全く表示されなくなりました。

```
[main] Version: 0.7.17-r1198-dirty
[main] CMD: /home/kojix2/C/bwa/bwa index -v0 poo.fa
[main] Real time: 0.105 sec; CPU: 0.004 sec
```

## これは、なんの役に立つの？

　ここでは、bwaをダウンロードして、少しだけソースコードを改変して、期待する挙動を作り出すことができました。ソースコードを自分で変更してソフトウェアをビルドできる環境を整えておくと、`printf` などを用いて、ソフトウェアの挙動を調査するのが可能になります。もしもなにか小さなバグを発見して、それを本家に変更したソースコードをプルリクエストという形で投げれば、自分のコードを世界中の人に使ってもらい、多くの人の作業を助けられる可能性もあります。

　**そしてなによりも、なんの役にも立たなくても、おもしろいじゃないですか!!**

## 余談：Heng Li氏のソフトウェア開発方針

余談ですが、Heng Li氏は、ソフトウェアのメンテナンスに関する考え方について自分のブログの記事を書いています。
そこでは、

 * ユーザーインターフェイスをなるべく簡素化する方針や、
 * プロジェクトを互いに独立させて外部ライブラリから独立させる方針、
 * 既存のコードベースに変更を加える場合は、コードを複製して開始する方針

が記載されていています。これはどのような状況でもベストの方法というわけではないと思いますが、彼がどのように考えてソフトウェアを作成しているかよくわかるため、とても参考になります。

https://lh3.github.io/2019/03/11/on-maintaining-bioinformatics-software

　つまり、実際にはverboseオプションを利用できるようにコードが整備されているのに、verboseオプションが利用できなかったのは、不必要にユーザーインターフェイスを複雑にせず簡単にする方針からきているものだと思います。また、プロジェクトが外部ライブラリから独立しているために `make` コマンド一つで容易にビルドできることを経験したと思います。

　このように、ソフトウェアのあり方は開発者の気持ちや姿勢が強く反映される面があります。

## まとめ

　今回は、オープンソースのバイオインフォマティックスのツールのbwaを少しだけ改造して、indexサブコマンドでも `bwa_verbose` の機能を使えるようにしてみました。R, Python, Perl, Rubyなどのスクリプト言語を一つ習得しており、なおかつ少しググるスキルがあれば、このぐらいの小さな改造を誰でも入れることができます。私自身、C言語のスキルはほとんどゼロに近いため、C言語のオプション引数の取得や、switch-case文の書き方など、この記事を通じて勉強することがいろいろありました。

　ともすればツールを使うだけになってしまいがち、なのですが、少しでもバイオインフォのソフトウェアで使われているコードを眺め、プロジェクトを自分で造らないまでも、コントリビューションできる人が増えてくれば、日本語の情報もだんだん充実してきて便利になってくるのではないかと思います。

　この記事は以上です。（記事に間違いを見つけた場合は、コメント欄でおしえてください！）