Last active
December 16, 2019 21:04
-
-
Save zacharystansell/d541ee2879b10472126719aa821e6513 to your computer and use it in GitHub Desktop.
Broccoli Diversity Poster
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
######################## | |
# 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 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment