added an example of calculating read depth

diff --git a/examples/Makefile b/examples/Makefile
index b6908faff0312ddf51e6a413ca03ad74508bcead..3fe3e5a63a8643c49e8bf75447d4a55fac847add 100644 (file)
--- a/examples/Makefile
+++ b/examples/Makefile
@@ -1,4 +1,4 @@
-all:../samtools ex1.glf ex1.pileup.gz ex1.bam.bai ex1.glfview.gz
+all:../libbam.a ../samtools ex1.glf ex1.pileup.gz ex1.bam.bai ex1.glfview.gz calDepth
                @echo; echo \# You can now launch the viewer with: \'samtools tview ex1.bam ex1.fa\'; echo;
 
 ex1.fa.fai:ex1.fa
@@ -17,5 +17,11 @@ ex1.glfview.gz:ex1.glf
 ../samtools:
                (cd ..; make samtools)
 
+../libbam.a:
+               (cd ..; make libbam.a)
+
+calDepth:../libbam.a calDepth.c
+               gcc -g -Wall -O2 -I.. calDepth.c -o $@ -lm -lz -L.. -lbam
+
 clean:
-               rm -f *.bam *.bai *.glf* *.fai *.pileup* *~
\ No newline at end of file
+               rm -fr *.bam *.bai *.glf* *.fai *.pileup* *~ calDepth *.dSYM
\ No newline at end of file

diff --git a/examples/calDepth.c b/examples/calDepth.c
new file mode 100644 (file)
index 0000000..11e596b
--- /dev/null
+++ b/examples/calDepth.c
@@ -0,0 +1,65 @@
+#include <stdlib.h>
+#include <string.h>
+#include <stdio.h>
+#include <unistd.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) {
+               sampileup(tmp.in, -1, pileup_func, &tmp);
+       } else {
+               int ref;
+               bam_index_t *idx;
+               bam_plbuf_t *buf;
+               idx = bam_index_load(argv[1]);
+               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);
+               if (ref < 0) {
+                       fprintf(stderr, "Invalid region %s\n", argv[2]);
+                       return 1;
+               }
+               buf = bam_plbuf_init(pileup_func, &tmp);
+               bam_fetch(tmp.in->x.bam, idx, ref, tmp.beg, tmp.end, buf, fetch_func);
+               bam_plbuf_push(0, buf);
+               bam_index_destroy(idx);
+               bam_plbuf_destroy(buf);
+       }
+       samclose(tmp.in);
+       return 0;
+}