bam2wig.cpp

   1 #include<cstdio>
   2 #include<cstring>
   3 #include<cstdlib>
   4 #include<cassert>
   5 #include<string>
   6
   7 #include<stdint.h>
   8
   9 #include "sam/bam.h"
  10 #include "sam/sam.h"
  11
  12 using namespace std;
  13
  14 samfile_t *bam_in;
  15 bam1_t *b;
  16
  17 int cur_tid; //current tid;
  18 float *wig_arr; // wiggle array
  19 FILE *fo;
  20
  21 void generateWiggle(int tid) {
  22         int chr_len = bam_in->header->target_len[tid];
  23         char *chr_name = bam_in->header->target_name[tid];
  24         int sp, ep;
  25
  26         sp = ep = -1;
  27         for (int i = 0; i < chr_len; i++) {
  28                 if (wig_arr[i] > 0) {
  29                         ep = i;
  30                 }
  31                 else {
  32                         if (sp < ep) {
  33                                 ++sp;
  34                                 fprintf(fo, "fixedStep chrom=%s start=%d step=1\n", chr_name, sp + 1);
  35                                 for (int j = sp; j <= ep; j++) fprintf(fo, "%.7g\n", wig_arr[j]);
  36                         }
  37                         sp = i;
  38                 }
  39         }
  40         if (sp < ep) {
  41                 ++sp;
  42                 fprintf(fo, "fixedStep chrom=%s start=%d step=1\n", chr_name, sp + 1);
  43                 for (int j = sp; j <= ep; j++) fprintf(fo, "%.7g\n", wig_arr[j]);
  44         }
  45 }
  46
  47 int main(int argc, char* argv[]) {
  48         int cnt = 0;
  49
  50         if (argc != 4) {
  51                 printf("Usage : rsem-bam2wig sorted_bam_input wig_output wiggle_name\n");
  52                 exit(-1);
  53         }
  54
  55         bam_in = samopen(argv[1], "rb", NULL);
  56         if (bam_in == 0) { fprintf(stderr, "Cannot open %s!\n", argv[1]); exit(-1); }
  57         //assert(bam_in != 0);
  58         b = bam_init1();
  59
  60         fo = fopen(argv[2], "w");
  61         fprintf(fo, "track type=wiggle_0 name=\"%s\" description=\"%s\" visibility=full\n", argv[3], argv[3]);
  62
  63         cur_tid = -1;
  64         wig_arr = NULL;
  65         while (samread(bam_in, b) >= 0) {
  66                 if (b->core.flag & 0x0004) continue;
  67
  68                 if (b->core.tid != cur_tid) {
  69                         if (cur_tid >= 0) generateWiggle(cur_tid);
  70                         cur_tid = b->core.tid;
  71                         size_t len = sizeof(float) * bam_in->header->target_len[cur_tid];
  72                         wig_arr = (float*)realloc(wig_arr, len);
  73                         memset(wig_arr, 0, len);
  74                 }
  75
  76                 uint8_t *p_tag = bam_aux_get(b, "ZW");
  77                 float w = (p_tag != NULL ? bam_aux2f(p_tag) : 1.0);
  78                 int pos = b->core.pos;
  79                 uint32_t *p = bam1_cigar(b);
  80
  81                 for (int i = 0; i < (int)b->core.n_cigar; i++, ++p) {
  82                         int op = *p & BAM_CIGAR_MASK;
  83                         int op_len = *p >> BAM_CIGAR_SHIFT;
  84
  85                         switch (op) {
  86                           //case BAM_CSOFT_CLIP : pos += op_len; break;
  87                         case BAM_CINS : pos += op_len; break;
  88                         case BAM_CMATCH :
  89                                 for (int j = 0; j < op_len; j++, ++pos) wig_arr[pos] += w;
  90                                 break;
  91                         case BAM_CREF_SKIP : pos += op_len; break;
  92                         default : assert(false);
  93                         }
  94                 }
  95
  96                 ++cnt;
  97                 if (cnt % 1000000 == 0) printf("%d FIN\n", cnt);
  98         }
  99         if (cur_tid >= 0) generateWiggle(cur_tid);
 100         free(wig_arr);
 101
 102         samclose(bam_in);
 103         bam_destroy1(b);
 104
 105         fclose(fo);
 106
 107         return 0;
 108 }