Mercurial > repos > ashvark > qiime_1_8_0
comparison bwa-0.6.2/bwase.c @ 0:dd1186b11b3b draft
Uploaded BWA
author | ashvark |
---|---|
date | Fri, 18 Jul 2014 07:55:14 -0400 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:dd1186b11b3b |
---|---|
1 #include <unistd.h> | |
2 #include <string.h> | |
3 #include <stdio.h> | |
4 #include <stdlib.h> | |
5 #include <math.h> | |
6 #include <time.h> | |
7 #include "stdaln.h" | |
8 #include "bwase.h" | |
9 #include "bwtaln.h" | |
10 #include "bntseq.h" | |
11 #include "utils.h" | |
12 #include "kstring.h" | |
13 | |
14 int g_log_n[256]; | |
15 char *bwa_rg_line, *bwa_rg_id; | |
16 | |
17 void bwa_print_sam_PG(); | |
18 | |
19 void bwa_aln2seq_core(int n_aln, const bwt_aln1_t *aln, bwa_seq_t *s, int set_main, int n_multi) | |
20 { | |
21 int i, cnt, best; | |
22 if (n_aln == 0) { | |
23 s->type = BWA_TYPE_NO_MATCH; | |
24 s->c1 = s->c2 = 0; | |
25 return; | |
26 } | |
27 | |
28 if (set_main) { | |
29 best = aln[0].score; | |
30 for (i = cnt = 0; i < n_aln; ++i) { | |
31 const bwt_aln1_t *p = aln + i; | |
32 if (p->score > best) break; | |
33 if (drand48() * (p->l - p->k + 1 + cnt) > (double)cnt) { | |
34 s->n_mm = p->n_mm; s->n_gapo = p->n_gapo; s->n_gape = p->n_gape; | |
35 s->score = p->score; | |
36 s->sa = p->k + (bwtint_t)((p->l - p->k + 1) * drand48()); | |
37 } | |
38 cnt += p->l - p->k + 1; | |
39 } | |
40 s->c1 = cnt; | |
41 for (; i < n_aln; ++i) cnt += aln[i].l - aln[i].k + 1; | |
42 s->c2 = cnt - s->c1; | |
43 s->type = s->c1 > 1? BWA_TYPE_REPEAT : BWA_TYPE_UNIQUE; | |
44 } | |
45 | |
46 if (n_multi) { | |
47 int k, rest, n_occ, z = 0; | |
48 for (k = n_occ = 0; k < n_aln; ++k) { | |
49 const bwt_aln1_t *q = aln + k; | |
50 n_occ += q->l - q->k + 1; | |
51 } | |
52 if (s->multi) free(s->multi); | |
53 if (n_occ > n_multi + 1) { // if there are too many hits, generate none of them | |
54 s->multi = 0; s->n_multi = 0; | |
55 return; | |
56 } | |
57 /* The following code is more flexible than what is required | |
58 * here. In principle, due to the requirement above, we can | |
59 * simply output all hits, but the following samples "rest" | |
60 * number of random hits. */ | |
61 rest = n_occ > n_multi + 1? n_multi + 1 : n_occ; // find one additional for ->sa | |
62 s->multi = calloc(rest, sizeof(bwt_multi1_t)); | |
63 for (k = 0; k < n_aln; ++k) { | |
64 const bwt_aln1_t *q = aln + k; | |
65 if (q->l - q->k + 1 <= rest) { | |
66 bwtint_t l; | |
67 for (l = q->k; l <= q->l; ++l) { | |
68 s->multi[z].pos = l; | |
69 s->multi[z].gap = q->n_gapo + q->n_gape; | |
70 s->multi[z++].mm = q->n_mm; | |
71 } | |
72 rest -= q->l - q->k + 1; | |
73 } else { // Random sampling (http://code.activestate.com/recipes/272884/). In fact, we never come here. | |
74 int j, i, k; | |
75 for (j = rest, i = q->l - q->k + 1, k = 0; j > 0; --j) { | |
76 double p = 1.0, x = drand48(); | |
77 while (x < p) p -= p * j / (i--); | |
78 s->multi[z].pos = q->l - i; | |
79 s->multi[z].gap = q->n_gapo + q->n_gape; | |
80 s->multi[z++].mm = q->n_mm; | |
81 } | |
82 rest = 0; | |
83 break; | |
84 } | |
85 } | |
86 s->n_multi = z; | |
87 } | |
88 } | |
89 | |
90 void bwa_aln2seq(int n_aln, const bwt_aln1_t *aln, bwa_seq_t *s) | |
91 { | |
92 bwa_aln2seq_core(n_aln, aln, s, 1, 0); | |
93 } | |
94 | |
95 int bwa_approx_mapQ(const bwa_seq_t *p, int mm) | |
96 { | |
97 int n; | |
98 if (p->c1 == 0) return 23; | |
99 if (p->c1 > 1) return 0; | |
100 if (p->n_mm == mm) return 25; | |
101 if (p->c2 == 0) return 37; | |
102 n = (p->c2 >= 255)? 255 : p->c2; | |
103 return (23 < g_log_n[n])? 0 : 23 - g_log_n[n]; | |
104 } | |
105 | |
106 bwtint_t bwa_sa2pos(const bntseq_t *bns, const bwt_t *bwt, bwtint_t sapos, int len, int *strand) | |
107 { | |
108 bwtint_t pos_f; | |
109 int is_rev; | |
110 pos_f = bns_depos(bns, bwt_sa(bwt, sapos), &is_rev); // pos_f | |
111 *strand = !is_rev; | |
112 /* NB: For gapped alignment, pacpos may not be correct, which will be fixed | |
113 * in bwa_refine_gapped_core(). This line also determines the way "x" is | |
114 * calculated in bwa_refine_gapped_core() when (ext < 0 && is_end == 0). */ | |
115 if (is_rev) pos_f = pos_f + 1 < len? 0 : pos_f - len + 1; // mapped to the forward strand | |
116 return pos_f; // FIXME: it is possible that pos_f < bns->anns[ref_id].offset | |
117 } | |
118 | |
119 /** | |
120 * Derive the actual position in the read from the given suffix array | |
121 * coordinates. Note that the position will be approximate based on | |
122 * whether indels appear in the read and whether calculations are | |
123 * performed from the start or end of the read. | |
124 */ | |
125 void bwa_cal_pac_pos_core(const bntseq_t *bns, const bwt_t *bwt, bwa_seq_t *seq, const int max_mm, const float fnr) | |
126 { | |
127 int max_diff, strand; | |
128 if (seq->type != BWA_TYPE_UNIQUE && seq->type != BWA_TYPE_REPEAT) return; | |
129 max_diff = fnr > 0.0? bwa_cal_maxdiff(seq->len, BWA_AVG_ERR, fnr) : max_mm; | |
130 seq->seQ = seq->mapQ = bwa_approx_mapQ(seq, max_diff); | |
131 seq->pos = bwa_sa2pos(bns, bwt, seq->sa, seq->len, &strand); | |
132 seq->strand = strand; | |
133 seq->seQ = seq->mapQ = bwa_approx_mapQ(seq, max_diff); | |
134 } | |
135 | |
136 void bwa_cal_pac_pos(const bntseq_t *bns, const char *prefix, int n_seqs, bwa_seq_t *seqs, int max_mm, float fnr) | |
137 { | |
138 int i, j, strand, n_multi; | |
139 char str[1024]; | |
140 bwt_t *bwt; | |
141 // load forward SA | |
142 strcpy(str, prefix); strcat(str, ".bwt"); bwt = bwt_restore_bwt(str); | |
143 strcpy(str, prefix); strcat(str, ".sa"); bwt_restore_sa(str, bwt); | |
144 for (i = 0; i != n_seqs; ++i) { | |
145 bwa_seq_t *p = &seqs[i]; | |
146 bwa_cal_pac_pos_core(bns, bwt, p, max_mm, fnr); | |
147 for (j = n_multi = 0; j < p->n_multi; ++j) { | |
148 bwt_multi1_t *q = p->multi + j; | |
149 q->pos = bwa_sa2pos(bns, bwt, q->pos, p->len, &strand); | |
150 q->strand = strand; | |
151 if (q->pos != p->pos) | |
152 p->multi[n_multi++] = *q; | |
153 } | |
154 p->n_multi = n_multi; | |
155 } | |
156 bwt_destroy(bwt); | |
157 } | |
158 | |
159 /* is_end_correct == 1 if (*pos+len) gives the correct coordinate on | |
160 * forward strand. This happens when p->pos is calculated by | |
161 * bwa_cal_pac_pos(). is_end_correct==0 if (*pos) gives the correct | |
162 * coordinate. This happens only for color-converted alignment. */ | |
163 bwa_cigar_t *bwa_refine_gapped_core(bwtint_t l_pac, const ubyte_t *pacseq, int len, const ubyte_t *seq, bwtint_t *_pos, | |
164 int ext, int *n_cigar, int is_end_correct) | |
165 { | |
166 bwa_cigar_t *cigar = 0; | |
167 ubyte_t *ref_seq; | |
168 int l = 0, path_len, ref_len; | |
169 AlnParam ap = aln_param_bwa; | |
170 path_t *path; | |
171 int64_t k, __pos = *_pos; | |
172 | |
173 ref_len = len + abs(ext); | |
174 if (ext > 0) { | |
175 ref_seq = (ubyte_t*)calloc(ref_len, 1); | |
176 for (k = __pos; k < __pos + ref_len && k < l_pac; ++k) | |
177 ref_seq[l++] = pacseq[k>>2] >> ((~k&3)<<1) & 3; | |
178 } else { | |
179 int64_t x = __pos + (is_end_correct? len : ref_len); | |
180 ref_seq = (ubyte_t*)calloc(ref_len, 1); | |
181 for (l = 0, k = x - ref_len > 0? x - ref_len : 0; k < x && k < l_pac; ++k) | |
182 ref_seq[l++] = pacseq[k>>2] >> ((~k&3)<<1) & 3; | |
183 } | |
184 path = (path_t*)calloc(l+len, sizeof(path_t)); | |
185 | |
186 aln_global_core(ref_seq, l, (ubyte_t*)seq, len, &ap, path, &path_len); | |
187 cigar = bwa_aln_path2cigar(path, path_len, n_cigar); | |
188 | |
189 if (ext < 0 && is_end_correct) { // fix coordinate for reads mapped to the forward strand | |
190 for (l = k = 0; k < *n_cigar; ++k) { | |
191 if (__cigar_op(cigar[k]) == FROM_D) l -= __cigar_len(cigar[k]); | |
192 else if (__cigar_op(cigar[k]) == FROM_I) l += __cigar_len(cigar[k]); | |
193 } | |
194 __pos += l; | |
195 } | |
196 | |
197 if (__cigar_op(cigar[0]) == FROM_D) { // deletion at the 5'-end | |
198 __pos += __cigar_len(cigar[0]); | |
199 for (k = 0; k < *n_cigar - 1; ++k) cigar[k] = cigar[k+1]; | |
200 --(*n_cigar); | |
201 } | |
202 if (__cigar_op(cigar[*n_cigar-1]) == FROM_D) --(*n_cigar); // deletion at the 3'-end | |
203 | |
204 // change "I" at either end of the read to S. just in case. This should rarely happen... | |
205 if (__cigar_op(cigar[*n_cigar-1]) == FROM_I) cigar[*n_cigar-1] = __cigar_create(3, (__cigar_len(cigar[*n_cigar-1]))); | |
206 if (__cigar_op(cigar[0]) == FROM_I) cigar[0] = __cigar_create(3, (__cigar_len(cigar[0]))); | |
207 | |
208 *_pos = (bwtint_t)__pos; | |
209 free(ref_seq); free(path); | |
210 return cigar; | |
211 } | |
212 | |
213 char *bwa_cal_md1(int n_cigar, bwa_cigar_t *cigar, int len, bwtint_t pos, ubyte_t *seq, | |
214 bwtint_t l_pac, ubyte_t *pacseq, kstring_t *str, int *_nm) | |
215 { | |
216 bwtint_t x, y; | |
217 int z, u, c, nm = 0; | |
218 str->l = 0; // reset | |
219 x = pos; y = 0; | |
220 if (cigar) { | |
221 int k, l; | |
222 for (k = u = 0; k < n_cigar; ++k) { | |
223 l = __cigar_len(cigar[k]); | |
224 if (__cigar_op(cigar[k]) == FROM_M) { | |
225 for (z = 0; z < l && x+z < l_pac; ++z) { | |
226 c = pacseq[(x+z)>>2] >> ((~(x+z)&3)<<1) & 3; | |
227 if (c > 3 || seq[y+z] > 3 || c != seq[y+z]) { | |
228 ksprintf(str, "%d", u); | |
229 kputc("ACGTN"[c], str); | |
230 ++nm; | |
231 u = 0; | |
232 } else ++u; | |
233 } | |
234 x += l; y += l; | |
235 } else if (__cigar_op(cigar[k]) == FROM_I || __cigar_op(cigar[k]) == FROM_S) { | |
236 y += l; | |
237 if (__cigar_op(cigar[k]) == FROM_I) nm += l; | |
238 } else if (__cigar_op(cigar[k]) == FROM_D) { | |
239 ksprintf(str, "%d", u); | |
240 kputc('^', str); | |
241 for (z = 0; z < l && x+z < l_pac; ++z) | |
242 kputc("ACGT"[pacseq[(x+z)>>2] >> ((~(x+z)&3)<<1) & 3], str); | |
243 u = 0; | |
244 x += l; nm += l; | |
245 } | |
246 } | |
247 } else { // no gaps | |
248 for (z = u = 0; z < (bwtint_t)len && x+z < l_pac; ++z) { | |
249 c = pacseq[(x+z)>>2] >> ((~(x+z)&3)<<1) & 3; | |
250 if (c > 3 || seq[y+z] > 3 || c != seq[y+z]) { | |
251 ksprintf(str, "%d", u); | |
252 kputc("ACGTN"[c], str); | |
253 ++nm; | |
254 u = 0; | |
255 } else ++u; | |
256 } | |
257 } | |
258 ksprintf(str, "%d", u); | |
259 *_nm = nm; | |
260 return strdup(str->s); | |
261 } | |
262 | |
263 void bwa_correct_trimmed(bwa_seq_t *s) | |
264 { | |
265 if (s->len == s->full_len) return; | |
266 if (s->strand == 0) { // forward | |
267 if (s->cigar && __cigar_op(s->cigar[s->n_cigar-1]) == FROM_S) { // the last is S | |
268 s->cigar[s->n_cigar-1] += s->full_len - s->len; | |
269 } else { | |
270 if (s->cigar == 0) { | |
271 s->n_cigar = 2; | |
272 s->cigar = calloc(s->n_cigar, sizeof(bwa_cigar_t)); | |
273 s->cigar[0] = __cigar_create(0, s->len); | |
274 } else { | |
275 ++s->n_cigar; | |
276 s->cigar = realloc(s->cigar, s->n_cigar * sizeof(bwa_cigar_t)); | |
277 } | |
278 s->cigar[s->n_cigar-1] = __cigar_create(3, (s->full_len - s->len)); | |
279 } | |
280 } else { // reverse | |
281 if (s->cigar && __cigar_op(s->cigar[0]) == FROM_S) { // the first is S | |
282 s->cigar[0] += s->full_len - s->len; | |
283 } else { | |
284 if (s->cigar == 0) { | |
285 s->n_cigar = 2; | |
286 s->cigar = calloc(s->n_cigar, sizeof(bwa_cigar_t)); | |
287 s->cigar[1] = __cigar_create(0, s->len); | |
288 } else { | |
289 ++s->n_cigar; | |
290 s->cigar = realloc(s->cigar, s->n_cigar * sizeof(bwa_cigar_t)); | |
291 memmove(s->cigar + 1, s->cigar, (s->n_cigar-1) * sizeof(bwa_cigar_t)); | |
292 } | |
293 s->cigar[0] = __cigar_create(3, (s->full_len - s->len)); | |
294 } | |
295 } | |
296 s->len = s->full_len; | |
297 } | |
298 | |
299 void bwa_refine_gapped(const bntseq_t *bns, int n_seqs, bwa_seq_t *seqs, ubyte_t *_pacseq, bntseq_t *ntbns) | |
300 { | |
301 ubyte_t *pacseq, *ntpac = 0; | |
302 int i, j; | |
303 kstring_t *str; | |
304 | |
305 if (ntbns) { // in color space | |
306 ntpac = (ubyte_t*)calloc(ntbns->l_pac/4+1, 1); | |
307 rewind(ntbns->fp_pac); | |
308 fread(ntpac, 1, ntbns->l_pac/4 + 1, ntbns->fp_pac); | |
309 } | |
310 | |
311 if (!_pacseq) { | |
312 pacseq = (ubyte_t*)calloc(bns->l_pac/4+1, 1); | |
313 rewind(bns->fp_pac); | |
314 fread(pacseq, 1, bns->l_pac/4+1, bns->fp_pac); | |
315 } else pacseq = _pacseq; | |
316 for (i = 0; i != n_seqs; ++i) { | |
317 bwa_seq_t *s = seqs + i; | |
318 seq_reverse(s->len, s->seq, 0); // IMPORTANT: s->seq is reversed here!!! | |
319 for (j = 0; j < s->n_multi; ++j) { | |
320 bwt_multi1_t *q = s->multi + j; | |
321 int n_cigar; | |
322 if (q->gap == 0) continue; | |
323 q->cigar = bwa_refine_gapped_core(bns->l_pac, pacseq, s->len, q->strand? s->rseq : s->seq, &q->pos, | |
324 (q->strand? 1 : -1) * q->gap, &n_cigar, 1); | |
325 q->n_cigar = n_cigar; | |
326 } | |
327 if (s->type == BWA_TYPE_NO_MATCH || s->type == BWA_TYPE_MATESW || s->n_gapo == 0) continue; | |
328 s->cigar = bwa_refine_gapped_core(bns->l_pac, pacseq, s->len, s->strand? s->rseq : s->seq, &s->pos, | |
329 (s->strand? 1 : -1) * (s->n_gapo + s->n_gape), &s->n_cigar, 1); | |
330 } | |
331 #if 0 | |
332 if (ntbns) { // in color space | |
333 for (i = 0; i < n_seqs; ++i) { | |
334 bwa_seq_t *s = seqs + i; | |
335 bwa_cs2nt_core(s, bns->l_pac, ntpac); | |
336 for (j = 0; j < s->n_multi; ++j) { | |
337 bwt_multi1_t *q = s->multi + j; | |
338 int n_cigar; | |
339 if (q->gap == 0) continue; | |
340 free(q->cigar); | |
341 q->cigar = bwa_refine_gapped_core(bns->l_pac, ntpac, s->len, q->strand? s->rseq : s->seq, &q->pos, | |
342 (q->strand? 1 : -1) * q->gap, &n_cigar, 0); | |
343 q->n_cigar = n_cigar; | |
344 } | |
345 if (s->type != BWA_TYPE_NO_MATCH && s->cigar) { // update cigar again | |
346 free(s->cigar); | |
347 s->cigar = bwa_refine_gapped_core(bns->l_pac, ntpac, s->len, s->strand? s->rseq : s->seq, &s->pos, | |
348 (s->strand? 1 : -1) * (s->n_gapo + s->n_gape), &s->n_cigar, 0); | |
349 } | |
350 } | |
351 } | |
352 #endif | |
353 // generate MD tag | |
354 str = (kstring_t*)calloc(1, sizeof(kstring_t)); | |
355 for (i = 0; i != n_seqs; ++i) { | |
356 bwa_seq_t *s = seqs + i; | |
357 if (s->type != BWA_TYPE_NO_MATCH) { | |
358 int nm; | |
359 s->md = bwa_cal_md1(s->n_cigar, s->cigar, s->len, s->pos, s->strand? s->rseq : s->seq, | |
360 bns->l_pac, ntbns? ntpac : pacseq, str, &nm); | |
361 s->nm = nm; | |
362 } | |
363 } | |
364 free(str->s); free(str); | |
365 | |
366 // correct for trimmed reads | |
367 if (!ntbns) // trimming is only enabled for Illumina reads | |
368 for (i = 0; i < n_seqs; ++i) bwa_correct_trimmed(seqs + i); | |
369 | |
370 if (!_pacseq) free(pacseq); | |
371 free(ntpac); | |
372 } | |
373 | |
374 int64_t pos_end(const bwa_seq_t *p) | |
375 { | |
376 if (p->cigar) { | |
377 int j; | |
378 int64_t x = p->pos; | |
379 for (j = 0; j != p->n_cigar; ++j) { | |
380 int op = __cigar_op(p->cigar[j]); | |
381 if (op == 0 || op == 2) x += __cigar_len(p->cigar[j]); | |
382 } | |
383 return x; | |
384 } else return p->pos + p->len; | |
385 } | |
386 | |
387 int64_t pos_end_multi(const bwt_multi1_t *p, int len) // analogy to pos_end() | |
388 { | |
389 if (p->cigar) { | |
390 int j; | |
391 int64_t x = p->pos; | |
392 for (j = 0; j != p->n_cigar; ++j) { | |
393 int op = __cigar_op(p->cigar[j]); | |
394 if (op == 0 || op == 2) x += __cigar_len(p->cigar[j]); | |
395 } | |
396 return x; | |
397 } else return p->pos + len; | |
398 } | |
399 | |
400 static int64_t pos_5(const bwa_seq_t *p) | |
401 { | |
402 if (p->type != BWA_TYPE_NO_MATCH) | |
403 return p->strand? pos_end(p) : p->pos; | |
404 return -1; | |
405 } | |
406 | |
407 void bwa_print_sam1(const bntseq_t *bns, bwa_seq_t *p, const bwa_seq_t *mate, int mode, int max_top2) | |
408 { | |
409 int j; | |
410 if (p->type != BWA_TYPE_NO_MATCH || (mate && mate->type != BWA_TYPE_NO_MATCH)) { | |
411 int seqid, nn, am = 0, flag = p->extra_flag; | |
412 char XT; | |
413 | |
414 if (p->type == BWA_TYPE_NO_MATCH) { | |
415 p->pos = mate->pos; | |
416 p->strand = mate->strand; | |
417 flag |= SAM_FSU; | |
418 j = 1; | |
419 } else j = pos_end(p) - p->pos; // j is the length of the reference in the alignment | |
420 | |
421 // get seqid | |
422 nn = bns_cnt_ambi(bns, p->pos, j, &seqid); | |
423 if (p->type != BWA_TYPE_NO_MATCH && p->pos + j - bns->anns[seqid].offset > bns->anns[seqid].len) | |
424 flag |= SAM_FSU; // flag UNMAP as this alignment bridges two adjacent reference sequences | |
425 | |
426 // update flag and print it | |
427 if (p->strand) flag |= SAM_FSR; | |
428 if (mate) { | |
429 if (mate->type != BWA_TYPE_NO_MATCH) { | |
430 if (mate->strand) flag |= SAM_FMR; | |
431 } else flag |= SAM_FMU; | |
432 } | |
433 err_printf("%s\t%d\t%s\t", p->name, flag, bns->anns[seqid].name); | |
434 err_printf("%d\t%d\t", (int)(p->pos - bns->anns[seqid].offset + 1), p->mapQ); | |
435 | |
436 // print CIGAR | |
437 if (p->cigar) { | |
438 for (j = 0; j != p->n_cigar; ++j) | |
439 err_printf("%d%c", __cigar_len(p->cigar[j]), "MIDS"[__cigar_op(p->cigar[j])]); | |
440 } else if (p->type == BWA_TYPE_NO_MATCH) err_printf("*"); | |
441 else err_printf("%dM", p->len); | |
442 | |
443 // print mate coordinate | |
444 if (mate && mate->type != BWA_TYPE_NO_MATCH) { | |
445 int m_seqid, m_is_N; | |
446 long long isize; | |
447 am = mate->seQ < p->seQ? mate->seQ : p->seQ; // smaller single-end mapping quality | |
448 // redundant calculation here, but should not matter too much | |
449 m_is_N = bns_cnt_ambi(bns, mate->pos, mate->len, &m_seqid); | |
450 err_printf("\t%s\t", (seqid == m_seqid)? "=" : bns->anns[m_seqid].name); | |
451 isize = (seqid == m_seqid)? pos_5(mate) - pos_5(p) : 0; | |
452 if (p->type == BWA_TYPE_NO_MATCH) isize = 0; | |
453 err_printf("%d\t%lld\t", (int)(mate->pos - bns->anns[m_seqid].offset + 1), isize); | |
454 } else if (mate) err_printf("\t=\t%d\t0\t", (int)(p->pos - bns->anns[seqid].offset + 1)); | |
455 else err_printf("\t*\t0\t0\t"); | |
456 | |
457 // print sequence and quality | |
458 if (p->strand == 0) | |
459 for (j = 0; j != p->full_len; ++j) putchar("ACGTN"[(int)p->seq[j]]); | |
460 else for (j = 0; j != p->full_len; ++j) putchar("TGCAN"[p->seq[p->full_len - 1 - j]]); | |
461 putchar('\t'); | |
462 if (p->qual) { | |
463 if (p->strand) seq_reverse(p->len, p->qual, 0); // reverse quality | |
464 err_printf("%s", p->qual); | |
465 } else err_printf("*"); | |
466 | |
467 if (bwa_rg_id) err_printf("\tRG:Z:%s", bwa_rg_id); | |
468 if (p->bc[0]) err_printf("\tBC:Z:%s", p->bc); | |
469 if (p->clip_len < p->full_len) err_printf("\tXC:i:%d", p->clip_len); | |
470 if (p->type != BWA_TYPE_NO_MATCH) { | |
471 int i; | |
472 // calculate XT tag | |
473 XT = "NURM"[p->type]; | |
474 if (nn > 10) XT = 'N'; | |
475 // print tags | |
476 err_printf("\tXT:A:%c\t%s:i:%d", XT, (mode & BWA_MODE_COMPREAD)? "NM" : "CM", p->nm); | |
477 if (nn) err_printf("\tXN:i:%d", nn); | |
478 if (mate) err_printf("\tSM:i:%d\tAM:i:%d", p->seQ, am); | |
479 if (p->type != BWA_TYPE_MATESW) { // X0 and X1 are not available for this type of alignment | |
480 err_printf("\tX0:i:%d", p->c1); | |
481 if (p->c1 <= max_top2) err_printf("\tX1:i:%d", p->c2); | |
482 } | |
483 err_printf("\tXM:i:%d\tXO:i:%d\tXG:i:%d", p->n_mm, p->n_gapo, p->n_gapo+p->n_gape); | |
484 if (p->md) err_printf("\tMD:Z:%s", p->md); | |
485 // print multiple hits | |
486 if (p->n_multi) { | |
487 err_printf("\tXA:Z:"); | |
488 for (i = 0; i < p->n_multi; ++i) { | |
489 bwt_multi1_t *q = p->multi + i; | |
490 int k; | |
491 j = pos_end_multi(q, p->len) - q->pos; | |
492 nn = bns_cnt_ambi(bns, q->pos, j, &seqid); | |
493 err_printf("%s,%c%d,", bns->anns[seqid].name, q->strand? '-' : '+', | |
494 (int)(q->pos - bns->anns[seqid].offset + 1)); | |
495 if (q->cigar) { | |
496 for (k = 0; k < q->n_cigar; ++k) | |
497 err_printf("%d%c", __cigar_len(q->cigar[k]), "MIDS"[__cigar_op(q->cigar[k])]); | |
498 } else err_printf("%dM", p->len); | |
499 err_printf(",%d;", q->gap + q->mm); | |
500 } | |
501 } | |
502 } | |
503 putchar('\n'); | |
504 } else { // this read has no match | |
505 ubyte_t *s = p->strand? p->rseq : p->seq; | |
506 int flag = p->extra_flag | SAM_FSU; | |
507 if (mate && mate->type == BWA_TYPE_NO_MATCH) flag |= SAM_FMU; | |
508 err_printf("%s\t%d\t*\t0\t0\t*\t*\t0\t0\t", p->name, flag); | |
509 for (j = 0; j != p->len; ++j) putchar("ACGTN"[(int)s[j]]); | |
510 putchar('\t'); | |
511 if (p->qual) { | |
512 if (p->strand) seq_reverse(p->len, p->qual, 0); // reverse quality | |
513 err_printf("%s", p->qual); | |
514 } else err_printf("*"); | |
515 if (bwa_rg_id) err_printf("\tRG:Z:%s", bwa_rg_id); | |
516 if (p->bc[0]) err_printf("\tBC:Z:%s", p->bc); | |
517 if (p->clip_len < p->full_len) err_printf("\tXC:i:%d", p->clip_len); | |
518 putchar('\n'); | |
519 } | |
520 } | |
521 | |
522 bntseq_t *bwa_open_nt(const char *prefix) | |
523 { | |
524 bntseq_t *ntbns; | |
525 char *str; | |
526 str = (char*)calloc(strlen(prefix) + 10, 1); | |
527 strcat(strcpy(str, prefix), ".nt"); | |
528 ntbns = bns_restore(str); | |
529 free(str); | |
530 return ntbns; | |
531 } | |
532 | |
533 void bwa_print_sam_SQ(const bntseq_t *bns) | |
534 { | |
535 int i; | |
536 for (i = 0; i < bns->n_seqs; ++i) | |
537 err_printf("@SQ\tSN:%s\tLN:%d\n", bns->anns[i].name, bns->anns[i].len); | |
538 if (bwa_rg_line) err_printf("%s\n", bwa_rg_line); | |
539 } | |
540 | |
541 void bwase_initialize() | |
542 { | |
543 int i; | |
544 for (i = 1; i != 256; ++i) g_log_n[i] = (int)(4.343 * log(i) + 0.5); | |
545 } | |
546 | |
547 char *bwa_escape(char *s) | |
548 { | |
549 char *p, *q; | |
550 for (p = q = s; *p; ++p) { | |
551 if (*p == '\\') { | |
552 ++p; | |
553 if (*p == 't') *q++ = '\t'; | |
554 else if (*p == 'n') *q++ = '\n'; | |
555 else if (*p == 'r') *q++ = '\r'; | |
556 else if (*p == '\\') *q++ = '\\'; | |
557 } else *q++ = *p; | |
558 } | |
559 *q = '\0'; | |
560 return s; | |
561 } | |
562 | |
563 int bwa_set_rg(const char *s) | |
564 { | |
565 char *p, *q, *r; | |
566 if (strstr(s, "@RG") != s) return -1; | |
567 if (bwa_rg_line) free(bwa_rg_line); | |
568 if (bwa_rg_id) free(bwa_rg_id); | |
569 bwa_rg_line = strdup(s); | |
570 bwa_rg_id = 0; | |
571 bwa_escape(bwa_rg_line); | |
572 p = strstr(bwa_rg_line, "\tID:"); | |
573 if (p == 0) return -1; | |
574 p += 4; | |
575 for (q = p; *q && *q != '\t' && *q != '\n'; ++q); | |
576 bwa_rg_id = calloc(q - p + 1, 1); | |
577 for (q = p, r = bwa_rg_id; *q && *q != '\t' && *q != '\n'; ++q) | |
578 *r++ = *q; | |
579 return 0; | |
580 } | |
581 | |
582 void bwa_sai2sam_se_core(const char *prefix, const char *fn_sa, const char *fn_fa, int n_occ) | |
583 { | |
584 extern bwa_seqio_t *bwa_open_reads(int mode, const char *fn_fa); | |
585 int i, n_seqs, tot_seqs = 0, m_aln; | |
586 bwt_aln1_t *aln = 0; | |
587 bwa_seq_t *seqs; | |
588 bwa_seqio_t *ks; | |
589 clock_t t; | |
590 bntseq_t *bns, *ntbns = 0; | |
591 FILE *fp_sa; | |
592 gap_opt_t opt; | |
593 | |
594 // initialization | |
595 bwase_initialize(); | |
596 bns = bns_restore(prefix); | |
597 srand48(bns->seed); | |
598 fp_sa = xopen(fn_sa, "r"); | |
599 | |
600 m_aln = 0; | |
601 fread(&opt, sizeof(gap_opt_t), 1, fp_sa); | |
602 if (!(opt.mode & BWA_MODE_COMPREAD)) // in color space; initialize ntpac | |
603 ntbns = bwa_open_nt(prefix); | |
604 bwa_print_sam_SQ(bns); | |
605 //bwa_print_sam_PG(); | |
606 // set ks | |
607 ks = bwa_open_reads(opt.mode, fn_fa); | |
608 // core loop | |
609 while ((seqs = bwa_read_seq(ks, 0x40000, &n_seqs, opt.mode, opt.trim_qual)) != 0) { | |
610 tot_seqs += n_seqs; | |
611 t = clock(); | |
612 | |
613 // read alignment | |
614 for (i = 0; i < n_seqs; ++i) { | |
615 bwa_seq_t *p = seqs + i; | |
616 int n_aln; | |
617 fread(&n_aln, 4, 1, fp_sa); | |
618 if (n_aln > m_aln) { | |
619 m_aln = n_aln; | |
620 aln = (bwt_aln1_t*)realloc(aln, sizeof(bwt_aln1_t) * m_aln); | |
621 } | |
622 fread(aln, sizeof(bwt_aln1_t), n_aln, fp_sa); | |
623 bwa_aln2seq_core(n_aln, aln, p, 1, n_occ); | |
624 } | |
625 | |
626 fprintf(stderr, "[bwa_aln_core] convert to sequence coordinate... "); | |
627 bwa_cal_pac_pos(bns, prefix, n_seqs, seqs, opt.max_diff, opt.fnr); // forward bwt will be destroyed here | |
628 fprintf(stderr, "%.2f sec\n", (float)(clock() - t) / CLOCKS_PER_SEC); t = clock(); | |
629 | |
630 fprintf(stderr, "[bwa_aln_core] refine gapped alignments... "); | |
631 bwa_refine_gapped(bns, n_seqs, seqs, 0, ntbns); | |
632 fprintf(stderr, "%.2f sec\n", (float)(clock() - t) / CLOCKS_PER_SEC); t = clock(); | |
633 | |
634 fprintf(stderr, "[bwa_aln_core] print alignments... "); | |
635 for (i = 0; i < n_seqs; ++i) | |
636 bwa_print_sam1(bns, seqs + i, 0, opt.mode, opt.max_top2); | |
637 fprintf(stderr, "%.2f sec\n", (float)(clock() - t) / CLOCKS_PER_SEC); t = clock(); | |
638 | |
639 bwa_free_read_seq(n_seqs, seqs); | |
640 fprintf(stderr, "[bwa_aln_core] %d sequences have been processed.\n", tot_seqs); | |
641 } | |
642 | |
643 // destroy | |
644 bwa_seq_close(ks); | |
645 if (ntbns) bns_destroy(ntbns); | |
646 bns_destroy(bns); | |
647 fclose(fp_sa); | |
648 free(aln); | |
649 } | |
650 | |
651 int bwa_sai2sam_se(int argc, char *argv[]) | |
652 { | |
653 extern char *bwa_infer_prefix(const char *hint); | |
654 int c, n_occ = 3; | |
655 char *prefix; | |
656 while ((c = getopt(argc, argv, "hn:f:r:")) >= 0) { | |
657 switch (c) { | |
658 case 'h': break; | |
659 case 'r': | |
660 if (bwa_set_rg(optarg) < 0) { | |
661 fprintf(stderr, "[%s] malformated @RG line\n", __func__); | |
662 return 1; | |
663 } | |
664 break; | |
665 case 'n': n_occ = atoi(optarg); break; | |
666 case 'f': xreopen(optarg, "w", stdout); break; | |
667 default: return 1; | |
668 } | |
669 } | |
670 | |
671 if (optind + 3 > argc) { | |
672 fprintf(stderr, "Usage: bwa samse [-n max_occ] [-f out.sam] [-r RG_line] <prefix> <in.sai> <in.fq>\n"); | |
673 return 1; | |
674 } | |
675 if ((prefix = bwa_infer_prefix(argv[optind])) == 0) { | |
676 fprintf(stderr, "[%s] fail to locate the index\n", __func__); | |
677 free(bwa_rg_line); free(bwa_rg_id); | |
678 return 0; | |
679 } | |
680 bwa_sai2sam_se_core(prefix, argv[optind+1], argv[optind+2], n_occ); | |
681 free(bwa_rg_line); free(bwa_rg_id); | |
682 return 0; | |
683 } |