R/plotContrasts.R

In variancePartition: Quantify and interpret divers of variation in multilevel gene expression experiments

#' Plot representation of contrast matrix
#' 
#' Plot contrast matrix to clarify interpretation of hypothesis tests with linear contrasts
#'
#' @param L contrast matrix
#' 
#' @return 
#' ggplot2 object
#'
#' @examples
#'
#' # load library
#' # library(variancePartition)
#'
#' # load simulated data:
#' # geneExpr: matrix of gene expression values
#' # info: information/metadata about each sample
#' data(varPartData)
#' 
#' # get contrast matrix testing if the coefficient for Batch2 is zero 
#' form <- ~ Batch + (1|Individual) + (1|Tissue) 
#' L1 = getContrast( geneExpr, form, info, "Batch3")
#' 
#' # get contrast matrix testing if the coefficient for Batch2 is different from Batch3 
#' form <- ~ Batch + (1|Individual) + (1|Tissue) 
#' L2 = getContrast( geneExpr, form, info, c("Batch2", "Batch3"))
#' 
#' # combine contrasts into single matrix
#' L_combined = cbind(L1, L2)
#' 
#' # plot contrasts
#' plotContrasts( L_combined )
#' 
#' @importFrom reshape2 melt
#' @import ggplot2
#' @export
plotContrasts = function( L ){

	if( is(L, "numeric") ){
		L = as.matrix(L, ncol=1)
	}
	if( is.null(colnames(L)) ){
		colnames(L) = paste0('L', seq_len(ncol(L)))
	}

	# check rownames of contrasts
	if( length(unique(colnames(L))) != ncol(L) ){
		stop(paste("Contrast names must be unique: ", paste(colnames(L), collapse=', ')))
	}

	# check that each contrast sum to zero
	tol = sqrt(.Machine$double.eps)
	sumZero = apply(L, 2, function(x){
		abs(sum(x)) < tol
		})

	# if( any(!sumZero) ){
	# 	stop('Each contrast must sum to 0. ', paste(names(sumZero[!sumZero]), collapse=', '), ' fails')
	# }

	df = melt(t(L))
	colnames(df)[1:2] = c("Var1", "Var2")
	df$Var1 = factor(df$Var1)

	if( identical(levels(df$Var1), "1") ){
		df$Var1 = factor(rep('', nrow(df)))
	}

	Var1 = Var2 = value = NULL

	h = length(unique(df$Var1))
	w = length(unique(df$Var2))

	ggplot(df, aes(Var2, y=Var1, fill=value)) + geom_tile(color="black") + theme_minimal() + theme(aspect.ratio=h/w,
		panel.grid.major =element_blank(),
		panel.grid.minor =element_blank(),plot.title = element_text(hjust = 0.5)) + scale_fill_gradient2(name="Contrast coef", limits=c(-1,1), low=alpha("blue",.8), mid="white", high=alpha("red", .8)) + xlab("Variable") + ylab("Contrasts") + ggtitle("Graphical representation of linear contrasts") + geom_text(aes(label=round(value, 2)), fontface = "bold")
}