0
|
1 #include <unistd.h>
|
|
2 #include <stdlib.h>
|
|
3 #include <string.h>
|
|
4 #include <stdio.h>
|
|
5 #include <math.h>
|
|
6 #include "bwt.h"
|
|
7 #include "bwtsw2.h"
|
|
8 #include "utils.h"
|
|
9
|
|
10 int bwa_bwtsw2(int argc, char *argv[])
|
|
11 {
|
|
12 extern char *bwa_infer_prefix(const char *hint);
|
|
13 bsw2opt_t *opt;
|
|
14 bwt_t *target;
|
|
15 char buf[1024], *prefix;
|
|
16 bntseq_t *bns;
|
|
17 int c;
|
|
18
|
|
19 opt = bsw2_init_opt();
|
|
20 srand48(11);
|
|
21 while ((c = getopt(argc, argv, "q:r:a:b:t:T:w:d:z:m:s:c:N:Hf:MI:S")) >= 0) {
|
|
22 switch (c) {
|
|
23 case 'q': opt->q = atoi(optarg); break;
|
|
24 case 'r': opt->r = atoi(optarg); break;
|
|
25 case 'a': opt->a = atoi(optarg); break;
|
|
26 case 'b': opt->b = atoi(optarg); break;
|
|
27 case 'w': opt->bw = atoi(optarg); break;
|
|
28 case 'T': opt->t = atoi(optarg); break;
|
|
29 case 't': opt->n_threads = atoi(optarg); break;
|
|
30 case 'z': opt->z = atoi(optarg); break;
|
|
31 case 's': opt->is = atoi(optarg); break;
|
|
32 case 'm': opt->mask_level = atof(optarg); break;
|
|
33 case 'c': opt->coef = atof(optarg); break;
|
|
34 case 'N': opt->t_seeds = atoi(optarg); break;
|
|
35 case 'M': opt->multi_2nd = 1; break;
|
|
36 case 'H': opt->hard_clip = 1; break;
|
|
37 case 'f': xreopen(optarg, "w", stdout); break;
|
|
38 case 'I': opt->max_ins = atoi(optarg); break;
|
|
39 case 'S': opt->skip_sw = 1; break;
|
|
40 }
|
|
41 }
|
|
42 opt->qr = opt->q + opt->r;
|
|
43
|
|
44 if (optind + 2 > argc) {
|
|
45 fprintf(stderr, "\n");
|
|
46 fprintf(stderr, "Usage: bwa bwasw [options] <target.prefix> <query.fa> [query2.fa]\n\n");
|
|
47 fprintf(stderr, "Options: -a INT score for a match [%d]\n", opt->a);
|
|
48 fprintf(stderr, " -b INT mismatch penalty [%d]\n", opt->b);
|
|
49 fprintf(stderr, " -q INT gap open penalty [%d]\n", opt->q);
|
|
50 fprintf(stderr, " -r INT gap extension penalty [%d]\n", opt->r);
|
|
51 fprintf(stderr, " -w INT band width [%d]\n", opt->bw);
|
|
52 fprintf(stderr, " -m FLOAT mask level [%.2f]\n", opt->mask_level);
|
|
53 fprintf(stderr, "\n");
|
|
54 fprintf(stderr, " -t INT number of threads [%d]\n", opt->n_threads);
|
|
55 fprintf(stderr, " -f FILE file to output results to instead of stdout\n");
|
|
56 fprintf(stderr, " -H in SAM output, use hard clipping instead of soft clipping\n");
|
|
57 fprintf(stderr, " -M mark multi-part alignments as secondary\n");
|
|
58 fprintf(stderr, " -S skip Smith-Waterman read pairing\n");
|
|
59 fprintf(stderr, " -I INT ignore pairs with insert >=INT for inferring the size distr [%d]\n", opt->max_ins);
|
|
60 fprintf(stderr, "\n");
|
|
61 fprintf(stderr, " -T INT score threshold divided by a [%d]\n", opt->t);
|
|
62 fprintf(stderr, " -c FLOAT coefficient of length-threshold adjustment [%.1f]\n", opt->coef);
|
|
63 fprintf(stderr, " -z INT Z-best [%d]\n", opt->z);
|
|
64 fprintf(stderr, " -s INT maximum seeding interval size [%d]\n", opt->is);
|
|
65 fprintf(stderr, " -N INT # seeds to trigger reverse alignment [%d]\n", opt->t_seeds);
|
|
66 fprintf(stderr, "\n");
|
|
67 fprintf(stderr, "Note: For long Illumina, 454 and Sanger reads, assembly contigs, fosmids and\n");
|
|
68 fprintf(stderr, " BACs, the default setting usually works well. For the current PacBio\n");
|
|
69 fprintf(stderr, " reads (end of 2010), '-b5 -q2 -r1 -z10' is recommended. One may also\n");
|
|
70 fprintf(stderr, " increase '-z' for better sensitivity.\n");
|
|
71 fprintf(stderr, "\n");
|
|
72
|
|
73 return 1;
|
|
74 }
|
|
75
|
|
76 // adjust opt for opt->a
|
|
77 opt->t *= opt->a;
|
|
78 opt->coef *= opt->a;
|
|
79
|
|
80 if ((prefix = bwa_infer_prefix(argv[optind])) == 0) {
|
|
81 fprintf(stderr, "[%s] fail to locate the index\n", __func__);
|
|
82 return 0;
|
|
83 }
|
|
84 strcpy(buf, prefix); target = bwt_restore_bwt(strcat(buf, ".bwt"));
|
|
85 strcpy(buf, prefix); bwt_restore_sa(strcat(buf, ".sa"), target);
|
|
86 bns = bns_restore(prefix);
|
|
87
|
|
88 bsw2_aln(opt, bns, target, argv[optind+1], optind+2 < argc? argv[optind+2] : 0);
|
|
89
|
|
90 bns_destroy(bns);
|
|
91 bwt_destroy(target);
|
|
92 free(opt); free(prefix);
|
|
93
|
|
94 return 0;
|
|
95 }
|