← Back to branch summary

~ubuntu-branches/ubuntu/quantal/genometools/quantal-backports

« back to all changes in this revision

Viewing changes to src/external/samtools-0.1.18/bam_sort.c

  • Committer: Package Import Robot
  • Author(s): Sascha Steinbiss
  • Date: 2012-07-09 14:10:23 UTC
  • Revision ID: package-import@ubuntu.com-20120709141023-juuu4spm6chqsf9o
Tags: upstream-1.4.1
Import upstream version 1.4.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/external/samtools-0.1.18/bam_sort.c
 
1
#include <stdlib.h>
 
2
#include <ctype.h>
 
3
#include <assert.h>
 
4
#include <errno.h>
 
5
#include <stdio.h>
 
6
#include <string.h>
 
7
#include <unistd.h>
 
8
#include "bam.h"
 
9
#include "ksort.h"
 
10
 
 
11
static int g_is_by_qname = 0;
 
12
 
 
13
static inline int strnum_cmp(const char *a, const char *b)
 
14
{
 
15
        char *pa, *pb;
 
16
        pa = (char*)a; pb = (char*)b;
 
17
        while (*pa && *pb) {
 
18
                if (isdigit(*pa) && isdigit(*pb)) {
 
19
                        long ai, bi;
 
20
                        ai = strtol(pa, &pa, 10);
 
21
                        bi = strtol(pb, &pb, 10);
 
22
                        if (ai != bi) return ai<bi? -1 : ai>bi? 1 : 0;
 
23
                } else {
 
24
                        if (*pa != *pb) break;
 
25
                        ++pa; ++pb;
 
26
                }
 
27
        }
 
28
        if (*pa == *pb)
 
29
                return (pa-a) < (pb-b)? -1 : (pa-a) > (pb-b)? 1 : 0;
 
30
        return *pa<*pb? -1 : *pa>*pb? 1 : 0;
 
31
}
 
32
 
 
33
#define HEAP_EMPTY 0xffffffffffffffffull
 
34
 
 
35
typedef struct {
 
36
        int i;
 
37
        uint64_t pos, idx;
 
38
        bam1_t *b;
 
39
} heap1_t;
 
40
 
 
41
#define __pos_cmp(a, b) ((a).pos > (b).pos || ((a).pos == (b).pos && ((a).i > (b).i || ((a).i == (b).i && (a).idx > (b).idx))))
 
42
 
 
43
static inline int heap_lt(const heap1_t a, const heap1_t b)
 
44
{
 
45
        if (g_is_by_qname) {
 
46
                int t;
 
47
                if (a.b == 0 || b.b == 0) return a.b == 0? 1 : 0;
 
48
                t = strnum_cmp(bam1_qname(a.b), bam1_qname(b.b));
 
49
                return (t > 0 || (t == 0 && __pos_cmp(a, b)));
 
50
        } else return __pos_cmp(a, b);
 
51
}
 
52
 
 
53
KSORT_INIT(heap, heap1_t, heap_lt)
 
54
 
 
55
static void swap_header_targets(bam_header_t *h1, bam_header_t *h2)
 
56
{
 
57
        bam_header_t t;
 
58
        t.n_targets = h1->n_targets, h1->n_targets = h2->n_targets, h2->n_targets = t.n_targets;
 
59
        t.target_name = h1->target_name, h1->target_name = h2->target_name, h2->target_name = t.target_name;
 
60
        t.target_len = h1->target_len, h1->target_len = h2->target_len, h2->target_len = t.target_len;
 
61
}
 
62
 
 
63
static void swap_header_text(bam_header_t *h1, bam_header_t *h2)
 
64
{
 
65
        int tempi;
 
66
        char *temps;
 
67
        tempi = h1->l_text, h1->l_text = h2->l_text, h2->l_text = tempi;
 
68
        temps = h1->text, h1->text = h2->text, h2->text = temps;
 
69
}
 
70
 
 
71
#define MERGE_RG     1
 
72
#define MERGE_UNCOMP 2
 
73
#define MERGE_LEVEL1 4
 
74
#define MERGE_FORCE  8
 
75
 
 
76
/*!
 
77
  @abstract    Merge multiple sorted BAM.
 
78
  @param  is_by_qname whether to sort by query name
 
79
  @param  out  output BAM file name
 
80
  @param  headers  name of SAM file from which to copy '@' header lines,
 
81
                   or NULL to copy them from the first file to be merged
 
82
  @param  n    number of files to be merged
 
83
  @param  fn   names of files to be merged
 
84
 
 
85
  @discussion Padding information may NOT correctly maintained. This
 
86
  function is NOT thread safe.
 
87
 */
 
88
int bam_merge_core(int by_qname, const char *out, const char *headers, int n, char * const *fn,
 
89
                                        int flag, const char *reg)
 
90
{
 
91
        bamFile fpout, *fp;
 
92
        heap1_t *heap;
 
93
        bam_header_t *hout = 0;
 
94
        bam_header_t *hheaders = NULL;
 
95
        int i, j, *RG_len = 0;
 
96
        uint64_t idx = 0;
 
97
        char **RG = 0;
 
98
        bam_iter_t *iter = 0;
 
99
 
 
100
        if (headers) {
 
101
                tamFile fpheaders = sam_open(headers);
 
102
                if (fpheaders == 0) {
 
103
                        const char *message = strerror(errno);
 
104
                        fprintf(stderr, "[bam_merge_core] cannot open '%s': %s\n", headers, message);
 
105
                        return -1;
 
106
                }
 
107
                hheaders = sam_header_read(fpheaders);
 
108
                sam_close(fpheaders);
 
109
        }
 
110
 
 
111
        g_is_by_qname = by_qname;
 
112
        fp = (bamFile*)calloc(n, sizeof(bamFile));
 
113
        heap = (heap1_t*)calloc(n, sizeof(heap1_t));
 
114
        iter = (bam_iter_t*)calloc(n, sizeof(bam_iter_t));
 
115
        // prepare RG tag
 
116
        if (flag & MERGE_RG) {
 
117
                RG = (char**)calloc(n, sizeof(void*));
 
118
                RG_len = (int*)calloc(n, sizeof(int));
 
119
                for (i = 0; i != n; ++i) {
 
120
                        int l = strlen(fn[i]);
 
121
                        const char *s = fn[i];
 
122
                        if (l > 4 && strcmp(s + l - 4, ".bam") == 0) l -= 4;
 
123
                        for (j = l - 1; j >= 0; --j) if (s[j] == '/') break;
 
124
                        ++j; l -= j;
 
125
                        RG[i] = calloc(l + 1, 1);
 
126
                        RG_len[i] = l;
 
127
                        strncpy(RG[i], s + j, l);
 
128
                }
 
129
        }
 
130
        // read the first
 
131
        for (i = 0; i != n; ++i) {
 
132
                bam_header_t *hin;
 
133
                fp[i] = bam_open(fn[i], "r");
 
134
                if (fp[i] == 0) {
 
135
                        int j;
 
136
                        fprintf(stderr, "[bam_merge_core] fail to open file %s\n", fn[i]);
 
137
                        for (j = 0; j < i; ++j) bam_close(fp[j]);
 
138
                        free(fp); free(heap);
 
139
                        // FIXME: possible memory leak
 
140
                        return -1;
 
141
                }
 
142
                hin = bam_header_read(fp[i]);
 
143
                if (i == 0) { // the first BAM
 
144
                        hout = hin;
 
145
                } else { // validate multiple baf
 
146
                        int min_n_targets = hout->n_targets;
 
147
                        if (hin->n_targets < min_n_targets) min_n_targets = hin->n_targets;
 
148
 
 
149
                        for (j = 0; j < min_n_targets; ++j)
 
150
                                if (strcmp(hout->target_name[j], hin->target_name[j]) != 0) {
 
151
                                        fprintf(stderr, "[bam_merge_core] different target sequence name: '%s' != '%s' in file '%s'\n",
 
152
                                                        hout->target_name[j], hin->target_name[j], fn[i]);
 
153
                                        return -1;
 
154
                                }
 
155
 
 
156
                        // If this input file has additional target reference sequences,
 
157
                        // add them to the headers to be output
 
158
                        if (hin->n_targets > hout->n_targets) {
 
159
                                swap_header_targets(hout, hin);
 
160
                                // FIXME Possibly we should also create @SQ text headers
 
161
                                // for the newly added reference sequences
 
162
                        }
 
163
 
 
164
                        bam_header_destroy(hin);
 
165
                }
 
166
        }
 
167
 
 
168
        if (hheaders) {
 
169
                // If the text headers to be swapped in include any @SQ headers,
 
170
                // check that they are consistent with the existing binary list
 
171
                // of reference information.
 
172
                if (hheaders->n_targets > 0) {
 
173
                        if (hout->n_targets != hheaders->n_targets) {
 
174
                                fprintf(stderr, "[bam_merge_core] number of @SQ headers in '%s' differs from number of target sequences\n", headers);
 
175
                                if (!reg) return -1;
 
176
                        }
 
177
                        for (j = 0; j < hout->n_targets; ++j)
 
178
                                if (strcmp(hout->target_name[j], hheaders->target_name[j]) != 0) {
 
179
                                        fprintf(stderr, "[bam_merge_core] @SQ header '%s' in '%s' differs from target sequence\n", hheaders->target_name[j], headers);
 
180
                                        if (!reg) return -1;
 
181
                                }
 
182
                }
 
183
 
 
184
                swap_header_text(hout, hheaders);
 
185
                bam_header_destroy(hheaders);
 
186
        }
 
187
 
 
188
        if (reg) {
 
189
                int tid, beg, end;
 
190
                if (bam_parse_region(hout, reg, &tid, &beg, &end) < 0) {
 
191
                        fprintf(stderr, "[%s] Malformated region string or undefined reference name\n", __func__);
 
192
                        return -1;
 
193
                }
 
194
                for (i = 0; i < n; ++i) {
 
195
                        bam_index_t *idx;
 
196
                        idx = bam_index_load(fn[i]);
 
197
                        iter[i] = bam_iter_query(idx, tid, beg, end);
 
198
                        bam_index_destroy(idx);
 
199
                }
 
200
        }
 
201
 
 
202
        for (i = 0; i < n; ++i) {
 
203
                heap1_t *h = heap + i;
 
204
                h->i = i;
 
205
                h->b = (bam1_t*)calloc(1, sizeof(bam1_t));
 
206
                if (bam_iter_read(fp[i], iter[i], h->b) >= 0) {
 
207
                        h->pos = ((uint64_t)h->b->core.tid<<32) | (uint32_t)((int32_t)h->b->core.pos+1)<<1 | bam1_strand(h->b);
 
208
                        h->idx = idx++;
 
209
                }
 
210
                else h->pos = HEAP_EMPTY;
 
211
        }
 
212
        if (flag & MERGE_UNCOMP) fpout = strcmp(out, "-")? bam_open(out, "wu") : bam_dopen(fileno(stdout), "wu");
 
213
        else if (flag & MERGE_LEVEL1) fpout = strcmp(out, "-")? bam_open(out, "w1") : bam_dopen(fileno(stdout), "w1");
 
214
        else fpout = strcmp(out, "-")? bam_open(out, "w") : bam_dopen(fileno(stdout), "w");
 
215
        if (fpout == 0) {
 
216
                fprintf(stderr, "[%s] fail to create the output file.\n", __func__);
 
217
                return -1;
 
218
        }
 
219
        bam_header_write(fpout, hout);
 
220
        bam_header_destroy(hout);
 
221
 
 
222
        ks_heapmake(heap, n, heap);
 
223
        while (heap->pos != HEAP_EMPTY) {
 
224
                bam1_t *b = heap->b;
 
225
                if (flag & MERGE_RG) {
 
226
                        uint8_t *rg = bam_aux_get(b, "RG");
 
227
                        if (rg) bam_aux_del(b, rg);
 
228
                        bam_aux_append(b, "RG", 'Z', RG_len[heap->i] + 1, (uint8_t*)RG[heap->i]);
 
229
                }
 
230
                bam_write1_core(fpout, &b->core, b->data_len, b->data);
 
231
                if ((j = bam_iter_read(fp[heap->i], iter[heap->i], b)) >= 0) {
 
232
                        heap->pos = ((uint64_t)b->core.tid<<32) | (uint32_t)((int)b->core.pos+1)<<1 | bam1_strand(b);
 
233
                        heap->idx = idx++;
 
234
                } else if (j == -1) {
 
235
                        heap->pos = HEAP_EMPTY;
 
236
                        free(heap->b->data); free(heap->b);
 
237
                        heap->b = 0;
 
238
                } else fprintf(stderr, "[bam_merge_core] '%s' is truncated. Continue anyway.\n", fn[heap->i]);
 
239
                ks_heapadjust(heap, 0, n, heap);
 
240
        }
 
241
 
 
242
        if (flag & MERGE_RG) {
 
243
                for (i = 0; i != n; ++i) free(RG[i]);
 
244
                free(RG); free(RG_len);
 
245
        }
 
246
        for (i = 0; i != n; ++i) {
 
247
                bam_iter_destroy(iter[i]);
 
248
                bam_close(fp[i]);
 
249
        }
 
250
        bam_close(fpout);
 
251
        free(fp); free(heap); free(iter);
 
252
        return 0;
 
253
}
 
254
 
 
255
int bam_merge(int argc, char *argv[])
 
256
{
 
257
        int c, is_by_qname = 0, flag = 0, ret = 0;
 
258
        char *fn_headers = NULL, *reg = 0;
 
259
 
 
260
        while ((c = getopt(argc, argv, "h:nru1R:f")) >= 0) {
 
261
                switch (c) {
 
262
                case 'r': flag |= MERGE_RG; break;
 
263
                case 'f': flag |= MERGE_FORCE; break;
 
264
                case 'h': fn_headers = strdup(optarg); break;
 
265
                case 'n': is_by_qname = 1; break;
 
266
                case '1': flag |= MERGE_LEVEL1; break;
 
267
                case 'u': flag |= MERGE_UNCOMP; break;
 
268
                case 'R': reg = strdup(optarg); break;
 
269
                }
 
270
        }
 
271
        if (optind + 2 >= argc) {
 
272
                fprintf(stderr, "\n");
 
273
                fprintf(stderr, "Usage:   samtools merge [-nr] [-h inh.sam] <out.bam> <in1.bam> <in2.bam> [...]\n\n");
 
274
                fprintf(stderr, "Options: -n       sort by read names\n");
 
275
                fprintf(stderr, "         -r       attach RG tag (inferred from file names)\n");
 
276
                fprintf(stderr, "         -u       uncompressed BAM output\n");
 
277
                fprintf(stderr, "         -f       overwrite the output BAM if exist\n");
 
278
                fprintf(stderr, "         -1       compress level 1\n");
 
279
                fprintf(stderr, "         -R STR   merge file in the specified region STR [all]\n");
 
280
                fprintf(stderr, "         -h FILE  copy the header in FILE to <out.bam> [in1.bam]\n\n");
 
281
                fprintf(stderr, "Note: Samtools' merge does not reconstruct the @RG dictionary in the header. Users\n");
 
282
                fprintf(stderr, "      must provide the correct header with -h, or uses Picard which properly maintains\n");
 
283
                fprintf(stderr, "      the header dictionary in merging.\n\n");
 
284
                return 1;
 
285
        }
 
286
        if (!(flag & MERGE_FORCE) && strcmp(argv[optind], "-")) {
 
287
                FILE *fp = fopen(argv[optind], "rb");
 
288
                if (fp != NULL) {
 
289
                        fclose(fp);
 
290
                        fprintf(stderr, "[%s] File '%s' exists. Please apply '-f' to overwrite. Abort.\n", __func__, argv[optind]);
 
291
                        return 1;
 
292
                }
 
293
        }
 
294
        if (bam_merge_core(is_by_qname, argv[optind], fn_headers, argc - optind - 1, argv + optind + 1, flag, reg) < 0) ret = 1;
 
295
        free(reg);
 
296
        free(fn_headers);
 
297
        return ret;
 
298
}
 
299
 
 
300
typedef bam1_t *bam1_p;
 
301
 
 
302
static inline int bam1_lt(const bam1_p a, const bam1_p b)
 
303
{
 
304
        if (g_is_by_qname) {
 
305
                int t = strnum_cmp(bam1_qname(a), bam1_qname(b));
 
306
                return (t < 0 || (t == 0 && (((uint64_t)a->core.tid<<32|(a->core.pos+1)) < ((uint64_t)b->core.tid<<32|(b->core.pos+1)))));
 
307
        } else return (((uint64_t)a->core.tid<<32|(a->core.pos+1)) < ((uint64_t)b->core.tid<<32|(b->core.pos+1)));
 
308
}
 
309
KSORT_INIT(sort, bam1_p, bam1_lt)
 
310
 
 
311
static void sort_blocks(int n, int k, bam1_p *buf, const char *prefix, const bam_header_t *h, int is_stdout)
 
312
{
 
313
        char *name, mode[3];
 
314
        int i;
 
315
        bamFile fp;
 
316
        ks_mergesort(sort, k, buf, 0);
 
317
        name = (char*)calloc(strlen(prefix) + 20, 1);
 
318
        if (n >= 0) {
 
319
                sprintf(name, "%s.%.4d.bam", prefix, n);
 
320
                strcpy(mode, "w1");
 
321
        } else {
 
322
                sprintf(name, "%s.bam", prefix);
 
323
                strcpy(mode, "w");
 
324
        }
 
325
        fp = is_stdout? bam_dopen(fileno(stdout), mode) : bam_open(name, mode);
 
326
        if (fp == 0) {
 
327
                fprintf(stderr, "[sort_blocks] fail to create file %s.\n", name);
 
328
                free(name);
 
329
                // FIXME: possible memory leak
 
330
                return;
 
331
        }
 
332
        free(name);
 
333
        bam_header_write(fp, h);
 
334
        for (i = 0; i < k; ++i)
 
335
                bam_write1_core(fp, &buf[i]->core, buf[i]->data_len, buf[i]->data);
 
336
        bam_close(fp);
 
337
}
 
338
 
 
339
/*!
 
340
  @abstract Sort an unsorted BAM file based on the chromosome order
 
341
  and the leftmost position of an alignment
 
342
 
 
343
  @param  is_by_qname whether to sort by query name
 
344
  @param  fn       name of the file to be sorted
 
345
  @param  prefix   prefix of the output and the temporary files; upon
 
346
                           sucessess, prefix.bam will be written.
 
347
  @param  max_mem  approxiate maximum memory (very inaccurate)
 
348
 
 
349
  @discussion It may create multiple temporary subalignment files
 
350
  and then merge them by calling bam_merge_core(). This function is
 
351
  NOT thread safe.
 
352
 */
 
353
void bam_sort_core_ext(int is_by_qname, const char *fn, const char *prefix, size_t max_mem, int is_stdout)
 
354
{
 
355
        int n, ret, k, i;
 
356
        size_t mem;
 
357
        bam_header_t *header;
 
358
        bamFile fp;
 
359
        bam1_t *b, **buf;
 
360
 
 
361
        g_is_by_qname = is_by_qname;
 
362
        n = k = 0; mem = 0;
 
363
        fp = strcmp(fn, "-")? bam_open(fn, "r") : bam_dopen(fileno(stdin), "r");
 
364
        if (fp == 0) {
 
365
                fprintf(stderr, "[bam_sort_core] fail to open file %s\n", fn);
 
366
                return;
 
367
        }
 
368
        header = bam_header_read(fp);
 
369
        buf = (bam1_t**)calloc(max_mem / BAM_CORE_SIZE, sizeof(bam1_t*));
 
370
        // write sub files
 
371
        for (;;) {
 
372
                if (buf[k] == 0) buf[k] = (bam1_t*)calloc(1, sizeof(bam1_t));
 
373
                b = buf[k];
 
374
                if ((ret = bam_read1(fp, b)) < 0) break;
 
375
                mem += ret;
 
376
                ++k;
 
377
                if (mem >= max_mem) {
 
378
                        sort_blocks(n++, k, buf, prefix, header, 0);
 
379
                        mem = 0; k = 0;
 
380
                }
 
381
        }
 
382
        if (ret != -1)
 
383
                fprintf(stderr, "[bam_sort_core] truncated file. Continue anyway.\n");
 
384
        if (n == 0) sort_blocks(-1, k, buf, prefix, header, is_stdout);
 
385
        else { // then merge
 
386
                char **fns, *fnout;
 
387
                fprintf(stderr, "[bam_sort_core] merging from %d files...\n", n+1);
 
388
                sort_blocks(n++, k, buf, prefix, header, 0);
 
389
                fnout = (char*)calloc(strlen(prefix) + 20, 1);
 
390
                if (is_stdout) sprintf(fnout, "-");
 
391
                else sprintf(fnout, "%s.bam", prefix);
 
392
                fns = (char**)calloc(n, sizeof(char*));
 
393
                for (i = 0; i < n; ++i) {
 
394
                        fns[i] = (char*)calloc(strlen(prefix) + 20, 1);
 
395
                        sprintf(fns[i], "%s.%.4d.bam", prefix, i);
 
396
                }
 
397
                bam_merge_core(is_by_qname, fnout, 0, n, fns, 0, 0);
 
398
                free(fnout);
 
399
                for (i = 0; i < n; ++i) {
 
400
                        unlink(fns[i]);
 
401
                        free(fns[i]);
 
402
                }
 
403
                free(fns);
 
404
        }
 
405
        for (k = 0; k < max_mem / BAM_CORE_SIZE; ++k) {
 
406
                if (buf[k]) {
 
407
                        free(buf[k]->data);
 
408
                        free(buf[k]);
 
409
                }
 
410
        }
 
411
        free(buf);
 
412
        bam_header_destroy(header);
 
413
        bam_close(fp);
 
414
}
 
415
 
 
416
void bam_sort_core(int is_by_qname, const char *fn, const char *prefix, size_t max_mem)
 
417
{
 
418
        bam_sort_core_ext(is_by_qname, fn, prefix, max_mem, 0);
 
419
}
 
420
 
 
421
int bam_sort(int argc, char *argv[])
 
422
{
 
423
        size_t max_mem = 500000000;
 
424
        int c, is_by_qname = 0, is_stdout = 0;
 
425
        while ((c = getopt(argc, argv, "nom:")) >= 0) {
 
426
                switch (c) {
 
427
                case 'o': is_stdout = 1; break;
 
428
                case 'n': is_by_qname = 1; break;
 
429
                case 'm': max_mem = atol(optarg); break;
 
430
                }
 
431
        }
 
432
        if (optind + 2 > argc) {
 
433
                fprintf(stderr, "Usage: samtools sort [-on] [-m <maxMem>] <in.bam> <out.prefix>\n");
 
434
                return 1;
 
435
        }
 
436
        bam_sort_core_ext(is_by_qname, argv[optind], argv[optind+1], max_mem, is_stdout);
 
437
        return 0;
 
438
}