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
  25
  26 load_readdepth_file <- function(filename) {
  27   data <- read.table(file = filename, sep = "\t", stringsAsFactors = FALSE)
  28   nrow <- dim(data)[1]
  29   readdepth <- list()
  30   for (i in 1:nrow) {
  31     readdepth[[data[i, 1]]] <- data[i, c(2, 3)]
  32   }
  33   readdepth
  34 }
  35
  36 build_t2gmap <- function(filename) {
  37   data <- read.table(file = filename, sep = "\t", stringsAsFactors = FALSE)
  38   t2gmap <- list()
  39
  40   nrow <- dim(data)[1]
  41   ncol <- dim(data)[2]
  42
  43   gene_id <- ""
  44   tids <- c()
  45   for (i in 1:nrow) {
  46     if (gene_id != data[i, ncol]) {
  47       if (gene_id != "") {
  48         t2gmap[[gene_id]] <- tids
  49       }
  50       gene_id <- data[i, ncol]
  51       tids <- c()
  52     }
  53     tids <- c(tids, data[i, 1])
  54   }
  55   if (gene_id != "") t2gmap[[gene_id]] <- tids
  56
  57   t2gmap
  58 }
  59
  60 generate_a_page <- function(tids, gene_id = NULL) {
  61   n <- length(tids)
  62   ncol <- ifelse(is_gene, floor(sqrt(n)), ncol_per_page)
  63   nrow <- ifelse(is_gene, ceiling(n / ncol), nrow_per_page)
  64
  65   par(mfrow = c(nrow, ncol), mar = c(2, 2, 2, 2))
  66   if (is_gene) par(oma = c(0, 0, 3, 0))
  67
  68   for (i in 1:n) {
  69     vec <- readdepth[[tids[i]]]
  70     if (is.null(vec)) exit_with_error(paste("Cannot find transcript", tids[i], sep = ""))
  71     if (is.na(vec[[2]])) wiggle <- rep(0, vec[[1]]) else wiggle <- as.numeric(unlist(strsplit(vec[[2]], split = " ")))
  72     len <- length(wiggle)
  73     if (!show_uniq) {
  74       plot(wiggle, type = "h")
  75     } else {
  76       vec <- readdepth_uniq[[tids[i]]]
  77       stopifnot(!is.null(vec))
  78       if (is.na(vec[[2]])) wiggle_uniq <- rep(0, vec[[1]]) else wiggle_uniq <- as.numeric(unlist(strsplit(vec[[2]], split = " ")))
  79       stopifnot(len == length(wiggle_uniq), len == sum(wiggle >= wiggle_uniq))
  80       heights <- rbind(wiggle_uniq, wiggle - wiggle_uniq)
  81       barplot(heights, space = 0, border = NA, names.arg = 1:len, col = c("black", "red"))
  82     }
  83     title(main = tids[i]) #, xlab = "Position in transcript", ylab = "Read depth")
  84   }
  85
  86   if (is_gene) mtext(gene_id, outer = TRUE, line = 1)
  87 }
  88
  89 readdepth <- load_readdepth_file(paste(sample_name, ".transcript.readdepth", sep = ""))
  90
  91 if (show_uniq) {
  92   readdepth_uniq <- load_readdepth_file(paste(sample_name, ".uniq.transcript.readdepth", sep = ""))
  93 }
  94
  95 ids <- scan(file = input_list, what = "", sep = "\n")
  96
  97 cat("Loading files is done!\n")
  98
  99 if (is_gene) {
 100   t2gmap <- build_t2gmap(paste(sample_name, ".isoforms.results", sep = ""))
 101   cat("Building transcript to gene map is done!\n")
 102 }
 103
 104 pdf(output_plot_file)
 105
 106 if (!is_gene) {
 107   n <- length(ids)
 108   ub <- (n - 1) %/% num_plots_per_page + 1
 109   for (i in 1:ub) {
 110     fr <- (i - 1) * num_plots_per_page + 1
 111     to <- min(i * num_plots_per_page, n)
 112     generate_a_page(ids[fr:to])
 113   }
 114 } else {
 115   for (gene_id in ids) {
 116     if (is.null(t2gmap[[gene_id]])) exit_with_error(paste("Cannot find gene", gene_id, sep = ""))
 117     generate_a_page(t2gmap[[gene_id]], gene_id)
 118   }
 119 }
 120
 121 cat("Plots are generated!\n)
 122
 123 dev.off.output <- dev.off()
 124
 125
 126
 127
 128
 129