From fcd4df2a734cbb89bb1936dbded8ce89aa48d16c Mon Sep 17 00:00:00 2001 From: Heng Li Date: Mon, 11 Mar 2024 17:19:13 -0400 Subject: [PATCH] r1190: output unadjusted dp_max to ms:i 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 --- align.c | 2 +- format.c | 2 +- minimap.h | 3 ++- 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/align.c b/align.c index 5a0af66..a1d90ea 100644 --- a/align.c +++ b/align.c @@ -295,7 +295,7 @@ static void mm_update_extra(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *ts 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); 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 } diff --git a/format.c b/format.c index 46521d4..d00a96d 100644 --- a/format.c +++ b/format.c @@ -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'; else type = r->inv? 'i' : 'S'; 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) mm_sprintf_lite(s, "\tts:A:%c", "?+-?"[r->p->trans_strand]); } diff --git a/minimap.h b/minimap.h index f186ddf..42f4867 100644 --- a/minimap.h +++ b/minimap.h @@ -5,7 +5,7 @@ #include #include -#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_DUAL (0x002LL) // skip pairs where query name is lexicographically larger than target name @@ -98,6 +98,7 @@ typedef struct { typedef struct { 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_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_cigar; // number of cigar operations in cigar[] uint32_t cigar[];