sam/bedcov.c

   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 }