r1256: working for one example, but still buggy

jump.c

@@ -1,3 +1,4 @@
+#include <stdio.h>
 #include "mmpriv.h"
 #include "kalloc.h"
 
@@ -22,17 +23,26 @@ static int32_t mm_jump_check(void *km, const mm_idx_t *mi, int32_t qlen, const u
 
 static uint8_t *mm_jump_get_qseq_seq(void *km, int32_t qlen, const uint8_t *qseq0, const mm_reg1_t *r, int32_t is_left, int32_t ql0, uint8_t *qseq)
 {
-	int32_t i, k;
+	extern unsigned char seq_nt4_table[256];
+	int32_t i, k = 0;
 	if (!r->rev) {
-		if (is_left) memcpy(qseq, qseq0, ql0);
-		else memcpy(qseq, &qseq0[qlen - ql0], ql0);
+		if (is_left)
+			for (i = 0; i < ql0; ++i)
+				qseq[k++] = seq_nt4_table[(uint8_t)qseq0[i]];
+		else
+			for (i = qlen - ql0; i < qlen; ++i)
+				qseq[k++] = seq_nt4_table[(uint8_t)qseq0[i]];
 	} else {
 		if (is_left)
-			for (i = qlen - 1, k = 0; i >= qlen - ql0; --i)
-				qseq[k++] = qseq0[i] >= 4? qseq0[i] : 3 - qseq0[i];
+			for (i = qlen - 1; i >= qlen - ql0; --i) {
+				uint8_t c = seq_nt4_table[(uint8_t)qseq0[i]];
+				qseq[k++] = c >= 4? c : 3 - c;
+			}
 		else
-			for (i = ql0 - 1, k = 0; i >= 0; --i)
-				qseq[k++] = qseq0[i] >= 4? qseq0[i] : 3 - qseq0[i];
+			for (i = ql0 - 1; i >= 0; --i) {
+				uint8_t c = seq_nt4_table[(uint8_t)qseq0[i]];
+				qseq[k++] = c >= 4? c : 3 - c;
+			}
 	}
 	return qseq;
 }
@@ -73,7 +83,7 @@ static void mm_jump_split_left(void *km, const mm_idx_t *mi, const mm_mapopt_t *
 		for (mm2 = 0; j < clip + ext; ++j)
 			if (qseq[j] != tseq[j] || qseq[j] > 3 || tseq[j] > 3)
 				++mm2;
-		if (mm1 == 0 && mm2 == 1)
+		if (mm1 == 0 && mm2 <= 1)
 			i0 = i, ++m; // i0 points to the rightmost i
 	}
 	kfree(km, tseq);
@@ -113,7 +123,7 @@ static void mm_jump_split_right(void *km, const mm_idx_t *mi, const mm_mapopt_t
 	for (i = 0; i < n; ++i) { // traverse possible jumps
 		const mm_idx_jjump1_t *ai = &a[i];
 		int32_t tlen, tl1, j, mm1, mm2;
-		assert(ai->off >= r->rs - extt && ai->off < r->rs + ext);
+		assert(ai->off >= r->re - ext && ai->off < r->re + extt);
 		if (ts_strand * ai->strand < 0) continue; // wrong strand
 		if (ai->off2 <= ai->off) continue; // wrong direction
 		if (ai->off2 + clip + ext > mi->seq[r->rid].len) continue; // not long enough
@@ -133,7 +143,7 @@ static void mm_jump_split_right(void *km, const mm_idx_t *mi, const mm_mapopt_t
 		for (mm1 = 0; j < clip + ext; ++j)
 			if (qseq[j] != tseq[j] || qseq[j] > 3 || tseq[j] > 3)
 				++mm1;
-		if (mm1 == 0 && mm2 == 1)
+		if (mm1 == 0 && mm2 <= 1)
 			i0 = i0 >= 0? i0 : i, ++m; // i0 points to the leftmost i
 	}
 	kfree(km, tseq);

main.c

@@ -82,6 +82,7 @@ static ko_longopt_t long_options[] = {
 	{ "spsc",           ko_required_argument, 357 },
 	{ "junc-pen",       ko_required_argument, 358 },
 	{ "pe-ind-chain",   ko_no_argument,       359 },
+	{ "jump-bed",       ko_required_argument, 360 },
 	{ "dbg-seed-occ",   ko_no_argument,       501 },
 	{ "help",           ko_no_argument,       'h' },
 	{ "max-intron-len", ko_required_argument, 'G' },
@@ -131,7 +132,7 @@ int main(int argc, char *argv[])
 	mm_mapopt_t opt;
 	mm_idxopt_t ipt;
 	int i, c, n_threads = 3, n_parts, old_best_n = -1;
-	char *fnw = 0, *rg = 0, *junc_bed = 0, *fn_spsc = 0, *s, *alt_list = 0;
+	char *fnw = 0, *rg = 0, *junc_bed = 0, *jump_bed = 0, *fn_spsc = 0, *s, *alt_list = 0;
 	FILE *fp_help = stderr;
 	mm_idx_reader_t *idx_rdr;
 	mm_idx_t *mi;
@@ -254,6 +255,7 @@ int main(int argc, char *argv[])
 		else if (c == 356) opt.rmq_inner_dist = mm_parse_num(o.arg); // --rmq-inner
 		else if (c == 357) fn_spsc = o.arg; // --spsc
 		else if (c == 359) opt.flag |= MM_F_PE_IND; // --pe-ind-chain
+		else if (c == 360) jump_bed = o.arg; // --jump-bed
 		else if (c == 501) mm_dbg_flag |= MM_DBG_SEED_FREQ; // --dbg-seed-occ
 		else if (c == 330) {
 			fprintf(stderr, "[WARNING] \033[1;31m --lj-min-ratio has been deprecated.\033[0m\n");
@@ -440,6 +442,11 @@ int main(int argc, char *argv[])
 			if (mi->I == 0 && mm_verbose >= 2)
 				fprintf(stderr, "[WARNING] failed to load the junction BED file\n");
 		}
+		if (jump_bed) {
+			mm_idx_bed_read2(mi, jump_bed, 1, 0, 1);
+			if (mi->J == 0 && mm_verbose >= 2)
+				fprintf(stderr, "[WARNING] failed to load the jump BED file\n");
+		}
 		if (fn_spsc) {
 			mm_idx_spsc_read(mi, fn_spsc, mm_max_spsc_bonus(&opt));
 			if (mi->spsc == 0 && mm_verbose >= 2)

map.c

@@ -359,6 +359,10 @@ void mm_map_frag_core(const mm_idx_t *mi, int n_segs, const int *qlens, const ch
 	kfree(b->km, u);
 	kfree(b->km, mini_pos);
 
+	if (mi->J && n_segs == 1 && is_splice)
+		for (i = 0; i < n_regs0; ++i)
+			mm_jump_split(b->km, mi, opt, qlens[0], (const uint8_t*)seqs[0], &regs0[i], 0);
+
 	if (b->km) {
 		km_stat(b->km, &kmst);
 		if (mm_dbg_flag & MM_DBG_PRINT_QNAME)

minimap.h

@@ -5,7 +5,7 @@
 #include <stdio.h>
 #include <sys/types.h>
 
-#define MM_VERSION "2.28-r1255-dirty"
+#define MM_VERSION "2.28-r1256-dirty"
 
 #define MM_F_NO_DIAG       (0x001LL) // no exact diagonal hit
 #define MM_F_NO_DUAL       (0x002LL) // skip pairs where query name is lexicographically larger than target name

mmpriv.h

@@ -84,6 +84,7 @@ const uint64_t *mm_idx_get(const mm_idx_t *mi, uint64_t minier, int *n);
 int32_t mm_idx_cal_max_occ(const mm_idx_t *mi, float f);
 int mm_idx_getseq2(const mm_idx_t *mi, int is_rev, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq);
 mm_reg1_t *mm_gen_regs(void *km, uint32_t hash, int qlen, int n_u, uint64_t *u, mm128_t *a, int is_qstrand);
+int mm_idx_bed_read2(mm_idx_t *mi, const char *fn, int read_junc, int for_score, int for_jump);
 const mm_idx_jjump1_t *mm_idx_jump_get(const mm_idx_t *db, int32_t cid, int32_t st, int32_t en, int32_t *n);
 
 mm128_t *mg_lchain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int max_iter, int min_cnt, int min_sc, float chn_pen_gap, float chn_pen_skip,
@@ -115,6 +116,8 @@ mm_seg_t *mm_seg_gen(void *km, uint32_t hash, int n_segs, const int *qlens, int
 void mm_seg_free(void *km, int n_segs, mm_seg_t *segs);
 void mm_pair(void *km, int max_gap_ref, int dp_bonus, int sub_diff, int match_sc, const int *qlens, int *n_regs, mm_reg1_t **regs);
 
+void mm_jump_split(void *km, const mm_idx_t *mi, const mm_mapopt_t *opt, int32_t qlen, const uint8_t *qseq, mm_reg1_t *r, int32_t ts_strand);
+
 FILE *mm_split_init(const char *prefix, const mm_idx_t *mi);
 mm_idx_t *mm_split_merge_prep(const char *prefix, int n_splits, FILE **fp, uint32_t *n_seq_part);
 int mm_split_merge(int n_segs, const char **fn, const mm_mapopt_t *opt, int n_split_idx);

options.c

@@ -63,7 +63,7 @@ void mm_mapopt_init(mm_mapopt_t *opt)
 	opt->pe_ori = 0; // FF
 	opt->pe_bonus = 33;
 
-	opt->jump_min_alen = 10;
+	opt->jump_min_alen = 5;
 }
 
 void mm_mapopt_update(mm_mapopt_t *opt, const mm_idx_t *mi)