|
0
|
1 #include <stdio.h>
|
|
|
2 #include "bwtsw2.h"
|
|
|
3
|
|
|
4 #ifdef USE_MALLOC_WRAPPERS
|
|
|
5 # include "malloc_wrap.h"
|
|
|
6 #endif
|
|
|
7
|
|
|
8 typedef struct {
|
|
|
9 uint32_t tbeg, tend;
|
|
|
10 int qbeg, qend;
|
|
|
11 uint32_t flag:1, idx:31;
|
|
|
12 int chain; // also reuse as a counter
|
|
|
13 } hsaip_t;
|
|
|
14
|
|
|
15 #define _hsaip_lt(a, b) ((a).qbeg < (b).qbeg)
|
|
|
16
|
|
|
17 #include "ksort.h"
|
|
|
18 KSORT_INIT(hsaip, hsaip_t, _hsaip_lt)
|
|
|
19
|
|
|
20 static int chaining(const bsw2opt_t *opt, int shift, int n, hsaip_t *z, hsaip_t *chain)
|
|
|
21 {
|
|
|
22 int j, k, m = 0;
|
|
|
23 ks_introsort(hsaip, n, z);
|
|
|
24 for (j = 0; j < n; ++j) {
|
|
|
25 hsaip_t *p = z + j;
|
|
|
26 for (k = m - 1; k >= 0; --k) {
|
|
|
27 hsaip_t *q = chain + k;
|
|
|
28 int x = p->qbeg - q->qbeg; // always positive
|
|
|
29 int y = p->tbeg - q->tbeg;
|
|
|
30 if (y > 0 && x < opt->max_chain_gap && y < opt->max_chain_gap && x - y <= opt->bw && y - x <= opt->bw) { // chained
|
|
|
31 if (p->qend > q->qend) q->qend = p->qend;
|
|
|
32 if (p->tend > q->tend) q->tend = p->tend;
|
|
|
33 ++q->chain;
|
|
|
34 p->chain = shift + k;
|
|
|
35 break;
|
|
|
36 } else if (q->chain > opt->t_seeds * 2) k = 0; // if the chain is strong enough, do not check the previous chains
|
|
|
37 }
|
|
|
38 if (k < 0) { // not added to any previous chains
|
|
|
39 chain[m] = *p;
|
|
|
40 chain[m].chain = 1;
|
|
|
41 chain[m].idx = p->chain = shift + m;
|
|
|
42 ++m;
|
|
|
43 }
|
|
|
44 }
|
|
|
45 return m;
|
|
|
46 }
|
|
|
47
|
|
|
48 void bsw2_chain_filter(const bsw2opt_t *opt, int len, bwtsw2_t *b[2])
|
|
|
49 {
|
|
|
50 hsaip_t *z[2], *chain[2];
|
|
|
51 int i, j, k, n[2], m[2], thres = opt->t_seeds * 2;
|
|
|
52 char *flag;
|
|
|
53 // initialization
|
|
|
54 n[0] = b[0]->n; n[1] = b[1]->n;
|
|
|
55 z[0] = calloc(n[0] + n[1], sizeof(hsaip_t));
|
|
|
56 z[1] = z[0] + n[0];
|
|
|
57 chain[0] = calloc(n[0] + n[1], sizeof(hsaip_t));
|
|
|
58 for (k = j = 0; k < 2; ++k) {
|
|
|
59 for (i = 0; i < b[k]->n; ++i) {
|
|
|
60 bsw2hit_t *p = b[k]->hits + i;
|
|
|
61 hsaip_t *q = z[k] + i;
|
|
|
62 q->flag = k; q->idx = i;
|
|
|
63 q->tbeg = p->k; q->tend = p->k + p->len;
|
|
|
64 q->chain = -1;
|
|
|
65 q->qbeg = p->beg; q->qend = p->end;
|
|
|
66 }
|
|
|
67 }
|
|
|
68 // chaining
|
|
|
69 m[0] = chaining(opt, 0, n[0], z[0], chain[0]);
|
|
|
70 chain[1] = chain[0] + m[0];
|
|
|
71 m[1] = chaining(opt, m[0], n[1], z[1], chain[1]);
|
|
|
72 // change query coordinate on the reverse strand
|
|
|
73 for (k = 0; k < m[1]; ++k) {
|
|
|
74 hsaip_t *p = chain[1] + k;
|
|
|
75 int tmp = p->qbeg;
|
|
|
76 p->qbeg = len - p->qend; p->qend = len - tmp;
|
|
|
77 }
|
|
|
78 //for (k = 0; k < m[0]; ++k) printf("%d, [%d,%d), [%d,%d)\n", chain[0][k].chain, chain[0][k].tbeg, chain[0][k].tend, chain[0][k].qbeg, chain[0][k].qend);
|
|
|
79 // filtering
|
|
|
80 flag = calloc(m[0] + m[1], 1);
|
|
|
81 ks_introsort(hsaip, m[0] + m[1], chain[0]);
|
|
|
82 for (k = 1; k < m[0] + m[1]; ++k) {
|
|
|
83 hsaip_t *p = chain[0] + k;
|
|
|
84 for (j = 0; j < k; ++j) {
|
|
|
85 hsaip_t *q = chain[0] + j;
|
|
|
86 if (flag[q->idx]) continue;
|
|
|
87 if (q->qend >= p->qend && q->chain > p->chain * thres && p->chain < thres) {
|
|
|
88 flag[p->idx] = 1;
|
|
|
89 break;
|
|
|
90 }
|
|
|
91 }
|
|
|
92 }
|
|
|
93 for (k = 0; k < n[0] + n[1]; ++k) {
|
|
|
94 hsaip_t *p = z[0] + k;
|
|
|
95 if (flag[p->chain])
|
|
|
96 b[p->flag]->hits[p->idx].G = 0;
|
|
|
97 }
|
|
|
98 free(flag);
|
|
|
99 // squeeze out filtered elements in b[2]
|
|
|
100 for (k = 0; k < 2; ++k) {
|
|
|
101 for (j = i = 0; j < n[k]; ++j) {
|
|
|
102 bsw2hit_t *p = b[k]->hits + j;
|
|
|
103 if (p->G) {
|
|
|
104 if (i != j) b[k]->hits[i++] = *p;
|
|
|
105 else ++i;
|
|
|
106 }
|
|
|
107 }
|
|
|
108 b[k]->n = i;
|
|
|
109 }
|
|
|
110 // free
|
|
|
111 free(z[0]); free(chain[0]);
|
|
|
112 }
|
