|
0
|
1 /* The MIT License
|
|
|
2
|
|
|
3 Copyright (c) 2008 Genome Research Ltd (GRL).
|
|
|
4
|
|
|
5 Permission is hereby granted, free of charge, to any person obtaining
|
|
|
6 a copy of this software and associated documentation files (the
|
|
|
7 "Software"), to deal in the Software without restriction, including
|
|
|
8 without limitation the rights to use, copy, modify, merge, publish,
|
|
|
9 distribute, sublicense, and/or sell copies of the Software, and to
|
|
|
10 permit persons to whom the Software is furnished to do so, subject to
|
|
|
11 the following conditions:
|
|
|
12
|
|
|
13 The above copyright notice and this permission notice shall be
|
|
|
14 included in all copies or substantial portions of the Software.
|
|
|
15
|
|
|
16 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
|
|
17 EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
|
|
|
18 MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
|
|
19 NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
|
|
|
20 BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
|
|
|
21 ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
|
|
|
22 CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
|
|
|
23 SOFTWARE.
|
|
|
24 */
|
|
|
25
|
|
|
26 /* Contact: Heng Li <lh3@sanger.ac.uk> */
|
|
|
27
|
|
|
28 #include <stdio.h>
|
|
|
29 #include <stdlib.h>
|
|
|
30 #include <string.h>
|
|
|
31 #include <unistd.h>
|
|
|
32 #include <time.h>
|
|
|
33 #include <zlib.h>
|
|
|
34 #include "bntseq.h"
|
|
|
35 #include "bwt.h"
|
|
|
36 #include "main.h"
|
|
|
37 #include "utils.h"
|
|
|
38
|
|
|
39 bwt_t *bwt_pac2bwt(const char *fn_pac, int use_is);
|
|
|
40 void bwa_pac_rev_core(const char *fn, const char *fn_rev);
|
|
|
41
|
|
|
42 int bwa_index(int argc, char *argv[])
|
|
|
43 {
|
|
|
44 char *prefix = 0, *str, *str2, *str3;
|
|
|
45 int c, algo_type = 3, is_color = 0;
|
|
|
46 clock_t t;
|
|
|
47
|
|
|
48 while ((c = getopt(argc, argv, "ca:p:")) >= 0) {
|
|
|
49 switch (c) {
|
|
|
50 case 'a':
|
|
|
51 if (strcmp(optarg, "div") == 0) algo_type = 1;
|
|
|
52 else if (strcmp(optarg, "bwtsw") == 0) algo_type = 2;
|
|
|
53 else if (strcmp(optarg, "is") == 0) algo_type = 3;
|
|
|
54 else err_fatal(__func__, "unknown algorithm: '%s'.", optarg);
|
|
|
55 break;
|
|
|
56 case 'p': prefix = strdup(optarg); break;
|
|
|
57 case 'c': is_color = 1; break;
|
|
|
58 default: return 1;
|
|
|
59 }
|
|
|
60 }
|
|
|
61
|
|
|
62 if (optind + 1 > argc) {
|
|
|
63 fprintf(stderr, "\n");
|
|
|
64 fprintf(stderr, "Usage: bwa index [-a bwtsw|div|is] [-c] <in.fasta>\n\n");
|
|
|
65 fprintf(stderr, "Options: -a STR BWT construction algorithm: bwtsw or is [is]\n");
|
|
|
66 fprintf(stderr, " -p STR prefix of the index [same as fasta name]\n");
|
|
|
67 fprintf(stderr, " -c build color-space index\n\n");
|
|
|
68 fprintf(stderr, "Warning: `-a bwtsw' does not work for short genomes, while `-a is' and\n");
|
|
|
69 fprintf(stderr, " `-a div' do not work not for long genomes. Please choose `-a'\n");
|
|
|
70 fprintf(stderr, " according to the length of the genome.\n\n");
|
|
|
71 return 1;
|
|
|
72 }
|
|
|
73 if (prefix == 0) prefix = strdup(argv[optind]);
|
|
|
74 str = (char*)calloc(strlen(prefix) + 10, 1);
|
|
|
75 str2 = (char*)calloc(strlen(prefix) + 10, 1);
|
|
|
76 str3 = (char*)calloc(strlen(prefix) + 10, 1);
|
|
|
77
|
|
|
78 if (is_color == 0) { // nucleotide indexing
|
|
|
79 gzFile fp = xzopen(argv[optind], "r");
|
|
|
80 t = clock();
|
|
|
81 fprintf(stderr, "[bwa_index] Pack FASTA... ");
|
|
|
82 bns_fasta2bntseq(fp, prefix);
|
|
|
83 fprintf(stderr, "%.2f sec\n", (float)(clock() - t) / CLOCKS_PER_SEC);
|
|
|
84 gzclose(fp);
|
|
|
85 } else { // color indexing
|
|
|
86 gzFile fp = xzopen(argv[optind], "r");
|
|
|
87 strcat(strcpy(str, prefix), ".nt");
|
|
|
88 t = clock();
|
|
|
89 fprintf(stderr, "[bwa_index] Pack nucleotide FASTA... ");
|
|
|
90 bns_fasta2bntseq(fp, str);
|
|
|
91 fprintf(stderr, "%.2f sec\n", (float)(clock() - t) / CLOCKS_PER_SEC);
|
|
|
92 gzclose(fp);
|
|
|
93 {
|
|
|
94 char *tmp_argv[3];
|
|
|
95 tmp_argv[0] = argv[0]; tmp_argv[1] = str; tmp_argv[2] = prefix;
|
|
|
96 t = clock();
|
|
|
97 fprintf(stderr, "[bwa_index] Convert nucleotide PAC to color PAC... ");
|
|
|
98 bwa_pac2cspac(3, tmp_argv);
|
|
|
99 fprintf(stderr, "%.2f sec\n", (float)(clock() - t) / CLOCKS_PER_SEC);
|
|
|
100 }
|
|
|
101 }
|
|
|
102 {
|
|
|
103 strcpy(str, prefix); strcat(str, ".pac");
|
|
|
104 strcpy(str2, prefix); strcat(str2, ".rpac");
|
|
|
105 t = clock();
|
|
|
106 fprintf(stderr, "[bwa_index] Reverse the packed sequence... ");
|
|
|
107 bwa_pac_rev_core(str, str2);
|
|
|
108 fprintf(stderr, "%.2f sec\n", (float)(clock() - t) / CLOCKS_PER_SEC);
|
|
|
109 }
|
|
|
110 {
|
|
|
111 strcpy(str, prefix); strcat(str, ".pac");
|
|
|
112 strcpy(str2, prefix); strcat(str2, ".bwt");
|
|
|
113 t = clock();
|
|
|
114 fprintf(stderr, "[bwa_index] Construct BWT for the packed sequence...\n");
|
|
|
115 if (algo_type == 2) bwt_bwtgen(str, str2);
|
|
|
116 else if (algo_type == 1 || algo_type == 3) {
|
|
|
117 bwt_t *bwt;
|
|
|
118 bwt = bwt_pac2bwt(str, algo_type == 3);
|
|
|
119 bwt_dump_bwt(str2, bwt);
|
|
|
120 bwt_destroy(bwt);
|
|
|
121 }
|
|
|
122 fprintf(stderr, "[bwa_index] %.2f seconds elapse.\n", (float)(clock() - t) / CLOCKS_PER_SEC);
|
|
|
123 }
|
|
|
124 {
|
|
|
125 strcpy(str, prefix); strcat(str, ".rpac");
|
|
|
126 strcpy(str2, prefix); strcat(str2, ".rbwt");
|
|
|
127 t = clock();
|
|
|
128 fprintf(stderr, "[bwa_index] Construct BWT for the reverse packed sequence...\n");
|
|
|
129 if (algo_type == 2) bwt_bwtgen(str, str2);
|
|
|
130 else if (algo_type == 1 || algo_type == 3) {
|
|
|
131 bwt_t *bwt;
|
|
|
132 bwt = bwt_pac2bwt(str, algo_type == 3);
|
|
|
133 bwt_dump_bwt(str2, bwt);
|
|
|
134 bwt_destroy(bwt);
|
|
|
135 }
|
|
|
136 fprintf(stderr, "[bwa_index] %.2f seconds elapse.\n", (float)(clock() - t) / CLOCKS_PER_SEC);
|
|
|
137 }
|
|
|
138 {
|
|
|
139 bwt_t *bwt;
|
|
|
140 strcpy(str, prefix); strcat(str, ".bwt");
|
|
|
141 t = clock();
|
|
|
142 fprintf(stderr, "[bwa_index] Update BWT... ");
|
|
|
143 bwt = bwt_restore_bwt(str);
|
|
|
144 bwt_bwtupdate_core(bwt);
|
|
|
145 bwt_dump_bwt(str, bwt);
|
|
|
146 bwt_destroy(bwt);
|
|
|
147 fprintf(stderr, "%.2f sec\n", (float)(clock() - t) / CLOCKS_PER_SEC);
|
|
|
148 }
|
|
|
149 {
|
|
|
150 bwt_t *bwt;
|
|
|
151 strcpy(str, prefix); strcat(str, ".rbwt");
|
|
|
152 t = clock();
|
|
|
153 fprintf(stderr, "[bwa_index] Update reverse BWT... ");
|
|
|
154 bwt = bwt_restore_bwt(str);
|
|
|
155 bwt_bwtupdate_core(bwt);
|
|
|
156 bwt_dump_bwt(str, bwt);
|
|
|
157 bwt_destroy(bwt);
|
|
|
158 fprintf(stderr, "%.2f sec\n", (float)(clock() - t) / CLOCKS_PER_SEC);
|
|
|
159 }
|
|
|
160 {
|
|
|
161 bwt_t *bwt;
|
|
|
162 strcpy(str, prefix); strcat(str, ".bwt");
|
|
|
163 strcpy(str3, prefix); strcat(str3, ".sa");
|
|
|
164 t = clock();
|
|
|
165 fprintf(stderr, "[bwa_index] Construct SA from BWT and Occ... ");
|
|
|
166 bwt = bwt_restore_bwt(str);
|
|
|
167 bwt_cal_sa(bwt, 32);
|
|
|
168 bwt_dump_sa(str3, bwt);
|
|
|
169 bwt_destroy(bwt);
|
|
|
170 fprintf(stderr, "%.2f sec\n", (float)(clock() - t) / CLOCKS_PER_SEC);
|
|
|
171 }
|
|
|
172 {
|
|
|
173 bwt_t *bwt;
|
|
|
174 strcpy(str, prefix); strcat(str, ".rbwt");
|
|
|
175 strcpy(str3, prefix); strcat(str3, ".rsa");
|
|
|
176 t = clock();
|
|
|
177 fprintf(stderr, "[bwa_index] Construct SA from reverse BWT and Occ... ");
|
|
|
178 bwt = bwt_restore_bwt(str);
|
|
|
179 bwt_cal_sa(bwt, 32);
|
|
|
180 bwt_dump_sa(str3, bwt);
|
|
|
181 bwt_destroy(bwt);
|
|
|
182 fprintf(stderr, "%.2f sec\n", (float)(clock() - t) / CLOCKS_PER_SEC);
|
|
|
183 }
|
|
|
184 free(str3); free(str2); free(str); free(prefix);
|
|
|
185 return 0;
|
|
|
186 }
|
