結論

公式リファレンス

https://www.htslib.org/doc/faidx.html

htslib公式サイトには次のように書かれています。

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

これはなに？

Fasta形式のリファレンスゲノムとかがあったときに、

samtools faidx hogehoge.fq

というコマンドを打つことがあると思います。すると、hogehoge.fai っていうファイルが生成されます。
この .fai っていうファイルは一体なにものなんでしょうか？

とりあえず中身を見てみようか

いろいろ考えるよりやった方が早いことが結構あります。
たとえばヒトゲノムのリファレンスにfaidxをつけてみましょう。

samtools faidx GRCh38.p13.genome.fa

GRCh38.p13.genome.fa.fai というファイルができあがります。ここで注目すべきことは、faidx サブコマンドには、 -@ などの並列化オプションが存在しないことです。faidx に限らず、インデックスを作るコマンドは並列化が効かないケースが多い印象がありますね。なんででしょうね。難しいことはわからない。

head GRCh38.p13.genome.fa.fai

さて、このインデックスファイルは中身はふつうのテキストファイルでした。

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

tail GRCh38.p13.genome.fa.fai

KI270748.1	93321	3320702127	60	61
KI270749.1	158759	3320797027	60	61
KI270750.1	148850	3320958455	60	61
KI270751.1	150742	3321109809	60	61
KI270752.1	27745	3321263087	60	61
KI270753.1	62944	3321291318	60	61
KI270754.1	40191	3321355335	60	61
KI270755.1	36723	3321396219	60	61
KI270756.1	79590	3321433578	60	61
KI270757.1	71251	3321514518	60	61

1列目は、染色体の名前、2列目は塩基配列の長さでしょうか。
なんでしょうね。この60とか、61とかって。

余談ですが、YouPlotというターミナル上にグラフ描出ソフトがございまして、

ターミナル上に棒グラフを描くことができます。

~~すきあらば自作ツールの宣伝~~

自分でインデックスを作成してみよう

こうやって見てみると思ったよりも簡単ですね。最後に本当に理解できていることを確認するために、自分でもインデックスを作るスクリプトを書いてみることにします。

私はRubyしか書けないのでRubyです。

完成イメージはこんな感じでしょうか？

cat genome.fa | ruby index.rb > genome.fa.fai # または
./index.rb genome.fa > genome.fa.fai

ちょっと富豪的プログラミングですが、以下のようなコードにしてみました。

#!/usr/bin/env ruby

name       = nil
index      = 0
bases      = []
width      = []
offset     = 0
line_width = 0

def puts_tab(*args)
  puts args.map(&:to_s).join("\t")
end

while l = ARGF.gets
  index += l.size

  if l.start_with?('>')
    puts_tab(name, bases.sum, offset, bases[0], width[0]) if name

    name = l[1..-1].split(' ')[0]
    offset = index
    bases = []
  else
    width << l.size
    bases << l.chomp.size
  end
end

puts_tab(name, bases.sum, offset, bases[0], width[0])

とりあえず GRCh38.p13.genome.fa に関しては samtools faidx 全く同じ結果が得られることを確認しました。

おわりに

samtools の faidx で求められるインデックスファイルは簡単な形式であることがわかりました。この記事は以上です。