diff --git a/hit.c b/hit.c index f5a50a2..1bd9683 100644 --- a/hit.c +++ b/hit.c @@ -418,16 +418,19 @@ static void mm_set_inv_mapq(void *km, int n_regs, mm_reg1_t *regs) kfree(km, aux); } -void mm_set_mapq(void *km, int n_regs, mm_reg1_t *regs, int min_chain_sc, int match_sc, int rep_len, int is_sr) +void mm_set_mapq2(void *km, int n_regs, mm_reg1_t *regs, int min_chain_sc, int match_sc, int rep_len, int is_sr, int is_splice) { static const float q_coef = 40.0f; int64_t sum_sc = 0; float uniq_ratio; - int i; + int i, n_2nd_splice = 0; if (n_regs == 0) return; - for (i = 0; i < n_regs; ++i) + for (i = 0; i < n_regs; ++i) { if (regs[i].parent == regs[i].id) sum_sc += regs[i].score; + else if (regs[i].is_spliced) + ++n_2nd_splice; + } uniq_ratio = (float)sum_sc / (sum_sc + rep_len); for (i = 0; i < n_regs; ++i) { mm_reg1_t *r = ®s[i]; @@ -440,13 +443,18 @@ void mm_set_mapq(void *km, int n_regs, mm_reg1_t *regs, int min_chain_sc, int ma pen_cm = pen_s1 < pen_cm? pen_s1 : pen_cm; subsc = r->subsc > min_chain_sc? r->subsc : min_chain_sc; if (r->p && r->p->dp_max2 > 0 && r->p->dp_max > 0) { - float identity = (float)r->mlen / r->blen; - float x = (float)r->p->dp_max2 * subsc / r->p->dp_max / r->score0; + float x, identity = (float)r->mlen / r->blen; + if (is_sr && is_splice) + x = (float)r->p->dp_max2 / r->p->dp_max; // ignore chaining score; for short RNA-seq reads, unspliced chaining score tends to be higher + else + x = (float)r->p->dp_max2 * subsc / r->p->dp_max / r->score0; mapq = (int)(identity * pen_cm * q_coef * (1.0f - x * x) * logf((float)r->p->dp_max / match_sc)); if (!is_sr) { int mapq_alt = (int)(6.02f * identity * identity * (r->p->dp_max - r->p->dp_max2) / match_sc + .499f); // BWA-MEM like mapQ, mostly for short reads mapq = mapq < mapq_alt? mapq : mapq_alt; // in case the long-read heuristic fails } + if (is_splice && is_sr && r->is_spliced && n_2nd_splice == 0) + mapq += 10; } else { float x = (float)subsc / r->score0; if (r->p) { diff --git a/map.c b/map.c index 88484b0..a62eab8 100644 --- a/map.c +++ b/map.c @@ -227,7 +227,7 @@ static mm_reg1_t *align_regs(const mm_mapopt_t *opt, const mm_idx_t *mi, void *k void mm_map_frag_core(const mm_idx_t *mi, int n_segs, const int *qlens, const char **seqs, int *n_regs, mm_reg1_t **regs, mm_tbuf_t *b, const mm_mapopt_t *opt, const char *qname) { int i, j, rep_len, qlen_sum, n_regs0, n_mini_pos; - int max_chain_gap_qry, max_chain_gap_ref, is_splice = !!(opt->flag & MM_F_SPLICE), is_sr = !!(opt->flag & MM_F_SR); + int max_chain_gap_qry, max_chain_gap_ref, is_splice = !!(opt->flag & MM_F_SPLICE), is_sr = !!(opt->flag & MM_F_SR), is_sr_rna = !!(opt->flag & MM_F_SR_RNA); uint32_t hash; int64_t n_a; uint64_t *u, *mini_pos; @@ -338,7 +338,7 @@ void mm_map_frag_core(const mm_idx_t *mi, int n_segs, const int *qlens, const ch if (n_segs == 1) { // uni-segment regs0 = align_regs(opt, mi, b->km, qlens[0], seqs[0], &n_regs0, regs0, a); regs0 = (mm_reg1_t*)realloc(regs0, sizeof(*regs0) * n_regs0); - mm_set_mapq(b->km, n_regs0, regs0, opt->min_chain_score, opt->a, rep_len, is_sr); + mm_set_mapq2(b->km, n_regs0, regs0, opt->min_chain_score, opt->a, rep_len, is_sr || is_sr_rna, is_splice); n_regs[0] = n_regs0, regs[0] = regs0; } else { // multi-segment mm_seg_t *seg; @@ -347,7 +347,7 @@ void mm_map_frag_core(const mm_idx_t *mi, int n_segs, const int *qlens, const ch for (i = 0; i < n_segs; ++i) { mm_set_parent(b->km, opt->mask_level, opt->mask_len, n_regs[i], regs[i], opt->a * 2 + opt->b, opt->flag&MM_F_HARD_MLEVEL, opt->alt_drop); // update mm_reg1_t::parent regs[i] = align_regs(opt, mi, b->km, qlens[i], seqs[i], &n_regs[i], regs[i], seg[i].a); - mm_set_mapq(b->km, n_regs[i], regs[i], opt->min_chain_score, opt->a, rep_len, is_sr); + mm_set_mapq2(b->km, n_regs[i], regs[i], opt->min_chain_score, opt->a, rep_len, is_sr || is_sr_rna, is_splice); } mm_seg_free(b->km, n_segs, seg); if (n_segs == 2 && opt->pe_ori >= 0 && (opt->flag&MM_F_CIGAR)) @@ -525,7 +525,7 @@ static void merge_hits(step_t *s) mm_select_sub(km, opt->pri_ratio, s->p->mi->k*2, opt->best_n, 0, opt->max_gap * 0.8, &s->n_reg[k], s->reg[k]); mm_set_sam_pri(s->n_reg[k], s->reg[k]); } - mm_set_mapq(km, s->n_reg[k], s->reg[k], opt->min_chain_score, opt->a, rep_len, !!(opt->flag & MM_F_SR)); + mm_set_mapq2(km, s->n_reg[k], s->reg[k], opt->min_chain_score, opt->a, rep_len, !!(opt->flag & (MM_F_SR|MM_F_SR_RNA)), !!(opt->flag & MM_F_SPLICE)); } if (s->n_seg[f] == 2 && opt->pe_ori >= 0 && (opt->flag&MM_F_CIGAR)) mm_pair(km, frag_gap_part[0], opt->pe_bonus, opt->a * 2 + opt->b, opt->a, qlens, &s->n_reg[k0], &s->reg[k0]); diff --git a/minimap.h b/minimap.h index dd86fd3..e39b4c1 100644 --- a/minimap.h +++ b/minimap.h @@ -5,7 +5,7 @@ #include #include -#define MM_VERSION "2.28-r1265-dirty" +#define MM_VERSION "2.28-r1266-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 diff --git a/mmpriv.h b/mmpriv.h index eec4882..c02a7f3 100644 --- a/mmpriv.h +++ b/mmpriv.h @@ -103,7 +103,7 @@ void mm_select_sub_multi(void *km, float pri_ratio, float pri1, float pri2, int int mm_filter_strand_retained(int n_regs, mm_reg1_t *r); void mm_filter_regs(const mm_mapopt_t *opt, int qlen, int *n_regs, mm_reg1_t *regs); void mm_hit_sort(void *km, int *n_regs, mm_reg1_t *r, float alt_diff_frac); -void mm_set_mapq(void *km, int n_regs, mm_reg1_t *regs, int min_chain_sc, int match_sc, int rep_len, int is_sr); +void mm_set_mapq2(void *km, int n_regs, mm_reg1_t *regs, int min_chain_sc, int match_sc, int rep_len, int is_sr, int is_splice); void mm_update_dp_max(int qlen, int n_regs, mm_reg1_t *regs, float frac, int a, int b); 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);