- int i, ret, last_tid, min_rpos = 0x7fffffff, capacity;
- bam1_t *b = bam_init1();
- bam1_core_t *c = &b->core;
- // squeeze out the empty cells at the beginning
- for (i = 0; i < buf->n; ++i)
- if (buf->buf[i].b) break;
- if (i < buf->n) { // squeeze
- if (i > 0) {
- memmove(buf->buf, buf->buf + i, sizeof(elem_t) * (buf->n - i));
- buf->n = buf->n - i;
- }
- } else buf->n = 0;
- // calculate min_rpos
- for (i = 0; i < buf->n; ++i) {
- elem_t *e = buf->buf + i;
- if (e->b && e->rpos >= 0 && e->rpos < min_rpos)
- min_rpos = buf->buf[i].rpos;
- }
- // fill the buffer
- buf->x = -1;
- last_tid = buf->n? buf->buf[0].b->core.tid : -1;
- capacity = buf->n + BLOCK_SIZE;
- while ((ret = samread(in, b)) >= 0) {
- elem_t *e;
- uint8_t *qual = bam1_qual(b);
- int is_mapped;
- if (last_tid < 0) last_tid = c->tid;
- if (c->tid != last_tid) {
- if (buf->x < 0) buf->x = buf->n;
- }
- if (buf->n >= buf->max) { // enlarge
- buf->max = buf->max? buf->max<<1 : 8;
- buf->buf = (elem_t*)realloc(buf->buf, sizeof(elem_t) * buf->max);
- }
- e = &buf->buf[buf->n++];
- e->b = bam_dup1(b);
- e->rpos = -1; e->score = 0;
- for (i = 0; i < c->l_qseq; ++i) e->score += qual[i] + 1;
- e->score = (double)e->score / sqrt(c->l_qseq + 1);
- is_mapped = (c->tid < 0 || c->tid >= in->header->n_targets || (c->flag&BAM_FUNMAP))? 0 : 1;
- if (!is_mapped) e->score = -1;
- if (is_mapped && (c->flag & BAM_FREVERSE)) {
- e->rpos = b->core.pos + bam_calend(&b->core, bam1_cigar(b));
- if (min_rpos > e->rpos) min_rpos = e->rpos;
- }
- if (buf->n >= capacity) {
- if (is_mapped && c->pos <= min_rpos) capacity += BLOCK_SIZE;
- else break;
- }
- }
- if (ret >= 0 && buf->x < 0) buf->x = buf->n;
- bam_destroy1(b);
- return buf->n;
+ khint_t k = kh_get(lib, aux, lib);
+ if (k == kh_end(aux)) {
+ int ret;
+ char *p = strdup(lib);
+ lib_aux_t *q;
+ k = kh_put(lib, aux, p, &ret);
+ q = &kh_val(aux, k);
+ q->left = kh_init(best);
+ q->rght = kh_init(best);
+ q->n_checked = q->n_removed = 0;
+ return q;
+ } else return &kh_val(aux, k);
+}
+
+static inline int sum_qual(const bam1_t *b)
+{
+ int i, q;
+ uint8_t *qual = bam1_qual(b);
+ for (i = q = 0; i < b->core.l_qseq; ++i) q += qual[i];
+ return q;
+}
+
+static inline elem_t *push_queue(queue_t *queue, const bam1_t *b, int endpos, int score)
+{
+ elem_t *p = kl_pushp(q, queue);
+ p->discarded = 0;
+ p->endpos = endpos; p->score = score;
+ if (p->b == 0) p->b = bam_init1();
+ bam_copy1(p->b, b);
+ return p;