]> git.donarmstrong.com Git - rsem.git/blobdiff - rsem-plot-model
Added check in rsem-plot-model for .stat output directory
[rsem.git] / rsem-plot-model
index 7466eea2398be524e3f073b1ac7d914b3c6a8d83..7ca89ec007e2889433ca4cb8111d5696344f47cc 100755 (executable)
@@ -2,15 +2,21 @@
 
 argv <- commandArgs(TRUE)
 if (length(argv) != 2) {
-  cat("Usage: rsem-plot-model sample_name outF\n")
+  cat("Usage: rsem-plot-model sample_name output_plot_file\n")
   q(status = 1)
 }
 
 strvec <- strsplit(argv[1], split = "/")[[1]]
 token <- strvec[length(strvec)]
 
-modelF <- paste(argv[1], ".stat/", token, ".model", sep = "")
-cntF <- paste(argv[1], ".stat/", token, ".cnt", sep = "")
+stat.dir <- paste(argv[1], ".stat", sep = "")
+if (!file.exists(stat.dir)) {
+  cat("Error: directory does not exist: ", stat.dir, "\n", sep = "")
+  cat(strwrap("This version of rsem-plot-model only works with the output of RSEM versions >= 1.1.8"), sep="\n")
+  q(status = 1)
+}
+modelF <- paste(stat.dir, "/", token, ".model", sep = "")
+cntF <- paste(stat.dir, "/", token, ".cnt", sep = "")
 
 pdf(argv[2])
 
@@ -82,10 +88,10 @@ if (model_type == 1 || model_type == 3) {
 
     if (sum(c(vecA, vecC, vecG, vecT)) < 1e-8) next
     x <- c(x, (i - 1))
-    peA <- c(peA, ifelse(sum(vecA) < 1e-8, NA, -10 * log(1.0 - vecA[1])))
-    peC <- c(peC, ifelse(sum(vecC) < 1e-8, NA, -10 * log(1.0 - vecC[2])))
-    peG <- c(peG, ifelse(sum(vecG) < 1e-8, NA, -10 * log(1.0 - vecG[3])))
-    peT <- c(peT, ifelse(sum(vecT) < 1e-8, NA, -10 * log(1.0 - vecT[4])))
+    peA <- c(peA, ifelse(sum(vecA) < 1e-8, NA, -10 * log10(1.0 - vecA[1])))
+    peC <- c(peC, ifelse(sum(vecC) < 1e-8, NA, -10 * log10(1.0 - vecC[2])))
+    peG <- c(peG, ifelse(sum(vecG) < 1e-8, NA, -10 * log10(1.0 - vecG[3])))
+    peT <- c(peT, ifelse(sum(vecT) < 1e-8, NA, -10 * log10(1.0 - vecT[4])))
   }
 
   matplot(x, cbind(peA, peC, peG, peT), type = "b", lty = 1:4, pch = 0:3, col = 1:4, main = "Phred Quality Score vs. Observed Quality", xlab = "Quality Score", ylab = "Observed Quality")