bam_lpileup.c

   1 #include <stdlib.h>
   2 #include <stdio.h>
   3 #include "bam.h"
   4 #include "ksort.h"
   5
   6 #define TV_GAP 2
   7
   8 typedef struct __freenode_t {
   9         uint32_t level:28, cnt:4;
  10         struct __freenode_t *next;
  11 } freenode_t, *freenode_p;
  12
  13 #define freenode_lt(a,b) ((a)->cnt < (b)->cnt || ((a)->cnt == (b)->cnt && (a)->level < (b)->level))
  14 KSORT_INIT(node, freenode_p, freenode_lt)
  15
  16 /* Memory pool, similar to the one in bam_pileup.c */
  17 typedef struct {
  18         int cnt, n, max;
  19         freenode_t **buf;
  20 } mempool_t;
  21
  22 static mempool_t *mp_init()
  23 {
  24         return (mempool_t*)calloc(1, sizeof(mempool_t));
  25 }
  26 static void mp_destroy(mempool_t *mp)
  27 {
  28         int k;
  29         for (k = 0; k < mp->n; ++k) free(mp->buf[k]);
  30         free(mp->buf); free(mp);
  31 }
  32 static inline freenode_t *mp_alloc(mempool_t *mp)
  33 {
  34         ++mp->cnt;
  35         if (mp->n == 0) return (freenode_t*)calloc(1, sizeof(freenode_t));
  36         else return mp->buf[--mp->n];
  37 }
  38 static inline void mp_free(mempool_t *mp, freenode_t *p)
  39 {
  40         --mp->cnt; p->next = 0; p->cnt = TV_GAP;
  41         if (mp->n == mp->max) {
  42                 mp->max = mp->max? mp->max<<1 : 256;
  43                 mp->buf = (freenode_t**)realloc(mp->buf, sizeof(freenode_t*) * mp->max);
  44         }
  45         mp->buf[mp->n++] = p;
  46 }
  47
  48 /* core part */
  49 struct __bam_lplbuf_t {
  50         int max, n_cur, n_pre;
  51         int max_level, *cur_level, *pre_level;
  52         mempool_t *mp;
  53         freenode_t **aux, *head, *tail;
  54         int n_nodes, m_aux;
  55         bam_pileup_f func;
  56         void *user_data;
  57         bam_plbuf_t *plbuf;
  58 };
  59
  60 void bam_lplbuf_reset(bam_lplbuf_t *buf)
  61 {
  62         freenode_t *p, *q;
  63         bam_plbuf_reset(buf->plbuf);
  64         for (p = buf->head; p->next;) {
  65                 q = p->next;
  66                 mp_free(buf->mp, p);
  67                 p = q;
  68         }
  69         buf->head = buf->tail;
  70         buf->max_level = 0;
  71         buf->n_cur = buf->n_pre = 0;
  72         buf->n_nodes = 0;
  73 }
  74
  75 static int tview_func(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *pl, void *data)
  76 {
  77         bam_lplbuf_t *tv = (bam_lplbuf_t*)data;
  78         freenode_t *p;
  79         int i, l, max_level;
  80         // allocate memory if necessary
  81         if (tv->max < n) { // enlarge
  82                 tv->max = n;
  83                 kroundup32(tv->max);
  84                 tv->cur_level = (int*)realloc(tv->cur_level, sizeof(int) * tv->max);
  85                 tv->pre_level = (int*)realloc(tv->pre_level, sizeof(int) * tv->max);
  86         }
  87         tv->n_cur = n;
  88         // update cnt
  89         for (p = tv->head; p->next; p = p->next)
  90                 if (p->cnt > 0) --p->cnt;
  91         // calculate cur_level[]
  92         max_level = 0;
  93         for (i = l = 0; i < n; ++i) {
  94                 const bam_pileup1_t *p = pl + i;
  95                 if (p->qpos == 0) {
  96                         if (tv->head->next && tv->head->cnt == 0) { // then take a free slot
  97                                 freenode_t *p = tv->head->next;
  98                                 tv->cur_level[i] = tv->head->level;
  99                                 mp_free(tv->mp, tv->head);
 100                                 tv->head = p;
 101                                 --tv->n_nodes;
 102                         } else tv->cur_level[i] = ++tv->max_level;
 103                 } else {
 104                         tv->cur_level[i] = tv->pre_level[l++];
 105                         if (p->qpos == p->b->core.l_qseq - 1) { // then return a free slot
 106                                 tv->tail->level = tv->cur_level[i];
 107                                 tv->tail->next = mp_alloc(tv->mp);
 108                                 tv->tail = tv->tail->next;
 109                                 ++tv->n_nodes;
 110                         }
 111                 }
 112                 if (tv->cur_level[i] > max_level) max_level = tv->cur_level[i];
 113                 ((bam_pileup1_t*)p)->level = tv->cur_level[i];
 114         }
 115         assert(l == tv->n_pre);
 116         tv->func(tid, pos, n, pl, tv->user_data);
 117         // sort the linked list
 118         if (tv->n_nodes) {
 119                 freenode_t *q;
 120                 if (tv->n_nodes + 1 > tv->m_aux) { // enlarge
 121                         tv->m_aux = tv->n_nodes + 1;
 122                         kroundup32(tv->m_aux);
 123                         tv->aux = (freenode_t**)realloc(tv->aux, sizeof(void*) * tv->m_aux);
 124                 }
 125                 for (p = tv->head, i = l = 0; p->next;) {
 126                         if (p->level > max_level) { // then discard this entry
 127                                 q = p->next;
 128                                 mp_free(tv->mp, p);
 129                                 p = q;
 130                         } else {
 131                                 tv->aux[i++] = p;
 132                                 p = p->next;
 133                         }
 134                 }
 135                 tv->aux[i] = tv->tail; // add a proper tail for the loop below
 136                 tv->n_nodes = i;
 137                 if (tv->n_nodes) {
 138                         ks_introsort(node, tv->n_nodes, tv->aux);
 139                         for (i = 0; i < tv->n_nodes; ++i) tv->aux[i]->next = tv->aux[i+1];
 140                         tv->head = tv->aux[0];
 141                 } else tv->head = tv->tail;
 142         }
 143         // clean up
 144         tv->max_level = max_level;
 145         memcpy(tv->pre_level, tv->cur_level, tv->n_cur * 4);
 146         // squeeze out terminated levels
 147         for (i = l = 0; i < n; ++i) {
 148                 const bam_pileup1_t *p = pl + i;
 149                 if (p->qpos != p->b->core.l_qseq - 1)
 150                         tv->pre_level[l++] = tv->pre_level[i];
 151         }
 152         tv->n_pre = l;
 153         return 0;
 154 }
 155
 156 bam_lplbuf_t *bam_lplbuf_init(bam_pileup_f func, void *data)
 157 {
 158         bam_lplbuf_t *tv;
 159         tv = (bam_lplbuf_t*)calloc(1, sizeof(bam_lplbuf_t));
 160         tv->mp = mp_init();
 161         tv->head = tv->tail = mp_alloc(tv->mp);
 162         tv->func = func;
 163         tv->user_data = data;
 164         tv->plbuf = bam_plbuf_init(tview_func, tv);
 165         return (bam_lplbuf_t*)tv;
 166 }
 167
 168 void bam_lplbuf_destroy(bam_lplbuf_t *tv)
 169 {
 170         mp_free(tv->mp, tv->head);
 171         mp_destroy(tv->mp);
 172         free(tv->cur_level); free(tv->pre_level);
 173         bam_plbuf_destroy(tv->plbuf);
 174         free(tv->aux);
 175         free(tv);
 176 }
 177
 178 int bam_lplbuf_push(const bam1_t *b, bam_lplbuf_t *tv)
 179 {
 180         return bam_plbuf_push(b, tv->plbuf);
 181 }
 182
 183 int bam_lpileup_file(bamFile fp, bam_pileup_f func, void *func_data)
 184 {
 185         bam_lplbuf_t *buf;
 186         int ret;
 187         bam1_t *b;
 188         b = (bam1_t*)calloc(1, sizeof(bam1_t));
 189         buf = bam_lplbuf_init(func, func_data);
 190         while ((ret = bam_read1(fp, b)) >= 0)
 191                 bam_lplbuf_push(b, buf);
 192         bam_lplbuf_push(0, buf);
 193         bam_lplbuf_destroy(buf);
 194         free(b->data); free(b);
 195         return 0;
 196 }