bam2wig.cpp

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