examples/calDepth.c

   1 #include <stdio.h>
   2 #include "sam.h"
   3
   4 typedef struct {
   5         int beg, end;
   6         samfile_t *in;
   7 } tmpstruct_t;
   8
   9 // callback for bam_fetch()
  10 static int fetch_func(const bam1_t *b, void *data)
  11 {
  12         bam_plbuf_t *buf = (bam_plbuf_t*)data;
  13         bam_plbuf_push(b, buf);
  14         return 0;
  15 }
  16 // callback for bam_plbuf_init()
  17 static int pileup_func(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *pl, void *data)
  18 {
  19         tmpstruct_t *tmp = (tmpstruct_t*)data;
  20         if ((int)pos >= tmp->beg && (int)pos < tmp->end)
  21                 printf("%s\t%d\t%d\n", tmp->in->header->target_name[tid], pos + 1, n);
  22         return 0;
  23 }
  24
  25 int main(int argc, char *argv[])
  26 {
  27         tmpstruct_t tmp;
  28         if (argc == 1) {
  29                 fprintf(stderr, "Usage: calDepth <in.bam> [region]\n");
  30                 return 1;
  31         }
  32         tmp.beg = 0; tmp.end = 0x7fffffff;
  33         tmp.in = samopen(argv[1], "rb", 0);
  34         if (tmp.in == 0) {
  35                 fprintf(stderr, "Fail to open BAM file %s\n", argv[1]);
  36                 return 1;
  37         }
  38         if (argc == 2) { // if a region is not specified
  39                 sampileup(tmp.in, -1, pileup_func, &tmp);
  40         } else {
  41                 int ref;
  42                 bam_index_t *idx;
  43                 bam_plbuf_t *buf;
  44                 idx = bam_index_load(argv[1]); // load BAM index
  45                 if (idx == 0) {
  46                         fprintf(stderr, "BAM indexing file is not available.\n");
  47                         return 1;
  48                 }
  49                 bam_parse_region(tmp.in->header, argv[2], &ref, &tmp.beg, &tmp.end); // parse the region
  50                 if (ref < 0) {
  51                         fprintf(stderr, "Invalid region %s\n", argv[2]);
  52                         return 1;
  53                 }
  54                 buf = bam_plbuf_init(pileup_func, &tmp); // initialize pileup
  55                 bam_fetch(tmp.in->x.bam, idx, ref, tmp.beg, tmp.end, buf, fetch_func);
  56                 bam_plbuf_push(0, buf); // finalize pileup
  57                 bam_index_destroy(idx);
  58                 bam_plbuf_destroy(buf);
  59         }
  60         samclose(tmp.in);
  61         return 0;
  62 }