From 14b853499f67333aa4bf43497e708cc95a2b4ed0 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Thu, 14 Sep 2017 22:44:10 -0400 Subject: [PATCH] r369: updated example with the latest API --- example.c | 67 +++++++++++++++++++++++++++---------------------------- main.c | 7 +++--- 2 files changed, 36 insertions(+), 38 deletions(-) diff --git a/example.c b/example.c index ccb1e68..3eac4c2 100644 --- a/example.c +++ b/example.c @@ -11,53 +11,52 @@ 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); + mopt.flag |= MM_F_CIGAR; // perform alignment if (argc < 3) { fprintf(stderr, "Usage: minimap2-lite \n"); return 1; } - + // open query file for reading; you may use your favorite FASTA/Q parser gzFile f = gzopen(argv[2], "r"); assert(f); kseq_t *ks = kseq_init(f); - // create index for target; we are creating one index for all target sequence - int n_threads = 4, w = 10, k = 15, is_hpc = 0; - mm_idx_t *mi = mm_idx_build(argv[1], w, k, is_hpc, n_threads); - assert(mi); - - // mapping - mm_mapopt_t opt; - mm_mapopt_init(&opt); // initialize mapping parameters - mm_mapopt_update(&opt, mi); // this sets the maximum minimizer occurrence; TODO: set a better default in mm_mapopt_init()! - opt.flag |= MM_F_CIGAR; // perform alignment - mm_tbuf_t *tbuf = mm_tbuf_init(); // thread buffer; for multi-threading, allocate one tbuf for each thread - while (kseq_read(ks) >= 0) { // each kseq_read() call reads one query sequence - mm_reg1_t *reg; - int j, i, n_reg; - // get all hits for the query - reg = mm_map(mi, ks->seq.l, ks->seq.s, &n_reg, tbuf, &opt, 0); - // traverse hits and print them out - for (j = 0; j < n_reg; ++j) { - mm_reg1_t *r = ®[j]; - assert(r->p); // with MM_F_CIGAR, this should not be NULL - 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\tcg:Z:", mi->seq[r->rid].name, mi->seq[r->rid].len, r->rs, r->re, - r->p->blen - r->p->n_ambi - r->p->n_diff, r->p->blen, r->mapq); - for (i = 0; i < r->p->n_cigar; ++i) // IMPORTANT: this gives the CIGAR in the aligned regions. NO soft/hard clippings! - printf("%d%c", r->p->cigar[i]>>4, "MIDSHN"[r->p->cigar[i]&0xf]); - putchar('\n'); - free(r->p); + // open index reader + mm_idx_reader_t *r = mm_idx_reader_open(argv[1], &iopt, 0); + mm_idx_t *mi; + while ((mi = mm_idx_reader_read(r, n_threads)) != 0) { // traverse each part of the index + mm_mapopt_update(&mopt, mi); // this sets the maximum minimizer occurrence; TODO: set a better default in mm_mapopt_init()! + mm_tbuf_t *tbuf = mm_tbuf_init(); // thread buffer; for multi-threading, allocate one tbuf for each thread + while (kseq_read(ks) >= 0) { // each kseq_read() call reads one query sequence + mm_reg1_t *reg; + int j, i, n_reg; + reg = mm_map(mi, ks->seq.l, ks->seq.s, &n_reg, tbuf, &mopt, 0); // get all hits for the query + for (j = 0; j < n_reg; ++j) { // traverse hits and print them out + mm_reg1_t *r = ®[j]; + assert(r->p); // with MM_F_CIGAR, this should not be NULL + 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\tcg:Z:", mi->seq[r->rid].name, mi->seq[r->rid].len, r->rs, r->re, + r->p->blen - r->p->n_ambi - r->p->n_diff, r->p->blen, r->mapq); + for (i = 0; i < r->p->n_cigar; ++i) // IMPORTANT: this gives the CIGAR in the aligned regions. NO soft/hard clippings! + printf("%d%c", r->p->cigar[i]>>4, "MIDSHN"[r->p->cigar[i]&0xf]); + putchar('\n'); + free(r->p); + } + free(reg); } - free(reg); + mm_tbuf_destroy(tbuf); + mm_idx_destroy(mi); } - mm_tbuf_destroy(tbuf); - - // deallocate index and close the query file - mm_idx_destroy(mi); - kseq_destroy(ks); + mm_idx_reader_close(r); // close the index reader + kseq_destroy(ks); // close the query file gzclose(f); return 0; } diff --git a/main.c b/main.c index 31dee73..853da1f 100644 --- a/main.c +++ b/main.c @@ -6,7 +6,7 @@ #include "mmpriv.h" #include "getopt.h" -#define MM_VERSION "2.1.1-r368-dirty" +#define MM_VERSION "2.1.1-r369-dirty" #ifdef __linux__ #include @@ -67,6 +67,7 @@ int main(int argc, char *argv[]) char *fnw = 0, *rg = 0, *s; FILE *fp_help = stderr; mm_idx_reader_t *idx_rdr; + mm_idx_t *mi; liftrlimit(); mm_realtime0 = realtime(); @@ -210,9 +211,7 @@ int main(int argc, char *argv[]) } if (opt.flag & MM_F_OUT_SAM) mm_write_sam_hdr_no_SQ(rg, MM_VERSION, argc, argv); - for (;;) { - mm_idx_t *mi; - if ((mi = mm_idx_reader_read(idx_rdr, n_threads)) == 0) break; + while ((mi = mm_idx_reader_read(idx_rdr, n_threads)) != 0) { if (mm_verbose >= 2 && idx_rdr->n_parts > 1 && (opt.flag&MM_F_OUT_SAM) && !(opt.flag&MM_F_NO_SAM_SQ)) fprintf(stderr, "[WARNING] \033[1;31mSAM output is malformated due to internal @SQ lines. Please add option --no-sam-sq or filter afterwards.\033[0m\n"); if (mm_verbose >= 3)