zzh
/
minimap2
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

r557: fixed another mapq underestimate 

When a chain is split during base-level alignment, its chaining score is
reduced. However, the chaining score of its suboptimal chain remains the same.
This leads to underestimated mapping quality.
This commit is contained in:
Heng Li 2017-11-07 23:20:49 -05:00
parent 65deedfa96
commit b24d68ae9f
4 changed files with 7 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
format.c
View File

@ -203,7 +203,6 @@ static void write_cs(void *km, kstring_t *s, const mm_idx_t *mi, const mm_bseq1_
static inline void write_tags(kstring_t *s, const mm_reg1_t *r)
{
int type = r->inv? 'I' : r->id == r->parent? 'P' : 'S';
if (r->iden_flt) mm_sprintf_lite(s, "\tom:i:%d", r->mapq);
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);
if (r->p->trans_strand == 1 || r->p->trans_strand == 2)
@ -336,9 +335,8 @@ void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int se
mm_sprintf_lite(s, "\t%s\t%d\t0\t*", mi->seq[this_rid].name, this_pos+1);
} else mm_sprintf_lite(s, "\t*\t0\t0\t*");
} else {
int mapq = !r->iden_flt? r->mapq : r->mapq < 3? r->mapq : 3;
this_rid = r->rid, this_pos = r->rs, this_rev = r->rev;
mm_sprintf_lite(s, "\t%s\t%d\t%d\t", mi->seq[r->rid].name, r->rs+1, mapq);
mm_sprintf_lite(s, "\t%s\t%d\t%d\t", mi->seq[r->rid].name, r->rs+1, r->mapq);
if ((opt_flag & MM_F_LONG_CIGAR) && r->p && r->p->n_cigar > max_bam_cigar_op - 2) {
int n_cigar = r->p->n_cigar;
if (r->qs != 0) ++n_cigar;

6
hit.c
View File

@ -76,7 +76,7 @@ mm_reg1_t *mm_gen_regs(void *km, uint32_t hash, int qlen, int n_u, uint64_t *u,
mm_reg1_t *ri = &r[i];
ri->id = i;
ri->parent = MM_PARENT_UNSET;
ri->score = z[i].x >> 32;
ri->score = ri->score0 = z[i].x >> 32;
ri->hash = (uint32_t)z[i].x;
ri->cnt = (int32_t)z[i].y;
ri->as = z[i].y >> 32;
@ -425,14 +425,14 @@ void mm_set_mapq(int n_regs, mm_reg1_t *regs, int min_chain_sc, int match_sc, in
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->score;
float 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
}
} else {
float x = (float)subsc / r->score;
float x = (float)subsc / r->score0;
if (r->p) {
float identity = (float)r->mlen / r->blen;
mapq = (int)(identity * pen_cm * q_coef * (1.0f - x) * logf((float)r->p->dp_max / match_sc));

2
main.c
View File

@ -6,7 +6,7 @@
#include "mmpriv.h"
#include "getopt.h"
#define MM_VERSION "2.4-r556-dirty"
#define MM_VERSION "2.4-r557-dirty"
#ifdef __linux__
#include <sys/resource.h>

4
minimap.h
View File

@ -64,7 +64,7 @@ typedef struct {
typedef struct {
int32_t id; // ID for internal uses (see also parent below)
uint32_t cnt:30, rev:1, seg_split:1; // number of minimizers; if on the reverse strand
uint32_t cnt:28, rev:1, seg_split:1, sam_pri:1, proper_frag:1; // number of minimizers; if on the reverse strand
uint32_t rid:31, inv:1; // reference index; if this is an alignment from inversion rescue
int32_t score; // DP alignment score
int32_t qs, qe, rs, re; // query start and end; reference start and end
@ -72,7 +72,7 @@ typedef struct {
int32_t as; // offset in the a[] array (for internal uses only)
int32_t mlen, blen; // seeded exact match length; seeded alignment block length
uint32_t mapq:8, split:2, n_sub:22; // mapQ; split pattern; number of suboptimal mappings
uint32_t sam_pri:1, proper_frag:1, iden_flt:1, pe_thru:1, dummy:28;
uint32_t pe_thru:1, score0:31;
uint32_t hash;
mm_extra_t *p;
} mm_reg1_t;