self-chain.c

// To compile:
//   gcc -g -O2 example.c libminimap2.a -lz

#include <stdlib.h>
#include <assert.h>
#include <stdio.h>
#include <zlib.h>
#include "minimap.h"
#include "kseq.h"
KSEQ_INIT(gzFile, gzread)

int main(int argc, char *argv[])
{
	mm_idxopt_t iopt;
	mm_mapopt_t mopt;
	int n_threads = 3;

	mm_verbose = 2; // disable message output to stderr
	mm_set_opt(0, &iopt, &mopt);
	iopt.w = 11, iopt.k = 19, iopt.flag |= MM_I_HPC;
	mopt.flag |= MM_F_NO_LJOIN | MM_F_ALL_CHAINS | MM_F_NO_DIAG | MM_F_FOR_ONLY;
	mopt.max_chain_skip = 50;
	mopt.mid_occ = 200;
	mopt.min_chain_score = 200;
	mopt.bw = 50;
	mopt.max_gap = 1000;

	if (argc < 2) {
		fprintf(stderr, "Usage: self-chain <in.fa>\n");
		return 1;
	}

	gzFile f = gzopen(argv[1], "r");
	assert(f);
	kseq_t *ks = kseq_init(f);

	mm_tbuf_t *tbuf = mm_tbuf_init();
	while (kseq_read(ks) >= 0) {
		mm_idx_t *mi = mm_idx_str(iopt.w, iopt.k, iopt.flag&MM_I_HPC, 10, 1, (const char**)&ks->seq.s, (const char**)&ks->name.s);
		int j, n_reg;
		mm_reg1_t *reg = mm_map(mi, ks->seq.l, ks->seq.s, &n_reg, tbuf, &mopt, ks->name.s);
		for (j = 0; j < n_reg; ++j) {
			mm_reg1_t *r = &reg[j];
			printf("%s\t%d\t%d\t%d\t%c\t", ks->name.s, ks->seq.l, r->qs, r->qe, "+-"[r->rev]);
			printf("%s\t%d\t%d\t%d\t%d\t%d\t%d\ts1:i:%d\tdv:f:%.4f\n", mi->seq[r->rid].name, mi->seq[r->rid].len, r->rs, r->re, r->mlen, r->blen, r->mapq, r->score, r->div);
		}
		free(reg);
		mm_idx_destroy(mi);
	}
	kseq_destroy(ks); // close the query file
	gzclose(f);
	return 0;
}