From: Bo Li <bli@cs.wisc.edu>
Date: Mon, 15 Apr 2013 16:26:02 +0000 (-0500)
Subject: Modified wiggle.cpp & wiggle.h for more accurate floating point number arithmetics... 
X-Git-Url: https://git.donarmstrong.com/?a=commitdiff_plain;h=672ce6113b1d773a4a6020b42cfb58a0a873fdec;p=rsem.git

Modified wiggle.cpp & wiggle.h for more accurate floating point number arithmetics. Modified rsem-gen-transcript-plots to output warning messages instead of exit the script when positions where multiread track has less coverage than uniq track are detected. Thanks Han Lin for pointing out the problem and suggesting possible fixes
---

diff --git a/rsem-gen-transcript-plots b/rsem-gen-transcript-plots
index 552632a..01bc9bd 100755
--- 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 4e68b44..00dcce8 100644
--- 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 ecce0fc..f94a0be 100644
--- 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;
 };