diff --git a/format.c b/format.c index 5a6d8ae..2e12669 100644 --- a/format.c +++ b/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; diff --git a/main.c b/main.c index 1e66df5..2ca8e69 100644 --- a/main.c +++ b/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"); diff --git a/map.c b/map.c index a62eab8..e2cb1c5 100644 --- a/map.c +++ b/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; diff --git a/minimap.h b/minimap.h index cdae5ec..b53a0ed 100644 --- a/minimap.h +++ b/minimap.h @@ -5,7 +5,7 @@ #include #include -#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 diff --git a/mmpriv.h b/mmpriv.h index c02a7f3..4b47e8e 100644 --- a/mmpriv.h +++ b/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);