Updated EBSeq from v1.1.5 to v1.1.6 and fixed a bug in 'rsem-generate-data-matrix...

[rsem.git] / rsem-gen-transcript-plots

diff --git a/rsem-gen-transcript-plots b/rsem-gen-transcript-plots
index 9a806d6c191710482b896105c3e965a49828af3e..01bc9bd6d60647b4cf7929ee2c9efad661a85e39 100755 (executable)
--- a/rsem-gen-transcript-plots
+++ b/rsem-gen-transcript-plots
@@ -67,7 +67,7 @@ generate_a_page <- function(tids, gene_id = NULL) {
 
   for (i in 1:n) {
     vec <- readdepth[[tids[i]]]
-    if (is.null(vec)) exit_with_error(paste("Cannot find transcript", tids[i], sep = ""))
+    if (is.null(vec)) exit_with_error(paste("Unknown transcript detected,", tids[i], "is not included in RSEM's indices."))
     if (is.na(vec[[2]])) wiggle <- rep(0, vec[[1]]) else wiggle <- as.numeric(unlist(strsplit(vec[[2]], split = " ")))
     len <- length(wiggle)
     if (!show_uniq) {
@@ -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")) 
     }
@@ -113,7 +117,7 @@ if (!is_gene) {
   }
 } else {
   for (gene_id in ids) {
-    if (is.null(t2gmap[[gene_id]])) exit_with_error(paste("Cannot find gene", gene_id, sep = ""))
+    if (is.null(t2gmap[[gene_id]])) exit_with_error(paste("Unknown gene detected,", gene_id, "is not included in RSEM's in indices."))
     generate_a_page(t2gmap[[gene_id]], gene_id)
   }
 }
@@ -121,9 +125,3 @@ if (!is_gene) {
 cat("Plots are generated!\n)
 
 dev.off.output <- dev.off()
-
-
-
-
-
-