|
0
|
1 #include <stdlib.h>
|
|
|
2 #include <stdio.h>
|
|
|
3 #include <math.h>
|
|
|
4 #ifdef HAVE_CONFIG_H
|
|
|
5 #include "config.h"
|
|
|
6 #endif
|
|
|
7 #ifdef HAVE_PTHREAD
|
|
|
8 #include <pthread.h>
|
|
|
9 #endif
|
|
|
10 #include "bntseq.h"
|
|
|
11 #include "bwt_lite.h"
|
|
|
12 #include "utils.h"
|
|
|
13 #include "bwtsw2.h"
|
|
|
14 #include "kstring.h"
|
|
|
15 #include "bwa.h"
|
|
|
16 #include "ksw.h"
|
|
|
17
|
|
|
18 #include "kseq.h"
|
|
|
19 KSEQ_DECLARE(gzFile)
|
|
|
20
|
|
|
21 #include "ksort.h"
|
|
|
22 #define __left_lt(a, b) ((a).end > (b).end)
|
|
|
23 KSORT_INIT(hit, bsw2hit_t, __left_lt)
|
|
|
24
|
|
|
25 #ifdef USE_MALLOC_WRAPPERS
|
|
|
26 # include "malloc_wrap.h"
|
|
|
27 #endif
|
|
|
28
|
|
|
29
|
|
|
30 extern unsigned char nst_nt4_table[256];
|
|
|
31
|
|
|
32 unsigned char nt_comp_table[256] = {
|
|
|
33 'N','N','N','N', 'N','N','N','N', 'N','N','N','N', 'N','N','N','N',
|
|
|
34 'N','N','N','N', 'N','N','N','N', 'N','N','N','N', 'N','N','N','N',
|
|
|
35 'N','N','N','N', 'N','N','N','N', 'N','N','N','N', 'N','N','N','N',
|
|
|
36 'N','N','N','N', 'N','N','N','N', 'N','N','N','N', 'N','N','N','N',
|
|
|
37 'N','T','V','G', 'H','N','N','C', 'D','N','N','M', 'N','K','N','N',
|
|
|
38 'N','N','Y','S', 'A','N','B','W', 'X','R','N','N', 'N','N','N','N',
|
|
|
39 'n','t','v','g', 'h','n','n','c', 'd','n','n','m', 'n','k','n','n',
|
|
|
40 'n','n','y','s', 'a','n','b','w', 'x','r','n','N', 'N','N','N','N',
|
|
|
41 'N','N','N','N', 'N','N','N','N', 'N','N','N','N', 'N','N','N','N',
|
|
|
42 'N','N','N','N', 'N','N','N','N', 'N','N','N','N', 'N','N','N','N',
|
|
|
43 'N','N','N','N', 'N','N','N','N', 'N','N','N','N', 'N','N','N','N',
|
|
|
44 'N','N','N','N', 'N','N','N','N', 'N','N','N','N', 'N','N','N','N',
|
|
|
45 'N','N','N','N', 'N','N','N','N', 'N','N','N','N', 'N','N','N','N',
|
|
|
46 'N','N','N','N', 'N','N','N','N', 'N','N','N','N', 'N','N','N','N',
|
|
|
47 'N','N','N','N', 'N','N','N','N', 'N','N','N','N', 'N','N','N','N',
|
|
|
48 'N','N','N','N', 'N','N','N','N', 'N','N','N','N', 'N','N','N','N'
|
|
|
49 };
|
|
|
50
|
|
|
51 extern int bsw2_resolve_duphits(const bntseq_t *bns, const bwt_t *bwt, bwtsw2_t *b, int IS);
|
|
|
52 extern int bsw2_resolve_query_overlaps(bwtsw2_t *b, float mask_level);
|
|
|
53
|
|
|
54 bsw2opt_t *bsw2_init_opt()
|
|
|
55 {
|
|
|
56 bsw2opt_t *o = (bsw2opt_t*)calloc(1, sizeof(bsw2opt_t));
|
|
|
57 o->a = 1; o->b = 3; o->q = 5; o->r = 2; o->t = 30;
|
|
|
58 o->bw = 50;
|
|
|
59 o->max_ins = 20000;
|
|
|
60 o->z = 1; o->is = 3; o->t_seeds = 5; o->hard_clip = 0; o->skip_sw = 0;
|
|
|
61 o->mask_level = 0.50f; o->coef = 5.5f;
|
|
|
62 o->qr = o->q + o->r; o->n_threads = 1; o->chunk_size = 10000000;
|
|
|
63 o->max_chain_gap = 10000;
|
|
|
64 o->cpy_cmt = 0;
|
|
|
65 return o;
|
|
|
66 }
|
|
|
67
|
|
|
68 void bsw2_destroy(bwtsw2_t *b)
|
|
|
69 {
|
|
|
70 int i;
|
|
|
71 if (b == 0) return;
|
|
|
72 if (b->aux)
|
|
|
73 for (i = 0; i < b->n; ++i) free(b->aux[i].cigar);
|
|
|
74 free(b->aux); free(b->hits);
|
|
|
75 free(b);
|
|
|
76 }
|
|
|
77
|
|
|
78 bwtsw2_t *bsw2_dup_no_cigar(const bwtsw2_t *b)
|
|
|
79 {
|
|
|
80 bwtsw2_t *p;
|
|
|
81 p = calloc(1, sizeof(bwtsw2_t));
|
|
|
82 p->max = p->n = b->n;
|
|
|
83 if (b->n) {
|
|
|
84 kroundup32(p->max);
|
|
|
85 p->hits = calloc(p->max, sizeof(bsw2hit_t));
|
|
|
86 memcpy(p->hits, b->hits, p->n * sizeof(bsw2hit_t));
|
|
|
87 }
|
|
|
88 return p;
|
|
|
89 }
|
|
|
90
|
|
|
91 #define __gen_ap(par, opt) do { \
|
|
|
92 int i; \
|
|
|
93 for (i = 0; i < 25; ++i) (par).matrix[i] = -(opt)->b; \
|
|
|
94 for (i = 0; i < 4; ++i) (par).matrix[i*5+i] = (opt)->a; \
|
|
|
95 (par).gap_open = (opt)->q; (par).gap_ext = (opt)->r; \
|
|
|
96 (par).gap_end = (opt)->r; \
|
|
|
97 (par).row = 5; (par).band_width = opt->bw; \
|
|
|
98 } while (0)
|
|
|
99
|
|
|
100 void bsw2_extend_left(const bsw2opt_t *opt, bwtsw2_t *b, uint8_t *_query, int lq, uint8_t *pac, bwtint_t l_pac, uint8_t *_mem)
|
|
|
101 {
|
|
|
102 int i;
|
|
|
103 bwtint_t k;
|
|
|
104 uint8_t *target = 0, *query;
|
|
|
105 int8_t mat[25];
|
|
|
106
|
|
|
107 bwa_fill_scmat(opt->a, opt->b, mat);
|
|
|
108 query = calloc(lq, 1);
|
|
|
109 // sort according to the descending order of query end
|
|
|
110 ks_introsort(hit, b->n, b->hits);
|
|
|
111 target = calloc(((lq + 1) / 2 * opt->a + opt->r) / opt->r + lq, 1);
|
|
|
112 // reverse _query
|
|
|
113 for (i = 0; i < lq; ++i) query[lq - i - 1] = _query[i];
|
|
|
114 // core loop
|
|
|
115 for (i = 0; i < b->n; ++i) {
|
|
|
116 bsw2hit_t *p = b->hits + i;
|
|
|
117 int lt = ((p->beg + 1) / 2 * opt->a + opt->r) / opt->r + lq;
|
|
|
118 int score, j, qle, tle;
|
|
|
119 p->n_seeds = 1;
|
|
|
120 if (p->l || p->k == 0) continue;
|
|
|
121 for (j = score = 0; j < i; ++j) {
|
|
|
122 bsw2hit_t *q = b->hits + j;
|
|
|
123 if (q->beg <= p->beg && q->k <= p->k && q->k + q->len >= p->k + p->len) {
|
|
|
124 if (q->n_seeds < (1<<13) - 2) ++q->n_seeds;
|
|
|
125 ++score;
|
|
|
126 }
|
|
|
127 }
|
|
|
128 if (score) continue;
|
|
|
129 if (lt > p->k) lt = p->k;
|
|
|
130 for (k = p->k - 1, j = 0; k > 0 && j < lt; --k) // FIXME: k=0 not considered!
|
|
|
131 target[j++] = pac[k>>2] >> (~k&3)*2 & 0x3;
|
|
|
132 lt = j;
|
|
|
133 score = ksw_extend(p->beg, &query[lq - p->beg], lt, target, 5, mat, opt->q, opt->r, opt->bw, 0, -1, p->G, &qle, &tle, 0, 0, 0);
|
|
|
134 if (score > p->G) { // extensible
|
|
|
135 p->G = score;
|
|
|
136 p->k -= tle;
|
|
|
137 p->len += tle;
|
|
|
138 p->beg -= qle;
|
|
|
139 }
|
|
|
140 }
|
|
|
141 free(query); free(target);
|
|
|
142 }
|
|
|
143
|
|
|
144 void bsw2_extend_rght(const bsw2opt_t *opt, bwtsw2_t *b, uint8_t *query, int lq, uint8_t *pac, bwtint_t l_pac, uint8_t *_mem)
|
|
|
145 {
|
|
|
146 int i;
|
|
|
147 bwtint_t k;
|
|
|
148 uint8_t *target;
|
|
|
149 int8_t mat[25];
|
|
|
150
|
|
|
151 bwa_fill_scmat(opt->a, opt->b, mat);
|
|
|
152 target = calloc(((lq + 1) / 2 * opt->a + opt->r) / opt->r + lq, 1);
|
|
|
153 for (i = 0; i < b->n; ++i) {
|
|
|
154 bsw2hit_t *p = b->hits + i;
|
|
|
155 int lt = ((lq - p->beg + 1) / 2 * opt->a + opt->r) / opt->r + lq;
|
|
|
156 int j, score, qle, tle;
|
|
|
157 if (p->l) continue;
|
|
|
158 for (k = p->k, j = 0; k < p->k + lt && k < l_pac; ++k)
|
|
|
159 target[j++] = pac[k>>2] >> (~k&3)*2 & 0x3;
|
|
|
160 lt = j;
|
|
|
161 score = ksw_extend(lq - p->beg, &query[p->beg], lt, target, 5, mat, opt->q, opt->r, opt->bw, 0, -1, 1, &qle, &tle, 0, 0, 0) - 1;
|
|
|
162 // if (score < p->G) fprintf(stderr, "[bsw2_extend_hits] %d < %d\n", score, p->G);
|
|
|
163 if (score >= p->G) {
|
|
|
164 p->G = score;
|
|
|
165 p->len = tle;
|
|
|
166 p->end = p->beg + qle;
|
|
|
167 }
|
|
|
168 }
|
|
|
169 free(target);
|
|
|
170 }
|
|
|
171
|
|
|
172 /* generate CIGAR array(s) in b->cigar[] */
|
|
|
173 static void gen_cigar(const bsw2opt_t *opt, int lq, uint8_t *seq[2], int64_t l_pac, const uint8_t *pac, bwtsw2_t *b, const char *name)
|
|
|
174 {
|
|
|
175 int i;
|
|
|
176 int8_t mat[25];
|
|
|
177
|
|
|
178 bwa_fill_scmat(opt->a, opt->b, mat);
|
|
|
179 for (i = 0; i < b->n; ++i) {
|
|
|
180 bsw2hit_t *p = b->hits + i;
|
|
|
181 bsw2aux_t *q = b->aux + i;
|
|
|
182 uint8_t *query;
|
|
|
183 int beg, end, score;
|
|
|
184 if (p->l) continue;
|
|
|
185 beg = (p->flag & 0x10)? lq - p->end : p->beg;
|
|
|
186 end = (p->flag & 0x10)? lq - p->beg : p->end;
|
|
|
187 query = seq[(p->flag & 0x10)? 1 : 0] + beg;
|
|
|
188 q->cigar = bwa_gen_cigar(mat, opt->q, opt->r, opt->bw, l_pac, pac, end - beg, query, p->k, p->k + p->len, &score, &q->n_cigar, &q->nm);
|
|
|
189 #if 0
|
|
|
190 if (name && score != p->G) { // debugging only
|
|
|
191 int j, glen = 0;
|
|
|
192 for (j = 0; j < q->n_cigar; ++j)
|
|
|
193 if ((q->cigar[j]&0xf) == 1 || (q->cigar[j]&0xf) == 2)
|
|
|
194 glen += q->cigar[j]>>4;
|
|
|
195 fprintf(stderr, "[E::%s] %s - unequal score: %d != %d; (qlen, aqlen, arlen, glen, bw) = (%d, %d, %d, %d, %d)\n",
|
|
|
196 __func__, name, score, p->G, lq, end - beg, p->len, glen, opt->bw);
|
|
|
197 }
|
|
|
198 #endif
|
|
|
199 if (q->cigar && (beg != 0 || end < lq)) { // write soft clipping
|
|
|
200 q->cigar = realloc(q->cigar, 4 * (q->n_cigar + 2));
|
|
|
201 if (beg != 0) {
|
|
|
202 memmove(q->cigar + 1, q->cigar, q->n_cigar * 4);
|
|
|
203 q->cigar[0] = beg<<4 | 4;
|
|
|
204 ++q->n_cigar;
|
|
|
205 }
|
|
|
206 if (end < lq) {
|
|
|
207 q->cigar[q->n_cigar] = (lq - end)<<4 | 4;
|
|
|
208 ++q->n_cigar;
|
|
|
209 }
|
|
|
210 }
|
|
|
211 }
|
|
|
212 }
|
|
|
213
|
|
|
214 /* this is for the debugging purpose only */
|
|
|
215 void bsw2_debug_hits(const bwtsw2_t *b)
|
|
|
216 {
|
|
|
217 int i;
|
|
|
218 printf("# raw hits: %d\n", b->n);
|
|
|
219 for (i = 0; i < b->n; ++i) {
|
|
|
220 bsw2hit_t *p = b->hits + i;
|
|
|
221 if (p->G > 0)
|
|
|
222 printf("G=%d, G2=%d, len=%d, [%d,%d), k=%lu, l=%lu, #seeds=%d, is_rev=%d\n", p->G, p->G2, p->len, p->beg, p->end, (long)p->k, (long)p->l, p->n_seeds, p->is_rev);
|
|
|
223 }
|
|
|
224 }
|
|
|
225
|
|
|
226 static void merge_hits(bwtsw2_t *b[2], int l, int is_reverse)
|
|
|
227 {
|
|
|
228 int i;
|
|
|
229 if (b[0]->n + b[1]->n > b[0]->max) {
|
|
|
230 b[0]->max = b[0]->n + b[1]->n;
|
|
|
231 b[0]->hits = realloc(b[0]->hits, b[0]->max * sizeof(bsw2hit_t));
|
|
|
232 }
|
|
|
233 for (i = 0; i < b[1]->n; ++i) {
|
|
|
234 bsw2hit_t *p = b[0]->hits + b[0]->n + i;
|
|
|
235 *p = b[1]->hits[i];
|
|
|
236 if (is_reverse) {
|
|
|
237 int x = p->beg;
|
|
|
238 p->beg = l - p->end;
|
|
|
239 p->end = l - x;
|
|
|
240 p->flag |= 0x10;
|
|
|
241 }
|
|
|
242 }
|
|
|
243 b[0]->n += b[1]->n;
|
|
|
244 bsw2_destroy(b[1]);
|
|
|
245 b[1] = 0;
|
|
|
246 }
|
|
|
247 /* seq[0] is the forward sequence and seq[1] is the reverse complement. */
|
|
|
248 static bwtsw2_t *bsw2_aln1_core(const bsw2opt_t *opt, const bntseq_t *bns, uint8_t *pac, const bwt_t *target,
|
|
|
249 int l, uint8_t *seq[2], bsw2global_t *pool)
|
|
|
250 {
|
|
|
251 extern void bsw2_chain_filter(const bsw2opt_t *opt, int len, bwtsw2_t *b[2]);
|
|
|
252 bwtsw2_t *b[2], **bb[2], **_b, *p;
|
|
|
253 int k, j;
|
|
|
254 bwtl_t *query;
|
|
|
255 query = bwtl_seq2bwtl(l, seq[0]);
|
|
|
256 _b = bsw2_core(bns, opt, query, target, pool);
|
|
|
257 bwtl_destroy(query);
|
|
|
258 for (k = 0; k < 2; ++k) {
|
|
|
259 bb[k] = calloc(2, sizeof(void*));
|
|
|
260 bb[k][0] = calloc(1, sizeof(bwtsw2_t));
|
|
|
261 bb[k][1] = calloc(1, sizeof(bwtsw2_t));
|
|
|
262 }
|
|
|
263 for (k = 0; k < 2; ++k) { // separate _b into bb[2] based on the strand
|
|
|
264 for (j = 0; j < _b[k]->n; ++j) {
|
|
|
265 bsw2hit_t *q;
|
|
|
266 p = bb[_b[k]->hits[j].is_rev][k];
|
|
|
267 if (p->n == p->max) {
|
|
|
268 p->max = p->max? p->max<<1 : 8;
|
|
|
269 p->hits = realloc(p->hits, p->max * sizeof(bsw2hit_t));
|
|
|
270 }
|
|
|
271 q = &p->hits[p->n++];
|
|
|
272 *q = _b[k]->hits[j];
|
|
|
273 if (_b[k]->hits[j].is_rev) {
|
|
|
274 int x = q->beg;
|
|
|
275 q->beg = l - q->end;
|
|
|
276 q->end = l - x;
|
|
|
277 }
|
|
|
278 }
|
|
|
279 }
|
|
|
280 b[0] = bb[0][1]; b[1] = bb[1][1]; // bb[*][1] are "narrow SA hits"
|
|
|
281 bsw2_chain_filter(opt, l, b); // NB: only unique seeds are chained
|
|
|
282 for (k = 0; k < 2; ++k) {
|
|
|
283 bsw2_extend_left(opt, bb[k][1], seq[k], l, pac, bns->l_pac, pool->aln_mem);
|
|
|
284 merge_hits(bb[k], l, 0); // bb[k][1] is merged to bb[k][0] here
|
|
|
285 bsw2_resolve_duphits(0, 0, bb[k][0], 0);
|
|
|
286 bsw2_extend_rght(opt, bb[k][0], seq[k], l, pac, bns->l_pac, pool->aln_mem);
|
|
|
287 bsw2_resolve_duphits(0, 0, bb[k][0], 0);
|
|
|
288 b[k] = bb[k][0];
|
|
|
289 free(bb[k]);
|
|
|
290 }
|
|
|
291 merge_hits(b, l, 1); // again, b[1] is merged to b[0]
|
|
|
292 bsw2_resolve_query_overlaps(b[0], opt->mask_level);
|
|
|
293 bsw2_destroy(_b[0]); bsw2_destroy(_b[1]); free(_b);
|
|
|
294 return b[0];
|
|
|
295 }
|
|
|
296
|
|
|
297 /* set ->flag to records the origin of the hit (to forward bwt or reverse bwt) */
|
|
|
298 static void flag_fr(bwtsw2_t *b[2])
|
|
|
299 {
|
|
|
300 int i, j;
|
|
|
301 for (i = 0; i < b[0]->n; ++i) {
|
|
|
302 bsw2hit_t *p = b[0]->hits + i;
|
|
|
303 p->flag |= 0x10000;
|
|
|
304 }
|
|
|
305 for (i = 0; i < b[1]->n; ++i) {
|
|
|
306 bsw2hit_t *p = b[1]->hits + i;
|
|
|
307 p->flag |= 0x20000;
|
|
|
308 }
|
|
|
309 for (i = 0; i < b[0]->n; ++i) {
|
|
|
310 bsw2hit_t *p = b[0]->hits + i;
|
|
|
311 for (j = 0; j < b[1]->n; ++j) {
|
|
|
312 bsw2hit_t *q = b[1]->hits + j;
|
|
|
313 if (q->beg == p->beg && q->end == p->end && q->k == p->k && q->len == p->len && q->G == p->G) {
|
|
|
314 q->flag |= 0x30000; p->flag |= 0x30000;
|
|
|
315 break;
|
|
|
316 }
|
|
|
317 }
|
|
|
318 }
|
|
|
319 }
|
|
|
320
|
|
|
321 typedef struct {
|
|
|
322 int n, max;
|
|
|
323 bsw2seq1_t *seq;
|
|
|
324 } bsw2seq_t;
|
|
|
325
|
|
|
326 static int fix_cigar(const bntseq_t *bns, bsw2hit_t *p, int n_cigar, uint32_t *cigar)
|
|
|
327 {
|
|
|
328 // FIXME: this routine does not work if the query bridge three reference sequences
|
|
|
329 int32_t coor, refl, lq;
|
|
|
330 int x, y, i, seqid;
|
|
|
331 bns_cnt_ambi(bns, p->k, p->len, &seqid);
|
|
|
332 coor = p->k - bns->anns[seqid].offset;
|
|
|
333 refl = bns->anns[seqid].len;
|
|
|
334 x = coor; y = 0;
|
|
|
335 // test if the alignment goes beyond the boundary
|
|
|
336 for (i = 0; i < n_cigar; ++i) {
|
|
|
337 int op = cigar[i]&0xf, ln = cigar[i]>>4;
|
|
|
338 if (op == 1 || op == 4 || op == 5) y += ln;
|
|
|
339 else if (op == 2) x += ln;
|
|
|
340 else x += ln, y += ln;
|
|
|
341 }
|
|
|
342 lq = y; // length of the query sequence
|
|
|
343 if (x > refl) { // then fix it
|
|
|
344 int j, nc, mq[2], nlen[2];
|
|
|
345 uint32_t *cn;
|
|
|
346 bwtint_t kk = 0;
|
|
|
347 nc = mq[0] = mq[1] = nlen[0] = nlen[1] = 0;
|
|
|
348 cn = calloc(n_cigar + 3, 4);
|
|
|
349 x = coor; y = 0;
|
|
|
350 for (i = j = 0; i < n_cigar; ++i) {
|
|
|
351 int op = cigar[i]&0xf, ln = cigar[i]>>4;
|
|
|
352 if (op == 4 || op == 5 || op == 1) { // ins or clipping
|
|
|
353 y += ln;
|
|
|
354 cn[j++] = cigar[i];
|
|
|
355 } else if (op == 2) { // del
|
|
|
356 if (x + ln >= refl && nc == 0) {
|
|
|
357 cn[j++] = (uint32_t)(lq - y)<<4 | 4;
|
|
|
358 nc = j;
|
|
|
359 cn[j++] = (uint32_t)y<<4 | 4;
|
|
|
360 kk = p->k + (x + ln - refl);
|
|
|
361 nlen[0] = x - coor;
|
|
|
362 nlen[1] = p->len - nlen[0] - ln;
|
|
|
363 } else cn[j++] = cigar[i];
|
|
|
364 x += ln;
|
|
|
365 } else if (op == 0) { // match
|
|
|
366 if (x + ln >= refl && nc == 0) {
|
|
|
367 // FIXME: not consider a special case where a split right between M and I
|
|
|
368 cn[j++] = (uint32_t)(refl - x)<<4 | 0; // write M
|
|
|
369 cn[j++] = (uint32_t)(lq - y - (refl - x))<<4 | 4; // write S
|
|
|
370 nc = j;
|
|
|
371 mq[0] += refl - x;
|
|
|
372 cn[j++] = (uint32_t)(y + (refl - x))<<4 | 4;
|
|
|
373 if (x + ln - refl) cn[j++] = (uint32_t)(x + ln - refl)<<4 | 0;
|
|
|
374 mq[1] += x + ln - refl;
|
|
|
375 kk = bns->anns[seqid].offset + refl;
|
|
|
376 nlen[0] = refl - coor;
|
|
|
377 nlen[1] = p->len - nlen[0];
|
|
|
378 } else {
|
|
|
379 cn[j++] = cigar[i];
|
|
|
380 mq[nc?1:0] += ln;
|
|
|
381 }
|
|
|
382 x += ln; y += ln;
|
|
|
383 }
|
|
|
384 }
|
|
|
385 if (mq[0] > mq[1]) { // then take the first alignment
|
|
|
386 n_cigar = nc;
|
|
|
387 memcpy(cigar, cn, 4 * nc);
|
|
|
388 p->len = nlen[0];
|
|
|
389 } else {
|
|
|
390 p->k = kk; p->len = nlen[1];
|
|
|
391 n_cigar = j - nc;
|
|
|
392 memcpy(cigar, cn + nc, 4 * (j - nc));
|
|
|
393 }
|
|
|
394 free(cn);
|
|
|
395 }
|
|
|
396 return n_cigar;
|
|
|
397 }
|
|
|
398
|
|
|
399 static void write_aux(const bsw2opt_t *opt, const bntseq_t *bns, int qlen, uint8_t *seq[2], const uint8_t *pac, bwtsw2_t *b, const char *name)
|
|
|
400 {
|
|
|
401 int i;
|
|
|
402 // allocate for b->aux
|
|
|
403 if (b->n<<1 < b->max) {
|
|
|
404 b->max = b->n;
|
|
|
405 kroundup32(b->max);
|
|
|
406 b->hits = realloc(b->hits, b->max * sizeof(bsw2hit_t));
|
|
|
407 }
|
|
|
408 b->aux = calloc(b->n, sizeof(bsw2aux_t));
|
|
|
409 // generate CIGAR
|
|
|
410 gen_cigar(opt, qlen, seq, bns->l_pac, pac, b, name);
|
|
|
411 // fix CIGAR, generate mapQ, and write chromosomal position
|
|
|
412 for (i = 0; i < b->n; ++i) {
|
|
|
413 bsw2hit_t *p = &b->hits[i];
|
|
|
414 bsw2aux_t *q = &b->aux[i];
|
|
|
415 q->flag = p->flag & 0xfe;
|
|
|
416 q->isize = 0;
|
|
|
417 if (p->l == 0) { // unique hit
|
|
|
418 float c = 1.0;
|
|
|
419 int subo;
|
|
|
420 // fix out-of-boundary CIGAR
|
|
|
421 q->n_cigar = fix_cigar(bns, p, q->n_cigar, q->cigar);
|
|
|
422 // compute mapQ
|
|
|
423 subo = p->G2 > opt->t? p->G2 : opt->t;
|
|
|
424 if (p->flag>>16 == 1 || p->flag>>16 == 2) c *= .5;
|
|
|
425 if (p->n_seeds < 2) c *= .2;
|
|
|
426 q->qual = (int)(c * (p->G - subo) * (250.0 / p->G + 0.03 / opt->a) + .499);
|
|
|
427 if (q->qual > 250) q->qual = 250;
|
|
|
428 if (q->qual < 0) q->qual = 0;
|
|
|
429 if (p->flag&1) q->qual = 0; // this is a random hit
|
|
|
430 q->pqual = q->qual; // set the paired qual as qual
|
|
|
431 // get the chromosomal position
|
|
|
432 q->nn = bns_cnt_ambi(bns, p->k, p->len, &q->chr);
|
|
|
433 q->pos = p->k - bns->anns[q->chr].offset;
|
|
|
434 } else q->qual = 0, q->n_cigar = 0, q->chr = q->pos = -1, q->nn = 0;
|
|
|
435 }
|
|
|
436 }
|
|
|
437
|
|
|
438 static void update_mate_aux(bwtsw2_t *b, const bwtsw2_t *m)
|
|
|
439 {
|
|
|
440 int i;
|
|
|
441 if (m == 0) return;
|
|
|
442 // update flag, mchr and mpos
|
|
|
443 for (i = 0; i < b->n; ++i) {
|
|
|
444 bsw2aux_t *q = &b->aux[i];
|
|
|
445 q->flag |= 1; // paired
|
|
|
446 if (m->n == 0) q->flag |= 8; // mate unmapped
|
|
|
447 if (m->n == 1) {
|
|
|
448 q->mchr = m->aux[0].chr;
|
|
|
449 q->mpos = m->aux[0].pos;
|
|
|
450 if (m->aux[0].flag&0x10) q->flag |= 0x20; // mate reverse strand
|
|
|
451 if (q->chr == q->mchr) { // set insert size
|
|
|
452 if (q->mpos + m->hits[0].len > q->pos)
|
|
|
453 q->isize = q->mpos + m->hits[0].len - q->pos;
|
|
|
454 else q->isize = q->mpos - q->pos - b->hits[0].len;
|
|
|
455 } else q->isize = 0;
|
|
|
456 } else q->mchr = q->mpos = -1;
|
|
|
457 }
|
|
|
458 // update mapping quality
|
|
|
459 if (b->n == 1 && m->n == 1) {
|
|
|
460 bsw2hit_t *p = &b->hits[0];
|
|
|
461 if (p->flag & BSW2_FLAG_MATESW) { // this alignment is found by Smith-Waterman
|
|
|
462 if (!(p->flag & BSW2_FLAG_TANDEM) && b->aux[0].pqual < 20)
|
|
|
463 b->aux[0].pqual = 20;
|
|
|
464 if (b->aux[0].pqual >= m->aux[0].qual) b->aux[0].pqual = m->aux[0].qual;
|
|
|
465 } else if ((p->flag & 2) && !(m->hits[0].flag & BSW2_FLAG_MATESW)) { // properly paired
|
|
|
466 if (!(p->flag & BSW2_FLAG_TANDEM)) { // pqual is bounded by [b->aux[0].qual,m->aux[0].qual]
|
|
|
467 b->aux[0].pqual += 20;
|
|
|
468 if (b->aux[0].pqual > m->aux[0].qual) b->aux[0].pqual = m->aux[0].qual;
|
|
|
469 if (b->aux[0].pqual < b->aux[0].qual) b->aux[0].pqual = b->aux[0].qual;
|
|
|
470 }
|
|
|
471 }
|
|
|
472 }
|
|
|
473 }
|
|
|
474
|
|
|
475 /* generate SAM lines for a sequence in ks with alignment stored in
|
|
|
476 * b. ks->name and ks->seq will be freed and set to NULL in the end. */
|
|
|
477 static void print_hits(const bntseq_t *bns, const bsw2opt_t *opt, bsw2seq1_t *ks, bwtsw2_t *b, int is_pe, bwtsw2_t *bmate)
|
|
|
478 {
|
|
|
479 int i, k;
|
|
|
480 kstring_t str;
|
|
|
481 memset(&str, 0, sizeof(kstring_t));
|
|
|
482 if (b == 0 || b->n == 0) { // no hits
|
|
|
483 ksprintf(&str, "%s\t4\t*\t0\t0\t*\t*\t0\t0\t", ks->name);
|
|
|
484 for (i = 0; i < ks->l; ++i) kputc(ks->seq[i], &str);
|
|
|
485 if (ks->qual) {
|
|
|
486 kputc('\t', &str);
|
|
|
487 for (i = 0; i < ks->l; ++i) kputc(ks->qual[i], &str);
|
|
|
488 } else kputs("\t*", &str);
|
|
|
489 kputc('\n', &str);
|
|
|
490 }
|
|
|
491 for (i = 0; b && i < b->n; ++i) {
|
|
|
492 bsw2hit_t *p = b->hits + i;
|
|
|
493 bsw2aux_t *q = b->aux + i;
|
|
|
494 int j, beg, end, type = 0;
|
|
|
495 // print mandatory fields before SEQ
|
|
|
496 if (q->cigar == 0) q->flag |= 0x4;
|
|
|
497 ksprintf(&str, "%s\t%d", ks->name, q->flag | (opt->multi_2nd && i? 0x100 : 0));
|
|
|
498 ksprintf(&str, "\t%s\t%ld", q->chr>=0? bns->anns[q->chr].name : "*", (long)q->pos + 1);
|
|
|
499 if (p->l == 0 && q->cigar) { // not a repetitive hit
|
|
|
500 ksprintf(&str, "\t%d\t", q->pqual);
|
|
|
501 for (k = 0; k < q->n_cigar; ++k)
|
|
|
502 ksprintf(&str, "%d%c", q->cigar[k]>>4, (opt->hard_clip? "MIDNHHP" : "MIDNSHP")[q->cigar[k]&0xf]);
|
|
|
503 } else ksprintf(&str, "\t0\t*");
|
|
|
504 if (!is_pe) kputs("\t*\t0\t0\t", &str);
|
|
|
505 else ksprintf(&str, "\t%s\t%d\t%d\t", q->mchr==q->chr? "=" : (q->mchr<0? "*" : bns->anns[q->mchr].name), q->mpos+1, q->isize);
|
|
|
506 // get the sequence begin and end
|
|
|
507 beg = 0; end = ks->l;
|
|
|
508 if (opt->hard_clip && q->cigar) {
|
|
|
509 if ((q->cigar[0]&0xf) == 4) beg += q->cigar[0]>>4;
|
|
|
510 if ((q->cigar[q->n_cigar-1]&0xf) == 4) end -= q->cigar[q->n_cigar-1]>>4;
|
|
|
511 }
|
|
|
512 for (j = beg; j < end; ++j) {
|
|
|
513 if (p->flag&0x10) kputc(nt_comp_table[(int)ks->seq[ks->l - 1 - j]], &str);
|
|
|
514 else kputc(ks->seq[j], &str);
|
|
|
515 }
|
|
|
516 // print base quality if present
|
|
|
517 if (ks->qual) {
|
|
|
518 kputc('\t', &str);
|
|
|
519 for (j = beg; j < end; ++j) {
|
|
|
520 if (p->flag&0x10) kputc(ks->qual[ks->l - 1 - j], &str);
|
|
|
521 else kputc(ks->qual[j], &str);
|
|
|
522 }
|
|
|
523 } else kputs("\t*", &str);
|
|
|
524 // print optional tags
|
|
|
525 ksprintf(&str, "\tAS:i:%d\tXS:i:%d\tXF:i:%d\tXE:i:%d\tNM:i:%d", p->G, p->G2, p->flag>>16, p->n_seeds, q->nm);
|
|
|
526 if (q->nn) ksprintf(&str, "\tXN:i:%d", q->nn);
|
|
|
527 if (p->l) ksprintf(&str, "\tXI:i:%d", p->l - p->k + 1);
|
|
|
528 if (p->flag&BSW2_FLAG_MATESW) type |= 1;
|
|
|
529 if (p->flag&BSW2_FLAG_TANDEM) type |= 2;
|
|
|
530 if (type) ksprintf(&str, "\tXT:i:%d", type);
|
|
|
531 if (opt->cpy_cmt && ks->comment) {
|
|
|
532 int l = strlen(ks->comment);
|
|
|
533 if (l >= 6 && ks->comment[2] == ':' && ks->comment[4] == ':') {
|
|
|
534 kputc('\t', &str); kputs(ks->comment, &str);
|
|
|
535 }
|
|
|
536 }
|
|
|
537 kputc('\n', &str);
|
|
|
538 }
|
|
|
539 ks->sam = str.s;
|
|
|
540 free(ks->seq); ks->seq = 0;
|
|
|
541 free(ks->qual); ks->qual = 0;
|
|
|
542 free(ks->name); ks->name = 0;
|
|
|
543 }
|
|
|
544
|
|
|
545 static void update_opt(bsw2opt_t *dst, const bsw2opt_t *src, int qlen)
|
|
|
546 {
|
|
|
547 double ll = log(qlen);
|
|
|
548 int i, k;
|
|
|
549 *dst = *src;
|
|
|
550 if (dst->t < ll * dst->coef) dst->t = (int)(ll * dst->coef + .499);
|
|
|
551 // set band width: the query length sets a boundary on the maximum band width
|
|
|
552 k = (qlen * dst->a - 2 * dst->q) / (2 * dst->r + dst->a);
|
|
|
553 i = (qlen * dst->a - dst->a - dst->t) / dst->r;
|
|
|
554 if (k > i) k = i;
|
|
|
555 if (k < 1) k = 1; // I do not know if k==0 causes troubles
|
|
|
556 dst->bw = src->bw < k? src->bw : k;
|
|
|
557 }
|
|
|
558
|
|
|
559 /* Core routine to align reads in _seq. It is separated from
|
|
|
560 * process_seqs() to realize multi-threading */
|
|
|
561 static void bsw2_aln_core(bsw2seq_t *_seq, const bsw2opt_t *_opt, const bntseq_t *bns, uint8_t *pac, const bwt_t *target, int is_pe)
|
|
|
562 {
|
|
|
563 int x;
|
|
|
564 bsw2opt_t opt;
|
|
|
565 bsw2global_t *pool = bsw2_global_init();
|
|
|
566 bwtsw2_t **buf;
|
|
|
567 buf = calloc(_seq->n, sizeof(void*));
|
|
|
568 for (x = 0; x < _seq->n; ++x) {
|
|
|
569 bsw2seq1_t *p = _seq->seq + x;
|
|
|
570 uint8_t *seq[2], *rseq[2];
|
|
|
571 int i, l, k;
|
|
|
572 bwtsw2_t *b[2];
|
|
|
573 l = p->l;
|
|
|
574 update_opt(&opt, _opt, p->l);
|
|
|
575 if (pool->max_l < l) { // then enlarge working space for aln_extend_core()
|
|
|
576 int tmp = ((l + 1) / 2 * opt.a + opt.r) / opt.r + l;
|
|
|
577 pool->max_l = l;
|
|
|
578 pool->aln_mem = realloc(pool->aln_mem, (tmp + 2) * 24);
|
|
|
579 }
|
|
|
580 // set seq[2] and rseq[2]
|
|
|
581 seq[0] = calloc(l * 4, 1);
|
|
|
582 seq[1] = seq[0] + l;
|
|
|
583 rseq[0] = seq[1] + l; rseq[1] = rseq[0] + l;
|
|
|
584 // convert sequences to 2-bit representation
|
|
|
585 for (i = k = 0; i < l; ++i) {
|
|
|
586 int c = nst_nt4_table[(int)p->seq[i]];
|
|
|
587 if (c >= 4) { c = (int)(drand48() * 4); ++k; } // FIXME: ambiguous bases are not properly handled
|
|
|
588 seq[0][i] = c;
|
|
|
589 seq[1][l-1-i] = 3 - c;
|
|
|
590 rseq[0][l-1-i] = 3 - c;
|
|
|
591 rseq[1][i] = c;
|
|
|
592 }
|
|
|
593 if (l - k < opt.t) { // too few unambiguous bases
|
|
|
594 buf[x] = calloc(1, sizeof(bwtsw2_t));
|
|
|
595 free(seq[0]); continue;
|
|
|
596 }
|
|
|
597 // alignment
|
|
|
598 b[0] = bsw2_aln1_core(&opt, bns, pac, target, l, seq, pool);
|
|
|
599 for (k = 0; k < b[0]->n; ++k)
|
|
|
600 if (b[0]->hits[k].n_seeds < opt.t_seeds) break;
|
|
|
601 if (k < b[0]->n) {
|
|
|
602 b[1] = bsw2_aln1_core(&opt, bns, pac, target, l, rseq, pool);
|
|
|
603 for (i = 0; i < b[1]->n; ++i) {
|
|
|
604 bsw2hit_t *p = &b[1]->hits[i];
|
|
|
605 int x = p->beg;
|
|
|
606 p->flag ^= 0x10, p->is_rev ^= 1; // flip the strand
|
|
|
607 p->beg = l - p->end;
|
|
|
608 p->end = l - x;
|
|
|
609 }
|
|
|
610 flag_fr(b);
|
|
|
611 merge_hits(b, l, 0);
|
|
|
612 bsw2_resolve_duphits(0, 0, b[0], 0);
|
|
|
613 bsw2_resolve_query_overlaps(b[0], opt.mask_level);
|
|
|
614 } else b[1] = 0;
|
|
|
615 // generate CIGAR and print SAM
|
|
|
616 buf[x] = bsw2_dup_no_cigar(b[0]);
|
|
|
617 // free
|
|
|
618 free(seq[0]);
|
|
|
619 bsw2_destroy(b[0]);
|
|
|
620 }
|
|
|
621 if (is_pe) bsw2_pair(&opt, bns->l_pac, pac, _seq->n, _seq->seq, buf);
|
|
|
622 for (x = 0; x < _seq->n; ++x) {
|
|
|
623 bsw2seq1_t *p = _seq->seq + x;
|
|
|
624 uint8_t *seq[2];
|
|
|
625 int i;
|
|
|
626 seq[0] = malloc(p->l * 2); seq[1] = seq[0] + p->l;
|
|
|
627 for (i = 0; i < p->l; ++i) {
|
|
|
628 int c = nst_nt4_table[(int)p->seq[i]];
|
|
|
629 if (c >= 4) c = (int)(drand48() * 4);
|
|
|
630 seq[0][i] = c;
|
|
|
631 seq[1][p->l-1-i] = 3 - c;
|
|
|
632 }
|
|
|
633 update_opt(&opt, _opt, p->l);
|
|
|
634 write_aux(&opt, bns, p->l, seq, pac, buf[x], _seq->seq[x].name);
|
|
|
635 free(seq[0]);
|
|
|
636 }
|
|
|
637 for (x = 0; x < _seq->n; ++x) {
|
|
|
638 if (is_pe) update_mate_aux(buf[x], buf[x^1]);
|
|
|
639 print_hits(bns, &opt, &_seq->seq[x], buf[x], is_pe, buf[x^1]);
|
|
|
640 }
|
|
|
641 for (x = 0; x < _seq->n; ++x) bsw2_destroy(buf[x]);
|
|
|
642 free(buf);
|
|
|
643 bsw2_global_destroy(pool);
|
|
|
644 }
|
|
|
645
|
|
|
646 #ifdef HAVE_PTHREAD
|
|
|
647 typedef struct {
|
|
|
648 int tid, is_pe;
|
|
|
649 bsw2seq_t *_seq;
|
|
|
650 const bsw2opt_t *_opt;
|
|
|
651 const bntseq_t *bns;
|
|
|
652 uint8_t *pac;
|
|
|
653 const bwt_t *target;
|
|
|
654 } thread_aux_t;
|
|
|
655
|
|
|
656 /* another interface to bsw2_aln_core() to facilitate pthread_create() */
|
|
|
657 static void *worker(void *data)
|
|
|
658 {
|
|
|
659 thread_aux_t *p = (thread_aux_t*)data;
|
|
|
660 bsw2_aln_core(p->_seq, p->_opt, p->bns, p->pac, p->target, p->is_pe);
|
|
|
661 return 0;
|
|
|
662 }
|
|
|
663 #endif
|
|
|
664
|
|
|
665 /* process sequences stored in _seq, generate SAM lines for these
|
|
|
666 * sequences and reset _seq afterwards. */
|
|
|
667 static void process_seqs(bsw2seq_t *_seq, const bsw2opt_t *opt, const bntseq_t *bns, uint8_t *pac, const bwt_t *target, int is_pe)
|
|
|
668 {
|
|
|
669 int i;
|
|
|
670 is_pe = is_pe? 1 : 0;
|
|
|
671
|
|
|
672 #ifdef HAVE_PTHREAD
|
|
|
673 if (opt->n_threads <= 1) {
|
|
|
674 bsw2_aln_core(_seq, opt, bns, pac, target, is_pe);
|
|
|
675 } else {
|
|
|
676 pthread_t *tid;
|
|
|
677 pthread_attr_t attr;
|
|
|
678 thread_aux_t *data;
|
|
|
679 int j;
|
|
|
680 pthread_attr_init(&attr);
|
|
|
681 pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
|
|
|
682 data = (thread_aux_t*)calloc(opt->n_threads, sizeof(thread_aux_t));
|
|
|
683 tid = (pthread_t*)calloc(opt->n_threads, sizeof(pthread_t));
|
|
|
684 for (j = 0; j < opt->n_threads; ++j) {
|
|
|
685 thread_aux_t *p = data + j;
|
|
|
686 p->tid = j; p->_opt = opt; p->bns = bns; p->is_pe = is_pe;
|
|
|
687 p->pac = pac; p->target = target;
|
|
|
688 p->_seq = calloc(1, sizeof(bsw2seq_t));
|
|
|
689 p->_seq->max = (_seq->n + opt->n_threads - 1) / opt->n_threads + 1;
|
|
|
690 p->_seq->n = 0;
|
|
|
691 p->_seq->seq = calloc(p->_seq->max, sizeof(bsw2seq1_t));
|
|
|
692 }
|
|
|
693 for (i = 0; i < _seq->n; ++i) { // assign sequences to each thread
|
|
|
694 bsw2seq_t *p = data[(i>>is_pe)%opt->n_threads]._seq;
|
|
|
695 p->seq[p->n++] = _seq->seq[i];
|
|
|
696 }
|
|
|
697 for (j = 0; j < opt->n_threads; ++j) pthread_create(&tid[j], &attr, worker, &data[j]);
|
|
|
698 for (j = 0; j < opt->n_threads; ++j) pthread_join(tid[j], 0);
|
|
|
699 for (j = 0; j < opt->n_threads; ++j) data[j]._seq->n = 0;
|
|
|
700 for (i = 0; i < _seq->n; ++i) { // copy the result from each thread back
|
|
|
701 bsw2seq_t *p = data[(i>>is_pe)%opt->n_threads]._seq;
|
|
|
702 _seq->seq[i] = p->seq[p->n++];
|
|
|
703 }
|
|
|
704 for (j = 0; j < opt->n_threads; ++j) {
|
|
|
705 thread_aux_t *p = data + j;
|
|
|
706 free(p->_seq->seq);
|
|
|
707 free(p->_seq);
|
|
|
708 }
|
|
|
709 free(data); free(tid);
|
|
|
710 }
|
|
|
711 #else
|
|
|
712 bsw2_aln_core(_seq, opt, bns, pac, target, is_pe);
|
|
|
713 #endif
|
|
|
714
|
|
|
715 // print and reset
|
|
|
716 for (i = 0; i < _seq->n; ++i) {
|
|
|
717 bsw2seq1_t *p = _seq->seq + i;
|
|
|
718 if (p->sam) err_printf("%s", p->sam);
|
|
|
719 free(p->name); free(p->seq); free(p->qual); free(p->sam);
|
|
|
720 p->tid = -1; p->l = 0;
|
|
|
721 p->name = p->seq = p->qual = p->sam = 0;
|
|
|
722 }
|
|
|
723 err_fflush(stdout);
|
|
|
724 _seq->n = 0;
|
|
|
725 }
|
|
|
726
|
|
|
727 void bsw2_aln(const bsw2opt_t *opt, const bntseq_t *bns, bwt_t * const target, const char *fn, const char *fn2)
|
|
|
728 {
|
|
|
729 gzFile fp, fp2;
|
|
|
730 kseq_t *ks, *ks2;
|
|
|
731 int l, is_pe = 0, i, n;
|
|
|
732 uint8_t *pac;
|
|
|
733 bsw2seq_t *_seq;
|
|
|
734 bseq1_t *bseq;
|
|
|
735
|
|
|
736 pac = calloc(bns->l_pac/4+1, 1);
|
|
|
737 for (l = 0; l < bns->n_seqs; ++l)
|
|
|
738 err_printf("@SQ\tSN:%s\tLN:%d\n", bns->anns[l].name, bns->anns[l].len);
|
|
|
739 err_fread_noeof(pac, 1, bns->l_pac/4+1, bns->fp_pac);
|
|
|
740 fp = xzopen(fn, "r");
|
|
|
741 ks = kseq_init(fp);
|
|
|
742 _seq = calloc(1, sizeof(bsw2seq_t));
|
|
|
743 if (fn2) {
|
|
|
744 fp2 = xzopen(fn2, "r");
|
|
|
745 ks2 = kseq_init(fp2);
|
|
|
746 is_pe = 1;
|
|
|
747 } else fp2 = 0, ks2 = 0, is_pe = 0;
|
|
|
748 while ((bseq = bseq_read(opt->chunk_size * opt->n_threads, &n, ks, ks2)) != 0) {
|
|
|
749 int size = 0;
|
|
|
750 if (n > _seq->max) {
|
|
|
751 _seq->max = n;
|
|
|
752 kroundup32(_seq->max);
|
|
|
753 _seq->seq = realloc(_seq->seq, _seq->max * sizeof(bsw2seq1_t));
|
|
|
754 }
|
|
|
755 _seq->n = n;
|
|
|
756 for (i = 0; i < n; ++i) {
|
|
|
757 bseq1_t *b = &bseq[i];
|
|
|
758 bsw2seq1_t *p = &_seq->seq[i];
|
|
|
759 p->tid = -1; p->l = b->l_seq;
|
|
|
760 p->name = b->name; p->seq = b->seq; p->qual = b->qual; p->comment = b->comment; p->sam = 0;
|
|
|
761 size += p->l;
|
|
|
762 }
|
|
|
763 fprintf(stderr, "[bsw2_aln] read %d sequences/pairs (%d bp) ...\n", n, size);
|
|
|
764 free(bseq);
|
|
|
765 process_seqs(_seq, opt, bns, pac, target, is_pe);
|
|
|
766 }
|
|
|
767 // free
|
|
|
768 free(pac);
|
|
|
769 free(_seq->seq); free(_seq);
|
|
|
770 kseq_destroy(ks);
|
|
|
771 err_gzclose(fp);
|
|
|
772 if (fn2) {
|
|
|
773 kseq_destroy(ks2);
|
|
|
774 err_gzclose(fp2);
|
|
|
775 }
|
|
|
776 }
|
