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
8 exit_with_error <- function(errmsg) {
9 cat(errmsg, "\n", sep = "", file = stderr())
10 quit(save = "no", status = 1)
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")
18 sample_name <- args[1]
22 output_plot_file <- args[5]
26 load_readdepth_file <- function(filename) {
27 data <- read.table(file = filename, sep = "\t", stringsAsFactors = FALSE)
31 readdepth[[data[i, 1]]] <- data[i, c(2, 3)]
36 build_t2gmap <- function(filename) {
37 data <- read.table(file = filename, sep = "\t", stringsAsFactors = FALSE)
46 if (gene_id != data[i, ncol]) {
48 t2gmap[[gene_id]] <- tids
50 gene_id <- data[i, ncol]
53 tids <- c(tids, data[i, 1])
55 if (gene_id != "") t2gmap[[gene_id]] <- tids
60 generate_a_page <- function(tids, gene_id = NULL) {
62 ncol <- ifelse(is_gene, floor(sqrt(n)), ncol_per_page)
63 nrow <- ifelse(is_gene, ceiling(n / ncol), nrow_per_page)
65 par(mfrow = c(nrow, ncol), mar = c(1, 1, 2, 1))
66 if (is_gene) par(oma = c(0, 0, 3, 0))
69 vec <- readdepth[[tids[i]]]
70 if (is.null(vec)) exit_with_error(paste("Cannot find transcript", tids[i], sep = ""))
71 wiggle <- ifelse(vec[2] == "NA", rep(0, vec[1]), as.numeric(unlist(strsplit(vec[2], split = " "))))
74 plot(wiggle, type = 'h')
77 vec <- readdepth_uniq[[tids[i]]]
78 stopifnot(!is.null(vec))
79 wiggle_uniq <- ifelse(vec[2] == "NA", rep(0, vec[1]), as.numeric(unlist(strsplit(vec[2], split = " "))))
80 stopifnot(len == length(wiggle_uniq), len == sum(wiggle >= wiggle_uniq))
81 heights <- rbind(wiggle_uniq, wiggle - wiggle_uniq)
82 barplot(heights, space = 0, border = NA, names.arg = 1:len)
84 title(main = tids[i], xlab = "Position in transcript", ylab = "Read depth")
87 if (is_gene) mtext(gene_id, outer = TRUE, line = 1)
90 readdepth <- load_readdepth_file(paste(sample_name, ".transcript.readdepth", sep = ""))
93 readdepth_uniq <- load_readdepth_file(paste(sample_name, ".uniq.transcript.readdepth", sep = ""))
96 ids <- scan(file = input_list, what = "", sep = "\n")
99 t2gmap <- build_t2gmap(paste(sample_name, ".isoforms.results", sep = ""))
102 pdf(output_plot_file)
106 ub <- (n - 1) %/% num_plots_per_page + 1
108 fr <- (i - 1) * num_plots_per_page + 1
109 to <- min(i * num_plots_per_page, n)
110 generate_a_page(ids[fr:to])
114 for (gene_id in ids) {
115 if (is.null(t2gmap[[gene_id]])) exit_with_error(paste("Cannot find gene", gene_id, sep = ""))
116 generate_a_page(t2gmap[[gene_id]], gene_id)
120 dev.off.output <- dev.off()