0
|
1 #include <stdlib.h>
|
|
2 #include <string.h>
|
|
3 #include <stdio.h>
|
|
4 #include <math.h>
|
|
5 #include "bwa.h"
|
|
6 #include "bwt.h"
|
|
7 #include "bwtgap.h"
|
|
8 #include "bntseq.h"
|
|
9
|
|
10 #ifndef kroundup32
|
|
11 #define kroundup32(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, ++(x))
|
|
12 #endif
|
|
13
|
|
14 extern unsigned char nst_nt4_table[256];
|
|
15 extern void seq_reverse(int len, uint8_t *seq, int is_comp);
|
|
16
|
|
17 bwa_opt_t bwa_def_opt = { 11, 4, -1, 1, 6, 32, 2, 0.04 };
|
|
18
|
|
19 struct bwa_idx_t {
|
|
20 bwt_t *bwt;
|
|
21 bntseq_t *bns;
|
|
22 uint8_t *pac;
|
|
23 };
|
|
24
|
|
25 struct bwa_buf_t {
|
|
26 int max_buf;
|
|
27 bwa_pestat_t pes;
|
|
28 gap_stack_t *stack;
|
|
29 gap_opt_t *opt;
|
|
30 int *diff_tab;
|
|
31 uint8_t *buf;
|
|
32 int *logn;
|
|
33 };
|
|
34
|
|
35 bwa_idx_t *bwa_idx_load(const char *prefix)
|
|
36 {
|
|
37 bwa_idx_t *p;
|
|
38 int l;
|
|
39 char *str;
|
|
40 l = strlen(prefix);
|
|
41 p = calloc(1, sizeof(bwa_idx_t));
|
|
42 str = malloc(l + 10);
|
|
43 strcpy(str, prefix);
|
|
44 p->bns = bns_restore(str);
|
|
45 strcpy(str + l, ".bwt");
|
|
46 p->bwt = bwt_restore_bwt(str);
|
|
47 str[l] = 0;
|
|
48 strcpy(str + l, ".sa");
|
|
49 bwt_restore_sa(str, p->bwt);
|
|
50 free(str);
|
|
51 p->pac = calloc(p->bns->l_pac/4+1, 1);
|
|
52 fread(p->pac, 1, p->bns->l_pac/4+1, p->bns->fp_pac);
|
|
53 fclose(p->bns->fp_pac);
|
|
54 p->bns->fp_pac = 0;
|
|
55 return p;
|
|
56 }
|
|
57
|
|
58 void bwa_idx_destroy(bwa_idx_t *p)
|
|
59 {
|
|
60 bns_destroy(p->bns);
|
|
61 bwt_destroy(p->bwt);
|
|
62 free(p->pac);
|
|
63 free(p);
|
|
64 }
|
|
65
|
|
66 bwa_buf_t *bwa_buf_init(const bwa_opt_t *opt, int max_score)
|
|
67 {
|
|
68 extern gap_opt_t *gap_init_opt(void);
|
|
69 extern int bwa_cal_maxdiff(int l, double err, double thres);
|
|
70 int i;
|
|
71 bwa_buf_t *p;
|
|
72 p = malloc(sizeof(bwa_buf_t));
|
|
73 p->stack = gap_init_stack2(max_score);
|
|
74 p->opt = gap_init_opt();
|
|
75 p->opt->s_gapo = opt->s_gapo;
|
|
76 p->opt->s_gape = opt->s_gape;
|
|
77 p->opt->max_diff = opt->max_diff;
|
|
78 p->opt->max_gapo = opt->max_gapo;
|
|
79 p->opt->max_gape = opt->max_gape;
|
|
80 p->opt->seed_len = opt->seed_len;
|
|
81 p->opt->max_seed_diff = opt->max_seed_diff;
|
|
82 p->opt->fnr = opt->fnr;
|
|
83 p->diff_tab = calloc(BWA_MAX_QUERY_LEN, sizeof(int));
|
|
84 for (i = 1; i < BWA_MAX_QUERY_LEN; ++i)
|
|
85 p->diff_tab[i] = bwa_cal_maxdiff(i, BWA_AVG_ERR, opt->fnr);
|
|
86 p->logn = calloc(256, sizeof(int));
|
|
87 for (i = 1; i != 256; ++i)
|
|
88 p->logn[i] = (int)(4.343 * log(i) + 0.499);
|
|
89 return p;
|
|
90 }
|
|
91
|
|
92 void bwa_buf_destroy(bwa_buf_t *p)
|
|
93 {
|
|
94 gap_destroy_stack(p->stack);
|
|
95 free(p->diff_tab); free(p->logn); free(p->opt);
|
|
96 free(p);
|
|
97 }
|
|
98
|
|
99 bwa_sai_t bwa_sai(const bwa_idx_t *idx, bwa_buf_t *buf, const char *seq)
|
|
100 {
|
|
101 extern int bwt_cal_width(const bwt_t *bwt, int len, const ubyte_t *str, bwt_width_t *width);
|
|
102 int i, seq_len, buf_len;
|
|
103 bwt_width_t *w, *seed_w;
|
|
104 uint8_t *s;
|
|
105 gap_opt_t opt2 = *buf->opt;
|
|
106 bwa_sai_t sai;
|
|
107
|
|
108 seq_len = strlen(seq);
|
|
109 // estimate the buffer length
|
|
110 buf_len = (buf->opt->seed_len + seq_len + 1) * sizeof(bwt_width_t) + seq_len;
|
|
111 if (buf_len > buf->max_buf) {
|
|
112 buf->max_buf = buf_len;
|
|
113 kroundup32(buf->max_buf);
|
|
114 buf->buf = realloc(buf->buf, buf->max_buf);
|
|
115 }
|
|
116 memset(buf->buf, 0, buf_len);
|
|
117 seed_w = (bwt_width_t*)buf->buf;
|
|
118 w = seed_w + buf->opt->seed_len;
|
|
119 s = (uint8_t*)(w + seq_len + 1);
|
|
120 if (opt2.fnr > 0.) opt2.max_diff = buf->diff_tab[seq_len];
|
|
121 // copy the sequence
|
|
122 for (i = 0; i < seq_len; ++i)
|
|
123 s[i] = nst_nt4_table[(int)seq[i]];
|
|
124 seq_reverse(seq_len, s, 0);
|
|
125 // mapping
|
|
126 bwt_cal_width(idx->bwt, seq_len, s, w);
|
|
127 if (opt2.seed_len >= seq_len) opt2.seed_len = 0x7fffffff;
|
|
128 if (seq_len > buf->opt->seed_len)
|
|
129 bwt_cal_width(idx->bwt, buf->opt->seed_len, s + (seq_len - buf->opt->seed_len), seed_w);
|
|
130 for (i = 0; i < seq_len; ++i) // complement; I forgot why...
|
|
131 s[i] = s[i] > 3? 4 : 3 - s[i];
|
|
132 sai.sai = (bwa_sai1_t*)bwt_match_gap(idx->bwt, seq_len, s, w, seq_len <= buf->opt->seed_len? 0 : seed_w, &opt2, &sai.n, buf->stack);
|
|
133 return sai;
|
|
134 }
|
|
135
|
|
136 static void compute_NM(const uint8_t *pac, uint64_t l_pac, uint8_t *seq, int64_t pos, int n_cigar, uint32_t *cigar, int *n_mm, int *n_gaps)
|
|
137 {
|
|
138 uint64_t x = pos, z;
|
|
139 int k, y = 0;
|
|
140 *n_mm = *n_gaps = 0;
|
|
141 for (k = 0; k < n_cigar; ++k) {
|
|
142 int l = cigar[k]>>4;
|
|
143 int op = cigar[k]&0xf;
|
|
144 if (op == 0) { // match/mismatch
|
|
145 for (z = 0; z < l && x + z < l_pac; ++z) {
|
|
146 int c = pac[(x+z)>>2] >> ((~(x+z)&3)<<1) & 3;
|
|
147 if (c > 3 || seq[y+z] > 3 || c != seq[y+z]) ++(*n_mm);
|
|
148 }
|
|
149 }
|
|
150 if (op == 1 || op == 2) (*n_gaps) += l;
|
|
151 if (op == 0 || op == 2) x += l;
|
|
152 if (op == 0 || op == 1 || op == 4) y += l;
|
|
153 }
|
|
154 }
|
|
155
|
|
156 void bwa_sa2aln(const bwa_idx_t *idx, bwa_buf_t *buf, const char *seq, uint64_t sa, int n_gaps, bwa_aln_t *aln)
|
|
157 {
|
|
158 extern bwtint_t bwa_sa2pos(const bntseq_t *bns, const bwt_t *bwt, bwtint_t sapos, int len, int *strand);
|
|
159 extern bwa_cigar_t *bwa_refine_gapped_core(bwtint_t l_pac, const ubyte_t *pacseq, int len, const uint8_t *seq, bwtint_t *_pos, int ext, int *n_cigar, int is_end_correct);
|
|
160 int strand, seq_len, i, n_gap, n_mm;
|
|
161 uint64_t pos3, pac_pos;
|
|
162 uint8_t *s[2];
|
|
163
|
|
164 memset(aln, 0, sizeof(bwa_aln_t));
|
|
165 seq_len = strlen(seq);
|
|
166 if (seq_len<<1 > buf->max_buf) {
|
|
167 buf->max_buf = seq_len<<1;
|
|
168 kroundup32(buf->max_buf);
|
|
169 buf->buf = realloc(buf->buf, buf->max_buf);
|
|
170 }
|
|
171 s[0] = buf->buf;
|
|
172 s[1] = s[0] + seq_len;
|
|
173 for (i = 0; i < seq_len; ++i)
|
|
174 s[0][i] = s[1][i] = nst_nt4_table[(int)seq[i]];
|
|
175 seq_reverse(seq_len, s[1], 1);
|
|
176 pac_pos = bwa_sa2pos(idx->bns, idx->bwt, sa, seq_len, &strand);
|
|
177 if (strand) aln->flag |= 16;
|
|
178 if (n_gaps) { // only for gapped alignment
|
|
179 int n_cigar;
|
|
180 bwa_cigar_t *cigar16;
|
|
181 cigar16 = bwa_refine_gapped_core(idx->bns->l_pac, idx->pac, seq_len, s[strand], &pac_pos, strand? n_gaps : -n_gaps, &n_cigar, 1);
|
|
182 aln->n_cigar = n_cigar;
|
|
183 aln->cigar = malloc(n_cigar * 4);
|
|
184 for (i = 0, pos3 = pac_pos; i < n_cigar; ++i) {
|
|
185 int op = cigar16[i]>>14;
|
|
186 int len = cigar16[i]&0x3fff;
|
|
187 if (op == 3) op = 4; // the 16-bit CIGAR is different from the 32-bit CIGAR
|
|
188 aln->cigar[i] = len<<4 | op;
|
|
189 if (op == 0 || op == 2) pos3 += len;
|
|
190 }
|
|
191 free(cigar16);
|
|
192 } else { // ungapped
|
|
193 aln->n_cigar = 1;
|
|
194 aln->cigar = malloc(4);
|
|
195 aln->cigar[0] = seq_len<<4 | 0;
|
|
196 pos3 = pac_pos + seq_len;
|
|
197 }
|
|
198 aln->n_n = bns_cnt_ambi(idx->bns, pac_pos, pos3 - pac_pos, &aln->ref_id);
|
|
199 aln->offset = pac_pos - idx->bns->anns[aln->ref_id].offset;
|
|
200 if (pos3 - idx->bns->anns[aln->ref_id].offset > idx->bns->anns[aln->ref_id].len) // read mapped beyond the end of a sequence
|
|
201 aln->flag |= 4; // read unmapped
|
|
202 compute_NM(idx->pac, idx->bns->l_pac, s[strand], pac_pos, aln->n_cigar, aln->cigar, &n_mm, &n_gap);
|
|
203 aln->n_mm = n_mm;
|
|
204 aln->n_gap = n_gap;
|
|
205 }
|
|
206
|
|
207 /************************
|
|
208 * Single-end alignment *
|
|
209 ************************/
|
|
210
|
|
211 bwa_one_t *bwa_se(const bwa_idx_t *idx, bwa_buf_t *buf, const char *seq, int gen_cigar)
|
|
212 {
|
|
213 bwa_one_t *one;
|
|
214 int best, cnt, i, seq_len;
|
|
215
|
|
216 seq_len = strlen(seq);
|
|
217 one = calloc(1, sizeof(bwa_one_t));
|
|
218 one->sai = bwa_sai(idx, buf, seq);
|
|
219 if (one->sai.n == 0) return one;
|
|
220 // count number of hits; randomly select one alignment
|
|
221 best = one->sai.sai[0].score;
|
|
222 for (i = cnt = 0; i < one->sai.n; ++i) {
|
|
223 bwa_sai1_t *p = &one->sai.sai[i];
|
|
224 if (p->score > best) break;
|
|
225 if (drand48() * (p->l - p->k + 1 + cnt) > (double)cnt) {
|
|
226 one->which = p;
|
|
227 one->sa = p->k + (bwtint_t)((p->l - p->k + 1) * drand48());
|
|
228 }
|
|
229 cnt += p->l - p->k + 1;
|
|
230 }
|
|
231 one->c1 = cnt;
|
|
232 for (; i < one->sai.n; ++i)
|
|
233 cnt += one->sai.sai[i].l - one->sai.sai[i].k + 1;
|
|
234 one->c2 = cnt - one->c1;
|
|
235 // estimate single-end mapping quality
|
|
236 one->mapQs = -1;
|
|
237 if (one->c1 == 0) one->mapQs = 23; // FIXME: is it possible?
|
|
238 else if (one->c1 > 1) one->mapQs = 0;
|
|
239 else {
|
|
240 int diff = one->which->n_mm + one->which->n_gapo + one->which->n_gape;
|
|
241 if (diff >= buf->diff_tab[seq_len]) one->mapQs = 25;
|
|
242 else if (one->c2 == 0) one->mapQs = 37;
|
|
243 }
|
|
244 if (one->mapQs < 0) {
|
|
245 cnt = (one->c2 >= 255)? 255 : one->c2;
|
|
246 one->mapQs = 23 < buf->logn[cnt]? 0 : 23 - buf->logn[cnt];
|
|
247 }
|
|
248 one->mapQ = one->mapQs;
|
|
249 // compute CIGAR on request
|
|
250 one->one.ref_id = -1;
|
|
251 if (gen_cigar) bwa_sa2aln(idx, buf, seq, one->sa, one->which->n_gapo + one->which->n_gape, &one->one);
|
|
252 return one;
|
|
253 }
|
|
254
|
|
255 void bwa_one_destroy(bwa_one_t *one)
|
|
256 {
|
|
257 free(one->sai.sai);
|
|
258 free(one->one.cigar);
|
|
259 free(one);
|
|
260 }
|
|
261
|
|
262 /************************
|
|
263 * Paired-end alignment *
|
|
264 ************************/
|
|
265
|
|
266 void bwa_pestat(bwa_buf_t *buf, int n, bwa_one_t **o[2])
|
|
267 {
|
|
268 }
|
|
269
|
|
270 void bwa_pe(const bwa_idx_t *idx, bwa_buf_t *buf, const char *seq[2], bwa_one_t *o[2])
|
|
271 {
|
|
272 }
|