From: Bo Li Date: Wed, 31 Jul 2013 07:19:15 +0000 (-0500) Subject: Install the latest version of EBSeq from Bioconductor and if fails, try to install... X-Git-Url: https://git.donarmstrong.com/?a=commitdiff_plain;h=1d786dca8b4c9c1e2176a28a669d2d16cd93e395;hp=649cbca29ee27c4c526d8500f125e3828d70e1cb;p=rsem.git Install the latest version of EBSeq from Bioconductor and if fails, try to install EBSeq v1.1.5 locally. Fixed a bug in 'rsem-gen-transcript-plots', which makes 'rsem-plot-transcript-wiggles' fail --- diff --git a/EBSeq/EBSeq_1.1.5.tar.gz b/EBSeq/EBSeq_1.1.5.tar.gz new file mode 100644 index 0000000..a5fb233 Binary files /dev/null and b/EBSeq/EBSeq_1.1.5.tar.gz differ diff --git a/EBSeq/EBSeq_1.1.6.tar.gz b/EBSeq/EBSeq_1.1.6.tar.gz deleted file mode 100644 index f2b0d5e..0000000 Binary files a/EBSeq/EBSeq_1.1.6.tar.gz and /dev/null differ diff --git a/EBSeq/install.R b/EBSeq/install.R new file mode 100755 index 0000000..2438adb --- /dev/null +++ b/EBSeq/install.R @@ -0,0 +1,19 @@ +#!/usr/bin/env Rscript + +result <- suppressWarnings(tryCatch({ + library("EBSeq", lib.loc = ".") + }, error = function(err) { + tryCatch({ + source("http://www.bioconductor.org/biocLite.R") + biocLite("EBSeq", lib = ".") + library("EBSeq", lib.loc = ".") + }, error = function(err) { + tryCatch({ + cat("Failed to install the latest version of EBSeq from Bioconductor! Try to install EBSeq v1.1.5 locally instead.\n") + install.packages(c("blockmodeling_0.1.8.tar.gz", "EBSeq_1.1.5.tar.gz"), lib = ".", repos = NULL) + library("EBSeq", lib.loc = ".") + }, error = function(err) { cat("Failed to install EBSeq v1.1.5 locally!\n") }) + }) + })) + + diff --git a/EBSeq/makefile b/EBSeq/makefile index ac97910..5bd86ff 100644 --- a/EBSeq/makefile +++ b/EBSeq/makefile @@ -1,16 +1,13 @@ CC = g++ -PROGRAMS = blockmodeling EBSeq rsem-for-ebseq-calculate-clustering-info +PROGRAMS = EBSeq rsem-for-ebseq-calculate-clustering-info all : $(PROGRAMS) -blockmodeling : blockmodeling_0.1.8.tar.gz - R CMD INSTALL -l "." blockmodeling_0.1.8.tar.gz - -EBSeq : blockmodeling EBSeq_1.1.6.tar.gz - R CMD INSTALL -l "." EBSeq_1.1.6.tar.gz +EBSeq : + ./install.R rsem-for-ebseq-calculate-clustering-info : calcClusteringInfo.cpp $(CC) -O3 -Wall calcClusteringInfo.cpp -o $@ clean : - rm -rf $(PROGRAMS) *~ + rm -rf blockmodeling $(PROGRAMS) *~ diff --git a/WHAT_IS_NEW b/WHAT_IS_NEW index 86e2194..4a777d5 100644 --- a/WHAT_IS_NEW +++ b/WHAT_IS_NEW @@ -1,3 +1,10 @@ +RSEM v1.2.6 + +- Install the latest version of EBSeq from Bioconductor and if fails, try to install EBSeq v1.1.5 locally +- Fixed a bug in 'rsem-gen-transcript-plots', which makes 'rsem-plot-transcript-wiggles' fail + +-------------------------------------------------------------------------------------------- + RSEM v1.2.5 - Updated EBSeq from v1.1.5 to v1.1.6 diff --git a/rsem-gen-transcript-plots b/rsem-gen-transcript-plots index 01bc9bd..4ea1d87 100755 --- a/rsem-gen-transcript-plots +++ b/rsem-gen-transcript-plots @@ -21,40 +21,42 @@ is_gene <- as.numeric(args[3]) show_uniq <- as.numeric(args[4]) output_plot_file <- args[5] - - 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 + readdepth <- split(data[, 2:3], data[, 1]) } 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() + tpos <- 1 # the position of transcript_id + gpos <- 2 # the position of gene_id + + data <- read.delim(file = filename, sep = "\t", stringsAsFactors = FALSE) + tmp <- aggregate(data[tpos], data[gpos], function(x) x) + t2gmap <- tmp[,2] + names(t2gmap) <- tmp[,1] + + t2gmap +} + +make_a_plot <- function(tid) { + vec <- readdepth[[tid]] + if (is.null(vec)) exit_with_error(paste("Unknown transcript detected,", tid, "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[[tid]] + 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: 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 = "") } - 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 = tid) #, xlab = "Position in transcript", ylab = "Read depth") } generate_a_page <- function(tids, gene_id = NULL) { @@ -64,30 +66,19 @@ generate_a_page <- function(tids, gene_id = NULL) { par(mfrow = c(nrow, ncol), mar = c(2, 2, 2, 2)) if (is_gene) par(oma = c(0, 0, 3, 0)) + sapply(tids, make_a_plot) + if (is_gene) mtext(gene_id, outer = TRUE, line = 1) +} - 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)) - if (len != sum(wiggle >= wiggle_uniq)) { - cat("Warning: transcript ", tids[i], " 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 = "") - } - - 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") - } +plot_a_transcript <- function(i) { + fr <- (i - 1) * num_plots_per_page + 1 + to <- min(i * num_plots_per_page, n) + generate_a_page(ids[fr:to]) +} - if (is_gene) mtext(gene_id, outer = TRUE, line = 1) +plot_a_gene <- function(gene_id) { + 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) } readdepth <- load_readdepth_file(paste(sample_name, ".transcript.readdepth", sep = "")) @@ -110,18 +101,14 @@ pdf(output_plot_file) if (!is_gene) { 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]) - } + tmp <- sapply(1:ub, plot_a_transcript) } 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) - } + tmp <- sapply(ids, plot_a_gene) } -cat("Plots are generated!\n) +cat("Plots are generated!\n") dev.off.output <- dev.off() + + +