#' Shape space in the ggplot style
#'
#' Make pretty plots from geometric morphometric ordinations
#' using the output of \code{\link[geomorph]{gm.prcomp}} (Adams et al. 2020).
#'
#' Beware that right now, the \code{save.as} feature only works with
#' \code{backtransform.examples = FALSE}.
#'
#' @references Geomorph: Software for geometric morphometric analyses.
#' R package version 3.2.1. D.C. Adams and M.L. Collyer and A. Kaliontzopoulou. 2020.
#' (\href{https://cran.r-project.org/package=geomorph}{Link})
#'
#' @source Dave Angelini \email{david.r.angelini@@gmail.com} [aut, cre]
#'
#' @param x An object of class \code{gm.prcomp}.
#' @param axis1 The PC axis to plot on the horizontal.
#' @param axis2 The PC axis to plot on the vertical.
#' @param main.title An optional title for the plot.
#' @param group A factor to group samples.
#' @param group.title An optional title for the group legend.
#' @param convex.hulls A logical factor specifying whether to plot convex hulls around groups.
#' @param label.groups A logical factor specifying whether to label groups.
#' @param labels An optional character vector of point labels.
#' @param include.legend A logical factor specifying whether to include a legend.
#' @param fixed.aspect A logical factor specifying whether to use a fixed aspect ratio.
#' @param xlim Set x scale limits
#' @param ylim Set y scale limits
#' @param color The color of points and hulls in each group.
#' @param pt.size The size of individual points.
#' @param pt.alpha The opacity of individual points.
#' @param hull.alpha The opacity of the convex hulls.
#' @param save.as An optional filename to export the plot.
#' @param height The height of the exported plot.
#' @param width The width of the exported plot.
#' @param scale A scale factor for elements of the exported plot.
#' @param units The units for height and width of the exported plot.
#' @param backtransform.examples A logical factor specifying whether to include example shapes at the extremes of each axis.
#' @param shape.method Method used to visualize shape differences, following \code{geomorph::plotRefToTarget}.
#' @param lm.labels A logical factor specifying whether to label landmarks in the example shapes.
#' @param ref.shape A reference of consensus shape.
#' @param bt.links A matrix with two columns indicating landmarks to connect by lines in the example shapes.
#' @param bt.shape.mag The desired magnification to be used on the example shape, in comparison to the reference.
#' @param bt.inset.size The relative size of the example shape insets.
#' @param bt.margin.factor A factor to control size of the margin where example shapes are plotted.
#' @param bt.legend.position A two-element vector specifying the justification and position of the legend in a back-transform plot.
#' @param ref.pt.size Scale factor for reference configuration points (single value or vector of values)
#' @param target.pt.size Scale factor for target configuration points (single value or vector of values)
#' @param viridis.color.option The color pallette to use from the \code{viridis} package.
#' @param group.label.size Font size for group labels.
#' @param point.label.size Font size for individual point labels.
#'
#' @export
#'
#' @examples
#' library(geomorph)
#' data(plethodon)
#'
#' # GPA
#' Y.gpa <- gpagen(plethodon$land)
#' # Create a gm.prcomp object
#' PCA <- gm.prcomp(Y.gpa$coords)
#'
#' # Morphospace plots
#' shape.space(PCA, group = plethodon$species, group.title = 'species')
#' shape.space(PCA, group = plethodon$species, group.title = 'species', convex.hulls = TRUE)
#' shape.space(PCA, group = plethodon$species, group.title = 'species', convex.hulls = TRUE,
#' axis1 = 1, axis2 = 3)
#' shape.space(PCA, group = plethodon$species, group.title = 'species', convex.hulls = TRUE,
#' viridis.color.option = 'plasma')
#' shape.space(PCA, group = plethodon$species, group.title = 'species', convex.hulls = TRUE,
#' color = c("royalblue","gray40"))
#'
#' # Labeling points
#' shape.space(PCA, group = plethodon$species, group.title = 'species', convex.hulls = TRUE,
#' labels = plethodon$site,
#' group.label.size = 8, point.label.size = 2,
#' color = c("royalblue","gray40"))
#'
#' # Plots with example shapes
#' shape.space(PCA, group = plethodon$species, group.title = 'species', convex.hulls = TRUE,
#' backtransform.examples = TRUE,
#' ref.shape = Y.gpa$consensus)
#'
#' # Custom links on the example shapes
#' pletho.links <- matrix(c(4,5,5,6,6,7,7,8,8,9,9,10,10,11,2,4,12,2,3,5),
#' ncol = 2, byrow = TRUE)
#' shape.space(PCA, group = plethodon$species, group.title = 'species', convex.hulls = TRUE,
#' backtransform.examples = TRUE,
#' ref.shape = Y.gpa$consensus,
#' bt.links = pletho.links)
#'
#' # No links on the example shapes
#' shape.space(PCA, group = plethodon$species, group.title = 'species', convex.hulls = TRUE,
#' backtransform.examples = TRUE,
#' ref.shape = Y.gpa$consensus,
#' bt.links = "none")
#'
#' # Example shapes displaysed as thin-plate spline (TPS) warp grids
#' shape.space(PCA, group = plethodon$species, group.title = 'species', convex.hulls = TRUE,
#' backtransform.examples = TRUE,
#' ref.shape = Y.gpa$consensus,
#' shape.method = "TPS",
#' lm.labels = FALSE)
#'
shape.space <- function (
x,
axis1 = 1,
axis2 = 2,
main.title = NULL,
group = NULL,
group.title = NULL,
convex.hulls = FALSE,
label.groups = TRUE,
labels = NULL,
include.legend = FALSE,
fixed.aspect = FALSE,
xlim = NULL,
ylim = NULL,
color = NULL,
pt.size = 1.5,
pt.alpha = 0.85,
hull.alpha = 0.2,
save.as = NULL,
height = 5,
width = 7,
scale = 1,
units = "in",
backtransform.examples = FALSE,
shape.method = "points",
lm.labels = TRUE,
ref.shape = NULL,
bt.links = NULL,
bt.shape.mag = 1,
bt.inset.size = 0.35,
bt.margin.factor = 2,
bt.legend.position = c(1,1),
ref.pt.size = 1,
target.pt.size = 0.75,
viridis.color.option = "viridis",
group.label.size = 2.5,
point.label.size = 2.5
) {
# Compatibility with umap::umap output
if (any(grepl("umap",class(x)))) {
class(x) <- "prcomp"
x$x <- x$layout
}
# Check that the input is a PCA object
if (!any(grepl("prcomp",class(x)))) {
stop("Error: 'x' must be an object of class 'prcomp' or 'gm.prcomp'.")
} else {
if (backtransform.examples & !any(grepl("gm.prcomp",class(x)))) {
warning("Warning: To include backtransformed example shapes, 'x' must be an object of class 'gm.prcomp'.")
backtransform.examples <- FALSE
}
}
# Load packages
if (!require(geomorph)) { stop("Package missing. First, try running `install.packages('geomorph')`")}
if (!require(ggplot2)) { stop("Package missing. First, try running `install.packages('ggplot2')`")}
if (!require(ggpubr)) { stop("Package missing. First, try running `install.packages('ggpubr')`")}
if (!require(grid)) { stop("Package missing. First, try running `install.packages('grid')`")}
if (!require(ggplotify)) { stop("Package missing. First, try running `install.packages('ggplotify')`")}
# Check that the requested axes are appropriate
if (!(is.numeric(axis1) & (axis1 %% 1 == 0) & (axis1 > 0) & (axis1 < dim(x$x)[2]))) {
warning("Axis 1 is out of bounds. Defaulting to axis1 = 1.")
axis1 <- 1
}
if (!(is.numeric(axis2) & (axis2 %% 1 == 0) & (axis2 > 0) & (axis2 < dim(x$x)[2]))) {
warning("Axis 2 is out of bounds. Defaulting to axis2 = 2.")
axis2 <- 2
}
pcx <- data.frame(
comp1 = x$x[,axis1],
comp2 = x$x[,axis2]
)
# Check the color options
if (is.null(color)) {
if (is.null(group)) {
color <- "black"
group.title <- NULL
convex.hulls <- FALSE
} else {
if(require("viridis")) {
if (!(viridis.color.option %in% c("magma", "A", "inferno", "B", "plasma", "C", "viridis", "D", "cividis", "E"))) {
warning("Warning: Color option is not recognized. See '?viridis::scale_color_viridis' for details. Defaulting to 'viridis'. ")
viridis.color.option <- "viridis"
}
color <- viridis::viridis(length(unique(group)), option = viridis.color.option)
}
else {
color <- sample(grDevices::colors()[grep('gr(a|e)y', grDevices::colors(), invert = TRUE)],
length(unique(group)) )
}
}
} # End if (is.null(color))
if (is.null(group)) {
group.title <- NULL
convex.hulls <- FALSE
}
pc.contributions <- pcvar(x, dimensions = dim(x$x)[2])[c(axis1,axis2)]
# The base PCA plot
base.plot <- ggplot(pcx, aes(x=comp1, y=comp2, color=group, fill = group)) +
theme_minimal() +
geom_point(size = pt.size, alpha = pt.alpha) +
scale_color_manual(values=color, name=group.title) +
# scale_color_viridis(discrete = TRUE, option = viridis.color.option, name=group.title) +
guides(fill = "none") +
labs(x = paste0('PC',axis1,' (',pc.contributions[1],')'),
y = paste0('PC',axis2,' (',pc.contributions[2],')')
)
# Add main title
if (!is.null(main.title)) {
base.plot <- base.plot + ggtitle(main.title)
}
# Fixed aspect?
if (fixed.aspect) {
base.plot <- base.plot + coord_fixed()
}
# Set scale limits
if (!is.null(xlim)) { base.plot <- base.plot + xlim(xlim) }
if (!is.null(ylim)) { base.plot <- base.plot + ylim(ylim) }
# Add convex hulls
if (convex.hulls) {
base.plot <- base.plot +
stat_chull(alpha = hull.alpha, geom = "polygon", show.legend = FALSE) +
scale_fill_manual(values=color)
if (label.groups) {
if (require("ggrepel") & require("magrittr") & require("dplyr")) {
pcx$group <- group
hull.labels <- pcx %>%
group_by(group) %>%
summarise(x = mean(comp1), y = mean(comp2), .groups = "drop" )
base.plot <- base.plot +
geom_text_repel(data=hull.labels,
aes(x=x, y=y, group=group, label=group ),
hjust=0, size = group.label.size, lineheight=0.825, alpha=0.7,
force = 1, max.iter = 10000)
}
} # End if (label.groups)
} # End if (convex.hulls)
if (!is.null(labels)) {
if (length(labels) != (dim(pcx)[1])) {
warning("Warning: Length of `labels` does not match specimen number. ")
} else {
if (require("ggrepel") & require("magrittr") & require("dplyr")) {
pt.labels <- data.frame(
x = pcx$comp1,
y = pcx$comp2,
label = labels
)
base.plot <- base.plot +
geom_text_repel(data=pt.labels,
aes(x=x, y=y, label=labels ),
hjust=0, size = point.label.size, lineheight=0.825, alpha=0.7,
force = 1, max.iter = 10000)
}
}
}
if (!include.legend) {
base.plot <- base.plot + theme(legend.position="none")
}
# Exit scenario: no backtransform & an output file is requested
if (!backtransform.examples & !is.null(save.as)) {
if (!(grepl("\\.pdf",save.as, ignore.case = TRUE) |
grepl("\\.png",save.as, ignore.case = TRUE) |
grepl("\\.jpg",save.as, ignore.case = TRUE) ) ) {
warning("Warning: Output file format is not recognized. Suported formats are pdf, png, and jpg. Defaulting to jpg. ")
save.as <- paste0(save.as,'.jpg')
}
ggsave(filename = save.as, plot = base.plot,
height = height, width = width, scale = scale, units = units)
return(base.plot)
}
# Add backtransform examples
if (backtransform.examples) {
# Check that a reference shape is provided. If not, return the base plot
if (is.null(ref.shape)) {
warning("Error: A reference shape is required to include back-transform examples. Try 'ref.shape = gpa$consensus'.")
return(base.plot)
}
# If no links are provided to connect landmarks, then outline them...
if (is.null(bt.links)) {
bt.links <- grDevices::chull(ref.shape)
bt.links <- matrix(c(bt.links,bt.links[-1],bt.links[1]), ncol=2, byrow = FALSE)
}
# Unless the input is "none", in which case don't show any links
else { if (any(bt.links == "none")) { bt.links <- NULL } }
# ggplotify::as.grob wants variables to global
ref.shape <<- ref.shape
x <<- x
shape.method <<- shape.method
lm.labels <<- lm.labels
bt.links <<- bt.links
bt.shape.mag <<- bt.shape.mag
ref.pt.size <<- ref.pt.size
target.pt.size <<- target.pt.size
axis1.shape.index <<- grep(paste0("comp",axis1,"$"),names(x$shapes))
axis2.shape.index <<- grep(paste0("comp",axis2,"$"),names(x$shapes))
# Create the extreme shape examples
axis1min <- as.grob( ~ plotRefToTarget(
ref.shape, x$shapes[[axis1.shape.index]][["min"]],
method = shape.method, links=bt.links, label = lm.labels, mag = bt.shape.mag,
gridPars = gridPar(pt.size = ref.pt.size, tar.pt.size = target.pt.size) ) )
axis1max <- as.grob( ~ plotRefToTarget(
ref.shape, x$shapes[[axis1.shape.index]][["max"]],
method = shape.method, links=bt.links, label = lm.labels, mag = bt.shape.mag,
gridPars = gridPar(pt.size = ref.pt.size, tar.pt.size = target.pt.size) ) )
axis2min <- as.grob( ~ plotRefToTarget(
ref.shape, x$shapes[[axis2.shape.index]][["min"]],
method = shape.method, links=bt.links, label = lm.labels, mag = bt.shape.mag,
gridPars = gridPar(pt.size = ref.pt.size, tar.pt.size = target.pt.size) ) )
axis2max <- as.grob( ~ plotRefToTarget(
ref.shape, x$shapes[[axis2.shape.index]][["max"]],
method = shape.method, links=bt.links, label = lm.labels, mag = bt.shape.mag,
gridPars = gridPar(pt.size = ref.pt.size, tar.pt.size = target.pt.size) ) )
# Clean-up
# Don't sure why this doesn't work, but right now there variable end up in .GlobalEnv
# rm(ref.shape, x, bt.links, bt.shape.mag)
# rm(axis1.shape.index, axis2.shape.index)
# Prepare to write out the plot
if (!is.null(save.as)) {
if (!(grepl("\\.pdf",save.as, ignore.case = TRUE) |
grepl("\\.png",save.as, ignore.case = TRUE) |
grepl("\\.jpg",save.as, ignore.case = TRUE) ) ) {
warning("Warning: Output file format is not recognized. Suported formats are pdf, png, and jpg. Defaulting to jpg. ")
save.as <- paste0(save.as,'.jpg')
}
if (grepl("\\.pdf",save.as, ignore.case = TRUE)) {
pdf(file = save.as, height = height, width = width)
}
if (grepl("\\.png",save.as, ignore.case = TRUE)) {
png(file = save.as, bg = "transparent", height = height, width = width, units = units)
}
if (grepl("\\.jpg",save.as, ignore.case = TRUE)) {
png(file = save.as, height = height, width = width, units = units)
}
}
# Use 'grid' tools to build the plot
if (bt.margin.factor < 0) { bt.margin.factor <- 1 }
grid.newpage()
grid.draw(
base.plot +
theme(legend.justification=bt.legend.position, legend.position=bt.legend.position) +
xlim(c(min(x$x[,axis1])*bt.margin.factor, max(x$x[,axis1])*bt.margin.factor)) +
ylim(c(min(x$x[,axis2])*bt.margin.factor, max(x$x[,axis2])*bt.margin.factor))
)
pushViewport(viewport(x=.15, y=.5, width=bt.inset.size, height=bt.inset.size))
grid.draw(axis1min)
upViewport()
pushViewport(viewport(x=.85, y=.5, width=bt.inset.size, height=bt.inset.size))
grid.draw(axis1max)
upViewport()
pushViewport(viewport(x=.5, y=.15, width=bt.inset.size, height=bt.inset.size))
grid.draw(axis2min)
upViewport()
pushViewport(viewport(x=.5, y=.85, width=bt.inset.size, height=bt.inset.size))
grid.draw(axis2max)
upViewport()
# Close the output file.
if (!is.null(save.as)) { dev.off() }
} # End back-transform steps
else {
return(base.plot)
}
} # End of function
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.