forked from genome/ibwa
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathdbset.c
339 lines (286 loc) · 9.31 KB
/
dbset.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
#include "dbset.h"
#include "bwaremap.h"
#include "kstring.h"
#include "utils.h"
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <limits.h>
extern char *bwa_rg_line;
extern char *bwa_rg_id;
bntseq_t *bwa_open_nt(const char *prefix);
/* binary search through sequences to find the one containing pos */
int coord2idx(const dbset_t *dbs, int64_t pos) {
int left = 0, right = dbs->count;
int mid = 0;
if (pos > dbs->l_pac)
return -1;
while (left < right) {
mid = (left + right) >> 1;
if (pos > dbs->db[mid]->offset) {
if (mid == dbs->count-1) break;
if (pos < dbs->db[mid+1]->offset) break;
left = mid + 1;
} else if (pos < dbs->db[mid]->offset) {
right = mid;
} else {
break;
}
}
assert(mid < dbs->count);
return mid;
}
static bwtdb_t *bwtdb_load(const char *prefix) {
bwtdb_t *db = calloc(1, sizeof(bwtdb_t));
db->prefix = prefix;
db->bwtcache = bwtcache_create();
return db;
}
static void bwtdb_load_sa(bwtdb_t *db, int which) {
char path[PATH_MAX];
const char *bwt_suffix = ".bwt";
const char *sa_suffix = ".sa";
if (db->bwt[which] != NULL)
return;
if (which == 1) {
bwt_suffix = ".rbwt";
sa_suffix = ".rsa";
}
strcpy(path, db->prefix); strcat(path, bwt_suffix);
db->bwt[which] = bwt_restore_bwt(path);
strcpy(path, db->prefix); strcat(path, sa_suffix);
bwt_restore_sa(path, db->bwt[which]);
}
static void bwtdb_unload_sa(bwtdb_t *db, int which) {
if (db->bwt[which])
bwt_destroy(db->bwt[which]);
db->bwt[which] = NULL;
}
static void bwtdb_destroy(bwtdb_t *db) {
bwtdb_unload_sa(db, 0);
bwtdb_unload_sa(db, 1);
bwtcache_destroy(db->bwtcache);
free(db);
}
static seq_t *seq_restore(const char *prefix, const char *extension, int remap) {
char path[PATH_MAX];
char rmpath[PATH_MAX];
strcat(strcpy(path, prefix), extension);
strcat(strcpy(rmpath, path), ".remap");
seq_t *s = calloc(1, sizeof(seq_t));
s->bns = bns_restore(path);
if (remap) {
s->remap = load_remappings(s, rmpath);
if (s->remap < 0) {
fprintf(stderr, "Fatal error loading sequence mappings from %s\n", rmpath);
exit(1);
} else if (s->remap) {
fprintf(stderr, " - Remapping enabled for sequence %s\n", prefix);
}
} else {
s->remap = 0;
}
return s;
}
static void seq_load_pac(seq_t *s) {
assert(s->data == NULL);
s->data = (ubyte_t*)calloc(s->bns->l_pac/4+1, 1);
rewind(s->bns->fp_pac);
fread(s->data, 1, s->bns->l_pac/4+1, s->bns->fp_pac);
}
static void seq_unload_pac(seq_t *s) {
if (s->data) free(s->data);
s->data = NULL;
}
static void seq_destroy(seq_t *s) {
int i;
if (!s) return;
seq_unload_pac(s);
if (s->mappings) {
for (i = 0; i < s->bns->n_seqs; ++i) {
if (s->mappings[i]) {
read_mapping_destroy(&s->mappings[i]->map);
free(s->mappings[i]);
}
}
free(s->mappings);
}
if (s->bns) bns_destroy(s->bns);
free(s);
}
dbset_t *dbset_restore(int count, const char **prefixes, int mode, int preload, int remap) {
int i;
dbset_t *dbs = calloc(1, sizeof(dbset_t));
dbs->count = count;
dbs->db = (bwtdb_t**)calloc(count, sizeof(bwtdb_t*));
dbs->bns = calloc(count, sizeof(seq_t*));
dbs->ntbns = calloc(count, sizeof(seq_t*));
dbs->preload = preload;
dbs->color_space = !(mode & BWA_MODE_COMPREAD);
for (i = 0; i < count; ++i) {
dbs->db[i] = bwtdb_load(prefixes[i]);
dbs->db[i]->offset = dbs->l_pac;
dbs->bns[i] = seq_restore(prefixes[i], "", remap);
dbs->db[i]->bns = dbs->bns[i];
dbs->l_pac += dbs->bns[i]->bns->l_pac;
/* TODO: get rid of these, we just need the length of the bwt */
bwtdb_load_sa(dbs->db[i], 0);
bwtdb_load_sa(dbs->db[i], 1);
dbs->total_bwt_seq_len[0] += dbs->db[i]->bwt[0]->seq_len;
dbs->total_bwt_seq_len[1] += dbs->db[i]->bwt[1]->seq_len;
if (dbs->color_space) {
dbs->ntbns[i] = seq_restore(prefixes[i], ".nt", remap);
dbs->db[i]->ntbns = dbs->ntbns[i];
dbs->color_space = 1;
} else if (preload) {
seq_load_pac(dbs->bns[i]);
bwtdb_load_sa(dbs->db[i], 0);
bwtdb_load_sa(dbs->db[i], 1);
}
}
return dbs;
}
void dbset_destroy(dbset_t *dbs) {
int i;
for (i = 0; i < dbs->count; ++i) {
bwtdb_destroy(dbs->db[i]);
seq_destroy(dbs->bns[i]);
seq_destroy(dbs->ntbns[i]);
}
free(dbs->db);
free(dbs->bns);
free(dbs->ntbns);
free(dbs);
}
void dbset_load_sa(dbset_t *dbs, int which) {
int i;
assert(which == 0 || which == 1);
if (dbs->preload) return;
for (i = 0; i < dbs->count; ++i) {
if (dbs->db[i]->bwt[which] == NULL)
bwtdb_load_sa(dbs->db[which], which);
}
}
void dbset_unload_sa(dbset_t *dbs, int which) {
int i;
if (dbs->preload) return;
for (i = 0; i < dbs->count; ++i)
bwtdb_unload_sa(dbs->db[i], which);
}
void dbset_load_pac(dbset_t *dbs) {
int i;
if (dbs->preload)
return;
for (i = 0; i < dbs->count; ++i)
if (dbs->bns[i]->data == NULL)
seq_load_pac(dbs->bns[i]);
}
void dbset_unload_pac(dbset_t *dbs) {
int i;
if (dbs->preload)
return;
for (i = 0; i < dbs->count; ++i)
seq_unload_pac(dbs->bns[i]);
}
void dbset_load_ntpac(dbset_t *dbs) {
int i;
for (i = 0; i < dbs->count; ++i)
if (dbs->ntbns[i]->data == NULL)
seq_load_pac(dbs->ntbns[i]);
}
void dbset_unload_ntpac(dbset_t *dbs) {
int i;
if (dbs->preload)
return;
for (i = 0; i < dbs->count; ++i)
seq_unload_pac(dbs->ntbns[i]);
}
uint64_t bwtdb_sa2seq(const bwtdb_t *db, int strand, uint32_t sa, uint32_t seq_len) {
if (strand) {
return db->offset + bwt_sa(db->bwt[0], sa);
} else {
return db->offset + db->bwt[1]->seq_len - (bwt_sa(db->bwt[1], sa) + seq_len);
}
}
int dbset_coor_pac2real(const dbset_t *dbs, int64_t pac_coor, int len, int32_t *real_seq,
const bntseq_t **bns, uint64_t *offset)
{
int idx = coord2idx(dbs, pac_coor);
*bns = dbs->bns[idx]->bns;
*offset = dbs->db[idx]->offset;
return bns_coor_pac2real(*bns, pac_coor - *offset, len, real_seq);
}
poslist_t bwtdb_cached_sa2seq(const bwtdb_t *db, const bwt_aln1_t* aln, uint32_t seq_len) {
return bwt_cached_sa(db->offset, db->bwtcache, (const bwt_t **const)db->bwt, aln, seq_len);
}
uint32_t dbset_extract_remapped(const dbset_t *dbs, seq_t **seqs, uint32_t dbidx, int32_t seqid, ubyte_t* ref_seq, uint64_t beg, uint32_t len) {
bwtdb_t *db = dbs->db[dbidx];
uint64_t seq_begin;
uint32_t total = 0;
int status;
bntann1_t *ann;
if (seqid < 0 || !db->bns->remap) {
return dbset_extract_sequence(dbs, seqs, ref_seq, beg, len);
}
ann = db->bns->bns->anns + seqid;
seq_begin = db->offset + ann->offset;
if (beg < seq_begin) {
uint64_t remapped_begin = bwa_remap_position_with_seqid(db->bns, dbs->db[0]->bns->bns, ann->offset, seqid, &status);
uint64_t sublen = seq_begin - beg;
uint64_t offset = remapped_begin - sublen;
if (sublen > remapped_begin || status == 0)
err_fatal(__func__, "request too far ahead of remapped region");
total += dbset_extract_sequence(dbs, seqs, &ref_seq[total], offset, sublen);
}
if (total < len) {
uint32_t sublen = len - total;
if (sublen > ann->len)
sublen = ann->len;
total += dbset_extract_sequence(dbs, seqs, &ref_seq[total], beg, sublen);
}
if (total < len) {
uint64_t remapped_end = bwa_remap_position_with_seqid(db->bns, dbs->db[0]->bns->bns, ann->offset + ann->len-1, seqid, &status)+1;
if (status == 0) {
err_fatal(__func__, "request too far ahead of remapped region");
}
total += dbset_extract_sequence(dbs, seqs, &ref_seq[total], remapped_end, len-total);
}
if (total != len) {
err_fatal(__func__, "logic error: got %lu bases instead of %lu\n", total, len);
}
return total;
}
uint32_t dbset_extract_sequence(const dbset_t *dbs, seq_t **seqs, ubyte_t* ref_seq, uint64_t beg, uint32_t len) {
uint32_t total = 0;
while (total < len) {
int64_t idx, pos;
seq_t *s;
bwtdb_t *db;
if (beg >= dbs->l_pac) break;
idx = coord2idx(dbs, beg);
s = seqs[idx];
db = dbs->db[idx];
pos = beg - db->offset;
/* TODO: this is wrong */
while (pos < s->bns->l_pac && total < len) {
ref_seq[total++] = bns_pac(s->data, pos);
++pos; /* bns_pac is a macro referencing idx more than once */
}
beg = pos + db->offset;
}
return total;
}
void dbset_print_sam_SQ(const dbset_t *dbs) {
int i, j;
for (i = 0; i < dbs->count; ++i) {
seq_t *s = dbs->bns[i];
for (j = 0; j < s->bns->n_seqs; ++j) {
if (!s->mappings || s->mappings[j] == NULL) {
printf("@SQ\tSN:%s\tLN:%d\n",
s->bns->anns[j].name, s->bns->anns[j].len);
}
}
}
if (bwa_rg_line) printf("%s\n", bwa_rg_line);
}