minimap2のリポジトリをながめる - その2

## はじめに

* この記事は素人の書いている個人用のメモです。minimap2の使い方を知りたい人や、minimap2のアルゴリズムの記事ではありませんのであらかじめご了承ください。また、内容の正しさについては70%ぐらいだと思っていてください。

* 本来であれば公開するクオリティの記事ではありませんが、日本語の情報は皆無なので公開しています。間違っている情報が公開されていると困る方は、ご自分でより正確な情報の記事を作成して頂き、それをインターネットに公開するという形で、私などが記事を書かなくても済みますように、建設的な解決策を取ってくださると大変助かります。

## 前回までのあらすじ

* [minimap2のリポジトリをながめる - その1](https://qiita.com/kojix2/items/e17bcc3f9595098c0ac8)

前回はminimap2のリポジトリの構成をざっくり見てきました。今回はminimap2の構造体をざっくり眺めていきたいと思います。

## minimap.h

* https://github.com/lh3/minimap2/blob/master/minimap.h

README.md 終末を読むと、`minimap.h` にC APIが書いてあるよ。`mmpriv.h` には内部向けのC APIが記載されているよ、と書かれていますので、これを見ることにします。

まず最初にダーっと定数が定義されています。

```c
#define MM_F_NO_DIAG       0x001 // no exact diagonal hit
#define MM_F_NO_DUAL       0x002 // skip pairs where query name is lexicographically larger than target name
#define MM_F_CIGAR         0x004
#define MM_F_OUT_SAM       0x008
#define MM_F_NO_QUAL       0x010
#define MM_F_OUT_CG        0x020
#define MM_F_OUT_CS        0x040
#define MM_F_SPLICE        0x080 // splice mode
#define MM_F_SPLICE_FOR    0x100 // match GT-AG
#define MM_F_SPLICE_REV    0x200 // match CT-AC, the reverse complement of GT-AG
#define MM_F_NO_LJOIN      0x400
#define MM_F_OUT_CS_LONG   0x800
#define MM_F_SR            0x1000
#define MM_F_FRAG_MODE     0x2000
#define MM_F_NO_PRINT_2ND  0x4000
#define MM_F_2_IO_THREADS  0x8000
#define MM_F_LONG_CIGAR    0x10000
#define MM_F_INDEPEND_SEG  0x20000
#define MM_F_SPLICE_FLANK  0x40000
#define MM_F_SOFTCLIP      0x80000
#define MM_F_FOR_ONLY      0x100000
#define MM_F_REV_ONLY      0x200000
#define MM_F_HEAP_SORT     0x400000
#define MM_F_ALL_CHAINS    0x800000
#define MM_F_OUT_MD        0x1000000
#define MM_F_COPY_COMMENT  0x2000000
#define MM_F_EQX           0x4000000 // use =/X instead of M
#define MM_F_PAF_NO_HIT    0x8000000 // output unmapped reads to PAF
#define MM_F_NO_END_FLT    0x10000000
#define MM_F_HARD_MLEVEL   0x20000000
#define MM_F_SAM_HIT_ONLY  0x40000000

#define MM_I_HPC          0x1
#define MM_I_NO_SEQ       0x2
#define MM_I_NO_NAME      0x4

#define MM_IDX_MAGIC   "MMI\2"

#define MM_MAX_SEG       255
```

はっきりとはわかりませんが、`MM` は minimap を、`F` は flagを意味しているのではないかな？と想像しました。実際に `example.c` には `MM_F_CIGAR` をFlagとして利用しているようなところがあるので、そう推測しました。

```c
// emulate 128-bit integers and arrays
typedef struct { uint64_t x, y; } mm128_t;
typedef struct { size_t n, m; mm128_t *a; } mm128_v
```

これは面白いですね。`uint128` というものは見たことがないですが、そういったものを実現しようとすると、こういう仕組みが必要になるんですね。普段RubyやPythonといった言語を触れていると、どうしても数値は数値だろ？という気持ちになってしまうのですが巨大な数を扱うためにはそれなりの工夫が必要になるということがわかります。

```c
// minimap2 index
typedef struct {
	char *name;      // name of the db sequence
	uint64_t offset; // offset in mm_idx_t::S
	uint32_t len;    // length
	uint32_t is_alt;
} mm_idx_seq_t;

typedef struct {
	int32_t b, w, k, flag;
	uint32_t n_seq;            // number of reference sequences
	int32_t index;
	int32_t n_alt;
	mm_idx_seq_t *seq;         // sequence name, length and offset
	uint32_t *S;               // 4-bit packed sequence
	struct mm_idx_bucket_s *B; // index (hidden)
	struct mm_idx_intv_s *I;   // intervals (hidden)
	void *km, *h;
} mm_idx_t;
```

この2つの構造体 `mm_idx_seq_t` および `mm_idx_t` はインデックスに関する構造体のようですね。C言語ですので、配列となるポインターと、その長さが別々に保存されているようです。全体的に略語が多いし、私の力でよくわからないところが多いですね。

```c
// minimap2 alignment
typedef struct {
	uint32_t capacity;                  // the capacity of cigar[]
	int32_t dp_score, dp_max, dp_max2;  // DP score; score of the max-scoring segment; score of the best alternate mappings
	uint32_t n_ambi:30, trans_strand:2; // number of ambiguous bases; transcript strand: 0 for unknown, 1 for +, 2 for -
	uint32_t n_cigar;                   // number of cigar operations in cigar[]
	uint32_t cigar[];
} mm_extra_t;

typedef struct {
	int32_t id;             // ID for internal uses (see also parent below)
	int32_t cnt;            // number of minimizers; if on the reverse strand
	int32_t rid;            // reference index; if this is an alignment from inversion rescue
	int32_t score;          // DP alignment score
	int32_t qs, qe, rs, re; // query start and end; reference start and end
	int32_t parent, subsc;  // parent==id if primary; best alternate mapping score
	int32_t as;             // offset in the a[] array (for internal uses only)
	int32_t mlen, blen;     // seeded exact match length; seeded alignment block length
	int32_t n_sub;          // number of suboptimal mappings
	int32_t score0;         // initial chaining score (before chain merging/spliting)
	uint32_t mapq:8, split:2, rev:1, inv:1, sam_pri:1, proper_frag:1, pe_thru:1, seg_split:1, seg_id:8, split_inv:1, is_alt:1, dummy:6;
	uint32_t hash;
	float div;
	mm_extra_t *p;
} mm_reg1_t;
```

　この `mm_extra_t` と `mm_reg1_t` は、アラインメントをあらわす構造体になります。アラインメントを表す構造体なので、`qs` は query start, `qe` は query end ということになりますが、いかんせん略語が多くてなかなかわかりにくいというのがあります。右側のコメントをみると、ちゃんと説明が書いてあります。一般的にRuby言語では、短すぎる命名は悪だとされ、もうすこし説明的な変数名が好まえるような気がします。その点で、C言語の流儀なのか、Heng Li氏の流儀なのかはわかりませんが、かなり略語を使っています。また、ここでもビットフィールドが大量に出てきます。`uint8` とかでそれぞれのフィールド作っておけばいいのに、とはならないんですね。ちゃんとメモリを細かく使うのが好きなようです…。

　こちらも、わかるメンバも結構あるのですが、何をやっているのかよくわからないメンバもあります。PythonバインディングのMappyでは、この`Reg1` をそのまま用いずに、`mm_hitpy_t` に移し替えて使っています。

```c:mappy.h
typedef struct {
	const char *ctg;
	int32_t ctg_start, ctg_end;
	int32_t qry_start, qry_end;
	int32_t blen, mlen, NM, ctg_len;
	uint8_t mapq, is_primary;
	int8_t strand, trans_strand;
	int32_t seg_id;
	int32_t n_cigar32;
	uint32_t *cigar32;
} mm_hitpy_t;

static inline void mm_reg2hitpy(const mm_idx_t *mi, mm_reg1_t *r, mm_hitpy_t *h)
{
	h->ctg = mi->seq[r->rid].name;
	h->ctg_len = mi->seq[r->rid].len;
	h->ctg_start = r->rs, h->ctg_end = r->re;
	h->qry_start = r->qs, h->qry_end = r->qe;
	h->strand = r->rev? -1 : 1;
	h->mapq = r->mapq;
	h->mlen = r->mlen;
	h->blen = r->blen;
	h->NM = r->blen - r->mlen + r->p->n_ambi;
	h->trans_strand = r->p->trans_strand == 1? 1 : r->p->trans_strand == 2? -1 : 0;
	h->is_primary = (r->id == r->parent);
	h->seg_id = r->seg_id;
	h->n_cigar32 = r->p->n_cigar;
	h->cigar32 = r->p->cigar;
}
```

こちらの `mm_hitpy_t` のメンバはわかりやすくて、その意味合は明瞭です。ほとんどの人はこちらを使うだけで十分そうです。

つぎに、オプションを表している構造体があります。オプションは2種類あって、`IdxOpt` と ` MapOpt` があります。Heng Li氏がREADME.mdで強調しているようにminimap2では、インデックスの作成にオプションが必要で、しかもバリエーションがあるようです。この点が、bwaのインデックスなどとは違うようです。

```c
// indexing and mapping options
typedef struct {
	short k, w, flag, bucket_bits;
	int64_t mini_batch_size;
	uint64_t batch_size;
} mm_idxopt_t;
```

インデックスのオプションです。何を意味しているのかはよくわかりません。
次にマッピングのオプションです。

```c
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;          // 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 mask_level;
	int mask_len;
	float pri_ratio;
	int best_n;      // top best_n chains are subjected to DP alignment

	int max_join_long, max_join_short;
	int min_join_flank_sc;
	float min_join_flank_ratio;

	float alt_drop;

	int a, b, q, e, q2, e2; // matching score, mismatch, gap-open and gap-ext penalties
	int sc_ambi; // score when one or both bases are "N"
	int noncan;      // cost of non-canonical splicing sites
	int junc_bonus;
	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 pe_ori, pe_bonus;

	float mid_occ_frac;  // only used by mm_mapopt_update(); see below
	int32_t min_mid_occ;
	int32_t mid_occ;     // ignore seeds with occurrences above this threshold
	int32_t max_occ;
	int64_t mini_batch_size; // size of a batch of query bases to process in parallel
	int64_t max_sw_mat;

	const char *split_prefix;
} mm_mapopt_t;
```

非常に多くのパラメーターがあることがわかります。しかし、ユーザーが実際にこれらのパラメータを扱う必要はなくて、いくつかのプリセットを指定すると、minimap側で、自動的にパラメータを指定してくれるようになっています。`options.c` を見るとそれがわかります。実際のところ、これらのオプションを適切に指定することはほとんどできないでしょう。Heng Li氏のプリセットを信用して利用する人が大半だと思います。また、非常に面倒くさいことに、presetは文字列で指定するくせに、0を投げ入れるとdefaultが指定されることになっています。

```c
// index reader
typedef struct {
	int is_idx, n_parts;
	int64_t idx_size;
	mm_idxopt_t opt;
	FILE *fp_out;
	union {
		struct mm_bseq_file_s *seq;
		FILE *idx;
	} fp;
} mm_idx_reader_t;
```

これはIndexReaderというやつで、メンバにファイルポインタを持ちます。

```c
// memory buffer for thread-local storage during mapping
typedef struct mm_tbuf_s mm_tbuf_t;
```

これが何かよくわからないのですが、スレッドローカルなメモリを確保するということだそうです。マルチコアでマッピングすることと何か関係があるのではないかと思います。

構造体はここで終わりで、以下関数が続きます。

関数は大量にあるのですが、なんとなく意味がわかるものがほとんどです、最も大切な関数はこれでしょう。

```c
/**
 * Align a query sequence against an index
 *
 * This function possibly finds multiple alignments of the query sequence.
 * The returned array and the mm_reg1_t::p field of each element are allocated
 * with malloc().
 *
 * @param mi         minimap2 index
 * @param l_seq      length of the query sequence
 * @param seq        the query sequence
 * @param n_regs     number of hits (out)
 * @param b          thread-local buffer; two mm_map() calls shall not use one buffer at the same time!
 * @param opt        mapping parameters
 * @param name       query name, used for all-vs-all overlapping and debugging
 *
 * @return an array of hits which need to be deallocated with free() together
 *         with mm_reg1_t::p of each element. The size is written to _n_regs_.
 */
mm_reg1_t *mm_map(const mm_idx_t *mi, int l_seq, const char *seq, int *n_regs, mm_tbuf_t *b, const mm_mapopt_t *opt, const char *name);
```

この関数でマッピングを行っています。ただし、Pythonバインディングのmappyでは、この`mm_map`関数を直接使用せずに `mm_map_aux`なる関数を使用しています。その理由は私にはまだよくわかっていません。

以上、ざっと `minimap.h` の構造体を眺めました。

この記事は以上です。