Logo Search packages:      
Sourcecode: samtools version File versions  Download package

calDepth.c

#include <stdio.h>
#include "sam.h"

typedef struct {
      int beg, end;
      samfile_t *in;
} tmpstruct_t;

// callback for bam_fetch()
static int fetch_func(const bam1_t *b, void *data)
{
      bam_plbuf_t *buf = (bam_plbuf_t*)data;
      bam_plbuf_push(b, buf);
      return 0;
}
// callback for bam_plbuf_init()
static int pileup_func(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *pl, void *data)
{
      tmpstruct_t *tmp = (tmpstruct_t*)data;
      if ((int)pos >= tmp->beg && (int)pos < tmp->end)
            printf("%s\t%d\t%d\n", tmp->in->header->target_name[tid], pos + 1, n);
      return 0;
}

int main(int argc, char *argv[])
{
      tmpstruct_t tmp;
      if (argc == 1) {
            fprintf(stderr, "Usage: calDepth <in.bam> [region]\n");
            return 1;
      }
      tmp.beg = 0; tmp.end = 0x7fffffff;
      tmp.in = samopen(argv[1], "rb", 0);
      if (tmp.in == 0) {
            fprintf(stderr, "Fail to open BAM file %s\n", argv[1]);
            return 1;
      }
      if (argc == 2) { // if a region is not specified
            sampileup(tmp.in, -1, pileup_func, &tmp);
      } else {
            int ref;
            bam_index_t *idx;
            bam_plbuf_t *buf;
            idx = bam_index_load(argv[1]); // load BAM index
            if (idx == 0) {
                  fprintf(stderr, "BAM indexing file is not available.\n");
                  return 1;
            }
            bam_parse_region(tmp.in->header, argv[2], &ref, &tmp.beg, &tmp.end); // parse the region
            if (ref < 0) {
                  fprintf(stderr, "Invalid region %s\n", argv[2]);
                  return 1;
            }
            buf = bam_plbuf_init(pileup_func, &tmp); // initialize pileup
            bam_fetch(tmp.in->x.bam, idx, ref, tmp.beg, tmp.end, buf, fetch_func);
            bam_plbuf_push(0, buf); // finalize pileup
            bam_index_destroy(idx);
            bam_plbuf_destroy(buf);
      }
      samclose(tmp.in);
      return 0;
}

Generated by  Doxygen 1.6.0   Back to index