X-Git-Url: https://git.donarmstrong.com/?p=rsem.git;a=blobdiff_plain;f=rsem-gen-transcript-plots;h=5b8cc98434d70083c1573868b6162b612ed1a572;hp=21d1db829174ca56e3c25a3f1afe59085db492af;hb=refs%2Fheads%2Fupstream;hpb=5a2752145fda053ac654baaf1102ac40961c65b6 diff --git a/rsem-gen-transcript-plots b/rsem-gen-transcript-plots index 21d1db8..5b8cc98 100755 --- a/rsem-gen-transcript-plots +++ b/rsem-gen-transcript-plots @@ -2,88 +2,89 @@ nrow_per_page <- 3 # if input_list is composed of transcript ids ncol_per_page <- 2 # if input_list is composed of transcript ids -num_plots_per_page <- nrow_per_page * ncol_per_page # if input_list is composed of transcript ids - +num_plots_per_page <- nrow_per_page * ncol_per_page # if input_list is composed of transcript/allele ids exit_with_error <- function(errmsg) { cat(errmsg, "\n", sep = "", file = stderr()) quit(save = "no", status = 1) } - args <- commandArgs(TRUE) -if (length(args) != 5) - exit_with_error("Usage: rsem-gen-transcript-plots sample_name input_list is_gene show_uniq output_plot_file") +if (length(args) != 6) + 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") sample_name <- args[1] input_list <- args[2] -is_gene <- as.numeric(args[3]) -show_uniq <- as.numeric(args[4]) -output_plot_file <- args[5] - - +alleleS <- as.numeric(args[3]) +id_type <- as.numeric(args[4]) +show_uniq <- as.numeric(args[5]) +output_plot_file <- args[6] +is_composite <- (!alleleS && (id_type == 2)) || (alleleS && (id_type > 0)) + 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 <- split(data[, 2:3], data[, 1]) +} + +build_map <- function(filename) { + value_pos <- 1 + if (!alleleS || (alleleS && id_type == 1)) { + key_pos <- 2 + } else { + key_pos <- 3 } - readdepth + + data <- read.delim(file = filename, sep = "\t", stringsAsFactors = FALSE) + tmp <- aggregate(data[value_pos], data[key_pos], function(x) x) + map <- tmp[,2] + names(map) <- tmp[,1] + + map } -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() +make_a_plot <- function(id) { + vec <- readdepth[[id]] + 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)) + 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[[id]] + 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: ", 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 = "") } - 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 = id) } -generate_a_page <- function(tids, gene_id = NULL) { - n <- length(tids) - ncol <- ifelse(is_gene, floor(sqrt(n)), ncol_per_page) - nrow <- ifelse(is_gene, ceiling(n / ncol), nrow_per_page) +generate_a_page <- function(ids, title = NULL) { + n <- length(ids) + ncol <- ifelse(is_composite, floor(sqrt(n)), ncol_per_page) + nrow <- ifelse(is_composite, ceiling(n / ncol), nrow_per_page) par(mfrow = c(nrow, ncol), mar = c(2, 2, 2, 2)) - if (is_gene) par(oma = c(0, 0, 3, 0)) - - for (i in 1:n) { - vec <- readdepth[[tids[i]]] - if (is.null(vec)) exit_with_error(paste("Unknown transcript detected,", tids[i], "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[[tids[i]]] - 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), len == sum(wiggle >= wiggle_uniq)) - heights <- rbind(wiggle_uniq, wiggle - wiggle_uniq) - barplot(heights, space = 0, border = NA, names.arg = 1:len, col = c("black", "red")) - } - title(main = tids[i]) #, xlab = "Position in transcript", ylab = "Read depth") - } + if (is_composite) par(oma = c(0, 0, 3, 0)) + sapply(ids, make_a_plot) + if (is_composite) mtext(title, outer = TRUE, line = 1) +} - if (is_gene) mtext(gene_id, outer = TRUE, line = 1) +plot_individual <- function(i) { + fr <- (i - 1) * num_plots_per_page + 1 + to <- min(i * num_plots_per_page, n) + generate_a_page(ids[fr:to]) +} + +# cid, composite id, can be either a gene id or transcript id (for allele-specific expression only) +plot_composite <- function(cid) { + 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)) + generate_a_page(map[[cid]], cid) } readdepth <- load_readdepth_file(paste(sample_name, ".transcript.readdepth", sep = "")) @@ -96,28 +97,22 @@ 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 = "")) +if (is_composite) { + file_name <- sprintf("%s.%s.results", sample_name, ifelse(alleleS, "alleles", "isoforms")) + map <- build_map(file_name) cat("Building transcript to gene map is done!\n") } pdf(output_plot_file) -if (!is_gene) { +if (!is_composite) { 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]) - } + dumbvar <- sapply(1:ub, plot_individual) } else { - for (gene_id in ids) { - 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) - } + dumbvar <- sapply(ids, plot_composite) } -cat("Plots are generated!\n) +cat("Plots are generated!\n") dev.off.output <- dev.off()