Broccoli Diversity Poster

broccoli_poster.R

########################
# Script to produce images with population structure of 97 broccoli accessions. 
# Data generated by Björkman Lab @ Cornell University
# Data collected by Zachary Stansell, Thomas Björkman, Colden Proe, and Teagan Zingg 
# Script by Zachary Stansell www.zacharystansell.com 
# With help from contributors (see additional functions below)
# Funded by the Eastern Broccoli Project (https://blogs.cornell.edu/easternbroccoliproject/)
# https://www.reddit.com/r/dataisbeautiful/comments/e9ymhg/i_grew_97_different_types_of_broccoli_this_summer/
########################

########################
# Packages
########################
library('png')

########################
# Read Data
########################
images <- read.csv("images.csv")
Q <- read.csv("structure.csv")
images <- merge(images, Q, by="geno", all.x=TRUE) 
images <- images[order(images$groupQ),]
clusters <- grep("Cluster",colnames(images))

########################
# Define getColors for popstructure plots
# Full credit to Roy Francis in his wonderful package popHelper (http://www.royfrancis.com/pophelper/articles/index.html)
# Modification and other mangling by me
########################
getColors <- getColours <- function(k) {
  if(length(k) > 1) stop("getColours: Input has more than one value. Argument k must be a single numeric or integer.")
  if(!is.integer(k) && !is.numeric(k) ) stop("getColours: Input is not an integer. Argument k must be a single numeric or integer.")
  k <- as.integer(k)
  # standard colours
  col1 <- c("#2121D9","#9999FF","#DF0101","#04B404","#FFFB23","#FF9326","#A945FF","#0089B2","#B26314","#610B5E","#FE2E9A","#BFF217")
  if(k <= 12) return(col1[1:k])
  if(k > 12) 
  {
    cr <- colorRampPalette(colors=c("#000040FF","#00004FFF","#000060FF","#000074FF","#000088FF","#00009DFF","#0000B2FF",
                                    "#0000C6FF","#000CD8FF","#0022E7FF","#0037F3FF","#004BFBFF","#005EFFFF","#0070FEFF",
                                    "#0081F8FF","#0091EEFF","#00A0E0FF","#00ADCFFF","#00BABCFF","#00C6A7FF","#01D092FF",
                                    "#02DA7EFF","#03E26AFF","#07E958FF","#0EF047FF","#1BF539FF","#31F92CFF","#54FC22FF",
                                    "#80FE1AFF","#ABFF13FF","#CEFF0EFF","#E4FE0AFF","#F1FB07FF","#F8F805FF","#FCF403FF",
                                    "#FDEE02FF","#FEE801FF","#FFE001FF","#FFD801FF","#FFCE00FF","#FFC300FF","#FFB800FF",
                                    "#FFAB00FF","#FF9D00FF","#FF8E00FF","#FF7E00FF","#FF6D00FF","#FF5B00FF","#FF4700FF",
                                    "#FF3300FF"),space="rgb")
    return(cr(k))
  }
}


########################
# Define addImg function. 
# Full credit to Marc Taylor (https://rdrr.io/github/marchtaylor/sinkr/src/R/addImg.R). 
# Very slight modification to his wonderful function; eg switching hortizontal and vertical axes
########################
addImg <- function(
  obj, # an image file imported as an array (e.g. png::readPNG, jpeg::readJPEG)
  x = NULL, # mid x coordinate for image
  y = NULL, # mid y coordinate for image
  HEIth = NULL, # HEIth of image (in x coordinate units)
  interpolate = TRUE # (passed to graphics::rasterImage) A logical vector (or scalar) indicating whether to apply linear interpolation to the image when drawing. 
){
  if (is.null(x) |
      is.null(y) |
      is.null(HEIth)) {
    stop("Must provide args 'x', 'y', and 'HEIth'")
  }
  USR <-
    par()$usr # A vector of the form c(x1, x2, y1, y2) giving the extremes of the user coordinates of the plotting region
  PIN <-
    par()$pin # The current plot dimensions, (HEIth, WEIght), in inches
  DIM <- dim(obj) # number of x-y pixels for the image
  ARp <- DIM[2] / DIM[1] # pixel aspect ratio (y/x)
  HEIi <-
    HEIth / (USR[4] - USR[3]) * PIN[1] # convert HEIth units to inches
  WEIi <- HEIi * ARp # WEIght in inches
  WEIu <- WEIi / PIN[1] * (USR[2] - USR[1]) # WEIght in units
  rasterImage(
    image = obj,
    ybottom = x - (HEIth / 2),
    ytop = x + (HEIth / 2),
    xleft = y - (WEIu / 2),
    xright = y + (WEIu / 2),
    interpolate = interpolate
  )
}

########################
# Corner text function 
# slightly modified from Benjamin Borgy (https://benborgy.wordpress.com/2014/08/01/add-text-in-the-corner-of-a-plot/)
########################
cornerText <- function(text, location = "topright", ...) {
  legend(location,
         legend = text,
         bty = "n",
         pch = NA,
         ...)
}

########################
# Set plotting parameters
########################
# Number of images
imgNo <- length(images$geno)

# set plot dimensions in inches ~ eg 4" x 4" cell
cellSize <- 4

# Automatically set grid size (eg. 10 x 10 for 97 images)
gridSize <- ceiling(sqrt(imgNo))

# Or manually define grid 
#gridSize <- 9

# Plotting dimisions and res
plotDim <- cellSize * gridSize
resolution <- 300


########################
# Plotting
########################
# Start plot;
START <- Sys.time()

png(
  "broccoli_fig_300.png",
  width = plotDim,
  height = plotDim,
  units = "in",
  res = resolution
)

par(mfrow = c(gridSize, gridSize), mar = c(0, 0, 0, 0))
for (i in 1:ifelse(imgNo == gridSize^2, gridSize^2,imgNo)) {
  plot(
    1:1,
    ty = "n",
    xlab = NULL,
    ylab = NULL,
    xaxt = 'n',
    yaxt = 'n',
    ann = FALSE
  )
  image <- as.character(images$filename[i])
  image <- readPNG(image)
  addImg(image, x = 1.05, y = 1.0, HEIth = 0.75)
  
  # Add name of accession
  cornerText(
    text = as.character(images$short.name[i]),
    location = "bottom",
    cex = 1.45,
    text.font = 2
  )
  
  # Add population structure barplots
  subplot(
    barplot(matrix(images[i, clusters]), col = getColors(8), yaxt = "n"),
    x = 1,
    y = 1,
    size=c(0.25,3.25),
    hadj=-7,
    type = "plt"
  )
}

dev.off()

END <- Sys.time()
END - START
#Time difference of 3.975424 mins