0
|
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 }
|