r1269: added --write-junc

in preparation for 2-pass

Author: lh3

format.c

@@ -253,6 +253,52 @@ static void write_cs_ds_core(kstring_t *s, const uint8_t *tseq, const uint8_t *q
 	assert(t_off == r->re - r->rs && q_off == r->qe - r->qs);
 }
 
+static inline void revcomp_splice(uint8_t s[2])
+{
+	uint8_t c = s[1] < 4? 3 - s[1] : 4;
+	s[1] = s[0] < 4? 3 - s[0] : 4;
+	s[0] = c;
+}
+
+void mm_write_junc(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r)
+{
+	int32_t i, t_off, swritten = 0;
+	s->l = 0;
+	if (!r->is_spliced || r->p == 0) return; // no junctions
+	if (r->p->trans_strand != 1 && r->p->trans_strand != 2) return; // no preferred strand
+	for (i = 0, t_off = r->rs; i < (int)r->p->n_cigar; ++i) {
+		int op = r->p->cigar[i]&0xf, len = r->p->cigar[i]>>4;
+		if (op == MM_CIGAR_MATCH || op == MM_CIGAR_EQ_MATCH || op == MM_CIGAR_X_MISMATCH || op == MM_CIGAR_DEL) {
+			t_off += len;
+		} else if (op == MM_CIGAR_N_SKIP) { // intron
+			uint8_t donor[2], acceptor[2];
+			int32_t score1 = 0, score2 = 0, rev;
+			assert(len >= 2);
+			rev = (r->p->trans_strand == 2) ^ r->rev;
+			if (!rev) {
+				mm_idx_getseq(mi, r->rid, t_off, t_off + 2, donor);
+				mm_idx_getseq(mi, r->rid, t_off + len - 2, t_off + len, acceptor);
+			} else {
+				mm_idx_getseq(mi, r->rid, t_off, t_off + 2, acceptor);
+				mm_idx_getseq(mi, r->rid, t_off + len - 2, t_off + len, donor);
+				revcomp_splice(donor);
+				revcomp_splice(acceptor);
+			}
+			//fprintf(stderr, "%c%c-%c%c\n", "ACGTN"[donor[0]], "ACGTN"[donor[1]], "ACGTN"[acceptor[0]], "ACGTN"[acceptor[1]]);
+			if (donor[0] == 2 && donor[1] == 3) score1 = 3;
+			else if (donor[0] == 2 && donor[1] == 1) score1 = 2;
+			else if (donor[0] == 0 && donor[1] == 3) score1 = 1;
+			if (acceptor[0] == 0 && acceptor[1] == 2) score2 = 3;
+			else if (acceptor[0] == 0 && acceptor[1] == 1) score2 = 1;
+			if (swritten) mm_sprintf_lite(s, "\n");
+			else swritten = 1;
+			mm_sprintf_lite(s, "%s\t%d\t%d\t%s\t%d\t%c", mi->seq[r->rid].name, t_off, t_off + len, t->name, score1 + score2, "+-"[rev]);
+			t_off += len;
+		}
+	}
+	assert(t_off == r->re);
+}
+
 static void write_MD_core(kstring_t *s, const uint8_t *tseq, const uint8_t *qseq, const mm_reg1_t *r, char *tmp, int write_tag)
 {
 	int i, q_off, t_off, l_MD = 0;

main.c

@@ -83,6 +83,7 @@ static ko_longopt_t long_options[] = {
 	{ "junc-pen",       ko_required_argument, 358 },
 	{ "pairing",        ko_required_argument, 359 },
 	{ "jump-min-match", ko_required_argument, 360 },
+	{ "write-junc",     ko_no_argument,       361 },
 	{ "dbg-seed-occ",   ko_no_argument,       501 },
 	{ "help",           ko_no_argument,       'h' },
 	{ "max-intron-len", ko_required_argument, 'G' },
@@ -255,6 +256,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 == 360) opt.jump_min_match = mm_parse_num(o.arg); // --jump-min-match
+		else if (c == 361) opt.flag |= MM_F_OUT_JUNC | MM_F_CIGAR; // --write-junc
 		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");

map.c

@@ -594,9 +594,16 @@ static void *worker_pipeline(void *shared, int step, void *in)
 							mm_err_fwrite(r->p, r->p->capacity, 4, p->fp_split);
 						}
 					}
+				} else if (p->opt->flag & MM_F_OUT_JUNC) { // extra logic for --write-junc
+					for (j = 0; j < s->n_reg[i]; ++j) {
+						const mm_reg1_t *r = &s->reg[i][j];
+						if (r->id != r->parent) continue;
+						mm_write_junc(&p->str, mi, t, r);
+						if (p->str.l > 0) mm_err_puts(p->str.s);
+					}
 				} else if (s->n_reg[i] > 0) { // the query has at least one hit
 					for (j = 0; j < s->n_reg[i]; ++j) {
-						mm_reg1_t *r = &s->reg[i][j];
+						const mm_reg1_t *r = &s->reg[i][j];
 						assert(!r->sam_pri || r->id == r->parent);
 						if ((p->opt->flag & MM_F_NO_PRINT_2ND) && r->id != r->parent)
 							continue;

minimap.h

@@ -5,7 +5,7 @@
 #include <stdio.h>
 #include <sys/types.h>
 
-#define MM_VERSION "2.28-r1268-dirty"
+#define MM_VERSION "2.28-r1269-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
@@ -47,6 +47,7 @@
 #define MM_F_OUT_DS        (0x2000000000LL)
 #define MM_F_WEAK_PAIRING  (0x4000000000LL)
 #define MM_F_SR_RNA        (0x8000000000LL)
+#define MM_F_OUT_JUNC      (0x10000000000LL)
 
 #define MM_I_HPC          0x1
 #define MM_I_NO_SEQ       0x2

mmpriv.h

@@ -78,6 +78,7 @@ void mm_write_paf4(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const
 void mm_write_sam(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, int n_regs, const mm_reg1_t *regs);
 void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int seg_idx, int reg_idx, int n_seg, const int *n_regs, const mm_reg1_t *const* regs, void *km, int64_t opt_flag);
 void mm_write_sam3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int seg_idx, int reg_idx, int n_seg, const int *n_regss, const mm_reg1_t *const* regss, void *km, int64_t opt_flag, int rep_len);
+void mm_write_junc(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r);
 
 void mm_idxopt_init(mm_idxopt_t *opt);
 const uint64_t *mm_idx_get(const mm_idx_t *mi, uint64_t minier, int *n);