【ゲノム】BGZFという圧縮形式を詳しくみる #bioinformatics

はじめに

ゲノム関連のファイルフォーマットでは、BGZF(Blocked GNU Zip Format)という圧縮形式をしばしば目にします。
たとえば、Bamや、Bcfも、BGZF圧縮されたファイルです。BGZFファイルは、ランダムアクセスをしたい大きなファイルを圧縮するときに用いられています。

ではBGZFとは一体どんなファイル形式なのでしょうか？

圧縮ファイルにランダムアクセスしたい

ゲノム関連ファイルはサイズが大きいので、圧縮して利用することが多いです。しかし、単純にファイルを圧縮してしまうと、ファイルの任意の位置を読み込むことができません。そこで、一定のサイズごとにファイルを切り分けて、分割して圧縮します。この分割・圧縮された単位のことをブロックと呼びます。ファイルのどの部分が、圧縮ファイルのどのブロックに対応するのか対応表を作成しておきます。この対応表をインデックスと呼びます。インデックスを参照すれば、圧縮ファイルのうち、必要な範囲だけを解凍することができます。このようにして圧縮ファイルへのランダムアクセスを実現するのがBGZF形式です。

画像出典：https://samtools.github.io/bcftools/bgzf-aes-encryption.pdf

bgzipコマンドの簡単な使い方

bgzip で圧縮

BGZFはバイオインフォマティックスのファイルフォーマットに特化した専用の圧縮形式だと誤解されがちです。私も一昨日まで誤解していました。たしかにBGZFで検索すると生命情報での利用事例ばかりがヒットして、他の世界ではこの圧縮形式はほとんど利用されていないようです。しかし、実際にはBGZFは全てのファイルを圧縮することができます。

bgzip file_name # file_name はどんなファイルでもOK

圧縮されたファイルは、file_name.gz となります。

gzip で解凍

BGZFファイルは gzipコマンドで解凍することができます。

gzip -d file_name.gz

これは、BGZFがgzipの一種だからです。後述のように、gzipの圧縮ファイルに追加のフィールド（自身のブロックのサイズ）を記録して、連結させたものがBGZFになります。

インデックスを作成する

bgzipはgziという形式のインデックスを作成することができます。圧縮するときに-i オプションを指定するか、圧縮後のファイルに-r オプションをつけてbgzipを実行します。

bgzip -i file_name
bgzip -r file_name.gz

実行すると file_name.gz.gzi が作成されます。内容は、後述するように、BGZFファイル内のブロック数、ブロックごとの、圧縮したときのファイルオフセット、解凍したときのファイルのオフセットが書き込まれています。

bgzip のインストール

Ubuntuではbgzipコマンドは tabix コマンドをインストールするとbgzipコマンドも使えるようになります。

apt install tabix`

gzipよりファイルサイズは大きい

BGZFはgzipよりも、数%〜20%程度ファイルサイズが大きくなります。だからランダムアクセスを必要としないファイルはgzipで圧縮しない方がよいでしょう。

BGZFの内容を詳しく見る

BGZFはGZIPが連結されたものです。GZIPはcatで結合することができます。だからGZIPを単純に結合するだけで一つのGZIPファイルになります。これはBGZFではなくもとからのGZIPの仕様です。

BGZFはGZIPの一種ですが、フィールドを拡張して、自身の圧縮後のファイルサイズをヘッダーに含めるようにしています。下の画像の赤い部分が、BGZFのヘッダーが一般的なGZIPファイルと異なるところです。

画像出典：https://samtools.github.io/bcftools/bgzf-aes-encryption.pdf

まず、Extra fieldがあることを示すために2bitのフラグ、つまり4を立てています。そして、Extra fieldとして、BLENを格納しています。BLENは(ブロック全体のファイルサイズ -1)です。

hts-specで公開されているSAMv1仕様書にBGZFの中身について詳しく書かれています。

ありがたいことに、GZIPの仕様を日本語に翻訳しているサイトがありました。これを参照すると以下のことがわかります。

ID1, ID2 はgzip フォーマットの識別子で、31, 139に決められています。
CM は圧縮方法で、8は"deflate"を意味します。
FLG は4ですが、これはbit2を意味すると思うので、Extra fieldがあるという意味です。
MTIME はタイムスタンプですが、手元のファイルで確認するとすべて０になっている。diffでファイルの差を確認する人に優しい。
XFL はlibdeflateのアルゴリズムに関するフラグで0が入っています。
OS オペレーティングシステムには255が入っており、これは「不明」を意味します。（改行コードの特定などに便利ということで存在するらしい）
XLEN はヘッダーの中の追加のフィールドのバイト数です。たとえば06であれば、6バイトを意味する。SI1 uint8 SI2 uint8 SLEN uint16 BSIZE uint16 とのことなので、合計すると6バイトで、一般的なBGZFファイルは6が多いはず。
- SI1 追加のフィールドの識別子その１ 66に決められています。
- SI2 追加のフィールドの識別子その２ 67に決められています。
- SLEN は追加のフィールドの長さを意味するらしい。ここでは2。2バイト、BSIZEがuint16だから、という理解でいいだろうか？
- BSIZE block size - 1。ブロック自身のサイズ(byte) 圧縮する前のサイズが64Kb以下であることが求められるため、圧縮後のデータのサイズは64Kbよりも小さくります。（圧縮する前のデータのサイズは後述のISIZEに格納される）　
CDATA は圧縮されたデータ。データサイズは（(BSIZE + 1) - XLEN(6) - 20) バイト
CRC32 CRCは日本語で巡回冗長検査というらしく、誤り検出符号だそうです。
ISIZE は圧縮される前のデータのサイズ。手元のファイルだと64kbよりわずかに小さい65280が格納されています。

以上のことは、バイナリファイルを閲覧できるソフト、たとえば od コマンド、hexylコマンド、ghex を利用することによって観察できます。(ghexを起動したところ）

GZI（インデックス）の内容を詳しく見る

インデックスgzi形式の中味は、bgzipコマンドのマニュアルに簡潔に書かれています。

ファイルの中身は、

uint64_t ブロックの数

のあとに以下の２つのペアが続く

uint64_t 圧縮ファイルのオフセット
uint64_t 非圧縮ファイルのオフセット

オフセットはファイル上の位置の座標の意味。

実際のファイルの中味を見るには、次のようなスクリプトを書きます。
（Rubyが得意なのでRubyで書きましたが、どの言語でも同様に書けると思います）

x = ARGF.read.unpack('Q*') # uint64

n = x.shift

# 一応チェック
raise "#{x.size} #{n * 2}" unless x.size == n * 2

n.times do |i|
  puts "#{i} 圧縮オフセット:#{x.shift} 非圧縮オフセット:#{x.shift}"
end

実行結果 ruby a.rb GRCh38.primary_assembly.genome.fa.gz などと実行すると次のような結果が出力されました。

0 圧縮オフセット:16783 非圧縮オフセット:65280
1 圧縮オフセット:36481 非圧縮オフセット:130560
2 圧縮オフセット:55815 非圧縮オフセット:195840
3 圧縮オフセット:61084 非圧縮オフセット:261120
4 圧縮オフセット:73972 非圧縮オフセット:326400

...
48270 圧縮オフセット:894444370 非圧縮オフセット:3151130880
48271 圧縮オフセット:894463973 非圧縮オフセット:3151196160
48272 圧縮オフセット:894478136 非圧縮オフセット:3151261440
48273 圧縮オフセット:894491446 非圧縮オフセット:3151326720
48274 圧縮オフセット:894503465 非圧縮オフセット:3151392000

この圧縮オフセットが正しいか、たとえばhexylなどのコマンドラインツールを使って、確認することができます。たとえば、オフセット36481から、本当にBGZFのブロックが始まっているか見てみましょう。

hexyl -s 36481 GRCh38.primary_assembly.genome.fa.gz -n 64

最初の 1f 8b 08 04 00 00 00 00 00 ff 06 00 42 43 02 00 がまさにBGZFファイルの開始パターンになっていることがわかります。

┌────────┬─────────────────────────┬─────────────────────────┬────────┬────────┐
│00008e81│ 1f 8b 08 04 00 00 00 00 ┊ 00 ff 06 00 42 43 02 00 │•×••0000┊0×•0BC•0│
│00008e91│ 85 4b 9d 9d 5b 9a dc 38 ┊ af 6c df 39 2b 7c 78 e0 │×K××[××8┊×l×9+|x×│
│00008ea1│ 04 38 ff b1 9c 4a ac 15 ┊ 94 b2 ec ee bf f7 71 b7 │•8×××J×•┊××××××q×│
│00008eb1│ ed 72 5e 24 8a 04 02 81 ┊ 00 28 55 ed ae da d5 a7 │×r^$×••×┊0(U×××××│
└────────┴─────────────────────────┴─────────────────────────┴────────┴────────┘

faidxを利用してBGZFを検索することを考える

jbrowse2などのゲノムブラウザなどのソフトウェアでは、リファレンスゲノムfai.gzを開くときに、fai と gzi の２つのインデックスファイルを必要とします。その理由を考えてみましょう。

ここで、faidxの中身を思い出してみよう。

faidxは

chr1    248956422   8   60  61
chr2    242193529   253105712   60  61
chr3    198295559   499335808   60  61
chr4    190214555   700936301   60  61
chr5    181538259   894321107   60  61
chr6    170805979   1078885012  60  61
chr7    159345973   1252537766  60  61
chr8    145138636   1414539514  60  61
chr9    138394717   1562097136  60  61
chr10   133797422   1702798442  60  61

こんな感じのファイルフォーマットで、

列	項目名	説明
1	名前	参照配列の名前
2	長さ	参照配列の全長 (塩基)
3	オフセット	この配列の最初の塩基の FASTA/FASTQ ファイル内のオフセット
4	ラインベース	各行の塩基数
5	線幅	改行を含む各行のバイト数
6	QUALOFFSET	FASTQ ファイル内のシーケンスの最初の品質のオフセット

3行目にオフセットがあります。

つまり、特定の染色体を開きたいと思ったときに、faiファイルを参照して目的の染色体の非圧縮オフセットを確認します。次にその非圧縮オフセットに対応する、BGZF圧縮ファイルのブロックのオフセットを取り出します。最後に、BGZFを開き、オフセットから必要なブロックだけ参照して解凍します。このようにすれば、圧縮ファイルのうち、必要な部分だけ解凍することが可能になります。

以上でだいたいBGZFファイルの仕組みがわかりました。

今後調べたいこと

しかし、まだわからない部分があります。

BamファイルやBcfファイルは単にSAMやVCFをBGZF形式に圧縮したものではなく、ファイルサイズがさらに小さくなるように、文字列ではなくバイナリとして表現しています。そのため、BamやBcfを解凍すると、何となく雰囲気はわかるけど文字化けしたような感じのファイルになります。さきほどのfaidxの場合は染色体ごとにインデックスを作成していますが、baiでは、染色体をビンで区切って、ビンをスパンに分割しているようです。検索にはR木というアルゴリズムを使っているらしいのですが、よくわかりません。それぞれのビンのサイズは同じでなくても良いらしくて、ビンのサイズを変えることによって効率化しうるということなのですが、実際のBamファイルが本当にそのような戦略をとっているのか、線形に区切っているだけなのかよくわかりませんでした。これらについては調べて自分が理解できるようなら記事にしたいと思います。

この記事は以上です。