12 KSTREAM_INIT(gzFile, gzread, 16384)
20 static int read_bam(void *data, bam1_t *b)
22 aux_t *aux = (aux_t*)data;
23 int ret = bam_iter_read(aux->fp, aux->iter, b);
24 if ((int)b->core.qual < aux->min_mapQ) b->core.flag |= BAM_FUNMAP;
28 int main_bedcov(int argc, char *argv[])
30 extern void bam_init_header_hash(bam_header_t*);
37 int *n_plp, dret, i, n, c, min_mapQ = 0;
39 const bam_pileup1_t **plp;
41 while ((c = getopt(argc, argv, "Q:")) >= 0) {
43 case 'Q': min_mapQ = atoi(optarg); break;
46 if (optind + 2 > argc) {
47 fprintf(stderr, "Usage: samtools bedcov <in.bed> <in1.bam> [...]\n");
50 memset(&str, 0, sizeof(kstring_t));
51 n = argc - optind - 1;
52 aux = calloc(n, sizeof(void*));
53 idx = calloc(n, sizeof(void*));
54 for (i = 0; i < n; ++i) {
55 aux[i] = calloc(1, sizeof(aux_t));
56 aux[i]->min_mapQ = min_mapQ;
57 aux[i]->fp = bam_open(argv[i+optind+1], "r");
58 idx[i] = bam_index_load(argv[i+optind+1]);
59 if (aux[i]->fp == 0 || idx[i] == 0) {
60 fprintf(stderr, "ERROR: fail to open index BAM file '%s'\n", argv[i+optind+1]);
63 bgzf_set_cache_size(aux[i]->fp, 20);
64 if (i == 0) h = bam_header_read(aux[0]->fp);
66 bam_init_header_hash(h);
69 fp = gzopen(argv[optind], "rb");
71 n_plp = calloc(n, sizeof(int));
72 plp = calloc(n, sizeof(void*));
73 while (ks_getuntil(ks, KS_SEP_LINE, &str, &dret) >= 0) {
75 int tid, beg, end, pos;
78 for (p = q = str.s; *p && *p != '\t'; ++p);
79 if (*p != '\t') goto bed_error;
80 *p = 0; tid = bam_get_tid(h, q); *p = '\t';
81 if (tid < 0) goto bed_error;
82 for (q = p = p + 1; isdigit(*p); ++p);
83 if (*p != '\t') goto bed_error;
84 *p = 0; beg = atoi(q); *p = '\t';
85 for (q = p = p + 1; isdigit(*p); ++p);
86 if (*p == '\t' || *p == 0) {
88 *p = 0; end = atoi(q); *p = c;
89 } else goto bed_error;
91 for (i = 0; i < n; ++i) {
92 if (aux[i]->iter) bam_iter_destroy(aux[i]->iter);
93 aux[i]->iter = bam_iter_query(idx[i], tid, beg, end);
95 mplp = bam_mplp_init(n, read_bam, (void**)aux);
96 bam_mplp_set_maxcnt(mplp, 64000);
97 memset(cnt, 0, 8 * n);
98 while (bam_mplp_auto(mplp, &tid, &pos, n_plp, plp) > 0)
99 if (pos >= beg && pos < end)
100 for (i = 0; i < n; ++i) cnt[i] += n_plp[i];
101 for (i = 0; i < n; ++i) {
106 bam_mplp_destroy(mplp);
110 fprintf(stderr, "Errors in BED line '%s'\n", str.s);
112 free(n_plp); free(plp);
117 for (i = 0; i < n; ++i) {
118 if (aux[i]->iter) bam_iter_destroy(aux[i]->iter);
119 bam_index_destroy(idx[i]);
120 bam_close(aux[i]->fp);
123 bam_header_destroy(h);
124 free(aux); free(idx);