zzh
/
fast-bwa
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

mark primary, instead of dropping secondary

This commit is contained in:
Heng Li 2013-02-12 15:34:44 -05:00
parent 2fc469d0c9
commit 22b79b3475
2 changed files with 16 additions and 10 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

24
bwamem.c
View File

@ -345,14 +345,18 @@ int mem_sort_and_dedup(int n, mem_alnreg_t *a)
return m;
}
int mem_choose_alnreg_se(const mem_opt_t *opt, int n, mem_alnreg_t *a) // IMPORTANT: must run mem_sort_and_dedup() before calling this function
void mem_mark_primary_se(const mem_opt_t *opt, int n, mem_alnreg_t *a) // IMPORTANT: must run mem_sort_and_dedup() before calling this function
{ // similar to the loop in mem_chain_flt()
int i, j, m, tmp;
if (n <= 1) return n;
for (i = 0; i < n; ++i) a[i].sub = 0;
int i, k, tmp;
kvec_t(int) z;
if (n == 0) return;
kv_init(z);
for (i = 0; i < n; ++i) a[i].sub = a[i].is_primary = 0;
tmp = opt->a + opt->b > opt->q + opt->r? opt->a + opt->b : opt->q + opt->r;
for (i = 1, m = 1; i < n; ++i) {
for (j = 0; j < m; ++j) {
kv_push(int, z, 0);
for (i = 1; i < n; ++i) {
for (k = 0; k < z.n; ++k) {
int j = z.a[k];
int b_max = a[j].qb > a[i].qb? a[j].qb : a[i].qb;
int e_min = a[j].qe < a[i].qe? a[j].qe : a[i].qe;
if (e_min > b_max) { // have overlap
@ -364,9 +368,10 @@ int mem_choose_alnreg_se(const mem_opt_t *opt, int n, mem_alnreg_t *a) // IMPORT
}
}
}
if (j == m) a[m++] = a[i];
if (k == z.n) kv_push(int, z, i);
}
return m;
for (k = 0; k < z.n; ++k) a[z.a[k]].is_primary = 1;
free(z.a);
}
/************************
@ -601,10 +606,11 @@ void mem_sam_se(const mem_opt_t *opt, const bntseq_t *bns, const uint8_t *pac, b
int k;
kstring_t str;
str.l = str.m = 0; str.s = 0;
a->n = mem_choose_alnreg_se(opt, a->n, a->a); // NOTE: mem_sort_and_dedup() called in worker1()
if (a->n > 0) {
mem_mark_primary_se(opt, a->n, a->a); // NOTE: mem_sort_and_dedup() called in worker1()
for (k = 0; k < a->n; ++k) {
bwahit_t h;
if (!a->a[k].is_primary) continue;
mem_alnreg2hit(&a->a[k], &h);
h.qual = approx_mapq_se(opt, &a->a[k]);
bwa_hit2sam(&str, opt->mat, opt->q, opt->r, opt->w, bns, pac, s, &h, opt->is_hard);

2
bwamem.h
View File

@ -34,7 +34,7 @@ typedef struct {
typedef struct {
int64_t rb, re;
int score, qb, qe, seedcov, sub, csub; // sub: suboptimal score; csub: suboptimal inside the chain
int sub_n;
int sub_n, is_primary;
} mem_alnreg_t;
typedef struct {