Modified wiggle.cpp & wiggle.h for more accurate floating point number arithmetics...

diff --git a/rsem-gen-transcript-plots b/rsem-gen-transcript-plots
index 552632a5534272dacba75f68accadbe24ffc7367..01bc9bd6d60647b4cf7929ee2c9efad661a85e39 100755 (executable)
--- a/rsem-gen-transcript-plots
+++ b/rsem-gen-transcript-plots
@@ -76,7 +76,11 @@ generate_a_page <- function(tids, gene_id = NULL) {
       vec <- readdepth_uniq[[tids[i]]]
       stopifnot(!is.null(vec))
       if (is.na(vec[[2]])) wiggle_uniq <- rep(0, vec[[1]]) else wiggle_uniq <- as.numeric(unlist(strsplit(vec[[2]], split = " ")))
-#      stopifnot(len == length(wiggle_uniq), len == sum(wiggle >= wiggle_uniq)) 
+      stopifnot(len == length(wiggle_uniq))
+      if (len != sum(wiggle >= wiggle_uniq)) {
+        cat("Warning: transcript ", tids[i], " has position(s) that read covarege with multireads is smaller than read covarge without multireads.\n", "         The 1-based position(s) is(are) : ", which(wiggle < wiggle_uniq), ".\n", "         This may be due to floating point arithmetics.\n", sep = "") 
+      }
+
       heights <- rbind(wiggle_uniq, wiggle - wiggle_uniq)      
       barplot(heights, space = 0, border = NA, names.arg = 1:len, col = c("black", "red")) 
     }

diff --git a/wiggle.cpp b/wiggle.cpp
index 4e68b4400f354c0160de91579c6d18906c641892..00dcce8d431edc2c38c0891d17b8085dbcbac25b 100644 (file)
--- a/wiggle.cpp
+++ b/wiggle.cpp
@@ -13,7 +13,7 @@
 bool no_fractional_weight = false;
 
 void add_bam_record_to_wiggle(const bam1_t *b, Wiggle& wiggle) {
-    float w;
+    double w;
 
     if (no_fractional_weight) w = 1.0;
     else {
@@ -104,14 +104,14 @@ void UCSCWiggleTrackWriter::process(const Wiggle& wiggle) {
     
     sp = ep = -1;
     for (size_t i = 0; i < wiggle.length; i++) {
-        if (wiggle.read_depth[i] > 0) {
+        if (wiggle.read_depth[i] >= 0.0095) {
             ep = i;
         }
         else {
             if (sp < ep) {
                 ++sp;
                 fprintf(fo, "fixedStep chrom=%s start=%d step=1\n", wiggle.name.c_str(), sp + 1);
-                for (int j = sp; j <= ep; j++) fprintf(fo, "%.7g\n", wiggle.read_depth[j]);
+                for (int j = sp; j <= ep; j++) fprintf(fo, "%.2f\n", wiggle.read_depth[j]);
             }
             sp = i;
         }
@@ -119,7 +119,7 @@ void UCSCWiggleTrackWriter::process(const Wiggle& wiggle) {
     if (sp < ep) {
         ++sp;
         fprintf(fo, "fixedStep chrom=%s start=%d step=1\n", wiggle.name.c_str(), sp + 1);
-        for (int j = sp; j <= ep; j++) fprintf(fo, "%.7g\n", wiggle.read_depth[j]);
+        for (int j = sp; j <= ep; j++) fprintf(fo, "%.2f\n", wiggle.read_depth[j]);
     }
 }
 

diff --git a/wiggle.h b/wiggle.h
index ecce0fcd62ceda827d9a9dd8e60ce0b59ae8115e..f94a0be39d573cd529474ad56748f4a089cf6154 100644 (file)
--- a/wiggle.h
+++ b/wiggle.h
@@ -10,7 +10,7 @@ extern bool no_fractional_weight; // if no_frac_weight == true, each alignment c
 
 struct Wiggle {
     std::string name;
-    std::vector<float> read_depth;
+    std::vector<double> read_depth;
     size_t length;
 };