結論
公式リファレンス
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 で求められるインデックスファイルは簡単な形式であることがわかりました。この記事は以上です。