X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=rsem-plot-model;h=e9225eddf41466424d265cc9b55ef7ebf96b7aa4;hb=58b9a4610246143498e712ab8b817a64204d9aa6;hp=aa2c724aec31b3b143d7d10a64d8d343a53893bb;hpb=029284bf74305b5540970f858d7194d05aac501b;p=rsem.git

diff --git a/rsem-plot-model b/rsem-plot-model
index aa2c724..e9225ed 100755
--- a/rsem-plot-model
+++ b/rsem-plot-model
@@ -23,7 +23,8 @@ std <- sqrt(weighted.mean((x - mean)^2, y))
 plot(x, y, type = "h", main = "Fragment Length Distribution", sub = paste("Mean = ", mean, ", Std = ", std), xlab = "Fragment Length", ylab = "Probability")  
 
 # mate length distribution
-bval <- as.numeric(readLines(con, n = 1)[1])
+if (model_type == 0 || model_type == 1) bval <- as.numeric(readLines(con, n = 1)[1]) else bval <- 1
+
 if (bval == 1) {
   list <- strsplit(readLines(con, n = 2), split = " ")
   vec <- as.numeric(list[[1]])
@@ -41,7 +42,9 @@ bval <- as.numeric(readLines(con, n = 1)[1])
 if (bval == 1) {
   bin_size <- as.numeric(readLines(con, n = 1)[1])
   y <- as.numeric(strsplit(readLines(con, n = 1), split = " ")[[1]])
-  barplot(y, space = 0, main = "Read Start Position Distribution", xlab = "Bin #", ylab = "Probability")
+  par(cex.axis = 0.7)
+  barplot(y, space = 0, names.arg = 1:bin_size, main = "Read Start Position Distribution", xlab = "Bin #", ylab = "Probability")
+  par(cex.axis = 1.0)
 }
 strvec <- readLines(con, n = 1)
 
@@ -68,15 +71,15 @@ if (model_type == 1 || model_type == 3) {
     vecG <- as.numeric(list[[3]])
     vecT <- as.numeric(list[[4]])
     if (sum(c(vecA, vecC, vecG, vecT)) < 1e-8) break
-    peA <- c(peA, ifelse(sum(vec) < 1e-8, NA, 1.0 - vecA[1]))
-    peC <- c(peC, ifelse(sum(vec) < 1e-8, NA, 1.0 - vecC[2]))
-    peG <- c(peG, ifelse(sum(vec) < 1e-8, NA, 1.0 - vecG[3]))
-    peT <- c(peT, ifelse(sum(vec) < 1e-8, NA, 1.0 - vecT[4]))
+    peA <- c(peA, ifelse(sum(vec) < 1e-8, NA, -10 * log(1.0 - vecA[1])))
+    peC <- c(peC, ifelse(sum(vec) < 1e-8, NA, -10 * log(1.0 - vecC[2])))
+    peG <- c(peG, ifelse(sum(vec) < 1e-8, NA, -10 * log(1.0 - vecG[3])))
+    peT <- c(peT, ifelse(sum(vec) < 1e-8, NA, -10 * log(1.0 - vecT[4])))
   }
 
   x <- 0 : (length(peA) - 1)
-  matplot(x, cbind(peA, peC, peG, peT), type = "b", lty = 1:4, pch = 0:3, main = "Quality Score v.s. Percentage Sequence Error", xlab = "Quality Score", ylab = "Percentage of Sequencing Error")
-  legend("topright", c("A", "C", "G", "T"), lty = 1:4, pch = 0:3)
+  matplot(x, cbind(peA, peC, peG, peT), type = "b", lty = 1:4, pch = 0:3, col = 1:4, main = "Quality Score vs. Observed Quality", xlab = "Quality Score", ylab = "Observed Quality")
+  legend("topleft", c("A", "C", "G", "T"), lty = 1:4, pch = 0:3, col = 1:4)
 } else {
   # Profile
   readLines(con, n = 1)
@@ -94,15 +97,15 @@ if (model_type == 1 || model_type == 3) {
     vecG <- as.numeric(list[[3]])
     vecT <- as.numeric(list[[4]])
     if (sum(c(vecA, vecC, vecG, vecT)) < 1e-8) break
-    peA <- c(peA, ifelse(sum(vec) < 1e-8, NA, 1.0 - vecA[1]))
-    peC <- c(peC, ifelse(sum(vec) < 1e-8, NA, 1.0 - vecC[2]))
-    peG <- c(peG, ifelse(sum(vec) < 1e-8, NA, 1.0 - vecG[3]))
-    peT <- c(peT, ifelse(sum(vec) < 1e-8, NA, 1.0 - vecT[4]))
+    peA <- c(peA, ifelse(sum(vec) < 1e-8, NA, (1.0 - vecA[1]) * 100))
+    peC <- c(peC, ifelse(sum(vec) < 1e-8, NA, (1.0 - vecC[2]) * 100))
+    peG <- c(peG, ifelse(sum(vec) < 1e-8, NA, (1.0 - vecG[3]) * 100))
+    peT <- c(peT, ifelse(sum(vec) < 1e-8, NA, (1.0 - vecT[4]) * 100))
   }
 
   x <- 1 : length(peA)
-  matplot(x, cbind(peA, peC, peG, peT), type = "b", lty = 1:4, pch = 0:3, main = "Position v.s. Percentage Sequence Error", xlab = "Position", ylab = "Percentage of Sequencing Error")
-  legend("topleft", c("A", "C", "G", "T"), lty = 1:4, pch = 0:3)       
+  matplot(x, cbind(peA, peC, peG, peT), type = "b", lty = 1:4, pch = 0:3, col = 1:4, main = "Position vs. Percentage Sequence Error", xlab = "Position", ylab = "Percentage of Sequencing Error")
+  legend("topleft", c("A", "C", "G", "T"), lty = 1:4, pch = 0:3, col = 1:4)       
 }
 
 dev.off()