Fixed a minor bug which only affects paired-end reads for reporting how many alignmen...

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

diff --git a/rsem-gen-transcript-plots b/rsem-gen-transcript-plots
index 6f7ae1038e9eded72c847b2904fcd34e3b6ab397..4ea1d87ba9ce14348dbb9d0b37173df66ee52149 100755 (executable)
--- a/rsem-gen-transcript-plots
+++ b/rsem-gen-transcript-plots
@@ -17,44 +17,46 @@ if (length(args) != 5)
 
 sample_name <- args[1]
 input_list <- args[2]
-is_gene <- args[3]
-show_uniq <- args[4]
+is_gene <- as.numeric(args[3])
+show_uniq <- as.numeric(args[4])
 output_plot_file <- args[5]
 
-
-
 load_readdepth_file <- function(filename) {
   data <- read.table(file = filename, sep = "\t", stringsAsFactors = FALSE)
-  nrow <- dim(data)[1]
-  readdepth <- list()
-  for (i in 1:nrow) {
-    readdepth[[data[i, 1]]] <- data[i, c(2, 3)]
-  }
-  readdepth
+  readdepth <- split(data[, 2:3], data[, 1])
 }
 
 build_t2gmap <- function(filename) {
-  data <- read.table(file = filename, sep = "\t", stringsAsFactors = FALSE)
-  t2gmap <- list()
-
-  nrow <- dim(data)[1]
-  ncol <- dim(data)[2]
-
-  gene_id <- ""
-  tids <- c()
-  for (i in 1:nrow) {
-    if (gene_id != data[i, ncol]) {
-      if (gene_id != "") { 
-        t2gmap[[gene_id]] <- tids
-      }
-      gene_id <- data[i, ncol]
-      tids <- c()
+  tpos <- 1   # the position of transcript_id
+  gpos <- 2   # the position of gene_id
+
+  data <- read.delim(file = filename, sep = "\t", stringsAsFactors = FALSE)
+  tmp <- aggregate(data[tpos], data[gpos], function(x) x)
+  t2gmap <- tmp[,2]
+  names(t2gmap) <- tmp[,1]
+
+  t2gmap
+}
+
+make_a_plot <- function(tid) {
+  vec <- readdepth[[tid]]      
+  if (is.null(vec)) exit_with_error(paste("Unknown transcript detected,", tid, "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) {
+    plot(wiggle, type = "h")
+  } else {
+    vec <- readdepth_uniq[[tid]]
+    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))
+    if (len != sum(wiggle >= wiggle_uniq)) {
+      cat("Warning: transcript ", tid, " 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 = "") 
     }
-    tids <- c(tids, data[i, 1])
+    heights <- rbind(wiggle_uniq, wiggle - wiggle_uniq)        
+    barplot(heights, space = 0, border = NA, names.arg = 1:len, col = c("black", "red")) 
   }
-  if (gene_id != "") t2gmap[[gene_id]] <- tids
-  
-  t2gmap
+  title(main = tid) #, xlab = "Position in transcript", ylab = "Read depth")
 }
 
 generate_a_page <- function(tids, gene_id = NULL) {
@@ -62,29 +64,21 @@ generate_a_page <- function(tids, gene_id = NULL) {
   ncol <- ifelse(is_gene, floor(sqrt(n)), ncol_per_page)
   nrow <- ifelse(is_gene, ceiling(n / ncol), nrow_per_page)
 
-  par(mfrow = c(nrow, ncol), mar = c(1, 1, 2, 1))
+  par(mfrow = c(nrow, ncol), mar = c(2, 2, 2, 2))
   if (is_gene) par(oma = c(0, 0, 3, 0)) 
+  sapply(tids, make_a_plot)
+  if (is_gene) mtext(gene_id, outer = TRUE, line = 1)
+}
 
-  for (i in 1:n) {
-    vec <- readdepth[[tids[i]]]
-    if (is.null(vec)) exit_with_error(paste("Cannot find transcript", tids[i], sep = ""))
-    wiggle <- ifelse(vec[2] == "NA", rep(0, vec[1]), as.numeric(unlist(strsplit(vec[2], split = " "))))
-    len <- length(wiggle)
-    if (!show_uniq) {
-      plot(wiggle, type = 'h')
-    }
-    else {
-      vec <- readdepth_uniq[[tids[i]]]
-      stopifnot(!is.null(vec))
-      wiggle_uniq <- ifelse(vec[2] == "NA", rep(0, vec[1]), as.numeric(unlist(strsplit(vec[2], split = " "))))
-      stopifnot(len == length(wiggle_uniq), len == sum(wiggle >= wiggle_uniq))
-      heights <- rbind(wiggle_uniq, wiggle - wiggle_uniq)      
-      barplot(heights, space = 0, border = NA, names.arg = 1:len) 
-    }
-    title(main = tids[i], xlab = "Position in transcript", ylab = "Read depth")
-  }
+plot_a_transcript <- function(i) {
+  fr <- (i - 1) * num_plots_per_page + 1
+  to <- min(i * num_plots_per_page, n)
+  generate_a_page(ids[fr:to])
+}
 
-  if (is_gene) mtext(gene_id, outer = TRUE, line = 1)
+plot_a_gene <- function(gene_id) {
+  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)
 }
 
 readdepth <- load_readdepth_file(paste(sample_name, ".transcript.readdepth", sep = ""))
@@ -95,8 +89,11 @@ if (show_uniq) {
 
 ids <- scan(file = input_list, what = "", sep = "\n")
 
+cat("Loading files is done!\n")
+
 if (is_gene) {
   t2gmap <- build_t2gmap(paste(sample_name, ".isoforms.results", sep = ""))
+  cat("Building transcript to gene map is done!\n")
 }
 
 pdf(output_plot_file)
@@ -104,21 +101,14 @@ pdf(output_plot_file)
 if (!is_gene) {        
   n <- length(ids)
   ub <- (n - 1) %/% num_plots_per_page + 1
-  for (i in 1:ub) {
-    fr <- (i - 1) * num_plots_per_page + 1
-    to <- min(i * num_plots_per_page, n)
-    generate_a_page(ids[fr:to])
-  }
-}
-else {
-  for (gene_id in ids) {
-    if (is.null(t2gmap[[gene_id]])) exit_with_error(paste("Cannot find gene", gene_id, sep = ""))
-    generate_a_page(t2gmap[[gene_id]], gene_id)
-  }
+  tmp <- sapply(1:ub, plot_a_transcript)
+} else {
+  tmp <- sapply(ids, plot_a_gene)
 }
 
-dev.off.output <- dev.off()
+cat("Plots are generated!\n")
 
+dev.off.output <- dev.off()