From 9710d1947970b4214e2189bd4db906dd73107512 Mon Sep 17 00:00:00 2001
From: Heng Li <lh3@live.co.uk>
Date: Sat, 13 Jun 2009 12:58:48 +0000
Subject: [PATCH] added an example of calculating read depth

---
 examples/Makefile   | 10 +++++--
 examples/calDepth.c | 65 +++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 73 insertions(+), 2 deletions(-)
 create mode 100644 examples/calDepth.c

diff --git a/examples/Makefile b/examples/Makefile
index b6908fa..3fe3e5a 100644
--- 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
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;
+}
-- 
2.39.2