8 typedef bam1_t *bam1_p;
10 KHASH_SET_INIT_STR(name)
11 KHASH_MAP_INIT_INT64(pos, bam1_p)
13 #define BUFFER_SIZE 0x40000
20 static inline void stack_insert(tmp_stack_t *stack, bam1_t *b)
22 if (stack->n == stack->max) {
23 stack->max = stack->max? stack->max<<1 : 0x10000;
24 stack->a = (bam1_t**)realloc(stack->a, sizeof(bam1_t*) * stack->max);
26 stack->a[stack->n++] = b;
29 static inline void dump_best(tmp_stack_t *stack, khash_t(pos) *best_hash, bamFile out)
32 for (i = 0; i != stack->n; ++i) {
33 bam_write1(out, stack->a[i]);
34 bam_destroy1(stack->a[i]);
37 if (kh_size(best_hash) > BUFFER_SIZE) kh_clear(pos, best_hash);
40 static void clear_del_set(khash_t(name) *del_set)
43 for (k = kh_begin(del_set); k < kh_end(del_set); ++k)
44 if (kh_exist(del_set, k))
45 free((char*)kh_key(del_set, k));
46 kh_clear(name, del_set);
49 void bam_rmdup_core(bamFile in, bamFile out)
53 int last_tid = -1, last_pos = -1;
54 uint64_t n_checked = 0, n_removed = 0;
57 khash_t(pos) *best_hash;
58 khash_t(name) *del_set;
60 best_hash = kh_init(pos);
61 del_set = kh_init(name);
63 memset(&stack, 0, sizeof(tmp_stack_t));
64 header = bam_header_read(in);
65 bam_header_write(out, header);
67 kh_resize(name, del_set, 4 * BUFFER_SIZE);
68 kh_resize(pos, best_hash, 3 * BUFFER_SIZE);
69 while (bam_read1(in, b) >= 0) {
70 bam1_core_t *c = &b->core;
71 if (c->tid != last_tid || last_pos != c->pos) {
72 dump_best(&stack, best_hash, out); // write the result
73 if (c->tid != last_tid) {
74 kh_clear(pos, best_hash);
75 if (kh_size(del_set)) { // check
76 fprintf(stderr, "[bam_rmdup_core] %llu unmatched pairs\n", (long long)kh_size(del_set));
77 clear_del_set(del_set);
80 fprintf(stderr, "[bam_rmdup_core] processing reference %s...\n", header->target_name[c->tid]);
83 if (!(c->flag&BAM_FPAIRED) || (c->flag&(BAM_FUNMAP|BAM_FMUNMAP)) || (c->mtid >= 0 && c->tid != c->mtid)) {
85 } else if (c->isize > 0) { // paired, head
86 uint64_t key = (uint64_t)c->pos<<32 | c->isize;
89 k = kh_put(pos, best_hash, key, &ret);
90 if (ret == 0) { // found in best_hash
91 bam1_t *p = kh_val(best_hash, k);
93 if (p->core.qual < c->qual) { // the current alignment is better
94 kh_put(name, del_set, strdup(bam1_qname(p)), &ret); // p will be removed
95 bam_copy1(p, b); // replaced as b
96 } else kh_put(name, del_set, strdup(bam1_qname(b)), &ret); // b will be removed
98 fprintf(stderr, "[bam_rmdup_core] inconsistent BAM file for pair '%s'. Continue anyway.\n", bam1_qname(b));
99 } else { // not found in best_hash
100 kh_val(best_hash, k) = bam_dup1(b);
101 stack_insert(&stack, kh_val(best_hash, k));
103 } else { // paired, tail
104 k = kh_get(name, del_set, bam1_qname(b));
105 if (k != kh_end(del_set)) {
106 free((char*)kh_key(del_set, k));
107 kh_del(name, del_set, k);
108 } else bam_write1(out, b);
112 dump_best(&stack, best_hash, out);
114 bam_header_destroy(header);
115 clear_del_set(del_set);
116 kh_destroy(name, del_set);
117 kh_destroy(pos, best_hash);
120 fprintf(stderr, "[bam_rmdup_core] %lld / %lld = %.4lf\n", (long long)n_removed, (long long)n_checked,
121 (double)n_removed/n_checked);
123 int bam_rmdup(int argc, char *argv[])
127 fprintf(stderr, "Usage: samtools rmdup <input.srt.bam> <output.bam>\n");
130 in = (strcmp(argv[1], "-") == 0)? bam_dopen(fileno(stdin), "r") : bam_open(argv[1], "r");
131 out = (strcmp(argv[2], "-") == 0)? bam_dopen(fileno(stdout), "w") : bam_open(argv[2], "w");
132 if (in == 0 || out == 0) {
133 fprintf(stderr, "[bam_rmdup] fail to read/write input files\n");
136 bam_rmdup_core(in, out);