R/plotExons.R

In hdeberg/limma: Linear Models for Microarray Data

plotExons <- function(fit, coef = ncol(fit), geneid = NULL, genecolname = "GeneID", exoncolname = NULL, rank = 1L, FDR = 0.05)
#	Plot log2 fold changes of the exons and mark the differentially expressed exons
#	Yifang Hu and Gordon Smyth.
#	Created 1 May 2014. Last modified 7 October 2014.
{

	# Check input
	if(!is(fit, "MArrayLM")) stop("fit is not MArrayLM object")
	if(is.null(fit$p.value)) stop("fit object should contain p value from eBayes")
	if(is.null(fit$method)) stop("fit object should have fitting method from least square or robust regression")

	genecolname <- as.character(genecolname)[1]
	if(! (genecolname %in% colnames(fit$genes))) stop(paste("genecolname", genecolname, "not found"))

	if( ! is.null(exoncolname)) {

		exoncolname <- as.character(exoncolname)[1]
		if( ! (exoncolname %in% colnames(fit$genes))) stop(paste("exoncolname", exoncolname, "not found"))

	}

	# Gene to plot using rank
	if (is.null(geneid)) {

		if (rank == 1L) igene <- which.min(fit$p.value[, coef])

		else {

			o.p <- order(fit$p.value[, coef])
			fit.o <- fit[o.p,]
			fit.uniq <- fit.o[!duplicated(fit.o$genes[, genecolname]),]
			igene <- match(fit.uniq$genes[rank, genecolname], fit$genes[, genecolname])

		}

	# Gene to plot using geneid
	} else {

		geneid <- as.character(geneid)
		igene <- match(geneid[1], as.character(fit$genes[, genecolname]))
		if (is.na(igene)) stop(paste("geneid", geneid, "not found"))

	}

	# Gene annotation
	ilab <- grep(paste0(genecolname, "|geneid|symbol"), colnames(fit$genes), ignore.case = TRUE)
	genecollab <- colnames(fit$genes)[ilab]
	genelab <- paste(sapply(fit$genes[igene,genecollab], as.character), collapse = ", ")

	# Get strand if possible
	strcol <- grepl("strand", colnames(fit$genes), ignore.case = TRUE)
	if(any(strcol)) genelab <- paste0(genelab, " (", as.character(fit$genes[igene, strcol])[1], ")")

	# Exon level DE results
	fdr <- p.adjust(fit$p.value[, coef], method = "BH")
	de <- data.frame(fit$genes, logFC = fit$coefficients[, coef], fdr)
	m <- which(de[, genecolname] %in% fit$genes[igene, genecolname])
	exon.de <- de[m,]

	# Add exon identifier
	if(is.null(exoncolname)) {

		exoncolname <- "ID"
		exon.de$ID <- 1:nrow(exon.de)

	}

	# Order exon level results by exon identifier
	exon.de <- exon.de[order(exon.de[,exoncolname]),]

	# Plot exons
	plot(exon.de$logFC, xlab = "", ylab = "Log2 Fold Change", main = genelab, type = "b", xaxt = "n")
	axis(1, at = 1:nrow(exon.de), labels = exon.de[, exoncolname], las = 2, cex.axis = 0.5)
	mtext(text = paste("Exon", exoncolname), side = 1, padj = 5.2)

	# Mark DE exons
	mark <- exon.de$fdr < FDR

	if (sum(mark) > 0) {

		col <- rep(NA, nrow(exon.de))
		col[mark & (exon.de$logFC > 0)] <- "red"
		col[mark & (exon.de$logFC < 0)] <- "blue"

		if (sum(mark) == 1) cex <- 1.5 else {

			abs.fdr <- abs(log10(exon.de$fdr[mark]))
			from <- range(abs.fdr)
			to <- c(1, 2)
			cex <- (abs.fdr - from[1])/diff(from) * diff(to) + to[1]

		}

		points((1:nrow(exon.de))[mark], exon.de$logFC[mark], col = col[mark], pch = 16, cex = cex)

	}

	abline(h = mean(exon.de$logFC), lty = 2)

}