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