From 672ce6113b1d773a4a6020b42cfb58a0a873fdec Mon Sep 17 00:00:00 2001
From: Bo Li <bli@cs.wisc.edu>
Date: Mon, 15 Apr 2013 11:26:02 -0500
Subject: [PATCH] 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

---
 rsem-gen-transcript-plots | 6 +++++-
 wiggle.cpp                | 8 ++++----
 wiggle.h                  | 2 +-
 3 files changed, 10 insertions(+), 6 deletions(-)

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;
 };
 
-- 
2.39.5