From da7109fd296b3a8ac8fc4ea633fdbbfc1209b74d Mon Sep 17 00:00:00 2001 From: Heng Li Date: Tue, 21 Apr 2020 12:37:35 -0400 Subject: [PATCH] r985: optionally report cs/cg on the query strand PAF only; not well tested --- align.c | 56 ++++++++++++++++++++++++++++++++++++++++--------------- format.c | 36 ++++++++++++++++++++++------------- hit.c | 18 +++++++++--------- index.c | 22 ++++++++++++++++++++++ main.c | 4 +++- map.c | 8 ++++++-- minimap.h | 2 ++ mmpriv.h | 5 +++-- options.c | 5 +++++ 9 files changed, 114 insertions(+), 42 deletions(-) diff --git a/align.c b/align.c index a7c4c3c..3d67eeb 100644 --- a/align.c +++ b/align.c @@ -533,8 +533,13 @@ static int mm_seed_ext_score(void *km, const mm_mapopt_t *opt, const mm_idx_t *m re = re + ext_len < (int32_t)mi->seq[rid].len? re + ext_len : mi->seq[rid].len; qe = qe + ext_len < qlen? qe + ext_len : qlen; tseq = (uint8_t*)kmalloc(km, re - rs); - mm_idx_getseq(mi, rid, rs, re, tseq); - qseq = qseq0[a->x>>63] + qs; + if (opt->flag & MM_F_QSTRAND) { + qseq = qseq0[0] + qs; + mm_idx_getseq2(mi, a->x>>63, rid, rs, re, tseq); + } else { + qseq = qseq0[a->x>>63] + qs; + mm_idx_getseq(mi, rid, rs, re, tseq); + } qp = ksw_ll_qinit(km, 2, qe - qs, qseq, 5, mat); score = ksw_ll_i16(qp, re - rs, tseq, opt->q, opt->e, &q_off, &t_off); kfree(km, tseq); @@ -688,8 +693,13 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int junc = (uint8_t*)kmalloc(km, re0 - rs0); if (qs > 0 && rs > 0) { // left extension; probably the condition can be changed to "qs > qs0 && rs > rs0" - qseq = &qseq0[rev][qs0]; - mm_idx_getseq(mi, rid, rs0, rs, tseq); + if (opt->flag & MM_F_QSTRAND) { + qseq = &qseq0[0][qs0]; + mm_idx_getseq2(mi, rev, rid, rs0, rs, tseq); + } else { + qseq = &qseq0[rev][qs0]; + mm_idx_getseq(mi, rid, rs0, rs, tseq); + } mm_idx_bed_junc(mi, rid, rs0, rs, junc); mm_seq_rev(qs - qs0, qseq); mm_seq_rev(rs - rs0, tseq); @@ -718,8 +728,13 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int if (a[as1+i].y & MM_SEED_LONG_JOIN) bw1 = qe - qs > re - rs? qe - qs : re - rs; // perform alignment - qseq = &qseq0[rev][qs]; - mm_idx_getseq(mi, rid, rs, re, tseq); + if (opt->flag & MM_F_QSTRAND) { + qseq = &qseq0[0][qs]; + mm_idx_getseq2(mi, rev, rid, rs, re, tseq); + } else { + qseq = &qseq0[rev][qs]; + mm_idx_getseq(mi, rid, rs, re, tseq); + } mm_idx_bed_junc(mi, rid, rs, re, junc); if (is_sr) { // perform ungapped alignment assert(qe - qs == re - rs); @@ -748,7 +763,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int re1 = rs + (ez->max_t + 1); qe1 = qs + (ez->max_q + 1); if (cnt1 - (j + 1) >= opt->min_cnt) { - mm_split_reg(r, r2, as1 + j + 1 - r->as, qlen, a); + mm_split_reg(r, r2, as1 + j + 1 - r->as, qlen, a, opt->flag&MM_F_QSTRAND); if (zdrop_code == 2) r2->split_inv = 1; } break; @@ -758,8 +773,13 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int } if (!dropped && qe < qe0 && re < re0) { // right extension - qseq = &qseq0[rev][qe]; - mm_idx_getseq(mi, rid, re, re0, tseq); + if (opt->flag & MM_F_QSTRAND) { + qseq = &qseq0[0][qe]; + mm_idx_getseq2(mi, rev, rid, re, re0, tseq); + } else { + qseq = &qseq0[rev][qe]; + mm_idx_getseq(mi, rid, re, re0, tseq); + } mm_idx_bed_junc(mi, rid, re, re0, junc); mm_align_pair(km, opt, qe0 - qe, qseq, re0 - re, tseq, junc, mat, bw, opt->end_bonus, opt->zdrop, extra_flag|KSW_EZ_EXTZ_ONLY, ez); if (ez->n_cigar > 0) { @@ -772,13 +792,19 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int assert(qe1 <= qlen); r->rs = rs1, r->re = re1; - if (rev) r->qs = qlen - qe1, r->qe = qlen - qs1; - else r->qs = qs1, r->qe = qe1; + if (!rev || (opt->flag & MM_F_QSTRAND)) r->qs = qs1, r->qe = qe1; + else r->qs = qlen - qe1, r->qe = qlen - qs1; assert(re1 - rs1 <= re0 - rs0); if (r->p) { - mm_idx_getseq(mi, rid, rs1, re1, tseq); - mm_update_extra(r, &qseq0[r->rev][qs1], tseq, mat, opt->q, opt->e, opt->flag & MM_F_EQX); + if (opt->flag & MM_F_QSTRAND) { + mm_idx_getseq2(mi, r->rev, rid, rs1, re1, tseq); + qseq = &qseq0[0][qs1]; + } else { + mm_idx_getseq(mi, rid, rs1, re1, tseq); + qseq = &qseq0[r->rev][qs1]; + } + mm_update_extra(r, qseq, tseq, mat, opt->q, opt->e, opt->flag & MM_F_EQX); if (rev && r->p->trans_strand) r->p->trans_strand ^= 3; // flip to the read strand } @@ -788,7 +814,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int } static int mm_align1_inv(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, uint8_t *qseq0[2], const mm_reg1_t *r1, const mm_reg1_t *r2, mm_reg1_t *r_inv, ksw_extz_t *ez) -{ +{ // NB: this doesn't work with the qstrand mode int tl, ql, score, ret = 0, q_off, t_off; uint8_t *tseq, *qseq; int8_t mat[25]; @@ -897,7 +923,7 @@ mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *m regs[i].p->trans_strand = opt->flag&MM_F_SPLICE_FOR? 1 : 2; } if (r2.cnt > 0) regs = mm_insert_reg(&r2, i, &n_regs, regs); - if (i > 0 && regs[i].split_inv) { + if (i > 0 && regs[i].split_inv && !(opt->flag & MM_F_NO_INV)) { if (mm_align1_inv(km, opt, mi, qlen, qseq0, ®s[i-1], ®s[i], &r2, &ez)) { regs = mm_insert_reg(&r2, i, &n_regs, regs); ++i; // skip the inserted INV alignment diff --git a/format.c b/format.c index 890335c..c3b8c09 100644 --- a/format.c +++ b/format.c @@ -217,7 +217,7 @@ static void write_MD_core(kstring_t *s, const uint8_t *tseq, const uint8_t *qseq assert(t_off == r->re - r->rs && q_off == r->qe - r->qs); } -static void write_cs_or_MD(void *km, kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, int no_iden, int is_MD, int write_tag) +static void write_cs_or_MD(void *km, kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, int no_iden, int is_MD, int write_tag, int is_qstrand) { extern unsigned char seq_nt4_table[256]; int i; @@ -227,14 +227,20 @@ static void write_cs_or_MD(void *km, kstring_t *s, const mm_idx_t *mi, const mm_ qseq = (uint8_t*)kmalloc(km, r->qe - r->qs); tseq = (uint8_t*)kmalloc(km, r->re - r->rs); tmp = (char*)kmalloc(km, r->re - r->rs > r->qe - r->qs? r->re - r->rs + 1 : r->qe - r->qs + 1); - mm_idx_getseq(mi, r->rid, r->rs, r->re, tseq); - if (!r->rev) { + if (is_qstrand) { + mm_idx_getseq2(mi, r->rev, r->rid, r->rs, r->re, tseq); for (i = r->qs; i < r->qe; ++i) qseq[i - r->qs] = seq_nt4_table[(uint8_t)t->seq[i]]; } else { - for (i = r->qs; i < r->qe; ++i) { - uint8_t c = seq_nt4_table[(uint8_t)t->seq[i]]; - qseq[r->qe - i - 1] = c >= 4? 4 : 3 - c; + mm_idx_getseq(mi, r->rid, r->rs, r->re, tseq); + if (!r->rev) { + for (i = r->qs; i < r->qe; ++i) + qseq[i - r->qs] = seq_nt4_table[(uint8_t)t->seq[i]]; + } else { + for (i = r->qs; i < r->qe; ++i) { + uint8_t c = seq_nt4_table[(uint8_t)t->seq[i]]; + qseq[r->qe - i - 1] = c >= 4? 4 : 3 - c; + } } } if (is_MD) write_MD_core(s, tseq, qseq, r, tmp, write_tag); @@ -242,14 +248,14 @@ static void write_cs_or_MD(void *km, kstring_t *s, const mm_idx_t *mi, const mm_ kfree(km, qseq); kfree(km, tseq); kfree(km, tmp); } -int mm_gen_cs_or_MD(void *km, char **buf, int *max_len, const mm_idx_t *mi, const mm_reg1_t *r, const char *seq, int is_MD, int no_iden) +int mm_gen_cs_or_MD(void *km, char **buf, int *max_len, const mm_idx_t *mi, const mm_reg1_t *r, const char *seq, int is_MD, int no_iden, int is_qstrand) { mm_bseq1_t t; kstring_t str; str.s = *buf, str.l = 0, str.m = *max_len; t.l_seq = strlen(seq); t.seq = (char*)seq; - write_cs_or_MD(km, &str, mi, &t, r, no_iden, is_MD, 0); + write_cs_or_MD(km, &str, mi, &t, r, no_iden, is_MD, 0, is_qstrand); *max_len = str.m; *buf = str.s; return str.l; @@ -257,12 +263,12 @@ int mm_gen_cs_or_MD(void *km, char **buf, int *max_len, const mm_idx_t *mi, cons int mm_gen_cs(void *km, char **buf, int *max_len, const mm_idx_t *mi, const mm_reg1_t *r, const char *seq, int no_iden) { - return mm_gen_cs_or_MD(km, buf, max_len, mi, r, seq, 0, no_iden); + return mm_gen_cs_or_MD(km, buf, max_len, mi, r, seq, 0, no_iden, 0); } int mm_gen_MD(void *km, char **buf, int *max_len, const mm_idx_t *mi, const mm_reg1_t *r, const char *seq) { - return mm_gen_cs_or_MD(km, buf, max_len, mi, r, seq, 1, 0); + return mm_gen_cs_or_MD(km, buf, max_len, mi, r, seq, 1, 0, 0); } double mm_event_identity(const mm_reg1_t *r) @@ -316,7 +322,11 @@ void mm_write_paf3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_sprintf_lite(s, "%s\t%d\t%d\t%d\t%c\t", t->name, t->l_seq, r->qs, r->qe, "+-"[r->rev]); if (mi->seq[r->rid].name) mm_sprintf_lite(s, "%s", mi->seq[r->rid].name); else mm_sprintf_lite(s, "%d", r->rid); - mm_sprintf_lite(s, "\t%d\t%d\t%d", mi->seq[r->rid].len, r->rs, r->re); + mm_sprintf_lite(s, "\t%d", mi->seq[r->rid].len); + if (opt_flag & MM_F_QSTRAND) + mm_sprintf_lite(s, "\t%d\t%d", mi->seq[r->rid].len - r->re, mi->seq[r->rid].len - r->rs); + else + mm_sprintf_lite(s, "\t%d\t%d", r->rs, r->re); mm_sprintf_lite(s, "\t%d\t%d", r->mlen, r->blen); mm_sprintf_lite(s, "\t%d", r->mapq); write_tags(s, r); @@ -328,7 +338,7 @@ void mm_write_paf3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_sprintf_lite(s, "%d%c", r->p->cigar[k]>>4, "MIDNSHP=XB"[r->p->cigar[k]&0xf]); } if (r->p && (opt_flag & (MM_F_OUT_CS|MM_F_OUT_MD))) - write_cs_or_MD(km, s, mi, t, r, !(opt_flag&MM_F_OUT_CS_LONG), opt_flag&MM_F_OUT_MD, 1); + write_cs_or_MD(km, s, mi, t, r, !(opt_flag&MM_F_OUT_CS_LONG), opt_flag&MM_F_OUT_MD, 1, opt_flag&MM_F_QSTRAND); if ((opt_flag & MM_F_COPY_COMMENT) && t->comment) mm_sprintf_lite(s, "\t%s", t->comment); } @@ -535,7 +545,7 @@ void mm_write_sam3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int se } } if (r->p && (opt_flag & (MM_F_OUT_CS|MM_F_OUT_MD))) - write_cs_or_MD(km, s, mi, t, r, !(opt_flag&MM_F_OUT_CS_LONG), opt_flag&MM_F_OUT_MD, 1); + write_cs_or_MD(km, s, mi, t, r, !(opt_flag&MM_F_OUT_CS_LONG), opt_flag&MM_F_OUT_MD, 1, 0); if (cigar_in_tag) write_sam_cigar(s, flag, 1, t->l_seq, r, opt_flag); } diff --git a/hit.c b/hit.c index 506464b..cd30ada 100644 --- a/hit.c +++ b/hit.c @@ -20,14 +20,14 @@ static inline void mm_cal_fuzzy_len(mm_reg1_t *r, const mm128_t *a) } } -static inline void mm_reg_set_coor(mm_reg1_t *r, int32_t qlen, const mm128_t *a) +static inline void mm_reg_set_coor(mm_reg1_t *r, int32_t qlen, const mm128_t *a, int is_qstrand) { // NB: r->as and r->cnt MUST BE set correctly for this function to work int32_t k = r->as, q_span = (int32_t)(a[k].y>>32&0xff); r->rev = a[k].x>>63; r->rid = a[k].x<<1>>33; r->rs = (int32_t)a[k].x + 1 > q_span? (int32_t)a[k].x + 1 - q_span : 0; // NB: target span may be shorter, so this test is necessary r->re = (int32_t)a[k + r->cnt - 1].x + 1; - if (!r->rev) { + if (!r->rev || is_qstrand) { r->qs = (int32_t)a[k].y + 1 - q_span; r->qe = (int32_t)a[k + r->cnt - 1].y + 1; } else { @@ -49,7 +49,7 @@ static inline uint64_t hash64(uint64_t key) return key; } -mm_reg1_t *mm_gen_regs(void *km, uint32_t hash, int qlen, int n_u, uint64_t *u, mm128_t *a) // convert chains to hits +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) // convert chains to hits { mm128_t *z, tmp; mm_reg1_t *r; @@ -81,7 +81,7 @@ mm_reg1_t *mm_gen_regs(void *km, uint32_t hash, int qlen, int n_u, uint64_t *u, ri->cnt = (int32_t)z[i].y; ri->as = z[i].y >> 32; ri->div = -1.0f; - mm_reg_set_coor(ri, qlen, a); + mm_reg_set_coor(ri, qlen, a, is_qstrand); } kfree(km, z); return r; @@ -103,7 +103,7 @@ static inline int mm_alt_score(int score, float alt_diff_frac) return score > 0? score : 1; } -void mm_split_reg(mm_reg1_t *r, mm_reg1_t *r2, int n, int qlen, mm128_t *a) +void mm_split_reg(mm_reg1_t *r, mm_reg1_t *r2, int n, int qlen, mm128_t *a, int is_qstrand) { if (n <= 0 || n >= r->cnt) return; *r2 = *r; @@ -115,10 +115,10 @@ void mm_split_reg(mm_reg1_t *r, mm_reg1_t *r2, int n, int qlen, mm128_t *a) r2->score = (int32_t)(r->score * ((float)r2->cnt / r->cnt) + .499); r2->as = r->as + n; if (r->parent == r->id) r2->parent = MM_PARENT_TMP_PRI; - mm_reg_set_coor(r2, qlen, a); + mm_reg_set_coor(r2, qlen, a, is_qstrand); r->cnt -= r2->cnt; r->score -= r2->score; - mm_reg_set_coor(r, qlen, a); + mm_reg_set_coor(r, qlen, a, is_qstrand); r->split |= 1, r2->split |= 2; } @@ -350,7 +350,7 @@ void mm_join_long(void *km, const mm_mapopt_t *opt, int qlen, int *n_regs_, mm_r // all conditions satisfied; join a[r1->as].y |= MM_SEED_LONG_JOIN; r0->cnt += r1->cnt, r0->score += r1->score; - mm_reg_set_coor(r0, qlen, a); + mm_reg_set_coor(r0, qlen, a, opt->flag&MM_F_QSTRAND); r1->cnt = 0; r1->parent = r0->id; ++n_drop; @@ -416,7 +416,7 @@ mm_seg_t *mm_seg_gen(void *km, uint32_t hash, int n_segs, const int *qlens, int } } for (s = 0; s < n_segs; ++s) { - regs[s] = mm_gen_regs(km, hash, qlens[s], seg[s].n_u, seg[s].u, seg[s].a); + regs[s] = mm_gen_regs(km, hash, qlens[s], seg[s].n_u, seg[s].u, seg[s].a, 0); n_regs[s] = seg[s].n_u; for (i = 0; i < n_regs[s]; ++i) { regs[s][i].seg_split = 1; diff --git a/index.c b/index.c index fd93dd3..3838382 100644 --- a/index.c +++ b/index.c @@ -161,6 +161,28 @@ int mm_idx_getseq(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, ui return en - st; } +int mm_idx_getseq_rev(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq) +{ + uint64_t i, st1, en1; + const mm_idx_seq_t *s; + if (rid >= mi->n_seq || st >= mi->seq[rid].len) return -1; + s = &mi->seq[rid]; + if (en > s->len) en = s->len; + st1 = s->offset + (s->len - en); + en1 = s->offset + (s->len - st); + for (i = st1; i < en1; ++i) { + uint8_t c = mm_seq4_get(mi->S, i); + seq[en1 - i - 1] = c < 4? 3 - c : c; + } + return en - st; +} + +int mm_idx_getseq2(const mm_idx_t *mi, int is_rev, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq) +{ + if (is_rev) return mm_idx_getseq_rev(mi, rid, st, en, seq); + else return mm_idx_getseq(mi, rid, st, en, seq); +} + int32_t mm_idx_cal_max_occ(const mm_idx_t *mi, float f) { int i; diff --git a/main.c b/main.c index 6780f41..e84a7f6 100644 --- a/main.c +++ b/main.c @@ -7,7 +7,7 @@ #include "mmpriv.h" #include "ketopt.h" -#define MM_VERSION "2.17-r974-dirty" +#define MM_VERSION "2.17-r985-qs-dirty" #ifdef __linux__ #include @@ -70,6 +70,7 @@ static ko_longopt_t long_options[] = { { "chain-gap-scale",ko_required_argument, 343 }, { "alt", ko_required_argument, 344 }, { "alt-drop", ko_required_argument, 345 }, + { "qstrand", ko_no_argument, 346 }, { "help", ko_no_argument, 'h' }, { "max-intron-len", ko_required_argument, 'G' }, { "version", ko_no_argument, 'V' }, @@ -217,6 +218,7 @@ int main(int argc, char *argv[]) else if (c == 343) opt.chain_gap_scale = atof(o.arg); // --chain-gap-scale else if (c == 344) alt_list = o.arg; // --alt else if (c == 345) opt.alt_drop = atof(o.arg); // --alt-drop + else if (c == 346) opt.flag |= MM_F_QSTRAND | MM_F_NO_INV; // --qstrand else if (c == 314) { // --frag yes_or_no(&opt, MM_F_FRAG_MODE, o.longidx, o.arg, 1); } else if (c == 315) { // --secondary diff --git a/map.c b/map.c index 01d323c..1cc4a05 100644 --- a/map.c +++ b/map.c @@ -232,9 +232,13 @@ static mm128_t *collect_seed_hits(void *km, const mm_mapopt_t *opt, int max_occ, if ((r[k]&1) == (q->q_pos&1)) { // forward strand p->x = (r[k]&0xffffffff00000000ULL) | rpos; p->y = (uint64_t)q->q_span << 32 | q->q_pos >> 1; - } else { // reverse strand + } else if (!(opt->flag & MM_F_QSTRAND)) { // reverse strand and not in the query-strand mode p->x = 1ULL<<63 | (r[k]&0xffffffff00000000ULL) | rpos; p->y = (uint64_t)q->q_span << 32 | (qlen - ((q->q_pos>>1) + 1 - q->q_span) - 1); + } else { // reverse strand; query-strand + int32_t len = mi->seq[r[k]>>32].len; + p->x = 1ULL<<63 | (r[k]&0xffffffff00000000ULL) | (len - (rpos + 1 - q->q_span) - 1); // coordinate only accurate for non-HPC seeds + p->y = (uint64_t)q->q_span << 32 | q->q_pos >> 1; } p->y |= (uint64_t)q->seg_id << MM_SEED_SEG_SHIFT; if (q->is_tandem) p->y |= MM_SEED_TANDEM; @@ -341,7 +345,7 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char ** b->frag_gap = max_chain_gap_ref; b->rep_len = rep_len; - regs0 = mm_gen_regs(b->km, hash, qlen_sum, n_regs0, u, a); + regs0 = mm_gen_regs(b->km, hash, qlen_sum, n_regs0, u, a, opt->flag&MM_F_QSTRAND); if (mi->n_alt) { mm_mark_alt(mi, n_regs0, regs0); mm_hit_sort(b->km, &n_regs0, regs0, opt->alt_drop); // this step can be merged into mm_gen_regs(); will do if this shows up in profile diff --git a/minimap.h b/minimap.h index 251879b..150f52c 100644 --- a/minimap.h +++ b/minimap.h @@ -36,6 +36,8 @@ #define MM_F_NO_END_FLT 0x10000000 #define MM_F_HARD_MLEVEL 0x20000000 #define MM_F_SAM_HIT_ONLY 0x40000000 +#define MM_F_QSTRAND (0x80000000LL) +#define MM_F_NO_INV (0x100000000LL) #define MM_I_HPC 0x1 #define MM_I_NO_SEQ 0x2 diff --git a/mmpriv.h b/mmpriv.h index ffba6eb..3e03d57 100644 --- a/mmpriv.h +++ b/mmpriv.h @@ -62,12 +62,13 @@ void mm_write_sam3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int se void mm_idxopt_init(mm_idxopt_t *opt); 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); mm128_t *mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int max_iter, int min_cnt, int min_sc, float gap_scale, int is_cdna, int n_segs, int64_t n, mm128_t *a, int *n_u_, uint64_t **_u, void *km); mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, const char *qstr, int *n_regs_, mm_reg1_t *regs, mm128_t *a); -mm_reg1_t *mm_gen_regs(void *km, uint32_t hash, int qlen, int n_u, uint64_t *u, mm128_t *a); +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); void mm_mark_alt(const mm_idx_t *mi, int n, mm_reg1_t *r); -void mm_split_reg(mm_reg1_t *r, mm_reg1_t *r2, int n, int qlen, mm128_t *a); +void mm_split_reg(mm_reg1_t *r, mm_reg1_t *r2, int n, int qlen, mm128_t *a, int is_qstrand); void mm_sync_regs(void *km, int n_regs, mm_reg1_t *regs); int mm_squeeze_a(void *km, int n_regs, mm_reg1_t *regs, mm128_t *a); int mm_set_sam_pri(int n, mm_reg1_t *r); diff --git a/options.c b/options.c index 8276b24..81f6f4c 100644 --- a/options.c +++ b/options.c @@ -191,5 +191,10 @@ int mm_check_opt(const mm_idxopt_t *io, const mm_mapopt_t *mo) fprintf(stderr, "[ERROR]\033[1;31m -X/-P and --secondary=no can't be applied at the same time\033[0m\n"); return -5; } + if ((mo->flag & MM_F_QSTRAND) && ((mo->flag & (MM_F_OUT_SAM|MM_F_SPLICE|MM_F_FRAG_MODE)) || (io->flag & MM_I_HPC))) { + if (mm_verbose >= 1) + fprintf(stderr, "[ERROR]\033[1;31m --qstrand doesn't work with -a, -H, --frag or --splice\033[0m\n"); + return -5; + } return 0; }