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