Rework single-end MAPQ computation

First, there was a bug in MAPQ computation due to initializing
best_score to -1000. In the first iteration, min_mapq_diff is updated thus:

min_mapq_diff = std::max(0, alignment.score - best_score);

This is just alignment.score + 1000, so most of the time, min_mapq_diff
ends up being something like 1600 or so for typical read lengths. (This is
clamped to 60 in the output, but still.)

Here, we keep track of the highest and second-highest alignment score in
best_score and second_best_score. Then, mapq is computed such that it is
* 60 for second_best_score=0 and
* 0 for second_best_score=best_score,
where anything in between is interpolated.

See #25

src/aln.cpp

@@ -92,12 +92,11 @@ static inline void align_SE(
     Nam n_max = nams[0];
 
     int best_edit_distance = std::numeric_limits<int>::max();
-    int best_score = -1000;
+    int best_score = 0;
+    int second_best_score = 0;
 
     Alignment best_alignment;
-    best_alignment.score = -100000;
     best_alignment.is_unaligned = true;
-    int min_mapq_diff = best_edit_distance;
     for (auto &nam : nams) {
         float score_dropoff = (float) nam.n_hits / n_max.n_hits;
         if (tries >= max_tries || (max_tries > 1 && best_edit_distance == 0) || score_dropoff < dropoff_threshold) {
@@ -109,25 +108,25 @@ static inline void align_SE(
         details.tried_alignment++;
         details.gapped += alignment.gapped;
 
-        int diff_to_best = std::abs(best_score - alignment.score);
-        min_mapq_diff = std::min(min_mapq_diff, diff_to_best);
-
         if (max_secondary > 0) {
             alignments.emplace_back(alignment);
         }
         if (alignment.score > best_score) {
-            min_mapq_diff = std::max(0, alignment.score - best_score); // new distance to next best match
+            second_best_score = best_score;
             best_score = alignment.score;
             best_alignment = std::move(alignment);
             if (max_secondary == 0) {
                 best_edit_distance = best_alignment.global_ed;
             }
+        } else if (alignment.score > second_best_score) {
+            second_best_score = alignment.score;
         }
         tries++;
     }
+    uint8_t mapq = 60.0 * (best_score - second_best_score) / best_score;
 
     if (max_secondary == 0) {
-        best_alignment.mapq = std::min(min_mapq_diff, 60);
+        best_alignment.mapq = mapq;
         sam.add(best_alignment, record, read.rc, true, details);
         return;
     }
@@ -146,7 +145,7 @@ static inline void align_SE(
         }
         bool is_primary = i == 0;
         if (is_primary) {
-            alignment.mapq = std::min(min_mapq_diff, 60);
+            alignment.mapq = mapq;
         } else {
             alignment.mapq = 255;
         }