0
|
1 #include <stdio.h>
|
|
2 #include <unistd.h>
|
|
3 #include <math.h>
|
|
4 #include <stdlib.h>
|
|
5 #include <string.h>
|
|
6 #include <time.h>
|
|
7 #include <stdint.h>
|
|
8 #ifdef HAVE_CONFIG_H
|
|
9 #include "config.h"
|
|
10 #endif
|
|
11 #include "bwtaln.h"
|
|
12 #include "bwtgap.h"
|
|
13 #include "utils.h"
|
|
14
|
|
15 #ifdef HAVE_PTHREAD
|
|
16 #include <pthread.h>
|
|
17 #endif
|
|
18
|
|
19 gap_opt_t *gap_init_opt()
|
|
20 {
|
|
21 gap_opt_t *o;
|
|
22 o = (gap_opt_t*)calloc(1, sizeof(gap_opt_t));
|
|
23 /* IMPORTANT: s_mm*10 should be about the average base error
|
|
24 rate. Voilating this requirement will break pairing! */
|
|
25 o->s_mm = 3; o->s_gapo = 11; o->s_gape = 4;
|
|
26 o->max_diff = -1; o->max_gapo = 1; o->max_gape = 6;
|
|
27 o->indel_end_skip = 5; o->max_del_occ = 10; o->max_entries = 2000000;
|
|
28 o->mode = BWA_MODE_GAPE | BWA_MODE_COMPREAD;
|
|
29 o->seed_len = 32; o->max_seed_diff = 2;
|
|
30 o->fnr = 0.04;
|
|
31 o->n_threads = 1;
|
|
32 o->max_top2 = 30;
|
|
33 o->trim_qual = 0;
|
|
34 return o;
|
|
35 }
|
|
36
|
|
37 int bwa_cal_maxdiff(int l, double err, double thres)
|
|
38 {
|
|
39 double elambda = exp(-l * err);
|
|
40 double sum, y = 1.0;
|
|
41 int k, x = 1;
|
|
42 for (k = 1, sum = elambda; k < 1000; ++k) {
|
|
43 y *= l * err;
|
|
44 x *= k;
|
|
45 sum += elambda * y / x;
|
|
46 if (1.0 - sum < thres) return k;
|
|
47 }
|
|
48 return 2;
|
|
49 }
|
|
50
|
|
51 // width must be filled as zero
|
|
52 int bwt_cal_width(const bwt_t *bwt, int len, const ubyte_t *str, bwt_width_t *width)
|
|
53 {
|
|
54 bwtint_t k, l, ok, ol;
|
|
55 int i, bid;
|
|
56 bid = 0;
|
|
57 k = 0; l = bwt->seq_len;
|
|
58 for (i = 0; i < len; ++i) {
|
|
59 ubyte_t c = str[i];
|
|
60 if (c < 4) {
|
|
61 bwt_2occ(bwt, k - 1, l, c, &ok, &ol);
|
|
62 k = bwt->L2[c] + ok + 1;
|
|
63 l = bwt->L2[c] + ol;
|
|
64 }
|
|
65 if (k > l || c > 3) { // then restart
|
|
66 k = 0;
|
|
67 l = bwt->seq_len;
|
|
68 ++bid;
|
|
69 }
|
|
70 width[i].w = l - k + 1;
|
|
71 width[i].bid = bid;
|
|
72 }
|
|
73 width[len].w = 0;
|
|
74 width[len].bid = ++bid;
|
|
75 return bid;
|
|
76 }
|
|
77
|
|
78 void bwa_cal_sa_reg_gap(int tid, bwt_t *const bwt, int n_seqs, bwa_seq_t *seqs, const gap_opt_t *opt)
|
|
79 {
|
|
80 int i, j, max_l = 0, max_len;
|
|
81 gap_stack_t *stack;
|
|
82 bwt_width_t *w, *seed_w;
|
|
83 gap_opt_t local_opt = *opt;
|
|
84
|
|
85 // initiate priority stack
|
|
86 for (i = max_len = 0; i != n_seqs; ++i)
|
|
87 if (seqs[i].len > max_len) max_len = seqs[i].len;
|
|
88 if (opt->fnr > 0.0) local_opt.max_diff = bwa_cal_maxdiff(max_len, BWA_AVG_ERR, opt->fnr);
|
|
89 if (local_opt.max_diff < local_opt.max_gapo) local_opt.max_gapo = local_opt.max_diff;
|
|
90 stack = gap_init_stack(local_opt.max_diff, local_opt.max_gapo, local_opt.max_gape, &local_opt);
|
|
91
|
|
92 seed_w = (bwt_width_t*)calloc(opt->seed_len+1, sizeof(bwt_width_t));
|
|
93 w = 0;
|
|
94 for (i = 0; i != n_seqs; ++i) {
|
|
95 bwa_seq_t *p = seqs + i;
|
|
96 #ifdef HAVE_PTHREAD
|
|
97 if (i % opt->n_threads != tid) continue;
|
|
98 #endif
|
|
99 p->sa = 0; p->type = BWA_TYPE_NO_MATCH; p->c1 = p->c2 = 0; p->n_aln = 0; p->aln = 0;
|
|
100 if (max_l < p->len) {
|
|
101 max_l = p->len;
|
|
102 w = (bwt_width_t*)realloc(w, (max_l + 1) * sizeof(bwt_width_t));
|
|
103 memset(w, 0, (max_l + 1) * sizeof(bwt_width_t));
|
|
104 }
|
|
105 bwt_cal_width(bwt, p->len, p->seq, w);
|
|
106 if (opt->fnr > 0.0) local_opt.max_diff = bwa_cal_maxdiff(p->len, BWA_AVG_ERR, opt->fnr);
|
|
107 local_opt.seed_len = opt->seed_len < p->len? opt->seed_len : 0x7fffffff;
|
|
108 if (p->len > opt->seed_len)
|
|
109 bwt_cal_width(bwt, opt->seed_len, p->seq + (p->len - opt->seed_len), seed_w);
|
|
110 // core function
|
|
111 for (j = 0; j < p->len; ++j) // we need to complement
|
|
112 p->seq[j] = p->seq[j] > 3? 4 : 3 - p->seq[j];
|
|
113 p->aln = bwt_match_gap(bwt, p->len, p->seq, w, p->len <= opt->seed_len? 0 : seed_w, &local_opt, &p->n_aln, stack);
|
|
114 // clean up the unused data in the record
|
|
115 free(p->name); free(p->seq); free(p->rseq); free(p->qual);
|
|
116 p->name = 0; p->seq = p->rseq = p->qual = 0;
|
|
117 }
|
|
118 free(seed_w); free(w);
|
|
119 gap_destroy_stack(stack);
|
|
120 }
|
|
121
|
|
122 #ifdef HAVE_PTHREAD
|
|
123 typedef struct {
|
|
124 int tid;
|
|
125 bwt_t *bwt;
|
|
126 int n_seqs;
|
|
127 bwa_seq_t *seqs;
|
|
128 const gap_opt_t *opt;
|
|
129 } thread_aux_t;
|
|
130
|
|
131 static void *worker(void *data)
|
|
132 {
|
|
133 thread_aux_t *d = (thread_aux_t*)data;
|
|
134 bwa_cal_sa_reg_gap(d->tid, d->bwt, d->n_seqs, d->seqs, d->opt);
|
|
135 return 0;
|
|
136 }
|
|
137 #endif
|
|
138
|
|
139 bwa_seqio_t *bwa_open_reads(int mode, const char *fn_fa)
|
|
140 {
|
|
141 bwa_seqio_t *ks;
|
|
142 if (mode & BWA_MODE_BAM) { // open BAM
|
|
143 int which = 0;
|
|
144 if (mode & BWA_MODE_BAM_SE) which |= 4;
|
|
145 if (mode & BWA_MODE_BAM_READ1) which |= 1;
|
|
146 if (mode & BWA_MODE_BAM_READ2) which |= 2;
|
|
147 if (which == 0) which = 7; // then read all reads
|
|
148 ks = bwa_bam_open(fn_fa, which);
|
|
149 } else ks = bwa_seq_open(fn_fa);
|
|
150 return ks;
|
|
151 }
|
|
152
|
|
153 void bwa_aln_core(const char *prefix, const char *fn_fa, const gap_opt_t *opt)
|
|
154 {
|
|
155 int i, n_seqs, tot_seqs = 0;
|
|
156 bwa_seq_t *seqs;
|
|
157 bwa_seqio_t *ks;
|
|
158 clock_t t;
|
|
159 bwt_t *bwt;
|
|
160
|
|
161 // initialization
|
|
162 ks = bwa_open_reads(opt->mode, fn_fa);
|
|
163
|
|
164 { // load BWT
|
|
165 char *str = (char*)calloc(strlen(prefix) + 10, 1);
|
|
166 strcpy(str, prefix); strcat(str, ".bwt"); bwt = bwt_restore_bwt(str);
|
|
167 free(str);
|
|
168 }
|
|
169
|
|
170 // core loop
|
|
171 err_fwrite(opt, sizeof(gap_opt_t), 1, stdout);
|
|
172 while ((seqs = bwa_read_seq(ks, 0x40000, &n_seqs, opt->mode, opt->trim_qual)) != 0) {
|
|
173 tot_seqs += n_seqs;
|
|
174 t = clock();
|
|
175
|
|
176 fprintf(stderr, "[bwa_aln_core] calculate SA coordinate... ");
|
|
177
|
|
178 #ifdef HAVE_PTHREAD
|
|
179 if (opt->n_threads <= 1) { // no multi-threading at all
|
|
180 bwa_cal_sa_reg_gap(0, bwt, n_seqs, seqs, opt);
|
|
181 } else {
|
|
182 pthread_t *tid;
|
|
183 pthread_attr_t attr;
|
|
184 thread_aux_t *data;
|
|
185 int j;
|
|
186 pthread_attr_init(&attr);
|
|
187 pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
|
|
188 data = (thread_aux_t*)calloc(opt->n_threads, sizeof(thread_aux_t));
|
|
189 tid = (pthread_t*)calloc(opt->n_threads, sizeof(pthread_t));
|
|
190 for (j = 0; j < opt->n_threads; ++j) {
|
|
191 data[j].tid = j; data[j].bwt = bwt;
|
|
192 data[j].n_seqs = n_seqs; data[j].seqs = seqs; data[j].opt = opt;
|
|
193 pthread_create(&tid[j], &attr, worker, data + j);
|
|
194 }
|
|
195 for (j = 0; j < opt->n_threads; ++j) pthread_join(tid[j], 0);
|
|
196 free(data); free(tid);
|
|
197 }
|
|
198 #else
|
|
199 bwa_cal_sa_reg_gap(0, bwt, n_seqs, seqs, opt);
|
|
200 #endif
|
|
201
|
|
202 fprintf(stderr, "%.2f sec\n", (float)(clock() - t) / CLOCKS_PER_SEC); t = clock();
|
|
203
|
|
204 t = clock();
|
|
205 fprintf(stderr, "[bwa_aln_core] write to the disk... ");
|
|
206 for (i = 0; i < n_seqs; ++i) {
|
|
207 bwa_seq_t *p = seqs + i;
|
|
208 err_fwrite(&p->n_aln, 4, 1, stdout);
|
|
209 if (p->n_aln) err_fwrite(p->aln, sizeof(bwt_aln1_t), p->n_aln, stdout);
|
|
210 }
|
|
211 fprintf(stderr, "%.2f sec\n", (float)(clock() - t) / CLOCKS_PER_SEC); t = clock();
|
|
212
|
|
213 bwa_free_read_seq(n_seqs, seqs);
|
|
214 fprintf(stderr, "[bwa_aln_core] %d sequences have been processed.\n", tot_seqs);
|
|
215 }
|
|
216
|
|
217 // destroy
|
|
218 bwt_destroy(bwt);
|
|
219 bwa_seq_close(ks);
|
|
220 }
|
|
221
|
|
222 char *bwa_infer_prefix(const char *hint)
|
|
223 {
|
|
224 char *prefix;
|
|
225 int l_hint;
|
|
226 FILE *fp;
|
|
227 l_hint = strlen(hint);
|
|
228 prefix = malloc(l_hint + 3 + 4 + 1);
|
|
229 strcpy(prefix, hint);
|
|
230 strcpy(prefix + l_hint, ".64.bwt");
|
|
231 if ((fp = fopen(prefix, "rb")) != 0) {
|
|
232 fclose(fp);
|
|
233 prefix[l_hint + 3] = 0;
|
|
234 return prefix;
|
|
235 } else {
|
|
236 strcpy(prefix + l_hint, ".bwt");
|
|
237 if ((fp = fopen(prefix, "rb")) == 0) {
|
|
238 free(prefix);
|
|
239 return 0;
|
|
240 } else {
|
|
241 fclose(fp);
|
|
242 prefix[l_hint] = 0;
|
|
243 return prefix;
|
|
244 }
|
|
245 }
|
|
246 }
|
|
247
|
|
248 int bwa_aln(int argc, char *argv[])
|
|
249 {
|
|
250 int c, opte = -1;
|
|
251 gap_opt_t *opt;
|
|
252 char *prefix;
|
|
253
|
|
254 opt = gap_init_opt();
|
|
255 while ((c = getopt(argc, argv, "n:o:e:i:d:l:k:cLR:m:t:NM:O:E:q:f:b012IYB:")) >= 0) {
|
|
256 switch (c) {
|
|
257 case 'n':
|
|
258 if (strstr(optarg, ".")) opt->fnr = atof(optarg), opt->max_diff = -1;
|
|
259 else opt->max_diff = atoi(optarg), opt->fnr = -1.0;
|
|
260 break;
|
|
261 case 'o': opt->max_gapo = atoi(optarg); break;
|
|
262 case 'e': opte = atoi(optarg); break;
|
|
263 case 'M': opt->s_mm = atoi(optarg); break;
|
|
264 case 'O': opt->s_gapo = atoi(optarg); break;
|
|
265 case 'E': opt->s_gape = atoi(optarg); break;
|
|
266 case 'd': opt->max_del_occ = atoi(optarg); break;
|
|
267 case 'i': opt->indel_end_skip = atoi(optarg); break;
|
|
268 case 'l': opt->seed_len = atoi(optarg); break;
|
|
269 case 'k': opt->max_seed_diff = atoi(optarg); break;
|
|
270 case 'm': opt->max_entries = atoi(optarg); break;
|
|
271 case 't': opt->n_threads = atoi(optarg); break;
|
|
272 case 'L': opt->mode |= BWA_MODE_LOGGAP; break;
|
|
273 case 'R': opt->max_top2 = atoi(optarg); break;
|
|
274 case 'q': opt->trim_qual = atoi(optarg); break;
|
|
275 case 'c': opt->mode &= ~BWA_MODE_COMPREAD; break;
|
|
276 case 'N': opt->mode |= BWA_MODE_NONSTOP; opt->max_top2 = 0x7fffffff; break;
|
|
277 case 'f': xreopen(optarg, "wb", stdout); break;
|
|
278 case 'b': opt->mode |= BWA_MODE_BAM; break;
|
|
279 case '0': opt->mode |= BWA_MODE_BAM_SE; break;
|
|
280 case '1': opt->mode |= BWA_MODE_BAM_READ1; break;
|
|
281 case '2': opt->mode |= BWA_MODE_BAM_READ2; break;
|
|
282 case 'I': opt->mode |= BWA_MODE_IL13; break;
|
|
283 case 'Y': opt->mode |= BWA_MODE_CFY; break;
|
|
284 case 'B': opt->mode |= atoi(optarg) << 24; break;
|
|
285 default: return 1;
|
|
286 }
|
|
287 }
|
|
288 if (opte > 0) {
|
|
289 opt->max_gape = opte;
|
|
290 opt->mode &= ~BWA_MODE_GAPE;
|
|
291 }
|
|
292
|
|
293 if (optind + 2 > argc) {
|
|
294 fprintf(stderr, "\n");
|
|
295 fprintf(stderr, "Usage: bwa aln [options] <prefix> <in.fq>\n\n");
|
|
296 fprintf(stderr, "Options: -n NUM max #diff (int) or missing prob under %.2f err rate (float) [%.2f]\n",
|
|
297 BWA_AVG_ERR, opt->fnr);
|
|
298 fprintf(stderr, " -o INT maximum number or fraction of gap opens [%d]\n", opt->max_gapo);
|
|
299 fprintf(stderr, " -e INT maximum number of gap extensions, -1 for disabling long gaps [-1]\n");
|
|
300 fprintf(stderr, " -i INT do not put an indel within INT bp towards the ends [%d]\n", opt->indel_end_skip);
|
|
301 fprintf(stderr, " -d INT maximum occurrences for extending a long deletion [%d]\n", opt->max_del_occ);
|
|
302 fprintf(stderr, " -l INT seed length [%d]\n", opt->seed_len);
|
|
303 fprintf(stderr, " -k INT maximum differences in the seed [%d]\n", opt->max_seed_diff);
|
|
304 fprintf(stderr, " -m INT maximum entries in the queue [%d]\n", opt->max_entries);
|
|
305 fprintf(stderr, " -t INT number of threads [%d]\n", opt->n_threads);
|
|
306 fprintf(stderr, " -M INT mismatch penalty [%d]\n", opt->s_mm);
|
|
307 fprintf(stderr, " -O INT gap open penalty [%d]\n", opt->s_gapo);
|
|
308 fprintf(stderr, " -E INT gap extension penalty [%d]\n", opt->s_gape);
|
|
309 fprintf(stderr, " -R INT stop searching when there are >INT equally best hits [%d]\n", opt->max_top2);
|
|
310 fprintf(stderr, " -q INT quality threshold for read trimming down to %dbp [%d]\n", BWA_MIN_RDLEN, opt->trim_qual);
|
|
311 fprintf(stderr, " -f FILE file to write output to instead of stdout\n");
|
|
312 fprintf(stderr, " -B INT length of barcode\n");
|
|
313 // fprintf(stderr, " -c input sequences are in the color space\n");
|
|
314 fprintf(stderr, " -L log-scaled gap penalty for long deletions\n");
|
|
315 fprintf(stderr, " -N non-iterative mode: search for all n-difference hits (slooow)\n");
|
|
316 fprintf(stderr, " -I the input is in the Illumina 1.3+ FASTQ-like format\n");
|
|
317 fprintf(stderr, " -b the input read file is in the BAM format\n");
|
|
318 fprintf(stderr, " -0 use single-end reads only (effective with -b)\n");
|
|
319 fprintf(stderr, " -1 use the 1st read in a pair (effective with -b)\n");
|
|
320 fprintf(stderr, " -2 use the 2nd read in a pair (effective with -b)\n");
|
|
321 fprintf(stderr, " -Y filter Casava-filtered sequences\n");
|
|
322 fprintf(stderr, "\n");
|
|
323 return 1;
|
|
324 }
|
|
325 if (opt->fnr > 0.0) {
|
|
326 int i, k;
|
|
327 for (i = 17, k = 0; i <= 250; ++i) {
|
|
328 int l = bwa_cal_maxdiff(i, BWA_AVG_ERR, opt->fnr);
|
|
329 if (l != k) fprintf(stderr, "[bwa_aln] %dbp reads: max_diff = %d\n", i, l);
|
|
330 k = l;
|
|
331 }
|
|
332 }
|
|
333 if ((prefix = bwa_infer_prefix(argv[optind])) == 0) {
|
|
334 fprintf(stderr, "[%s] fail to locate the index\n", __func__);
|
|
335 free(opt);
|
|
336 return 0;
|
|
337 }
|
|
338 bwa_aln_core(prefix, argv[optind+1], opt);
|
|
339 free(opt); free(prefix);
|
|
340 return 0;
|
|
341 }
|
|
342
|
|
343 /* rgoya: Temporary clone of aln_path2cigar to accomodate for bwa_cigar_t,
|
|
344 __cigar_op and __cigar_len while keeping stdaln stand alone */
|
|
345 bwa_cigar_t *bwa_aln_path2cigar(const path_t *path, int path_len, int *n_cigar)
|
|
346 {
|
|
347 uint32_t *cigar32;
|
|
348 bwa_cigar_t *cigar;
|
|
349 int i;
|
|
350 cigar32 = aln_path2cigar32((path_t*) path, path_len, n_cigar);
|
|
351 cigar = (bwa_cigar_t*)cigar32;
|
|
352 for (i = 0; i < *n_cigar; ++i)
|
|
353 cigar[i] = __cigar_create( (cigar32[i]&0xf), (cigar32[i]>>4) );
|
|
354 return cigar;
|
|
355 }
|
|
356
|