Minipileupのコードを少しだけ眺めてみた　BAMファイルからバリアントを数えるシンプルなツール

## Minipileupのコードを少しだけ眺める

ゴールデンウィーク中に[Minipileup](https://github.com/lh3/minipileup)のコードを少しだけ眺めた。  
Minipileupは、BAMファイルのバリアントの数を数えてくれるシンプルなツールである。

```bash
minipileup -f ref.fa in.bam
```

出力例：

```
chr1	31	C	C,G	0/1:15,1
chr1	37	G	G,T	0/1:21,1
chr1	38	G	G,T	0/1:21,1
chr1	39	G	G,T	0/1:21,1
chr1	41	T	T,A	0/1:24,1
chr1	42	C	C,G	0/1:24,1
chr1	46	C	C,G	0/1:25,1
chr1	49	G	G,T	0/1:25,2
chr1	50	G	G,T	0/1:26,1
chr1	51	T	T,A	0/1:26,1
```

このツールは、以前は[htsbox](https://github.com/lh3/htsbox)の中に含まれていたが、最近新しいリポジトリとして公開された。古いツールが新たに公開されると、新しいツールの伏線になっていることが多い。例えば、[miniprot](https://github.com/lh3/miniprot)と[pangene](https://github.com/lh3/pangene)のように。しかし今回は、必ずしもそうではなく、ロングリードのRNA-Seqバリアントコーラーである[longcallR](https://github.com/huangnengCSU/longcallR)と関係しているのかもしれない。

### 何を知りたかったのか？

C言語のコードを素人が読むのは手間がかかる。今回知りたかったのは、どうやってBAMファイルのリードの特定の位置の塩基を取り出しているかということだ。リードには、挿入（insertion）や欠損（deletion）が含まれているので、単純にインデックスから塩基を見ることはできない。これらの情報は、BAMファイルではCIGARという形式で記録されている。本格的なバリアントコーラーなら再アラインメントをしているかもしれないが、MinipileupはBAMから直接情報を読んでいるはずだ。（CIGARとよく似た概念で、[MDタグ](https://qiita.com/kojix2/items/4342e8927047e68a4373)というものがあり、これはリファレンスとの関係を表現している。Minipileupは引数でリファレンスファイルを指定するため、このタグは考慮しなくても良さそうだ）

似たようなツールで普及しているものに`bcftools mpileup`がある。こちらもオリジナルの作者はlh3さんのはずだが、現在は別の方がメンテナンスしているようだ。

### ツールの構成

まずは、構造体や関数に注目して、ツールの構成を読んでいく。

構造体には次のようなものがある。

- `aux_t`
- `allele_t`
- `paux_t`

名前が少し分かりにくいが、`aux`という名前を見るとBAMファイルのレコードについている`aux`を思い出してしまうが関係ない。auxiliaryは「補助」を意味する単語だそうだ。`aux_t` にはリードのフィルタリング用のオプションが詰められているし、`paux_t` は pileup の auxiliary という意味だと思うが、カウントデータなどデータ一式をまとめて保存するための構造体のようである。

関数には次のようなものがある。

- `read_bam`
  - BAMの読み込みというよりは、リードのフィルタリングを担当する。マッピングクオリティ、リードの長さ、BEDファイルとのオーバーラップなどで、pileupするリードをフィルタリングする。この関数は`bam_mplp_init`関数の引数として使われる。
- `pileup2allele`
  - `bam_pileup1_t`の構造体を、より扱いやすい構造体に変換するための関数。ここでも塩基のクオリティなどによるフィルタリングが行われる。
- `print_allele`
- `main`

まず、メイン関数が呼ばれる。メイン関数内で、コマンドラインオプションの処理が行われる。コマンドラインオプションは、リードのフィルタリング、SNPのフィルタリング、出力（VCF形式かどうか）に関するものなどがある。

### BAMのpileup

結論からいうと、Minipileupは自分でBAMのパイルアップをしているわけではなく、htslibに依存している。

関数：

- `bam_mplp_init`
- `bam_mplp_auto`

および構造体：

- `bam_pileup1_t`

がポイントである。しかし、htslibのpileup関連の関数は`sam.h`を読んでもドキュメントがなく、おそらくはAPIの改良を計画しているうちに現在に至ってしまったのではないかと推測する(根拠なし)。そのため、完全な使い方がわかりにくいが、ChatGPTの助けを得ながら大体の感じを書いて[Gistにアップ](https://gist.github.com/kojix2/75f16c8c5e4503e9ba9311dfe7cb9983)した。

`bam_mplp_init`は関数を引数に取る関数である。ここでは`read_bam`を引数に取っている。それから`bam_mplp_auto`関数で、`tid`や`pos`が渡される。ここで、`bam_pileup1_t`には`qpos`というメンバがあり、これがポイントである。`qpos`はクエリーポジションを意味しているようだ。以下のような感じである。

```c
// bam_pileup1_t *p
uint8_t *seq = bam_get_seq(p->b);
putchar(seq_nt16_str[bam_seqi(seq, p->qpos)]);
```

こうして、BAMの特定のカラムの塩基を参照できるようだ。挿入や欠損の場合は`p->indel`などが使えるらしい。

### Pysam

原理はわかったが、これをC言語で書くのは大変である。Heng Liさんはライブラリのオリジナルの作者なのでコードが頭に入っているかもしれないが、一般的な人がC言語で`bam_mplp_auto`を使うのは難しいはずだ。

こういうときはPythonが便利だ。Pythonは[Pysam](https://github.com/pysam-developers/pysam)で[pileup](https://pysam.readthedocs.io/en/latest/api.html#sam-bam-cram-files)に対応している。

```python
import pysam
samfile = pysam.AlignmentFile("ex1.bam", "rb")
for pileupcolumn in samfile.pileup("chr1", 100, 120):
    print("\nCoverage at base %s = %s" % (pileupcolumn.pos, pileupcolumn.n))
    for pileupread in pileupcolumn.pileups:
        if not pileupread.is_del and not pileupread.is_refskip:
            # query position is None if is_del or is_refskip is set.
            print('\tBase in read %s = %s' %
                  (pileupread.alignment.query_name,
                   pileupread.alignment.query_sequence[pileupread.query_position]))

samfile.close()
```

私はPythonは苦手なので[ruby-htslib](https://github.com/kojix2/ruby-htslib)や[hts.cr](https://github.com/kojix2/hts.cr) を細々とメンテナンスしているが、pileupに関係するAPIはまだ提供できていない。[hts-nim](https://github.com/brentp/hts-nim)でもpileupに関する関数は提供していないので、なかなか難しい部分があるようだ。Rubyでは2つの関数を引数に取るメソッドは原則作らず、2つに分割した方が良い。また、Crystalではクロージャーを定義するのが少々面倒である。

### GATKなどの他のツールはどうだろう

結論めいたものを書ければ良いのだが、現時点では勉強不足で明確なことは言えない。確かにhtsjdkなどがあるけれども、htslibに依存していないような雰囲気を感じる。独自のpileupを実装しているのだろうか？どうなんだろうか。冒頭で挙げたlongcallRもhtslibのpileupに依存していないように見える。ロングリードの時代になって、htslibのpileupに頼るべきか、頼らないべきかが新しいツールの一つのポイントになっているのかもしれない。知らない。

### その他

- 全体的に`del`をSNPとして扱うことについて例外的に扱っているようだ。この辺りがあまり腑に落ちていない。
- C言語なので`realloc`などがところどころに見られる。すごいけど大変だなと思う。
- `count_allele`の実際に塩基を数えているところはうまく効率化されているようだが、いまいちよくわからない。

正直わからないところだらけです。それでもコードを眺めているとわかることも。そんな中途半端な知識をメモしました。この記事は以上です。