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

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

ゴールデンウィーク中にMinipileupのコードを少しだけ眺めた。
Minipileupは、BAMファイルのバリアントの数を数えてくれるシンプルなツールである。

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の中に含まれていたが、最近新しいリポジトリとして公開された。古いツールが新たに公開されると、新しいツールの伏線になっていることが多い。例えば、miniprotとpangeneのように。しかし今回は、必ずしもそうではなく、ロングリードのRNA-SeqバリアントコーラーであるlongcallRと関係しているのかもしれない。

何を知りたかったのか？

C言語のコードを素人が読むのは手間がかかる。今回知りたかったのは、どうやってBAMファイルのリードの特定の位置の塩基を取り出しているかということだ。リードには、挿入（insertion）や欠損（deletion）が含まれているので、単純にインデックスから塩基を見ることはできない。これらの情報は、BAMファイルではCIGARという形式で記録されている。本格的なバリアントコーラーなら再アラインメントをしているかもしれないが、MinipileupはBAMから直接情報を読んでいるはずだ。（CIGARとよく似た概念で、MDタグというものがあり、これはリファレンスとの関係を表現している。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にアップした。

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

// 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でpileupに対応している。

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やhts.cr を細々とメンテナンスしているが、pileupに関係するAPIはまだ提供できていない。hts-nimでもpileupに関する関数は提供していないので、なかなか難しい部分があるようだ。Rubyでは2つの関数を引数に取るメソッドは原則作らず、2つに分割した方が良い。また、Crystalではクロージャーを定義するのが少々面倒である。

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

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

その他

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

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

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