Mercurial > repos > siyuan > prada
comparison pyPRADA_1.2/tools/bwa-0.5.7-mh/bwtsw2_main.c @ 0:acc2ca1a3ba4
Uploaded
author | siyuan |
---|---|
date | Thu, 20 Feb 2014 00:44:58 -0500 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:acc2ca1a3ba4 |
---|---|
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 | |
9 int bwa_bwtsw2(int argc, char *argv[]) | |
10 { | |
11 bsw2opt_t *opt; | |
12 bwt_t *target[2]; | |
13 char buf[1024]; | |
14 bntseq_t *bns; | |
15 int c; | |
16 | |
17 opt = bsw2_init_opt(); | |
18 srand48(11); | |
19 while ((c = getopt(argc, argv, "q:r:a:b:t:T:w:d:z:m:y:s:c:N:H:f:")) >= 0) { | |
20 switch (c) { | |
21 case 'q': opt->q = atoi(optarg); break; | |
22 case 'r': opt->r = atoi(optarg); break; | |
23 case 'a': opt->a = atoi(optarg); break; | |
24 case 'b': opt->b = atoi(optarg); break; | |
25 case 'w': opt->bw = atoi(optarg); break; | |
26 case 'T': opt->t = atoi(optarg); break; | |
27 case 't': opt->n_threads = atoi(optarg); break; | |
28 case 'z': opt->z = atoi(optarg); break; | |
29 case 'y': opt->yita = atof(optarg); break; | |
30 case 's': opt->is = atoi(optarg); break; | |
31 case 'm': opt->mask_level = atof(optarg); break; | |
32 case 'c': opt->coef = atof(optarg); break; | |
33 case 'N': opt->t_seeds = atoi(optarg); break; | |
34 case 'H': opt->hard_clip = 1; break; | |
35 case 'f': freopen(optarg, "w", stdout); | |
36 } | |
37 } | |
38 opt->qr = opt->q + opt->r; | |
39 | |
40 if (optind + 2 > argc) { | |
41 fprintf(stderr, "\n"); | |
42 fprintf(stderr, "Usage: bwa bwasw [options] <target.prefix> <query.fa>\n\n"); | |
43 fprintf(stderr, "Options: -a INT score for a match [%d]\n", opt->a); | |
44 fprintf(stderr, " -b INT mismatch penalty [%d]\n", opt->b); | |
45 fprintf(stderr, " -q INT gap open penalty [%d]\n", opt->q); | |
46 fprintf(stderr, " -r INT gap extension penalty [%d]\n", opt->r); | |
47 // fprintf(stderr, " -y FLOAT error recurrence coef. (4..16) [%.1f]\n", opt->yita); | |
48 fprintf(stderr, "\n"); | |
49 fprintf(stderr, " -t INT nmber of threads [%d]\n", opt->n_threads); | |
50 fprintf(stderr, " -s INT size of a chunk of reads [%d]\n", opt->chunk_size); | |
51 fprintf(stderr, "\n"); | |
52 fprintf(stderr, " -w INT band width [%d]\n", opt->bw); | |
53 fprintf(stderr, " -m FLOAT mask level [%.2f]\n", opt->mask_level); | |
54 fprintf(stderr, "\n"); | |
55 fprintf(stderr, " -T INT score threshold divided by a [%d]\n", opt->t); | |
56 fprintf(stderr, " -s INT maximum seeding interval size [%d]\n", opt->is); | |
57 fprintf(stderr, " -z INT Z-best [%d]\n", opt->z); | |
58 fprintf(stderr, " -N INT # seeds to trigger reverse alignment [%d]\n", opt->t_seeds); | |
59 fprintf(stderr, " -c FLOAT coefficient of length-threshold adjustment [%.1f]\n", opt->coef); | |
60 fprintf(stderr, " -H in SAM output, use hard clipping rather than soft\n"); | |
61 fprintf(stderr, " -f FILE file to output results to instead of stdout\n"); | |
62 fprintf(stderr, "\n"); | |
63 | |
64 { | |
65 double c, theta, eps, delta; | |
66 c = opt->a / log(opt->yita); | |
67 theta = exp(-opt->b / c) / opt->yita; | |
68 eps = exp(-opt->q / c); | |
69 delta = exp(-opt->r / c); | |
70 fprintf(stderr, "mismatch: %lf, gap_open: %lf, gap_ext: %lf\n\n", | |
71 theta, eps, delta); | |
72 } | |
73 return 1; | |
74 } | |
75 | |
76 // adjust opt for opt->a | |
77 opt->t *= opt->a; | |
78 opt->coef *= opt->a; | |
79 | |
80 strcpy(buf, argv[optind]); target[0] = bwt_restore_bwt(strcat(buf, ".bwt")); | |
81 strcpy(buf, argv[optind]); bwt_restore_sa(strcat(buf, ".sa"), target[0]); | |
82 strcpy(buf, argv[optind]); target[1] = bwt_restore_bwt(strcat(buf, ".rbwt")); | |
83 strcpy(buf, argv[optind]); bwt_restore_sa(strcat(buf, ".rsa"), target[1]); | |
84 bns = bns_restore(argv[optind]); | |
85 | |
86 bsw2_aln(opt, bns, target, argv[optind+1]); | |
87 | |
88 bns_destroy(bns); | |
89 bwt_destroy(target[0]); bwt_destroy(target[1]); | |
90 free(opt); | |
91 | |
92 return 0; | |
93 } |