From 19c39e704f3aa6191fb93c563f403ba3771e8293 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Wed, 4 Oct 2017 13:37:40 -0400 Subject: [PATCH] r464: fixed a bug in pairing, due to randomization --- main.c | 2 +- pe.c | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/main.c b/main.c index 56d2afc..84054c3 100644 --- a/main.c +++ b/main.c @@ -6,7 +6,7 @@ #include "mmpriv.h" #include "getopt.h" -#define MM_VERSION "2.2-r463-dirty" +#define MM_VERSION "2.2-r464-dirty" #ifdef __linux__ #include diff --git a/pe.c b/pe.c index c3174bc..ffef9ec 100644 --- a/pe.c +++ b/pe.c @@ -99,7 +99,7 @@ void mm_pair(void *km, int max_gap_ref, int pe_bonus, int sub_diff, int match_sc q = a[j].r; if (r->rid != q->rid || r->rs - q->re > max_gap_ref) break; if (r->p->dp_max + q->p->dp_max < dp_thres) continue; - score = (int64_t)(r->p->dp_max + q->p->dp_max) | (r->hash + q->hash); + score = (int64_t)(r->p->dp_max + q->p->dp_max) << 32 | (r->hash + q->hash); if (score > max) max = score, max_idx[a[j].s] = j, max_idx[a[i].s] = i; kv_push(uint64_t, km, sc, score); @@ -126,11 +126,11 @@ void mm_pair(void *km, int max_gap_ref, int pe_bonus, int sub_diff, int match_sc } mapq_pe = r[0]->mapq > r[1]->mapq? r[0]->mapq : r[1]->mapq; for (i = 0; i < sc.n; ++i) - if (sc.a[i] + sub_diff >= max) + if ((sc.a[i]>>32) + sub_diff >= max>>32) ++n_sub; if (sc.n > 1) { int mapq_pe_alt; - mapq_pe_alt = (int)(6.02f * (max - sc.a[sc.n - 2]) / match_sc - 4.343f * logf(n_sub)); // n_sub > 0 because it counts the optimal, too + mapq_pe_alt = (int)(6.02f * ((max>>32) - (sc.a[sc.n - 2]>>32)) / match_sc - 4.343f * logf(n_sub)); // n_sub > 0 because it counts the optimal, too mapq_pe = mapq_pe < mapq_pe_alt? mapq_pe : mapq_pe_alt; } if (r[0]->mapq < mapq_pe) r[0]->mapq = (r[0]->mapq + mapq_pe) / 2;