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

r711: assign proper mapq to primary inversions

This commit is contained in:
Heng Li 2018-02-15 14:34:59 -05:00
parent a0d62519c1
commit 8fc5f8dc90
4 changed files with 32 additions and 5 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

29
hit.c
View File

@ -408,7 +408,33 @@ void mm_seg_free(void *km, int n_segs, mm_seg_t *segs)
kfree(km, segs);
}
void mm_set_mapq(int n_regs, mm_reg1_t *regs, int min_chain_sc, int match_sc, int rep_len, int is_sr)
static void mm_set_inv_mapq(void *km, int n_regs, mm_reg1_t *regs)
{
int i, n_aux;
uint64_t *aux;
if (n_regs < 3) return;
for (i = 0; i < n_regs; ++i)
if (regs[i].inv) break;
if (i == n_regs) return; // no inversion hits
aux = (uint64_t*)kmalloc(km, n_regs * 8);
for (i = n_aux = 0; i < n_regs; ++i)
if (regs[i].parent == i || regs[i].parent < 0)
aux[n_aux++] = (uint64_t)regs[i].as << 32 | i;
radix_sort_64(aux, aux + n_aux);
for (i = 1; i < n_aux - 1; ++i) {
mm_reg1_t *inv = &regs[(int32_t)aux[i]];
if (inv->inv) {
mm_reg1_t *l = &regs[(int32_t)aux[i-1]];
mm_reg1_t *r = &regs[(int32_t)aux[i+1]];
inv->mapq = l->mapq < r->mapq? l->mapq : r->mapq;
}
}
kfree(km, aux);
}
void mm_set_mapq(void *km, int n_regs, mm_reg1_t *regs, int min_chain_sc, int match_sc, int rep_len, int is_sr)
{
static const float q_coef = 40.0f;
int64_t sum_sc = 0;
@ -451,4 +477,5 @@ void mm_set_mapq(int n_regs, mm_reg1_t *regs, int min_chain_sc, int match_sc, in
if (r->p && r->p->dp_max > r->p->dp_max2 && r->mapq == 0) r->mapq = 1;
} else r->mapq = 0;
}
mm_set_inv_mapq(km, n_regs, regs);
}

2
main.c
View File

@ -6,7 +6,7 @@
#include "mmpriv.h"
#include "getopt.h"
#define MM_VERSION "2.8-r710-dirty"
#define MM_VERSION "2.8-r711-dirty"
#ifdef __linux__
#include <sys/resource.h>

4
map.c
View File

@ -340,7 +340,7 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
if (n_segs == 1) { // uni-segment
regs0 = align_regs(opt, mi, b->km, qlens[0], seqs[0], &n_regs0, regs0, a);
mm_set_mapq(n_regs0, regs0, opt->min_chain_score, opt->a, rep_len, is_sr);
mm_set_mapq(b->km, n_regs0, regs0, opt->min_chain_score, opt->a, rep_len, is_sr);
n_regs[0] = n_regs0, regs[0] = regs0;
} else { // multi-segment
mm_seg_t *seg;
@ -349,7 +349,7 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
for (i = 0; i < n_segs; ++i) {
mm_set_parent(b->km, opt->mask_level, n_regs[i], regs[i], opt->a * 2 + opt->b); // update mm_reg1_t::parent
regs[i] = align_regs(opt, mi, b->km, qlens[i], seqs[i], &n_regs[i], regs[i], seg[i].a);
mm_set_mapq(n_regs[i], regs[i], opt->min_chain_score, opt->a, rep_len, is_sr);
mm_set_mapq(b->km, n_regs[i], regs[i], opt->min_chain_score, opt->a, rep_len, is_sr);
}
mm_seg_free(b->km, n_segs, seg);
if (n_segs == 2 && opt->pe_ori >= 0 && (opt->flag&MM_F_CIGAR))

2
mmpriv.h
View File

@ -78,7 +78,7 @@ void mm_select_sub_multi(void *km, float pri_ratio, float pri1, float pri2, int
void mm_filter_regs(void *km, const mm_mapopt_t *opt, int qlen, 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, int min_chain_sc, int match_sc, int rep_len, int is_sr);
void mm_set_mapq(void *km, int n_regs, mm_reg1_t *regs, int min_chain_sc, int match_sc, int rep_len, int is_sr);
void mm_est_err(const mm_idx_t *mi, int qlen, int n_regs, mm_reg1_t *regs, const mm128_t *a, int32_t n, const uint64_t *mini_pos);