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

r1190: output unadjusted dp_max to ms:i

Browse Source 
This was an oversight affecting v2.22+. The latest minimap2 ranks hits and
estimates mapping quality with an adjusted alignment score (see the minimap2
update paper). This score however is not calculated when there is only one hit.
As a result, the ms:i tag varies depends on other sequences in the reference
genome, which is confusing. This change lets minimap2 to output the unadjusted
score at ms:i. At present, the adjusted score is not outputted.

Resolves #1146
This commit is contained in:
Heng Li
2024-03-11 17:19:13 -04:00
parent 0efc886012
commit fcd4df2a73
3 changed files with 4 additions and 3 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
align.c
View File

@@ -295,7 +295,7 @@ static void mm_update_extra(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *ts
toff += len; toff += len;
} }
} }
p->dp_max = (int32_t)(max + .499); p->dp_max = p->dp_max0 = (int32_t)(max + .499);
assert(qoff == r->qe - r->qs && toff == r->re - r->rs); assert(qoff == r->qe - r->qs && toff == r->re - r->rs);
if (is_eqx) mm_update_cigar_eqx(r, qseq, tseq); // NB: it has to be called here as changes to qseq and tseq are not returned if (is_eqx) mm_update_cigar_eqx(r, qseq, tseq); // NB: it has to be called here as changes to qseq and tseq are not returned
} }

2
format.c
View File

@@ -344,7 +344,7 @@ static inline void write_tags(kstring_t *s, const mm_reg1_t *r)
if (r->id == r->parent) type = r->inv? 'I' : 'P'; if (r->id == r->parent) type = r->inv? 'I' : 'P';
else type = r->inv? 'i' : 'S'; else type = r->inv? 'i' : 'S';
if (r->p) { if (r->p) {
mm_sprintf_lite(s, "\tNM:i:%d\tms:i:%d\tAS:i:%d\tnn:i:%d", r->blen - r->mlen + r->p->n_ambi, r->p->dp_max, r->p->dp_score, r->p->n_ambi); mm_sprintf_lite(s, "\tNM:i:%d\tms:i:%d\tAS:i:%d\tnn:i:%d", r->blen - r->mlen + r->p->n_ambi, r->p->dp_max0, r->p->dp_score, r->p->n_ambi);
if (r->p->trans_strand == 1 || r->p->trans_strand == 2) if (r->p->trans_strand == 1 || r->p->trans_strand == 2)
mm_sprintf_lite(s, "\tts:A:%c", "?+-?"[r->p->trans_strand]); mm_sprintf_lite(s, "\tts:A:%c", "?+-?"[r->p->trans_strand]);
} }

3
minimap.h
View File

@@ -5,7 +5,7 @@
#include <stdio.h> #include <stdio.h>
#include <sys/types.h> #include <sys/types.h>
#define MM_VERSION "2.26-r1188-dirty" #define MM_VERSION "2.26-r1190-dirty"
#define MM_F_NO_DIAG (0x001LL) // no exact diagonal hit #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 #define MM_F_NO_DUAL (0x002LL) // skip pairs where query name is lexicographically larger than target name
@@ -98,6 +98,7 @@ typedef struct {
typedef struct { typedef struct {
uint32_t capacity; // the capacity of cigar[] uint32_t capacity; // the capacity of cigar[]
int32_t dp_score, dp_max, dp_max2; // DP score; score of the max-scoring segment; score of the best alternate mappings int32_t dp_score, dp_max, dp_max2; // DP score; score of the max-scoring segment; score of the best alternate mappings
int32_t dp_max0; // DP score before mm_update_dp_max() adjustment
uint32_t n_ambi:30, trans_strand:2; // number of ambiguous bases; transcript strand: 0 for unknown, 1 for +, 2 for - uint32_t n_ambi:30, trans_strand:2; // number of ambiguous bases; transcript strand: 0 for unknown, 1 for +, 2 for -
uint32_t n_cigar; // number of cigar operations in cigar[] uint32_t n_cigar; // number of cigar operations in cigar[]
uint32_t cigar[]; uint32_t cigar[];