fork/minimap2
1
0
Fork 0
mirror of https://github.com/lh3/minimap2.git synced 2026-09-15 13:07:55 +08:00
Code Issues Packages Projects Releases Wiki Activity

r1265: prefer spliced alignment

Browse Source 
mapq needs to be elevated
This commit is contained in:
Heng Li
2025-04-07 00:31:04 -04:00
parent fbb9c0fcba
commit 75619c7b51
4 changed files with 17 additions and 11 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5
README.md
View File

@@ -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**.

15
align.c
View File

@@ -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 = &regs[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, &regs[i], &r2, n_a, a, &ez, opt->flag);
if (opt->flag&MM_F_SPLICE)

4
jump.c
View File

@@ -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;

4
minimap.h
View File

@@ -5,7 +5,7 @@
#include <stdio.h>
#include <sys/types.h>
#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;