bedcov.c

   1 #include <zlib.h>
   2 #include <stdio.h>
   3 #include <ctype.h>
   4 #include <stdlib.h>
   5 #include <string.h>
   6 #include "kstring.h"
   7 #include "bgzf.h"
   8 #include "bam.h"
   9
  10 #include "kseq.h"
  11 KSTREAM_INIT(gzFile, gzread, 16384)
  12
  13 typedef struct {
  14         bamFile fp;
  15         bam_iter_t iter;
  16         int min_mapQ;
  17 } aux_t;
  18
  19 static int read_bam(void *data, bam1_t *b)
  20 {
  21         aux_t *aux = (aux_t*)data;
  22         int ret = bam_iter_read(aux->fp, aux->iter, b);
  23         if ((int)b->core.qual < aux->min_mapQ) b->core.flag |= BAM_FUNMAP;
  24         return ret;
  25 }
  26
  27 int main_bedcov(int argc, char *argv[])
  28 {
  29         extern void bam_init_header_hash(bam_header_t*);
  30         gzFile fp;
  31         kstring_t str;
  32         kstream_t *ks;
  33         bam_index_t **idx;
  34         bam_header_t *h = 0;
  35         aux_t **aux;
  36         int *n_plp, dret, i, n, c, min_mapQ = 0;
  37         int64_t *cnt;
  38         const bam_pileup1_t **plp;
  39
  40         while ((c = getopt(argc, argv, "Q:")) >= 0) {
  41                 switch (c) {
  42                 case 'Q': min_mapQ = atoi(optarg); break;
  43                 }
  44         }
  45         if (optind + 2 > argc) {
  46                 fprintf(stderr, "Usage: samtools bedcov <in.bed> <in1.bam> [...]\n");
  47                 return 1;
  48         }
  49         memset(&str, 0, sizeof(kstring_t));
  50         n = argc - optind - 1;
  51         aux = calloc(n, sizeof(void*));
  52         idx = calloc(n, sizeof(void*));
  53         for (i = 0; i < n; ++i) {
  54                 aux[i] = calloc(1, sizeof(aux_t));
  55                 aux[i]->min_mapQ = min_mapQ;
  56                 aux[i]->fp = bam_open(argv[i+optind+1], "r");
  57                 idx[i] = bam_index_load(argv[i+optind+1]);
  58                 if (aux[i]->fp == 0 || idx[i] == 0) {
  59                         fprintf(stderr, "ERROR: fail to open index BAM file '%s'\n", argv[i+optind+1]);
  60                         return 2;
  61                 }
  62                 bgzf_set_cache_size(aux[i]->fp, 20);
  63                 if (i == 0) h = bam_header_read(aux[0]->fp);
  64         }
  65         bam_init_header_hash(h);
  66         cnt = calloc(n, 8);
  67
  68         fp = gzopen(argv[optind], "rb");
  69         ks = ks_init(fp);
  70         n_plp = calloc(n, sizeof(int));
  71         plp = calloc(n, sizeof(void*));
  72         while (ks_getuntil(ks, KS_SEP_LINE, &str, &dret) >= 0) {
  73                 char *p, *q;
  74                 int tid, beg, end, pos;
  75                 bam_mplp_t mplp;
  76
  77                 for (p = q = str.s; *p && *p != '\t'; ++p);
  78                 if (*p != '\t') goto bed_error;
  79                 *p = 0; tid = bam_get_tid(h, q); *p = '\t';
  80                 if (tid < 0) goto bed_error;
  81                 for (q = p = p + 1; isdigit(*p); ++p);
  82                 if (*p != '\t') goto bed_error;
  83                 *p = 0; beg = atoi(q); *p = '\t';
  84                 for (q = p = p + 1; isdigit(*p); ++p);
  85                 if (*p == '\t' || *p == 0) {
  86                         int c = *p;
  87                         *p = 0; end = atoi(q); *p = c;
  88                 } else goto bed_error;
  89
  90                 for (i = 0; i < n; ++i) {
  91                         if (aux[i]->iter) bam_iter_destroy(aux[i]->iter);
  92                         aux[i]->iter = bam_iter_query(idx[i], tid, beg, end);
  93                 }
  94                 mplp = bam_mplp_init(n, read_bam, (void**)aux);
  95                 bam_mplp_set_maxcnt(mplp, 64000);
  96                 memset(cnt, 0, 8 * n);
  97                 while (bam_mplp_auto(mplp, &tid, &pos, n_plp, plp) > 0)
  98                         if (pos >= beg && pos < end)
  99                                 for (i = 0; i < n; ++i) cnt[i] += n_plp[i];
 100                 for (i = 0; i < n; ++i) {
 101                         kputc('\t', &str);
 102                         kputl(cnt[i], &str);
 103                 }
 104                 puts(str.s);
 105                 bam_mplp_destroy(mplp);
 106                 continue;
 107
 108 bed_error:
 109                 fprintf(stderr, "Errors in BED line '%s'\n", str.s);
 110         }
 111         free(n_plp); free(plp);
 112         ks_destroy(ks);
 113         gzclose(fp);
 114
 115         free(cnt);
 116         for (i = 0; i < n; ++i) {
 117                 if (aux[i]->iter) bam_iter_destroy(aux[i]->iter);
 118                 bam_index_destroy(idx[i]);
 119                 bam_close(aux[i]->fp);
 120                 free(aux[i]);
 121         }
 122         bam_header_destroy(h);
 123         free(aux); free(idx);
 124         free(str.s);
 125         return 0;
 126 }