3 argv <- commandArgs(TRUE)
4 if (length(argv) != 2) {
5 cat("Usage: rsem-plot-model modelF outF\n")
9 con <- file(argv[1], open = "r")
12 # model type and forward probability
13 model_type <- as.numeric(readLines(con, n = 4)[1])
15 # fragment length distribution
16 strvec <- readLines(con, n = 3)
17 vec <- as.numeric(strsplit(strvec[1], split = " ")[[1]])
18 maxL <- vec[2] # maxL used for Profile
19 x <- (vec[1] + 1) : vec[2]
20 y <- as.numeric(strsplit(strvec[2], split = " ")[[1]])
21 mean <- weighted.mean(x, y)
22 std <- sqrt(weighted.mean((x - mean)^2, y))
23 plot(x, y, type = "h", main = "Fragment Length Distribution", sub = paste("Mean = ", mean, ", Std = ", std), xlab = "Fragment Length", ylab = "Probability")
25 # mate length distribution
26 bval <- as.numeric(readLines(con, n = 1)[1])
28 list <- strsplit(readLines(con, n = 2), split = " ")
29 vec <- as.numeric(list[[1]])
31 x <- (vec[1] + 1) : vec[2]
32 y <- as.numeric(list[[2]])
33 mean <- weighted.mean(x, y)
34 std <- sqrt(weighted.mean((x - mean)^2, y))
35 plot(x, y, type = "h", main = "Mate Length Distribution", sub = paste("Mean = ", mean, ", Std = ", std), xlab = "Mate Length", ylab = "Probability")
37 strvec <- readLines(con, n = 1)
40 bval <- as.numeric(readLines(con, n = 1)[1])
42 bin_size <- as.numeric(readLines(con, n = 1)[1])
43 y <- as.numeric(strsplit(readLines(con, n = 1), split = " ")[[1]])
45 barplot(y, space = 0, names.arg = 1:bin_size, main = "Read Start Position Distribution", xlab = "Bin #", ylab = "Probability")
47 strvec <- readLines(con, n = 1)
49 # plot sequencing errors
50 if (model_type == 1 || model_type == 3) {
52 N <- as.numeric(readLines(con, n = 1)[1])
53 readLines(con, n = N + 1)
54 readLines(con, n = 1) # for the blank line
59 peA <- c() # probability of sequencing error given reference base is A
65 strvec <- readLines(con, n = 6)
66 list <- strsplit(strvec[1:4], split = " ")
67 vecA <- as.numeric(list[[1]])
68 vecC <- as.numeric(list[[2]])
69 vecG <- as.numeric(list[[3]])
70 vecT <- as.numeric(list[[4]])
71 if (sum(c(vecA, vecC, vecG, vecT)) < 1e-8) break
72 peA <- c(peA, ifelse(sum(vec) < 1e-8, NA, -10 * log(1.0 - vecA[1])))
73 peC <- c(peC, ifelse(sum(vec) < 1e-8, NA, -10 * log(1.0 - vecC[2])))
74 peG <- c(peG, ifelse(sum(vec) < 1e-8, NA, -10 * log(1.0 - vecG[3])))
75 peT <- c(peT, ifelse(sum(vec) < 1e-8, NA, -10 * log(1.0 - vecT[4])))
78 x <- 0 : (length(peA) - 1)
79 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")
80 legend("topleft", c("A", "C", "G", "T"), lty = 1:4, pch = 0:3, col = 1:4)
85 peA <- c() # probability of sequencing error given reference base is A
91 strvec <- readLines(con, n = 6)
92 list <- strsplit(strvec[1:4], split = " ")
93 vecA <- as.numeric(list[[1]])
94 vecC <- as.numeric(list[[2]])
95 vecG <- as.numeric(list[[3]])
96 vecT <- as.numeric(list[[4]])
97 if (sum(c(vecA, vecC, vecG, vecT)) < 1e-8) break
98 peA <- c(peA, ifelse(sum(vec) < 1e-8, NA, 1.0 - vecA[1]))
99 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]))
101 peT <- c(peT, ifelse(sum(vec) < 1e-8, NA, 1.0 - vecT[4]))
105 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")
106 legend("topleft", c("A", "C", "G", "T"), lty = 1:4, pch = 0:3, col = 1:4)