1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
|
#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;
}
|
