sam/bam2depth.c

   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 }