17 int cur_tid; //current tid;
18 float *wig_arr; // wiggle array
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];
27 for (int i = 0; i < chr_len; i++) {
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]);
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]);
47 int main(int argc, char* argv[]) {
51 printf("Usage : rsem-bam2wig sorted_bam_input wig_output wiggle_name\n");
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);
60 fo = fopen(argv[2], "w");
61 fprintf(fo, "track type=wiggle_0 name=\"%s\" description=\"%s\" visibility=full\n", argv[3], argv[3]);
65 while (samread(bam_in, b) >= 0) {
66 if (b->core.flag & 0x0004) continue;
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);
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);
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;
86 //case BAM_CSOFT_CLIP : pos += op_len; break;
87 case BAM_CINS : pos += op_len; break;
89 for (int j = 0; j < op_len; j++, ++pos) wig_arr[pos] += w;
91 case BAM_CREF_SKIP : pos += op_len; break;
92 default : assert(false);
97 if (cnt % 1000000 == 0) printf("%d FIN\n", cnt);
99 if (cur_tid >= 0) generateWiggle(cur_tid);