Mercurial > repos > siyuan > prada
view pyPRADA_1.2/tools/samtools-0.1.16/examples/calDepth.c @ 3:f17965495ec9 draft default tip
Uploaded
author | siyuan |
---|---|
date | Tue, 11 Mar 2014 12:14:01 -0400 |
parents | acc2ca1a3ba4 |
children |
line wrap: on
line source
#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; }