- barcodeplot() is now able to plot different weights for different
genes, as specified by the new argument 'gene.weights'.
- new function plotExons() to plot log-fold-changes by exon for a
given gene.
- plotSplice() now plots strand information if this can be found.
- first argument of plotMA() and plotMA3by2() is now 'object' instead
of 'MA'.
- mdplot() now accepts arguments 'xlab' and 'ylab'.
- all default methods for S3 generic functions are now registered in
NAMESPACE.
git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/limma@95275 bc3139a8-67e5-0310-9ffc-ced21a209358
Gordon Smyth authored on 10/10/2014 07:17:04
| 1 | 1 |
new file mode 100644 |
| ... | ... |
@@ -0,0 +1,100 @@ |
| 1 |
+plotExons <- function(fit, coef = ncol(fit), geneid = NULL, genecolname = "GeneID", exoncolname = NULL, rank = 1L, FDR = 0.05) |
|
| 2 |
+# Plot log2 fold changes of the exons and mark the differentially expressed exons |
|
| 3 |
+# Yifang Hu and Gordon Smyth. |
|
| 4 |
+# Created 1 May 2014. Last modified 7 October 2014. |
|
| 5 |
+{
|
|
| 6 |
+ |
|
| 7 |
+ # Check input |
|
| 8 |
+ if(!is(fit, "MArrayLM")) stop("fit is not MArrayLM object")
|
|
| 9 |
+ if(is.null(fit$p.value)) stop("fit object should contain p value from eBayes")
|
|
| 10 |
+ if(is.null(fit$method)) stop("fit object should have fitting method from least square or robust regression")
|
|
| 11 |
+ |
|
| 12 |
+ genecolname <- as.character(genecolname)[1] |
|
| 13 |
+ if(! (genecolname %in% colnames(fit$genes))) stop(paste("genecolname", genecolname, "not found"))
|
|
| 14 |
+ |
|
| 15 |
+ if( ! is.null(exoncolname)) {
|
|
| 16 |
+ |
|
| 17 |
+ exoncolname <- as.character(exoncolname)[1] |
|
| 18 |
+ if( ! (exoncolname %in% colnames(fit$genes))) stop(paste("exoncolname", exoncolname, "not found"))
|
|
| 19 |
+ |
|
| 20 |
+ } |
|
| 21 |
+ |
|
| 22 |
+ # Gene to plot using rank |
|
| 23 |
+ if (is.null(geneid)) {
|
|
| 24 |
+ |
|
| 25 |
+ if (rank == 1L) igene <- which.min(fit$p.value[, coef]) |
|
| 26 |
+ |
|
| 27 |
+ else {
|
|
| 28 |
+ |
|
| 29 |
+ o.p <- order(fit$p.value[, coef]) |
|
| 30 |
+ fit.o <- fit[o.p,] |
|
| 31 |
+ fit.uniq <- fit.o[!duplicated(fit.o$genes[, genecolname]),] |
|
| 32 |
+ igene <- match(fit.uniq$genes[rank, genecolname], fit$genes[, genecolname]) |
|
| 33 |
+ |
|
| 34 |
+ } |
|
| 35 |
+ |
|
| 36 |
+ # Gene to plot using geneid |
|
| 37 |
+ } else {
|
|
| 38 |
+ |
|
| 39 |
+ geneid <- as.character(geneid) |
|
| 40 |
+ igene <- match(geneid[1], as.character(fit$genes[, genecolname])) |
|
| 41 |
+ if (is.na(igene)) stop(paste("geneid", geneid, "not found"))
|
|
| 42 |
+ |
|
| 43 |
+ } |
|
| 44 |
+ |
|
| 45 |
+ # Gene annotation |
|
| 46 |
+ ilab <- grep(paste0(genecolname, "|geneid|symbol"), colnames(fit$genes), ignore.case = TRUE) |
|
| 47 |
+ genecollab <- colnames(fit$genes)[ilab] |
|
| 48 |
+ genelab <- paste(sapply(fit$genes[igene,genecollab], as.character), collapse = ", ") |
|
| 49 |
+ |
|
| 50 |
+ # Get strand if possible |
|
| 51 |
+ strcol <- grepl("strand", colnames(fit$genes), ignore.case = TRUE)
|
|
| 52 |
+ if(any(strcol)) genelab <- paste0(genelab, " (", as.character(fit$genes[igene, strcol])[1], ")")
|
|
| 53 |
+ |
|
| 54 |
+ # Exon level DE results |
|
| 55 |
+ fdr <- p.adjust(fit$p.value[, coef], method = "BH") |
|
| 56 |
+ de <- data.frame(fit$genes, logFC = fit$coefficients[, coef], fdr) |
|
| 57 |
+ m <- which(de[, genecolname] %in% fit$genes[igene, genecolname]) |
|
| 58 |
+ exon.de <- de[m,] |
|
| 59 |
+ |
|
| 60 |
+ # Add exon identifier |
|
| 61 |
+ if(is.null(exoncolname)) {
|
|
| 62 |
+ |
|
| 63 |
+ exoncolname <- "ID" |
|
| 64 |
+ exon.de$ID <- 1:nrow(exon.de) |
|
| 65 |
+ |
|
| 66 |
+ } |
|
| 67 |
+ |
|
| 68 |
+ # Order exon level results by exon identifier |
|
| 69 |
+ exon.de <- exon.de[order(exon.de[,exoncolname]),] |
|
| 70 |
+ |
|
| 71 |
+ # Plot exons |
|
| 72 |
+ plot(exon.de$logFC, xlab = "", ylab = "Log2 Fold Change", main = genelab, type = "b", xaxt = "n") |
|
| 73 |
+ axis(1, at = 1:nrow(exon.de), labels = exon.de[, exoncolname], las = 2, cex.axis = 0.5) |
|
| 74 |
+ mtext(text = paste("Exon", exoncolname), side = 1, padj = 5.2)
|
|
| 75 |
+ |
|
| 76 |
+ # Mark DE exons |
|
| 77 |
+ mark <- exon.de$fdr < FDR |
|
| 78 |
+ |
|
| 79 |
+ if (sum(mark) > 0) {
|
|
| 80 |
+ |
|
| 81 |
+ col <- rep(NA, nrow(exon.de)) |
|
| 82 |
+ col[mark & (exon.de$logFC > 0)] <- "red" |
|
| 83 |
+ col[mark & (exon.de$logFC < 0)] <- "blue" |
|
| 84 |
+ |
|
| 85 |
+ if (sum(mark) == 1) cex <- 1.5 else {
|
|
| 86 |
+ |
|
| 87 |
+ abs.fdr <- abs(log10(exon.de$fdr[mark])) |
|
| 88 |
+ from <- range(abs.fdr) |
|
| 89 |
+ to <- c(1, 2) |
|
| 90 |
+ cex <- (abs.fdr - from[1])/diff(from) * diff(to) + to[1] |
|
| 91 |
+ |
|
| 92 |
+ } |
|
| 93 |
+ |
|
| 94 |
+ points((1:nrow(exon.de))[mark], exon.de$logFC[mark], col = col[mark], pch = 16, cex = cex) |
|
| 95 |
+ |
|
| 96 |
+ } |
|
| 97 |
+ |
|
| 98 |
+ abline(h = mean(exon.de$logFC), lty = 2) |
|
| 99 |
+ |
|
| 100 |
+} |