|
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 "utils.h"
|
|
|
37
|
|
|
38 #ifdef _DIVBWT
|
|
|
39 #include "divsufsort.h"
|
|
|
40 #endif
|
|
|
41
|
|
|
42 #ifdef USE_MALLOC_WRAPPERS
|
|
|
43 # include "malloc_wrap.h"
|
|
|
44 #endif
|
|
|
45
|
|
|
46
|
|
|
47 int is_bwt(ubyte_t *T, int n);
|
|
|
48
|
|
|
49 int64_t bwa_seq_len(const char *fn_pac)
|
|
|
50 {
|
|
|
51 FILE *fp;
|
|
|
52 int64_t pac_len;
|
|
|
53 ubyte_t c;
|
|
|
54 fp = xopen(fn_pac, "rb");
|
|
|
55 err_fseek(fp, -1, SEEK_END);
|
|
|
56 pac_len = err_ftell(fp);
|
|
|
57 err_fread_noeof(&c, 1, 1, fp);
|
|
|
58 err_fclose(fp);
|
|
|
59 return (pac_len - 1) * 4 + (int)c;
|
|
|
60 }
|
|
|
61
|
|
|
62 bwt_t *bwt_pac2bwt(const char *fn_pac, int use_is)
|
|
|
63 {
|
|
|
64 bwt_t *bwt;
|
|
|
65 ubyte_t *buf, *buf2;
|
|
|
66 int i, pac_size;
|
|
|
67 FILE *fp;
|
|
|
68
|
|
|
69 // initialization
|
|
|
70 bwt = (bwt_t*)calloc(1, sizeof(bwt_t));
|
|
|
71 bwt->seq_len = bwa_seq_len(fn_pac);
|
|
|
72 bwt->bwt_size = (bwt->seq_len + 15) >> 4;
|
|
|
73 fp = xopen(fn_pac, "rb");
|
|
|
74
|
|
|
75 // prepare sequence
|
|
|
76 pac_size = (bwt->seq_len>>2) + ((bwt->seq_len&3) == 0? 0 : 1);
|
|
|
77 buf2 = (ubyte_t*)calloc(pac_size, 1);
|
|
|
78 err_fread_noeof(buf2, 1, pac_size, fp);
|
|
|
79 err_fclose(fp);
|
|
|
80 memset(bwt->L2, 0, 5 * 4);
|
|
|
81 buf = (ubyte_t*)calloc(bwt->seq_len + 1, 1);
|
|
|
82 for (i = 0; i < bwt->seq_len; ++i) {
|
|
|
83 buf[i] = buf2[i>>2] >> ((3 - (i&3)) << 1) & 3;
|
|
|
84 ++bwt->L2[1+buf[i]];
|
|
|
85 }
|
|
|
86 for (i = 2; i <= 4; ++i) bwt->L2[i] += bwt->L2[i-1];
|
|
|
87 free(buf2);
|
|
|
88
|
|
|
89 // Burrows-Wheeler Transform
|
|
|
90 if (use_is) {
|
|
|
91 bwt->primary = is_bwt(buf, bwt->seq_len);
|
|
|
92 } else {
|
|
|
93 #ifdef _DIVBWT
|
|
|
94 bwt->primary = divbwt(buf, buf, 0, bwt->seq_len);
|
|
|
95 #else
|
|
|
96 err_fatal_simple("libdivsufsort is not compiled in.");
|
|
|
97 #endif
|
|
|
98 }
|
|
|
99 bwt->bwt = (u_int32_t*)calloc(bwt->bwt_size, 4);
|
|
|
100 for (i = 0; i < bwt->seq_len; ++i)
|
|
|
101 bwt->bwt[i>>4] |= buf[i] << ((15 - (i&15)) << 1);
|
|
|
102 free(buf);
|
|
|
103 return bwt;
|
|
|
104 }
|
|
|
105
|
|
|
106 int bwa_pac2bwt(int argc, char *argv[]) // the "pac2bwt" command; IMPORTANT: bwt generated at this step CANNOT be used with BWA. bwtupdate is required!
|
|
|
107 {
|
|
|
108 bwt_t *bwt;
|
|
|
109 int c, use_is = 1;
|
|
|
110 while ((c = getopt(argc, argv, "d")) >= 0) {
|
|
|
111 switch (c) {
|
|
|
112 case 'd': use_is = 0; break;
|
|
|
113 default: return 1;
|
|
|
114 }
|
|
|
115 }
|
|
|
116 if (optind + 2 > argc) {
|
|
|
117 fprintf(stderr, "Usage: bwa pac2bwt [-d] <in.pac> <out.bwt>\n");
|
|
|
118 return 1;
|
|
|
119 }
|
|
|
120 bwt = bwt_pac2bwt(argv[optind], use_is);
|
|
|
121 bwt_dump_bwt(argv[optind+1], bwt);
|
|
|
122 bwt_destroy(bwt);
|
|
|
123 return 0;
|
|
|
124 }
|
|
|
125
|
|
|
126 #define bwt_B00(b, k) ((b)->bwt[(k)>>4]>>((~(k)&0xf)<<1)&3)
|
|
|
127
|
|
|
128 void bwt_bwtupdate_core(bwt_t *bwt)
|
|
|
129 {
|
|
|
130 bwtint_t i, k, c[4], n_occ;
|
|
|
131 uint32_t *buf;
|
|
|
132
|
|
|
133 n_occ = (bwt->seq_len + OCC_INTERVAL - 1) / OCC_INTERVAL + 1;
|
|
|
134 bwt->bwt_size += n_occ * sizeof(bwtint_t); // the new size
|
|
|
135 buf = (uint32_t*)calloc(bwt->bwt_size, 4); // will be the new bwt
|
|
|
136 c[0] = c[1] = c[2] = c[3] = 0;
|
|
|
137 for (i = k = 0; i < bwt->seq_len; ++i) {
|
|
|
138 if (i % OCC_INTERVAL == 0) {
|
|
|
139 memcpy(buf + k, c, sizeof(bwtint_t) * 4);
|
|
|
140 k += sizeof(bwtint_t); // in fact: sizeof(bwtint_t)=4*(sizeof(bwtint_t)/4)
|
|
|
141 }
|
|
|
142 if (i % 16 == 0) buf[k++] = bwt->bwt[i/16]; // 16 == sizeof(uint32_t)/2
|
|
|
143 ++c[bwt_B00(bwt, i)];
|
|
|
144 }
|
|
|
145 // the last element
|
|
|
146 memcpy(buf + k, c, sizeof(bwtint_t) * 4);
|
|
|
147 xassert(k + sizeof(bwtint_t) == bwt->bwt_size, "inconsistent bwt_size");
|
|
|
148 // update bwt
|
|
|
149 free(bwt->bwt); bwt->bwt = buf;
|
|
|
150 }
|
|
|
151
|
|
|
152 int bwa_bwtupdate(int argc, char *argv[]) // the "bwtupdate" command
|
|
|
153 {
|
|
|
154 bwt_t *bwt;
|
|
|
155 if (argc < 2) {
|
|
|
156 fprintf(stderr, "Usage: bwa bwtupdate <the.bwt>\n");
|
|
|
157 return 1;
|
|
|
158 }
|
|
|
159 bwt = bwt_restore_bwt(argv[1]);
|
|
|
160 bwt_bwtupdate_core(bwt);
|
|
|
161 bwt_dump_bwt(argv[1], bwt);
|
|
|
162 bwt_destroy(bwt);
|
|
|
163 return 0;
|
|
|
164 }
|
|
|
165
|
|
|
166 int bwa_bwt2sa(int argc, char *argv[]) // the "bwt2sa" command
|
|
|
167 {
|
|
|
168 bwt_t *bwt;
|
|
|
169 int c, sa_intv = 32;
|
|
|
170 while ((c = getopt(argc, argv, "i:")) >= 0) {
|
|
|
171 switch (c) {
|
|
|
172 case 'i': sa_intv = atoi(optarg); break;
|
|
|
173 default: return 1;
|
|
|
174 }
|
|
|
175 }
|
|
|
176 if (optind + 2 > argc) {
|
|
|
177 fprintf(stderr, "Usage: bwa bwt2sa [-i %d] <in.bwt> <out.sa>\n", sa_intv);
|
|
|
178 return 1;
|
|
|
179 }
|
|
|
180 bwt = bwt_restore_bwt(argv[optind]);
|
|
|
181 bwt_cal_sa(bwt, sa_intv);
|
|
|
182 bwt_dump_sa(argv[optind+1], bwt);
|
|
|
183 bwt_destroy(bwt);
|
|
|
184 return 0;
|
|
|
185 }
|
|
|
186
|
|
|
187 int bwa_index(int argc, char *argv[]) // the "index" command
|
|
|
188 {
|
|
|
189 extern void bwa_pac_rev_core(const char *fn, const char *fn_rev);
|
|
|
190
|
|
|
191 char *prefix = 0, *str, *str2, *str3;
|
|
|
192 int c, algo_type = 0, is_64 = 0;
|
|
|
193 clock_t t;
|
|
|
194 int64_t l_pac;
|
|
|
195
|
|
|
196 while ((c = getopt(argc, argv, "6a:p:")) >= 0) {
|
|
|
197 switch (c) {
|
|
|
198 case 'a': // if -a is not set, algo_type will be determined later
|
|
|
199 if (strcmp(optarg, "div") == 0) algo_type = 1;
|
|
|
200 else if (strcmp(optarg, "bwtsw") == 0) algo_type = 2;
|
|
|
201 else if (strcmp(optarg, "is") == 0) algo_type = 3;
|
|
|
202 else err_fatal(__func__, "unknown algorithm: '%s'.", optarg);
|
|
|
203 break;
|
|
|
204 case 'p': prefix = strdup(optarg); break;
|
|
|
205 case '6': is_64 = 1; break;
|
|
|
206 default: return 1;
|
|
|
207 }
|
|
|
208 }
|
|
|
209
|
|
|
210 if (optind + 1 > argc) {
|
|
|
211 fprintf(stderr, "\n");
|
|
|
212 fprintf(stderr, "Usage: bwa index [-a bwtsw|is] [-c] <in.fasta>\n\n");
|
|
|
213 fprintf(stderr, "Options: -a STR BWT construction algorithm: bwtsw or is [auto]\n");
|
|
|
214 fprintf(stderr, " -p STR prefix of the index [same as fasta name]\n");
|
|
|
215 fprintf(stderr, " -6 index files named as <in.fasta>.64.* instead of <in.fasta>.* \n");
|
|
|
216 fprintf(stderr, "\n");
|
|
|
217 fprintf(stderr, "Warning: `-a bwtsw' does not work for short genomes, while `-a is' and\n");
|
|
|
218 fprintf(stderr, " `-a div' do not work not for long genomes. Please choose `-a'\n");
|
|
|
219 fprintf(stderr, " according to the length of the genome.\n\n");
|
|
|
220 return 1;
|
|
|
221 }
|
|
|
222 if (prefix == 0) {
|
|
|
223 prefix = malloc(strlen(argv[optind]) + 4);
|
|
|
224 strcpy(prefix, argv[optind]);
|
|
|
225 if (is_64) strcat(prefix, ".64");
|
|
|
226 }
|
|
|
227 str = (char*)calloc(strlen(prefix) + 10, 1);
|
|
|
228 str2 = (char*)calloc(strlen(prefix) + 10, 1);
|
|
|
229 str3 = (char*)calloc(strlen(prefix) + 10, 1);
|
|
|
230
|
|
|
231 { // nucleotide indexing
|
|
|
232 gzFile fp = xzopen(argv[optind], "r");
|
|
|
233 t = clock();
|
|
|
234 fprintf(stderr, "[bwa_index] Pack FASTA... ");
|
|
|
235 l_pac = bns_fasta2bntseq(fp, prefix, 0);
|
|
|
236 fprintf(stderr, "%.2f sec\n", (float)(clock() - t) / CLOCKS_PER_SEC);
|
|
|
237 err_gzclose(fp);
|
|
|
238 }
|
|
|
239 if (algo_type == 0) algo_type = l_pac > 50000000? 2 : 3; // set the algorithm for generating BWT
|
|
|
240 {
|
|
|
241 strcpy(str, prefix); strcat(str, ".pac");
|
|
|
242 strcpy(str2, prefix); strcat(str2, ".bwt");
|
|
|
243 t = clock();
|
|
|
244 fprintf(stderr, "[bwa_index] Construct BWT for the packed sequence...\n");
|
|
|
245 if (algo_type == 2) bwt_bwtgen(str, str2);
|
|
|
246 else if (algo_type == 1 || algo_type == 3) {
|
|
|
247 bwt_t *bwt;
|
|
|
248 bwt = bwt_pac2bwt(str, algo_type == 3);
|
|
|
249 bwt_dump_bwt(str2, bwt);
|
|
|
250 bwt_destroy(bwt);
|
|
|
251 }
|
|
|
252 fprintf(stderr, "[bwa_index] %.2f seconds elapse.\n", (float)(clock() - t) / CLOCKS_PER_SEC);
|
|
|
253 }
|
|
|
254 {
|
|
|
255 bwt_t *bwt;
|
|
|
256 strcpy(str, prefix); strcat(str, ".bwt");
|
|
|
257 t = clock();
|
|
|
258 fprintf(stderr, "[bwa_index] Update BWT... ");
|
|
|
259 bwt = bwt_restore_bwt(str);
|
|
|
260 bwt_bwtupdate_core(bwt);
|
|
|
261 bwt_dump_bwt(str, bwt);
|
|
|
262 bwt_destroy(bwt);
|
|
|
263 fprintf(stderr, "%.2f sec\n", (float)(clock() - t) / CLOCKS_PER_SEC);
|
|
|
264 }
|
|
|
265 {
|
|
|
266 gzFile fp = xzopen(argv[optind], "r");
|
|
|
267 t = clock();
|
|
|
268 fprintf(stderr, "[bwa_index] Pack forward-only FASTA... ");
|
|
|
269 l_pac = bns_fasta2bntseq(fp, prefix, 1);
|
|
|
270 fprintf(stderr, "%.2f sec\n", (float)(clock() - t) / CLOCKS_PER_SEC);
|
|
|
271 err_gzclose(fp);
|
|
|
272 }
|
|
|
273 {
|
|
|
274 bwt_t *bwt;
|
|
|
275 strcpy(str, prefix); strcat(str, ".bwt");
|
|
|
276 strcpy(str3, prefix); strcat(str3, ".sa");
|
|
|
277 t = clock();
|
|
|
278 fprintf(stderr, "[bwa_index] Construct SA from BWT and Occ... ");
|
|
|
279 bwt = bwt_restore_bwt(str);
|
|
|
280 bwt_cal_sa(bwt, 32);
|
|
|
281 bwt_dump_sa(str3, bwt);
|
|
|
282 bwt_destroy(bwt);
|
|
|
283 fprintf(stderr, "%.2f sec\n", (float)(clock() - t) / CLOCKS_PER_SEC);
|
|
|
284 }
|
|
|
285 free(str3); free(str2); free(str); free(prefix);
|
|
|
286 return 0;
|
|
|
287 }
|
