|
0
|
1 #include <zlib.h>
|
|
|
2 #include <stdio.h>
|
|
|
3 #include <ctype.h>
|
|
|
4 #include <stdlib.h>
|
|
|
5 #include <string.h>
|
|
|
6 #include <unistd.h>
|
|
|
7 #include "kstring.h"
|
|
|
8 #include "bgzf.h"
|
|
|
9 #include "bam.h"
|
|
|
10
|
|
|
11 #include "kseq.h"
|
|
|
12 KSTREAM_INIT(gzFile, gzread, 16384)
|
|
|
13
|
|
|
14 typedef struct {
|
|
|
15 bamFile fp;
|
|
|
16 bam_iter_t iter;
|
|
|
17 int min_mapQ;
|
|
|
18 } aux_t;
|
|
|
19
|
|
|
20 static int read_bam(void *data, bam1_t *b)
|
|
|
21 {
|
|
|
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;
|
|
|
25 return ret;
|
|
|
26 }
|
|
|
27
|
|
|
28 int main_bedcov(int argc, char *argv[])
|
|
|
29 {
|
|
|
30 extern void bam_init_header_hash(bam_header_t*);
|
|
|
31 gzFile fp;
|
|
|
32 kstring_t str;
|
|
|
33 kstream_t *ks;
|
|
|
34 bam_index_t **idx;
|
|
|
35 bam_header_t *h = 0;
|
|
|
36 aux_t **aux;
|
|
|
37 int *n_plp, dret, i, n, c, min_mapQ = 0;
|
|
|
38 int64_t *cnt;
|
|
|
39 const bam_pileup1_t **plp;
|
|
|
40
|
|
|
41 while ((c = getopt(argc, argv, "Q:")) >= 0) {
|
|
|
42 switch (c) {
|
|
|
43 case 'Q': min_mapQ = atoi(optarg); break;
|
|
|
44 }
|
|
|
45 }
|
|
|
46 if (optind + 2 > argc) {
|
|
|
47 fprintf(stderr, "Usage: samtools bedcov <in.bed> <in1.bam> [...]\n");
|
|
|
48 return 1;
|
|
|
49 }
|
|
|
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]);
|
|
|
61 return 2;
|
|
|
62 }
|
|
|
63 bgzf_set_cache_size(aux[i]->fp, 20);
|
|
|
64 if (i == 0) h = bam_header_read(aux[0]->fp);
|
|
|
65 }
|
|
|
66 bam_init_header_hash(h);
|
|
|
67 cnt = calloc(n, 8);
|
|
|
68
|
|
|
69 fp = gzopen(argv[optind], "rb");
|
|
|
70 ks = ks_init(fp);
|
|
|
71 n_plp = calloc(n, sizeof(int));
|
|
|
72 plp = calloc(n, sizeof(void*));
|
|
|
73 while (ks_getuntil(ks, KS_SEP_LINE, &str, &dret) >= 0) {
|
|
|
74 char *p, *q;
|
|
|
75 int tid, beg, end, pos;
|
|
|
76 bam_mplp_t mplp;
|
|
|
77
|
|
|
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) {
|
|
|
87 int c = *p;
|
|
|
88 *p = 0; end = atoi(q); *p = c;
|
|
|
89 } else goto bed_error;
|
|
|
90
|
|
|
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);
|
|
|
94 }
|
|
|
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) {
|
|
|
102 kputc('\t', &str);
|
|
|
103 kputl(cnt[i], &str);
|
|
|
104 }
|
|
|
105 puts(str.s);
|
|
|
106 bam_mplp_destroy(mplp);
|
|
|
107 continue;
|
|
|
108
|
|
|
109 bed_error:
|
|
|
110 fprintf(stderr, "Errors in BED line '%s'\n", str.s);
|
|
|
111 }
|
|
|
112 free(n_plp); free(plp);
|
|
|
113 ks_destroy(ks);
|
|
|
114 gzclose(fp);
|
|
|
115
|
|
|
116 free(cnt);
|
|
|
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);
|
|
|
121 free(aux[i]);
|
|
|
122 }
|
|
|
123 bam_header_destroy(h);
|
|
|
124 free(aux); free(idx);
|
|
|
125 free(str.s);
|
|
|
126 return 0;
|
|
|
127 }
|
