From 75619c7b511b2f47799f291ee934d8a9269ba4e0 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Mon, 7 Apr 2025 00:31:04 -0400 Subject: [PATCH] r1265: prefer spliced alignment mapq needs to be elevated --- README.md | 5 ++--- align.c | 15 +++++++++------ jump.c | 4 ++++ minimap.h | 4 ++-- 4 files changed, 17 insertions(+), 11 deletions(-) diff --git a/README.md b/README.md index c5e7c62..600ae8c 100644 --- a/README.md +++ b/README.md @@ -174,9 +174,8 @@ or the last exons. Minimap2 rates an alignment by the score of the max-scoring sub-segment, *excluding* introns, and marks the best alignment as primary in SAM. When a -spliced gene also has unspliced pseudogenes, minimap2 does not intentionally -prefer spliced alignment, though in practice it more often marks the spliced -alignment as the primary. By default, minimap2 outputs up to five secondary +spliced gene also has unspliced pseudogenes, minimap2 slightly prefers +the spliced alignment. By default, minimap2 outputs up to five secondary alignments (i.e. likely pseudogenes in the context of RNA-seq mapping). This can be tuned with option **-N**. diff --git a/align.c b/align.c index 97c29f6..56345a5 100644 --- a/align.c +++ b/align.c @@ -258,7 +258,7 @@ static void mm_update_extra(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *ts if (p == 0) return; mm_fix_cigar(r, qseq, tseq, &qshift, &tshift); qseq += qshift, tseq += tshift; // qseq and tseq may be shifted due to the removal of leading I/D - r->blen = r->mlen = 0; + r->blen = r->mlen = 0, r->is_spliced = 0; for (k = 0; k < p->n_cigar; ++k) { uint32_t op = p->cigar[k]&0xf, len = p->cigar[k]>>4; if (op == MM_CIGAR_MATCH) { @@ -292,7 +292,7 @@ static void mm_update_extra(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *ts if (s < 0) s = 0; toff += len; } else if (op == MM_CIGAR_N_SKIP) { - toff += len; + r->is_spliced = 1, toff += len; } } p->dp_max = p->dp_max0 = (int32_t)(max + .499); @@ -852,7 +852,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int 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, !(opt->flag & MM_F_SR)); + mm_update_extra(r, qseq, tseq, mat, opt->q, opt->e, opt->flag & MM_F_EQX, !(is_sr || is_sr_rna)); if (rev && r->p->trans_strand) r->p->trans_strand ^= 3; // flip to the read strand } @@ -911,7 +911,7 @@ static int mm_align1_inv(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, i } r_inv->rs = r1->re + t_off; r_inv->re = r_inv->rs + ez->max_t + 1; - mm_update_extra(r_inv, &qseq[q_off], &tseq[t_off], mat, opt->q, opt->e, opt->flag & MM_F_EQX, !(opt->flag & MM_F_SR)); + mm_update_extra(r_inv, &qseq[q_off], &tseq[t_off], mat, opt->q, opt->e, opt->flag & MM_F_EQX, !(opt->flag & (MM_F_SR|MM_F_SR_RNA))); ret = 1; end_align1_inv: kfree(km, tseq); @@ -1014,7 +1014,7 @@ mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *m for (i = 0; i < n_regs; ++i) { mm_reg1_t r2; if ((opt->flag&MM_F_SPLICE) && (opt->flag&MM_F_SPLICE_FOR) && (opt->flag&MM_F_SPLICE_REV)) { // then do two rounds of alignments for both strands - mm_reg1_t s[2], s2[2]; + mm_reg1_t s[2], s2[2], *r; int which, trans_strand; s[0] = s[1] = regs[i]; mm_align1(km, opt, mi, qlen, qseq0, &s[0], &s2[0], n_a, a, &ez, MM_F_SPLICE_FOR); @@ -1029,7 +1029,10 @@ mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *m regs[i] = s[1], r2 = s2[1]; free(s[0].p); } - regs[i].p->trans_strand = trans_strand; + r = ®s[i]; + r->p->trans_strand = trans_strand; + if (r->is_spliced && (trans_strand == 1 || trans_strand == 2)) + r->p->dp_max += (opt->a + opt->b) + ((opt->a + opt->b) >> 1); } else { // one round of alignment mm_align1(km, opt, mi, qlen, qseq0, ®s[i], &r2, n_a, a, &ez, opt->flag); if (opt->flag&MM_F_SPLICE) diff --git a/jump.c b/jump.c index c169cfb..5a953e4 100644 --- a/jump.c +++ b/jump.c @@ -101,6 +101,8 @@ static void mm_jump_split_left(void *km, const mm_idx_t *mi, const mm_mapopt_t * else r->qe = qlen; r->blen += clip, r->mlen += clip - mm0; r->p->dp_max0 += (clip - mm0) * opt->a - mm0 * opt->b; + r->p->dp_max += (clip - mm0) * opt->a - mm0 * opt->b; + if (!r->is_spliced) r->is_spliced = 1, r->p->dp_max += (opt->a + opt->b) + ((opt->a + opt->b) >> 1); } else if (m > 0 && a[i0].off > r->rs) { // trim by l; l is always positive r->p->cigar[0] -= l << 4 | MM_CIGAR_MATCH; r->rs += l; @@ -164,6 +166,8 @@ static void mm_jump_split_right(void *km, const mm_idx_t *mi, const mm_mapopt_t else r->qs = 0; r->blen += clip, r->mlen += clip - mm0; r->p->dp_max0 += (clip - mm0) * opt->a - mm0 * opt->b; + r->p->dp_max += (clip - mm0) * opt->a - mm0 * opt->b; + if (!r->is_spliced) r->is_spliced = 1, r->p->dp_max += (opt->a + opt->b) + ((opt->a + opt->b) >> 1); } else if (m > 0 && r->re > a[i0].off) { // trim by l; l is always positive r->p->cigar[r->p->n_cigar - 1] -= l << 4 | MM_CIGAR_MATCH; r->re -= l; diff --git a/minimap.h b/minimap.h index 3116f9b..dd86fd3 100644 --- a/minimap.h +++ b/minimap.h @@ -5,7 +5,7 @@ #include #include -#define MM_VERSION "2.28-r1264-dirty" +#define MM_VERSION "2.28-r1265-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 @@ -119,7 +119,7 @@ typedef struct { int32_t mlen, blen; // seeded exact match length; seeded alignment block length int32_t n_sub; // number of suboptimal mappings int32_t score0; // initial chaining score (before chain merging/spliting) - uint32_t mapq:8, split:2, rev:1, inv:1, sam_pri:1, proper_frag:1, pe_thru:1, seg_split:1, seg_id:8, split_inv:1, is_alt:1, strand_retained:1, dummy:5; + uint32_t mapq:8, split:2, rev:1, inv:1, sam_pri:1, proper_frag:1, pe_thru:1, seg_split:1, seg_id:8, split_inv:1, is_alt:1, strand_retained:1, is_spliced:1, is_spliced2:1, dummy:3; uint32_t hash; float div; mm_extra_t *p;