|
0
|
1 /* This program demonstrates how to generate pileup from multiple BAMs
|
|
|
2 * simutaneously, to achieve random access and to use the BED interface.
|
|
|
3 * To compile this program separately, you may:
|
|
|
4 *
|
|
|
5 * gcc -g -O2 -Wall -o bam2depth -D_MAIN_BAM2DEPTH bam2depth.c -L. -lbam -lz
|
|
|
6 */
|
|
|
7 #include <stdlib.h>
|
|
|
8 #include <string.h>
|
|
|
9 #include <stdio.h>
|
|
|
10 #include <unistd.h>
|
|
|
11 #include "bam.h"
|
|
|
12
|
|
|
13 typedef struct { // auxiliary data structure
|
|
|
14 bamFile fp; // the file handler
|
|
|
15 bam_iter_t iter; // NULL if a region not specified
|
|
|
16 int min_mapQ, min_len; // mapQ filter; length filter
|
|
|
17 } aux_t;
|
|
|
18
|
|
|
19 void *bed_read(const char *fn); // read a BED or position list file
|
|
|
20 void bed_destroy(void *_h); // destroy the BED data structure
|
|
|
21 int bed_overlap(const void *_h, const char *chr, int beg, int end); // test if chr:beg-end overlaps
|
|
|
22
|
|
|
23 // This function reads a BAM alignment from one BAM file.
|
|
|
24 static int read_bam(void *data, bam1_t *b) // read level filters better go here to avoid pileup
|
|
|
25 {
|
|
|
26 aux_t *aux = (aux_t*)data; // data in fact is a pointer to an auxiliary structure
|
|
|
27 int ret = aux->iter? bam_iter_read(aux->fp, aux->iter, b) : bam_read1(aux->fp, b);
|
|
|
28 if (!(b->core.flag&BAM_FUNMAP)) {
|
|
|
29 if ((int)b->core.qual < aux->min_mapQ) b->core.flag |= BAM_FUNMAP;
|
|
|
30 else if (aux->min_len && bam_cigar2qlen(&b->core, bam1_cigar(b)) < aux->min_len) b->core.flag |= BAM_FUNMAP;
|
|
|
31 }
|
|
|
32 return ret;
|
|
|
33 }
|
|
|
34
|
|
|
35 int read_file_list(const char *file_list,int *n,char **argv[]);
|
|
|
36
|
|
|
37 #ifdef _MAIN_BAM2DEPTH
|
|
|
38 int main(int argc, char *argv[])
|
|
|
39 #else
|
|
|
40 int main_depth(int argc, char *argv[])
|
|
|
41 #endif
|
|
|
42 {
|
|
|
43 int i, n, tid, beg, end, pos, *n_plp, baseQ = 0, mapQ = 0, min_len = 0, nfiles;
|
|
|
44 const bam_pileup1_t **plp;
|
|
|
45 char *reg = 0; // specified region
|
|
|
46 void *bed = 0; // BED data structure
|
|
|
47 char *file_list = NULL, **fn = NULL;
|
|
|
48 bam_header_t *h = 0; // BAM header of the 1st input
|
|
|
49 aux_t **data;
|
|
|
50 bam_mplp_t mplp;
|
|
|
51
|
|
|
52 // parse the command line
|
|
|
53 while ((n = getopt(argc, argv, "r:b:q:Q:l:f:")) >= 0) {
|
|
|
54 switch (n) {
|
|
|
55 case 'l': min_len = atoi(optarg); break; // minimum query length
|
|
|
56 case 'r': reg = strdup(optarg); break; // parsing a region requires a BAM header
|
|
|
57 case 'b': bed = bed_read(optarg); break; // BED or position list file can be parsed now
|
|
|
58 case 'q': baseQ = atoi(optarg); break; // base quality threshold
|
|
|
59 case 'Q': mapQ = atoi(optarg); break; // mapping quality threshold
|
|
|
60 case 'f': file_list = optarg; break;
|
|
|
61 }
|
|
|
62 }
|
|
|
63 if (optind == argc && !file_list) {
|
|
|
64 fprintf(stderr, "\n");
|
|
|
65 fprintf(stderr, "Usage: samtools depth [options] in1.bam [in2.bam [...]]\n");
|
|
|
66 fprintf(stderr, "Options:\n");
|
|
|
67 fprintf(stderr, " -b <bed> list of positions or regions\n");
|
|
|
68 fprintf(stderr, " -f <list> list of input BAM filenames, one per line [null]\n");
|
|
|
69 fprintf(stderr, " -l <int> minQLen\n");
|
|
|
70 fprintf(stderr, " -q <int> base quality threshold\n");
|
|
|
71 fprintf(stderr, " -Q <int> mapping quality threshold\n");
|
|
|
72 fprintf(stderr, " -r <chr:from-to> region\n");
|
|
|
73 fprintf(stderr, "\n");
|
|
|
74 return 1;
|
|
|
75 }
|
|
|
76
|
|
|
77 // initialize the auxiliary data structures
|
|
|
78 if (file_list)
|
|
|
79 {
|
|
|
80 if ( read_file_list(file_list,&nfiles,&fn) ) return 1;
|
|
|
81 n = nfiles;
|
|
|
82 argv = fn;
|
|
|
83 optind = 0;
|
|
|
84 }
|
|
|
85 else
|
|
|
86 n = argc - optind; // the number of BAMs on the command line
|
|
|
87 data = calloc(n, sizeof(void*)); // data[i] for the i-th input
|
|
|
88 beg = 0; end = 1<<30; tid = -1; // set the default region
|
|
|
89 for (i = 0; i < n; ++i) {
|
|
|
90 bam_header_t *htmp;
|
|
|
91 data[i] = calloc(1, sizeof(aux_t));
|
|
|
92 data[i]->fp = bam_open(argv[optind+i], "r"); // open BAM
|
|
|
93 data[i]->min_mapQ = mapQ; // set the mapQ filter
|
|
|
94 data[i]->min_len = min_len; // set the qlen filter
|
|
|
95 htmp = bam_header_read(data[i]->fp); // read the BAM header
|
|
|
96 if (i == 0) {
|
|
|
97 h = htmp; // keep the header of the 1st BAM
|
|
|
98 if (reg) bam_parse_region(h, reg, &tid, &beg, &end); // also parse the region
|
|
|
99 } else bam_header_destroy(htmp); // if not the 1st BAM, trash the header
|
|
|
100 if (tid >= 0) { // if a region is specified and parsed successfully
|
|
|
101 bam_index_t *idx = bam_index_load(argv[optind+i]); // load the index
|
|
|
102 data[i]->iter = bam_iter_query(idx, tid, beg, end); // set the iterator
|
|
|
103 bam_index_destroy(idx); // the index is not needed any more; phase out of the memory
|
|
|
104 }
|
|
|
105 }
|
|
|
106
|
|
|
107 // the core multi-pileup loop
|
|
|
108 mplp = bam_mplp_init(n, read_bam, (void**)data); // initialization
|
|
|
109 n_plp = calloc(n, sizeof(int)); // n_plp[i] is the number of covering reads from the i-th BAM
|
|
|
110 plp = calloc(n, sizeof(void*)); // plp[i] points to the array of covering reads (internal in mplp)
|
|
|
111 while (bam_mplp_auto(mplp, &tid, &pos, n_plp, plp) > 0) { // come to the next covered position
|
|
|
112 if (pos < beg || pos >= end) continue; // out of range; skip
|
|
|
113 if (bed && bed_overlap(bed, h->target_name[tid], pos, pos + 1) == 0) continue; // not in BED; skip
|
|
|
114 fputs(h->target_name[tid], stdout); printf("\t%d", pos+1); // a customized printf() would be faster
|
|
|
115 for (i = 0; i < n; ++i) { // base level filters have to go here
|
|
|
116 int j, m = 0;
|
|
|
117 for (j = 0; j < n_plp[i]; ++j) {
|
|
|
118 const bam_pileup1_t *p = plp[i] + j; // DON'T modfity plp[][] unless you really know
|
|
|
119 if (p->is_del || p->is_refskip) ++m; // having dels or refskips at tid:pos
|
|
|
120 else if (bam1_qual(p->b)[p->qpos] < baseQ) ++m; // low base quality
|
|
|
121 }
|
|
|
122 printf("\t%d", n_plp[i] - m); // this the depth to output
|
|
|
123 }
|
|
|
124 putchar('\n');
|
|
|
125 }
|
|
|
126 free(n_plp); free(plp);
|
|
|
127 bam_mplp_destroy(mplp);
|
|
|
128
|
|
|
129 bam_header_destroy(h);
|
|
|
130 for (i = 0; i < n; ++i) {
|
|
|
131 bam_close(data[i]->fp);
|
|
|
132 if (data[i]->iter) bam_iter_destroy(data[i]->iter);
|
|
|
133 free(data[i]);
|
|
|
134 }
|
|
|
135 free(data); free(reg);
|
|
|
136 if (bed) bed_destroy(bed);
|
|
|
137 if ( file_list )
|
|
|
138 {
|
|
|
139 for (i=0; i<n; i++) free(fn[i]);
|
|
|
140 free(fn);
|
|
|
141 }
|
|
|
142 return 0;
|
|
|
143 }
|
