はじめに
minimap2はロングリードシーケンス用のマッピングツールです。
minimap2では、リファレンスゲノムに対するインデックスを作成することができます。しかし、インデックスファイルを、どのfasta・オプションで作成したか、からなくなってしまう場合があるかもしれません。
今回はMinimpa2のインデックスファイルについて調べました。
Minimap2の2つのオプション
minimap2 でインデックスを作成するときに -x オプションでプリセットの数値を利用することができます。
Minimap2の内部には2つのオプション構造体があります。
mm_idxopt_t インデックスのオプション
typedef struct {
short k, w, flag, bucket_bits;
int64_t mini_batch_size;
uint64_t batch_size;
} mm_idxopt_t;
mm_mapopt_t マッピングのオプション
typedef struct {
int64_t flag; // see MM_F_* macros
int seed;
int sdust_thres; // score threshold for SDUST; 0 to disable
int max_qlen; // max query length
int bw, bw_long; // bandwidth
int max_gap, max_gap_ref; // break a chain if there are no minimizers in a max_gap window
int max_frag_len;
int max_chain_skip, max_chain_iter;
int min_cnt; // min number of minimizers on each chain
int min_chain_score; // min chaining score
float chain_gap_scale;
float chain_skip_scale;
int rmq_size_cap, rmq_inner_dist;
int rmq_rescue_size;
float rmq_rescue_ratio;
float mask_level;
int mask_len;
float pri_ratio;
int best_n; // top best_n chains are subjected to DP alignment
float alt_drop;
int a, b, q, e, q2, e2; // matching score, mismatch, gap-open and gap-ext penalties
int transition; // transition mismatch score (A:G, C:T)
int sc_ambi; // score when one or both bases are "N"
int noncan; // cost of non-canonical splicing sites
int junc_bonus, junc_pen;
int zdrop, zdrop_inv; // break alignment if alignment score drops too fast along the diagonal
int end_bonus;
int min_dp_max; // drop an alignment if the score of the max scoring segment is below this threshold
int min_ksw_len;
int anchor_ext_len, anchor_ext_shift;
float max_clip_ratio; // drop an alignment if BOTH ends are clipped above this ratio
int rank_min_len;
float rank_frac;
int pe_ori, pe_bonus;
float mid_occ_frac; // only used by mm_mapopt_update(); see below
float q_occ_frac;
int32_t min_mid_occ, max_mid_occ;
int32_t mid_occ; // ignore seeds with occurrences above this threshold
int32_t max_occ, max_max_occ, occ_dist;
int64_t mini_batch_size; // size of a batch of query bases to process in parallel
int64_t max_sw_mat;
int64_t cap_kalloc;
const char *split_prefix;
} mm_mapopt_t;
これはとても難しいので理解する必要はありません。ただ、2つのオプションがあることを知っていればいいと思います。
ではプリセットでどのようにオプションが設定されるのか？
これは option.c に記載されています。
mmiファイルの作成に関係しているのは主に、k, w, flag の3点ですね。これに加えて、染色体名の一覧を出力できれば、mmi ファイルの概要がわかるというものです。ChatGPTに表を生成してもらいました
| フィールド名 | データ型 | 説明 | デフォルト値・備考 |
|---|---|---|---|
k |
short |
k-merのサイズ（文字数）。シーケンスをk-merに分割する際の長さを指定します。 | 通常 15（長いほど精度が向上、短いと柔軟性が増す） |
w |
short |
シードウィンドウサイズ（window size）。minimizer（代表k-mer）を選ぶウィンドウの長さ。 | 通常 10（小さいと感度向上、大きいと速度向上） |
flag |
short |
オプションフラグ。インデックスやマッピング時の動作を制御するフラグビット。 | 例: MM_I_HPC（重複k-merの圧縮）など |
bucket_bits |
short |
バケットビット数。インデックス作成時のバケットの分割数（2^bucket_bits）。 |
通常 14（2^14 = 16384 バケット） |
よって、mmiファイルを参照して、これらのオプションがどうなっているかを見れば、どのプリセットを指定したか推測することができるわけです。
Minimap2 インデックスファイル仕様 (MMI)
ChatGPTに index.c を見せて、MMIファイルの仕様書を書いてもらいました。
ファイルフォーマット概要
| オフセット | サイズ | データ型 | 内容 |
|---|---|---|---|
| 0x00 | 4バイト | char[4] | マジックナンバー (MMI\2) |
| 0x04 | 4バイト | uint32_t | k-mer サイズ (w) |
| 0x08 | 4バイト | uint32_t | minimizer サイズ (k) |
| 0x0C | 4バイト | uint32_t | バケット数のビット数 (b) |
| 0x10 | 4バイト | uint32_t | シーケンス数 (n_seq) |
| 0x14 | 4バイト | uint32_t | フラグ (flag) |
シーケンス情報 (シーケンスリスト)
続く部分では、各シーケンスの名前と長さが格納されています。
| 内容 | サイズ | 説明 |
|---|---|---|
| 名前長さ | 1バイト | シーケンス名の長さ (0の場合は無名) |
| シーケンス名 | 可変 (名前長さ) | NULL 終端なし |
| シーケンス長 | 4バイト | シーケンスの長さ |
バケット情報
インデックスはバケット（B）単位で分割され、各バケットに minimizer 情報が格納されます。
| 内容 | サイズ | 説明 |
|---|---|---|
重複 minimizer 数 (n) |
4バイト | minimizer が複数回出現した位置の数 |
位置情報 (p) |
8バイト × n 個 | minimizer が複数回出現した位置情報 |
ハッシュサイズ (size) |
4バイト | minimizer のハッシュテーブルのサイズ |
ハッシュデータ (h) |
16バイト × size 個 | minimizer とその位置情報のペア |
シーケンスデータ (S)
flag に MM_I_NO_SEQ が設定されていない場合、シーケンスデータそのものもインデックスに含まれます。
| 内容 | サイズ | 説明 |
|---|---|---|
シーケンスデータ (S) |
4バイト × ((全長 + 7) / 8) 個 | 圧縮されたシーケンス情報 (2ビット/塩基) |
ハッシュテーブル (minimizer の索引)
バケットごとのハッシュテーブルには minimizer とその出現位置が格納されています。
- キー: minimizer 値 (64ビット)
-
値:
-
奇数キー: 単一位置 (
uint64_t) -
偶数キー:
p配列へのオフセットと出現回数を 64ビットで表現- 上位32ビット:
p配列のオフセット - 下位32ビット: 出現回数
- 上位32ビット:
-
奇数キー: 単一位置 (
フラグ (flag) の詳細
| ビット位置 | フラグ名 | 説明 |
|---|---|---|
| 0 | MM_I_NO_SEQ |
シーケンスデータを保存しない |
| 1 | MM_I_HPC |
Homopolymer-compressed モード |
データの読み出しフロー
-
マジックナンバーの確認:
MMI\2か確認 -
メタデータの読み込み:
w,k,b,n_seq,flag - シーケンスリストの読み込み: 各シーケンスの名前と長さ
- バケット情報の読み込み: minimizer 情報、ハッシュテーブル
- シーケンスデータの読み込み: 必要に応じて
補足
- minimizer の出現回数が 1 の場合は直接格納し、複数回出現する場合は位置配列
pを参照することで高速化しています。 - minimizer のソートには Radix Sort を使用しています。
上記はChatGPTの出力ですので、実際に正しいかどうかはソースコードを確認してください。
インデックスのメタデータを表示するコマンドラインツールを作成する
minimap2はC言語のライブラリとして利用することもできます。これを利用して、minimap2のインデックスのメタデータを表示するコマンドラインツールを作成しました（正確にはChatGPTによって書いてもらいました。）
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "minimap.h"
#include "mmpriv.h"
void print_mmi_metadata(const char *mmi_file) {
FILE *fp = fopen(mmi_file, "rb");
if (fp == NULL) {
fprintf(stderr, "Error: Cannot open .mmi file '%s'\n", mmi_file);
exit(1);
}
// Load the .mmi file
mm_idx_t *mi = mm_idx_load(fp);
fclose(fp);
if (mi == NULL) {
fprintf(stderr, "Error: Failed to load .mmi file '%s'\n", mmi_file);
exit(1);
}
// Print general index parameters
printf("MMI Metadata for '%s':\n", mmi_file);
printf(" k-mer size: %d\n", mi->k);
printf(" Seed width: %d\n", mi->w);
printf(" Number of sequences: %d\n", mi->n_seq);
// Print information for each reference sequence
printf("Reference Sequences:\n");
for (int i = 0; i < mi->n_seq; ++i) {
printf(" [%d] Name: %s, Length: %d\n", i + 1, mi->seq[i].name, mi->seq[i].len);
}
// Free the index structure
mm_idx_destroy(mi);
}
int main(int argc, char *argv[]) {
if (argc != 2) {
fprintf(stderr, "Usage: %s <index.mmi>\n", argv[0]);
return 1;
}
const char *mmi_file = argv[1];
print_mmi_metadata(mmi_file);
return 0;
}
さて、このコマンドは、Minimpa2をライブラリとして使用しているため信頼性は高いと思いますが、全てのインデックスを読み込むので少し遅いです。実際にはヘッダーの部分（がどこかはminimap2では明示されていませんが）だけ読み込めば、大まかななことがわかるので、minimap2をライブラリとして使わず、Crystalで簡単なツールを作成しました。
if ARGV.size != 1
puts "Usage: crystal #{__FILE__} <index.mmi>"
exit 1
end
File.open(ARGV[0], "rb") do |file|
# Read and verify the magic number
magic = file.read_string(4)
puts "Magic number (raw): #{magic.inspect}"
if magic != "MMI\2"
raise "Invalid .mmi file: Magic number mismatch"
end
# Read the header values
seed_width = file.read_bytes(Int32)
kmer_size = file.read_bytes(Int32)
bucket_bits = file.read_bytes(Int32)
num_sequences = file.read_bytes(Int32)
flags = file.read_bytes(Int32)
puts "MMI Header:"
puts " k-mer size: #{kmer_size}"
puts " Seed width: #{seed_width}"
puts " Bucket bits: #{bucket_bits}"
puts " Number of sequences: #{num_sequences}"
puts " Flags: 0x#{flags.to_s(16)}"
puts
# Read sequence information
puts "Sequences:"
num_digits = num_sequences.to_s.size
puts " Idx\tName\tLength"
num_sequences.times do |i|
name_length = file.read_bytes(UInt8)
name = file.read_string(name_length)
seq_length = file.read_bytes(Int32)
puts " [#{(i + 1).to_s.rjust(num_digits)}]\t#{name}\t#{seq_length}"
end
end
これはかなり高速です。私はRubyやCrystalが好きなのですが、もちろんPythonであっても全く同様のコードを簡単に作成することができるでしょう。
終わりに
仕様書のないバイナリファイルは怖いですが、AIの進展により、単にC言語のコードを読ませることで、バイナリファイルの仕様書を書いてもらえることがわかりました。AIは間違えることもあるので、怪しいと思ったらソースコードに直接当たることも必要ですが、それでも以前に比べると、少ないスキルでも、できることのリーチが大きく伸びたと感じます。
この記事は以上です。