From e9dc1ce2b659377f390cb93f6aa208ffcb18f034 Mon Sep 17 00:00:00 2001
From: Heng Li <lh3@me.com>
Date: Wed, 26 Jul 2017 11:34:14 -0400
Subject: [PATCH] r205: when computing mapq, consider min_chain_sc

Not doing this was a mistake.
---
 hit.c    | 9 +++++----
 main.c   | 2 +-
 map.c    | 2 +-
 mmpriv.h | 2 +-
 4 files changed, 8 insertions(+), 7 deletions(-)

diff --git a/hit.c b/hit.c
index dddba12..477dedb 100644
--- a/hit.c
+++ b/hit.c
@@ -288,17 +288,18 @@ void mm_join_long(void *km, const mm_mapopt_t *opt, int qlen, int *n_regs_, mm_r
 	}
 }
 
-void mm_set_mapq(int n_regs, mm_reg1_t *regs)
+void mm_set_mapq(int n_regs, mm_reg1_t *regs, int min_chain_sc)
 {
 	int i;
 	for (i = 0; i < n_regs; ++i) {
 		mm_reg1_t *r = &regs[i];
 		if (r->parent == r->id) {
-			int mapq;
+			int mapq, subsc;
+			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->p->blen - r->p->n_diff - r->p->n_ambi) / (r->p->blen - r->p->n_ambi);
-				mapq = (int)(identity * 30.0 * (1. - (float)r->p->dp_max2 * r->subsc / r->p->dp_max / r->score) * logf(r->score));
-			} else mapq = (int)(30.0 * (1. - (float)r->subsc / r->score) * logf(r->score));
+				mapq = (int)(identity * 30.0 * (1. - (float)r->p->dp_max2 * subsc / r->p->dp_max / r->score) * logf(r->score));
+			} else mapq = (int)(30.0 * (1. - (float)subsc / r->score) * logf(r->score));
 			mapq = mapq > 0? mapq : 0;
 			r->mapq = mapq < 60? mapq : 60;
 		} else r->mapq = 0;
diff --git a/main.c b/main.c
index 2ee1a36..175d085 100644
--- a/main.c
+++ b/main.c
@@ -8,7 +8,7 @@
 #include "minimap.h"
 #include "mmpriv.h"
 
-#define MM_VERSION "2.0-r203-dirty"
+#define MM_VERSION "2.0-r205-dirty"
 
 void liftrlimit()
 {
diff --git a/map.c b/map.c
index cad4dd6..94320ab 100644
--- a/map.c
+++ b/map.c
@@ -255,7 +255,7 @@ mm_reg1_t *mm_map_frag(const mm_mapopt_t *opt, const mm_idx_t *mi, mm_tbuf_t *b,
 			mm_set_sam_pri(*n_regs, regs);
 		}
 	}
-	mm_set_mapq(*n_regs, regs);
+	mm_set_mapq(*n_regs, regs, opt->min_chain_score);
 
 	// free
 	kfree(b->km, a);
diff --git a/mmpriv.h b/mmpriv.h
index 35ffa1c..6135172 100644
--- a/mmpriv.h
+++ b/mmpriv.h
@@ -49,7 +49,7 @@ void mm_select_sub(void *km, float mask_level, float pri_ratio, int min_diff, in
 void mm_filter_regs(void *km, const mm_mapopt_t *opt, int *n_regs, mm_reg1_t *regs);
 void mm_join_long(void *km, const mm_mapopt_t *opt, int qlen, int *n_regs, mm_reg1_t *regs, mm128_t *a);
 void mm_hit_sort_by_dp(void *km, int *n_regs, mm_reg1_t *r);
-void mm_set_mapq(int n_regs, mm_reg1_t *regs);
+void mm_set_mapq(int n_regs, mm_reg1_t *regs, int min_chain_sc);
 
 #ifdef __cplusplus
 }