diff --git a/Makefile b/Makefile index 4118616..90b1cc4 100644 --- a/Makefile +++ b/Makefile @@ -8,6 +8,10 @@ PROG= minimap2 PROG_EXTRA= sdust minimap2-lite LIBS= -lm -lz -lpthread +ifneq ($(aarch64),) + arm_neon=1 +endif + ifeq ($(arm_neon),) # if arm_neon is not defined ifeq ($(sse2only),) # if sse2only is not defined OBJS+=ksw2_extz2_sse41.o ksw2_extd2_sse41.o ksw2_exts2_sse41.o ksw2_extz2_sse2.o ksw2_extd2_sse2.o ksw2_exts2_sse2.o ksw2_dispatch.o diff --git a/index.c b/index.c index 0d8a2ae..c9bd01f 100644 --- a/index.c +++ b/index.c @@ -358,7 +358,7 @@ static void *worker_pipeline(void *shared, int step, void *in) for (i = 0; i < s->n_seq; ++i) { mm_bseq1_t *t = &s->seq[i]; if (t->l_seq > 0) - mm_sketch(0, t->seq, t->l_seq, p->mi->w, p->mi->k, t->rid, p->mi->flag&MM_I_HPC, &s->a); + mm_sketch2(0, t->seq, t->l_seq, p->mi->w, p->mi->k, t->rid, p->mi->flag&MM_I_HPC, p->mi->flag&MM_I_SYNCMER, &s->a); else if (mm_verbose >= 2) fprintf(stderr, "[WARNING] the length database sequence '%s' is 0\n", t->name); free(t->seq); free(t->name); @@ -446,7 +446,7 @@ mm_idx_t *mm_idx_str(int w, int k, int is_hpc, int bucket_bits, int n, const cha sum_len += p->len; if (p->len > 0) { a.n = 0; - mm_sketch(0, s, p->len, w, k, i, is_hpc, &a); + mm_sketch2(0, s, p->len, w, k, i, is_hpc, 0, &a); // TODO: mm_idx_str() doesn't support syncmer mm_idx_add(mi, a.n, a.a); } } diff --git a/lchain.c b/lchain.c index d004157..2da543c 100644 --- a/lchain.c +++ b/lchain.c @@ -251,7 +251,7 @@ mm128_t *mg_lchain_rmq(int max_dist, int max_dist_inner, int bw, int max_chn_ski int64_t n, mm128_t *a, int *n_u_, uint64_t **_u, void *km) { int32_t *f,*t, *v, n_u, n_v, mmax_f = 0, max_rmq_size = 0, max_drop = bw; - int64_t *p, i, i0, st = 0, st_inner = 0, n_iter = 0; + int64_t *p, i, i0, st = 0, st_inner = 0; uint64_t *u; lc_elem_t *root = 0, *root_inner = 0; void *mem_mp = 0; @@ -345,7 +345,6 @@ mm128_t *mg_lchain_rmq(int max_dist, int max_dist_inner, int bw, int max_chn_ski } if (!krmq_itr_prev(lc_elem, &itr)) break; } - n_iter += n_rmq_iter; } } } diff --git a/main.c b/main.c index 0be9933..af5ffa5 100644 --- a/main.c +++ b/main.c @@ -7,7 +7,7 @@ #include "mmpriv.h" #include "ketopt.h" -#define MM_VERSION "2.24-r1122" +#define MM_VERSION "2.24-r1149-dirty" #ifdef __linux__ #include @@ -121,7 +121,7 @@ static inline void yes_or_no(mm_mapopt_t *opt, int64_t flag, int long_idx, const int main(int argc, char *argv[]) { - const char *opt_str = "2aSDw:k:K:t:r:f:Vv:g:G:I:d:XT:s:x:Hcp:M:n:z:A:B:O:E:m:N:Qu:R:hF:LC:yYPo:e:U:"; + const char *opt_str = "2aSDw:k:K:t:r:f:Vv:g:G:I:d:XT:s:x:Hcp:M:n:z:A:B:O:E:m:N:Qu:R:hF:LC:yYPo:e:U:j:"; ketopt_t o = KETOPT_INIT; mm_mapopt_t opt; mm_idxopt_t ipt; @@ -153,7 +153,8 @@ int main(int argc, char *argv[]) o = KETOPT_INIT; while ((c = ketopt(&o, argc, argv, 1, opt_str, long_options)) >= 0) { - if (c == 'w') ipt.w = atoi(o.arg); + if (c == 'w') ipt.w = atoi(o.arg), ipt.flag &= ~MM_I_SYNCMER; + else if (c == 'j') ipt.w = atoi(o.arg), ipt.flag |= MM_I_SYNCMER; else if (c == 'k') ipt.k = atoi(o.arg); else if (c == 'H') ipt.flag |= MM_I_HPC; else if (c == 'd') fnw = o.arg; // the above are indexing related options, except -I @@ -322,6 +323,7 @@ int main(int argc, char *argv[]) fprintf(fp_help, " -H use homopolymer-compressed k-mer (preferrable for PacBio)\n"); fprintf(fp_help, " -k INT k-mer size (no larger than 28) [%d]\n", ipt.k); fprintf(fp_help, " -w INT minimizer window size [%d]\n", ipt.w); + fprintf(fp_help, " -j INT syncmer submer size (overriding -w) []\n"); fprintf(fp_help, " -I NUM split index for every ~NUM input bases [4G]\n"); fprintf(fp_help, " -d FILE dump index to FILE []\n"); fprintf(fp_help, " Mapping:\n"); diff --git a/map.c b/map.c index 5311468..2342c9e 100644 --- a/map.c +++ b/map.c @@ -67,7 +67,7 @@ static void collect_minimizers(void *km, const mm_mapopt_t *opt, const mm_idx_t mv->n = 0; for (i = n = 0; i < n_segs; ++i) { size_t j; - mm_sketch(km, seqs[i], qlens[i], mi->w, mi->k, i, mi->flag&MM_I_HPC, mv); + mm_sketch2(km, seqs[i], qlens[i], mi->w, mi->k, i, mi->flag&MM_I_HPC, mi->flag&MM_I_SYNCMER, mv); for (j = n; j < mv->n; ++j) mv->a[j].y += sum << 1; if (opt->sdust_thres > 0) // mask low-complexity minimizers diff --git a/minimap.h b/minimap.h index 13e12e0..3fb6900 100644 --- a/minimap.h +++ b/minimap.h @@ -44,6 +44,7 @@ #define MM_I_HPC 0x1 #define MM_I_NO_SEQ 0x2 #define MM_I_NO_NAME 0x4 +#define MM_I_SYNCMER 0x8 #define MM_IDX_MAGIC "MMI\2" diff --git a/minimap2.1 b/minimap2.1 index 8a759db..08d9c94 100644 --- a/minimap2.1 +++ b/minimap2.1 @@ -79,6 +79,19 @@ Minimizer k-mer length [15] .BI -w \ INT Minimizer window size [10]. A minimizer is the smallest k-mer in a window of w consecutive k-mers. +.TP +.BI -j \ INT +Syncmer submer size [10]. Option +.B -j +and +.B -w +will override each: if +.B -w +is applied after +.BR -j , +.B -j +will have no effect, and vice versa. + .TP .B -H Use homopolymer-compressed (HPC) minimizers. An HPC sequence is constructed by diff --git a/mmpriv.h b/mmpriv.h index 2f5034b..7b51b98 100644 --- a/mmpriv.h +++ b/mmpriv.h @@ -60,6 +60,8 @@ void radix_sort_64(uint64_t *beg, uint64_t *end); uint32_t ks_ksmall_uint32_t(size_t n, uint32_t arr[], size_t kk); void mm_sketch(void *km, const char *str, int len, int w, int k, uint32_t rid, int is_hpc, mm128_v *p); +void mm_sketch_syncmer(void *km, const char *str, int len, int smer, int k, uint32_t rid, int is_hpc, mm128_v *p); +void mm_sketch2(void *km, const char *str, int len, int w, int k, uint32_t rid, int is_hpc, int is_syncmer, mm128_v *p); mm_seed_t *mm_collect_matches(void *km, int *_n_m, int qlen, int max_occ, int max_max_occ, int dist, const mm_idx_t *mi, const mm128_v *mv, int64_t *n_a, int *rep_len, int *n_mini_pos, uint64_t **mini_pos); void mm_seed_mz_flt(void *km, mm128_v *mv, int32_t q_occ_max, float q_occ_frac); diff --git a/sketch.c b/sketch.c index f830693..42ab4f7 100644 --- a/sketch.c +++ b/sketch.c @@ -141,3 +141,58 @@ void mm_sketch(void *km, const char *str, int len, int w, int k, uint32_t rid, i if (min.x != UINT64_MAX) kv_push(mm128_t, km, *p, min); } + +void mm_sketch_syncmer(void *km, const char *str, int len, int smer, int k, uint32_t rid, int is_hpc, mm128_v *p) +{ + uint64_t shift1 = 2 * (k - 1), mask = (1ULL<<2*k) - 1, smask = (1ULL<<2*smer) - 1, kmer[2] = {0,0}; + int i, j, l, buf_pos, min_pos, kmer_span = 0; + tiny_queue_t tq; + + assert(len > 0 && (smer > 0 && smer <= k) && (k > 0 && k <= 28)); // 56 bits for k-mer; could use long k-mers, but 28 enough in practice + memset(&tq, 0, sizeof(tiny_queue_t)); + kv_resize(mm128_t, km, *p, p->n + len/(k - smer)); + + for (i = l = buf_pos = min_pos = 0; i < len; ++i) { + int c = seq_nt4_table[(uint8_t)str[i]]; + if (c < 4) { // not an ambiguous base + int z; + if (is_hpc) { + int skip_len = 1; + if (i + 1 < len && seq_nt4_table[(uint8_t)str[i + 1]] == c) { + for (skip_len = 2; i + skip_len < len; ++skip_len) + if (seq_nt4_table[(uint8_t)str[i + skip_len]] != c) + break; + i += skip_len - 1; // put $i at the end of the current homopolymer run + } + tq_push(&tq, skip_len); + kmer_span += skip_len; + if (tq.count > k) kmer_span -= tq_shift(&tq); + } else kmer_span = l + 1 < k? l + 1 : k; + kmer[0] = (kmer[0] << 2 | c) & mask; // forward k-mer + kmer[1] = (kmer[1] >> 2) | (3ULL^c) << shift1; // reverse k-mer + if (kmer[0] == kmer[1]) continue; // skip "symmetric k-mers" as we don't know it strand + z = kmer[0] < kmer[1]? 0 : 1; // strand + ++l; + if (l >= k && kmer_span < 256) { + uint64_t x, min = UINT64_MAX; + x = hash64(kmer[z], mask); + for (j = 0; j < k - smer; ++j) { + uint64_t y = x >> (j + j) & smask; + min = min < y? min : y; + } + if ((x & smask) == min) { + mm128_t t; + t.x = x << 8 | kmer_span; + t.y = (uint64_t)rid<<32 | (uint32_t)i<<1 | z; + kv_push(mm128_t, km, *p, t); + } + } + } else l = 0, tq.count = tq.front = 0, kmer_span = 0; + } +} + +void mm_sketch2(void *km, const char *str, int len, int w, int k, uint32_t rid, int is_hpc, int is_syncmer, mm128_v *p) +{ + if (is_syncmer) mm_sketch_syncmer(km, str, len, w, k, rid, is_hpc, p); + else mm_sketch(km, str, len, w, k, rid, is_hpc, p); +}