#' Plot phenotypic effect for two markers as a bar plot
#'
#' This internal function is called by
#' \code{\link{plot_effects}} to generate a
#' bar plot showing mean trait values for all
#' combinations of genotypes for two markers.
#' This function also indicates the predicted
#' additive effects of the two markers with
#' a dashed line, as well as the error of the
#' predicted additive effect with an orange box.
#' The true effect of the marker combination is
#' shown with a gray bar.
#' The effects can centered on the reference
#' allele to better show relative effects of
#' each allele and allele combination.
#'
#' @param phenoV A vector of trait values
#' @param marker1_vals A vector of genotype values
#' for marker1
#' @param marker2_vals A vector of genotype values
#' for marker2.
#' @param pheno_name A string indicating the name of
#' the trait being plotted.
#' @param marker1_label A string indicating the name
#' of marker1
#' @param marker2_label A string indicating the name
#' of marker2
#' @param ymin A numeric value indicating the minimum
#' y value for the plot. If NULL, it will be calculated
#' based on phenoV and error bars.
#' @param ymax A numeric value indicating the maximum
#' y value for the plot. If NULL, it will be calculated
#' based on phenoV and error bars.
#' @param error_bars A string indicating the type of error
#' bars to draw. Can be "sd" for standard deviation, "se"
#' for standard error, or "none".
#' @param ref_centered Whether to center the effects
#' on the reference genotype.
#' @return None
#'
#' @importFrom graphics mtext
#' @keywords internal
plot_bars <- function(phenoV, marker1_vals, marker2_vals, pheno_name, marker1_label,
marker2_label, ymin = NULL, ymax = NULL, error_bars, ref_centered){
oldPar <- par(no.readonly = TRUE)
on.exit(oldPar)
error_bar_width = 0.1
addline_width = 0.5
addline_offset = 0.55
error_bar_lwd = 2
text_offset = 0.1
if(is.null(marker2_vals)){
genotypes <- sort(unique(marker1_vals[which(!is.na(marker1_vals))]))
pheno_vals <- lapply(genotypes, function(x) phenoV[which(marker1_vals == x)])
pheno_means <- sapply(pheno_vals, function(x) mean(x, na.rm = TRUE))
if(error_bars == "sd"){
pheno_error <- sapply(pheno_vals, function(x) sd(x, na.rm = TRUE))
}
if(error_bars == "se"){
pheno_error <- sapply(pheno_vals, function(x) sd(x, na.rm = TRUE)/sqrt(length(x)))
}
if(error_bars == "none"){
pheno_error <- rep(0, length(pheno_vals))
}
# pheno_error[which(is.na(pheno_error))] <- 0
if(is.null(ymin)){ymin <- min(c(pheno_means - pheno_error, 0), na.rm = TRUE)}
if(is.null(ymax)){ymax <- max(pheno_means + pheno_error, na.rm = TRUE)}
plot_height = ymax - ymin
a <- barplot(pheno_means, ylim = c(ymin, ymax*1.1))
segments(a, pheno_means-pheno_error, a, pheno_means+pheno_error, lwd = error_bar_lwd)
abline(h = 0)
#lower bar
segments((a+error_bar_width), (pheno_means-pheno_error), (a-error_bar_width), pheno_means-pheno_error, lwd = error_bar_lwd)
#upper bar
segments(a+error_bar_width, pheno_means+pheno_error, a-error_bar_width, pheno_means+pheno_error, lwd = error_bar_lwd)
par(xpd = TRUE)
text(x = a[,1], y = ymin-(plot_height*0.05), labels = genotypes)
mtext(pheno_name, side = 3, line = 1.5)
} else { #instead if there are two markers
#========================================================================
# internal functions
#========================================================================
plot_bars_int <- function(pheno_vals, pheno_error, pheno_name, marker1, marker2, geno_n){
if(length(pheno_vals) == 4){
base_val = pheno_vals[1]
}else{
base_val = 0
}
pred_add <- base_val + ((pheno_vals[2]-base_val) + (pheno_vals[3]-base_val))
if(is.na(pred_add)){pred_add <- 0}
pred_error <- pheno_error[2] + pheno_error[3]
if(is.null(ymin)){
if(error_bars != "none"){
ymin <- min(c((pheno_vals-pheno_error), pred_add - pred_error, 0), na.rm = TRUE)
}else{
ymin <- min(c(pheno_vals, pred_add, 0), na.rm = TRUE)
}
}
if(is.null(ymax)){
if(error_bars != "none"){
ymax <- max(c((pheno_vals+pheno_error), pred_add + pred_error), na.rm = TRUE)
}else{
ymax <- max(c(pheno_vals, pred_add), na.rm = TRUE)
}
}
full_range = ymax-ymin
geno_n_y <- ymax*1.1
ymax <- ymax*1.1
label_y <- min(pheno_vals)-full_range*0.1
a <- barplot(pheno_vals, ylim = c(ymin, ymax), axes = FALSE,
xlim = c(0, length(pheno_vals)*1.5), cex.names = 1.5, main = pheno_name, names.arg = NA)
axis(2, cex.axis = 1.5)
abline(h = 0)
par(xpd = TRUE)
text_cex = 1.5
text(x = 0, y = ymin-(full_range*0.1), labels = marker1_label, adj = 1, cex = text_cex)
text(x = 0, y = ymin-(full_range*0.22), labels = marker2_label, adj = 1, cex = text_cex)
text(x = c(a), y = ymin-(full_range*0.1), labels = c(min_geno, max_geno, min_geno, max_geno), cex = text_cex)
text(x = c(a), y = ymin-(full_range*0.22), labels = c(min_geno, min_geno, max_geno, max_geno), cex = text_cex)
if(error_bars != "none"){
segments(a, pheno_vals-pheno_error, a, pheno_vals+pheno_error,
lwd = error_bar_lwd)
#lower bar
segments((a+error_bar_width), (pheno_vals-pheno_error),
(a-error_bar_width), pheno_vals-pheno_error, lwd = error_bar_lwd)
#upper bar
segments(a+error_bar_width, pheno_vals+pheno_error, a-error_bar_width,
pheno_vals+pheno_error, lwd = error_bar_lwd)
}
segments(a[length(a)]+addline_width, pred_add, a[length(a)]-addline_width,
pred_add, lty = 2, lwd = 2)
poly_x <- c(a[length(a)]-addline_width, a[length(a)]+addline_width)
poly_y <- c(pred_add-pred_error, pred_add+pred_error)
polygon(x = c(poly_x, rev(poly_x)), y = rep(poly_y, each = 2),
col = rgb(253/256,192/256,134/256, alpha = 0.5))
arrows(a[length(a)]+addline_width+addline_offset, pred_add,
a[length(a)]+addline_offset, pred_add, lty = 1, lwd = 2, length = 0.1)
text(x = a[length(a)]+addline_width+addline_offset+text_offset, y = pred_add,
labels = "Additive\nPrediction", adj = 0)
text(x = a[,1], y = rep(geno_n_y, length(a[,1])), labels = geno_n)
par(xpd = FALSE)
mtext(pheno_name, side = 2, line = 2.5)
}
#========================================================================
#the genotype combinations used are based on the
#maximum genotype present. We want to show the
#maximum effect for the "mutant" phenotype.
min_geno <- min(c(marker1_vals, marker2_vals), na.rm = TRUE)
max_geno <- max(c(marker1_vals, marker2_vals), na.rm = TRUE)
w_w <- intersect(which(marker1_vals == min_geno), which(marker2_vals == min_geno))
w_m <- intersect(which(marker1_vals == max_geno), which(marker2_vals == min_geno))
m_w <- intersect(which(marker1_vals == min_geno), which(marker2_vals == max_geno))
m_m <- intersect(which(marker1_vals == max_geno), which(marker2_vals == max_geno))
pheno_list <- vector(mode = "list", length = 4)
names(pheno_list) <- c(paste(min_geno, min_geno, sep = "/"),
paste(min_geno, max_geno, sep = "/"),
paste(max_geno, min_geno, sep = "/"),
paste(max_geno, max_geno, sep = "/"))
pheno_list[[1]] <- phenoV[w_w]
pheno_list[[2]] <- phenoV[w_m]
pheno_list[[3]] <- phenoV[m_w]
pheno_list[[4]] <- phenoV[m_m]
geno_n <- sapply(pheno_list, length)
pheno_means <- sapply(pheno_list, function(x) mean(x, na.rm = TRUE))
if(error_bars == "se"){
pheno_error <- sapply(pheno_list, function(x) sd(x, na.rm = TRUE)/sqrt(length(x)))
}
if(error_bars == "sd"){
pheno_error <- sapply(pheno_list, function(x) sd(x, na.rm = TRUE))
}
if(error_bars == "none"){
pheno_error <- rep(0, length(pheno_list))
}
pheno_error[which(is.na(pheno_error))] <- 0
if(!ref_centered){
plot_bars_int(pheno_vals = pheno_means, pheno_error = pheno_error,
pheno_name = pheno_name, marker1 = marker1_vals, marker2 = marker2_vals,
geno_n = geno_n)
}else{
#also make a plot using B6 as the reference point
#center on B6 genotype
ref_pheno <- pheno_means-pheno_means[1]
plot_bars_int(pheno_vals = ref_pheno, pheno_error, pheno_name = pheno_name, marker1 = marker1_vals, marker2 = marker2_vals, geno_n = geno_n)
}
}
} #end case for two markers
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.