install the rsem_perl_utils.pm to /usr/bin even though this is wrong

[rsem.git] / bam2wig.cpp

diff --git a/bam2wig.cpp b/bam2wig.cpp
index e46fa67bc3812c8ca6978aeb8c6764680c3549b9..d205e8f7f4e9fde1cacef0ace99b626b18cc01d4 100644 (file)
--- a/bam2wig.cpp
+++ b/bam2wig.cpp
@@ -1,108 +1,27 @@
-#include<cstdio>
-#include<cstring>
-#include<cstdlib>
-#include<cassert>
-#include<string>
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
 
-#include<stdint.h>
-
-#include "sam/bam.h"
-#include "sam/sam.h"
+#include "wiggle.h"
 
 using namespace std;
 
-samfile_t *bam_in;
-bam1_t *b;
-
-int cur_tid; //current tid;
-float *wig_arr; // wiggle array
-FILE *fo;
-
-void generateWiggle(int tid) {
-       int chr_len = bam_in->header->target_len[tid];
-       char *chr_name = bam_in->header->target_name[tid];
-       int sp, ep;
-
-       sp = ep = -1;
-       for (int i = 0; i < chr_len; i++) {
-               if (wig_arr[i] > 0) {
-                       ep = i;
-               }
-               else {
-                       if (sp < ep) {
-                               ++sp;
-                               fprintf(fo, "fixedStep chrom=%s start=%d step=1\n", chr_name, sp + 1);
-                               for (int j = sp; j <= ep; j++) fprintf(fo, "%.7g\n", wig_arr[j]);
-                       }
-                       sp = i;
-               }
-       }
-       if (sp < ep) {
-               ++sp;
-               fprintf(fo, "fixedStep chrom=%s start=%d step=1\n", chr_name, sp + 1);
-               for (int j = sp; j <= ep; j++) fprintf(fo, "%.7g\n", wig_arr[j]);
-       }
+void printUsage() {
+  printf("Usage: rsem-bam2wig sorted_bam_input wig_output wiggle_name [--no-fractional-weight]\n");
+  printf("sorted_bam_input\t: Input BAM format file, must be sorted\n");
+  printf("wig_output\t\t: Output wiggle file's name, e.g. output.wig\n");
+  printf("wiggle_name\t\t: the name of this wiggle plot\n");
+  printf("--no-fractional-weight\t: If this is set, RSEM will not look for \"ZW\" tag and each alignment appeared in the BAM file has weight 1. Set this if your BAM file is not generated by RSEM. Please note that this option must be at the end of the command line.\n");
+  exit(-1);
 }
 
 int main(int argc, char* argv[]) {
-       int cnt = 0;
-
-       if (argc != 4) {
-               printf("Usage : rsem-bam2wig sorted_bam_input wig_output wiggle_name\n");
-               exit(-1);
-       }
-
-       bam_in = samopen(argv[1], "rb", NULL);
-       if (bam_in == 0) { fprintf(stderr, "Cannot open %s!\n", argv[1]); exit(-1); }
-       //assert(bam_in != 0);
-       b = bam_init1();
-
-       fo = fopen(argv[2], "w");
-       fprintf(fo, "track type=wiggle_0 name=\"%s\" description=\"%s\" visibility=full\n", argv[3], argv[3]);
-
-       cur_tid = -1;
-       wig_arr = NULL;
-       while (samread(bam_in, b) >= 0) {
-               if (b->core.flag & 0x0004) continue;
-               
-               if (b->core.tid != cur_tid) {
-                       if (cur_tid >= 0) generateWiggle(cur_tid);
-                       cur_tid = b->core.tid;
-                       size_t len = sizeof(float) * bam_in->header->target_len[cur_tid];
-                       wig_arr = (float*)realloc(wig_arr, len);
-                       memset(wig_arr, 0, len);
-               }
-
-               uint8_t *p_tag = bam_aux_get(b, "ZW");
-               float w = (p_tag != NULL ? bam_aux2f(p_tag) : 1.0);
-               int pos = b->core.pos;
-               uint32_t *p = bam1_cigar(b);
-
-               for (int i = 0; i < (int)b->core.n_cigar; i++, ++p) {
-                       int op = *p & BAM_CIGAR_MASK;
-                       int op_len = *p >> BAM_CIGAR_SHIFT;
-
-                       switch (op) {
-                         //case BAM_CSOFT_CLIP : pos += op_len; break;
-                       case BAM_CINS : pos += op_len; break;
-                       case BAM_CMATCH :
-                               for (int j = 0; j < op_len; j++, ++pos) wig_arr[pos] += w;
-                               break;
-                       case BAM_CREF_SKIP : pos += op_len; break;
-                       default : assert(false);
-                       }
-               }
-
-               ++cnt;
-               if (cnt % 1000000 == 0) printf("%d FIN\n", cnt);
-       }
-       if (cur_tid >= 0) generateWiggle(cur_tid);
-       free(wig_arr);
-
-       samclose(bam_in);
-       bam_destroy1(b);
+       if (argc < 4 || argc > 5) { printf("Number of arguments is not correct!\n"); printUsage(); }
+       if (argc == 5 && strcmp(argv[4], "--no-fractional-weight")) { printf("Cannot recognize option %s!\n", argv[4]); printUsage(); }
 
-       fclose(fo);
+       no_fractional_weight = (argc == 5 && !strcmp(argv[4], "--no-fractional-weight"));
+       UCSCWiggleTrackWriter track_writer(argv[2], argv[3]);
+       build_wiggles(argv[1], track_writer);
 
        return 0;
 }