rsem-gen-transcript-plots

   1 #!/usr/bin/env Rscript
   2
   3 nrow_per_page <- 3 # if input_list is composed of transcript ids
   4 ncol_per_page <- 2 # if input_list is composed of transcript ids
   5 num_plots_per_page <- nrow_per_page * ncol_per_page # if input_list is composed of transcript ids
   6
   7
   8 exit_with_error <- function(errmsg) {
   9   cat(errmsg, "\n", sep = "", file = stderr())
  10   quit(save = "no", status = 1)
  11 }
  12
  13
  14 args <- commandArgs(TRUE)
  15 if (length(args) != 5)
  16   exit_with_error("Usage: rsem-gen-transcript-plots sample_name input_list is_gene show_uniq output_plot_file")
  17
  18 sample_name <- args[1]
  19 input_list <- args[2]
  20 is_gene <- as.numeric(args[3])
  21 show_uniq <- as.numeric(args[4])
  22 output_plot_file <- args[5]
  23
  24 load_readdepth_file <- function(filename) {
  25   data <- read.table(file = filename, sep = "\t", stringsAsFactors = FALSE)
  26   readdepth <- split(data[, 2:3], data[, 1])
  27 }
  28
  29 build_t2gmap <- function(filename) {
  30   tpos <- 1   # the position of transcript_id
  31   gpos <- 2   # the position of gene_id
  32
  33   data <- read.delim(file = filename, sep = "\t", stringsAsFactors = FALSE)
  34   tmp <- aggregate(data[tpos], data[gpos], function(x) x)
  35   t2gmap <- tmp[,2]
  36   names(t2gmap) <- tmp[,1]
  37
  38   t2gmap
  39 }
  40
  41 make_a_plot <- function(tid) {
  42   vec <- readdepth[[tid]]
  43   if (is.null(vec)) exit_with_error(paste("Unknown transcript detected,", tid, "is not included in RSEM's indices."))
  44   if (is.na(vec[[2]])) wiggle <- rep(0, vec[[1]]) else wiggle <- as.numeric(unlist(strsplit(vec[[2]], split = " ")))
  45   len <- length(wiggle)
  46   if (!show_uniq) {
  47     plot(wiggle, type = "h")
  48   } else {
  49     vec <- readdepth_uniq[[tid]]
  50     stopifnot(!is.null(vec))
  51     if (is.na(vec[[2]])) wiggle_uniq <- rep(0, vec[[1]]) else wiggle_uniq <- as.numeric(unlist(strsplit(vec[[2]], split = " ")))
  52     stopifnot(len == length(wiggle_uniq))
  53     if (len != sum(wiggle >= wiggle_uniq)) {
  54       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 = "")
  55     }
  56     heights <- rbind(wiggle_uniq, wiggle - wiggle_uniq)
  57     barplot(heights, space = 0, border = NA, names.arg = 1:len, col = c("black", "red"))
  58   }
  59   title(main = tid) #, xlab = "Position in transcript", ylab = "Read depth")
  60 }
  61
  62 generate_a_page <- function(tids, gene_id = NULL) {
  63   n <- length(tids)
  64   ncol <- ifelse(is_gene, floor(sqrt(n)), ncol_per_page)
  65   nrow <- ifelse(is_gene, ceiling(n / ncol), nrow_per_page)
  66
  67   par(mfrow = c(nrow, ncol), mar = c(2, 2, 2, 2))
  68   if (is_gene) par(oma = c(0, 0, 3, 0))
  69   sapply(tids, make_a_plot)
  70   if (is_gene) mtext(gene_id, outer = TRUE, line = 1)
  71 }
  72
  73 plot_a_transcript <- function(i) {
  74   fr <- (i - 1) * num_plots_per_page + 1
  75   to <- min(i * num_plots_per_page, n)
  76   generate_a_page(ids[fr:to])
  77 }
  78
  79 plot_a_gene <- function(gene_id) {
  80   if (is.null(t2gmap[[gene_id]])) exit_with_error(paste("Unknown gene detected,", gene_id, "is not included in RSEM's in indices."))
  81   generate_a_page(t2gmap[[gene_id]], gene_id)
  82 }
  83
  84 readdepth <- load_readdepth_file(paste(sample_name, ".transcript.readdepth", sep = ""))
  85
  86 if (show_uniq) {
  87   readdepth_uniq <- load_readdepth_file(paste(sample_name, ".uniq.transcript.readdepth", sep = ""))
  88 }
  89
  90 ids <- scan(file = input_list, what = "", sep = "\n")
  91
  92 cat("Loading files is done!\n")
  93
  94 if (is_gene) {
  95   t2gmap <- build_t2gmap(paste(sample_name, ".isoforms.results", sep = ""))
  96   cat("Building transcript to gene map is done!\n")
  97 }
  98
  99 pdf(output_plot_file)
 100
 101 if (!is_gene) {
 102   n <- length(ids)
 103   ub <- (n - 1) %/% num_plots_per_page + 1
 104   tmp <- sapply(1:ub, plot_a_transcript)
 105 } else {
 106   tmp <- sapply(ids, plot_a_gene)
 107 }
 108
 109 cat("Plots are generated!\n")
 110
 111 dev.off.output <- dev.off()
 112
 113
 114