See http://htslib.org/ for the new 1.x releases of SAMtools, BCFtools, and HTSlib. This website contains information pertaining to the old 0.1.19 samtools release, and so is useful but somewhat out of date. As time permits, this information will be updated for the new samtools/bcftools versions and moved to the new website.


view plainprint?
  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,  
  50.                          &tmp.beg, &tmp.end); // parse the region  
  51.         if (ref < 0) {  
  52.             fprintf(stderr, "Invalid region %s\n", argv[2]);  
  53.             return 1;  
  54.         }  
  55.         buf = bam_plbuf_init(pileup_func, &tmp); // initialize pileup  
  56.         bam_fetch(tmp.in->x.bam, idx, ref, tmp.beg, tmp.end, buf, fetch_func);  
  57.         bam_plbuf_push(0, buf); // finalize pileup  
  58.         bam_index_destroy(idx);  
  59.         bam_plbuf_destroy(buf);  
  60.     }  
  61.     samclose(tmp.in);  
  62.     return 0;  
  63. }