Nothing
#' 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
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.