Working with multiple Image Files....

I showcased this script, analysing multiple images, during our Image Analysis Workshop today. It was a very interesting morning. Dr Joaquín de Navascués supplied the Drosophila embryo images that we are analysing here.

In my previous two scripts, I have looked at analysing an image step by step. However, what we really want to do is to analyse multiple images and run the same analysis across all of them. One good way to do this in R is to bring all the images into one object - a list.

This allows us to write functions and then use the lapply( ) function to apply it to each element of the list.

The script below downloads 48 images from a folder on Github. These images are turned grey to allow masking. Then the regions of fluorescence are identified and the number of regions found with each staining is determined. We can look at the distribution of cells through the different layers.

The layers detail the staining of cells in a Drosophila embryo with different methods. The include fluorescence markers and nuclear stains.

Then I have extracted some data into a data frame and drawn graphs with ggplot2. Here are the graphs:

The graphs imply that the cell number as determined by the nuclear staining (blue graph) are similar in all the layers. They layers are generated as the fly embryo is imaged horizontally. The numbers of cells detected with Stain number 1 also stays constant while the number of cells detected with Stain number 2 diminish as we transit down the embryo.

Here is an example of the images:

Three of the images from the stack of 48. Left hand panel (green) is stain 1, middle panel is nuclei staining, right hand panel (white) is stain 2.

Here is the script. It's a bit long - 207 lines ! It might be easier to cut and paste from Github.
SCRIPT START
# install.packages("EBImage", "ggplot2", "RCurl")  # if required. 

library(EBImage)
library(ggplot2)
library(RCurl)

# this is a link to the list of files in our folider
link <- "https://raw.githubusercontent.com/brennanpincardiff/RforBiochemists/master/data/images/imageSeq/filenameslist"
# the download.file() function downloads and saves the file with the name given
download.file(url=link, destfile="filenameslist", mode="wb")
load("filenameslist")

# how many files are in this list?
length(filelist)
# 48

# this creates the first function called downloadImages
downloadImages <- function(x){
  # assemble URLs
  url <- paste0("https://raw.githubusercontent.com/brennanpincardiff/RforBiochemists/master/data/images/imageSeq/", x)
  cat("Image downloaded") # so that I know something is happening. 
  img <- readImage(url)
}

# read in all the files into a list using the lapply( ) function
# lapply - applies a function to a list or a vector. 
# in this case, it applies the function readImage() to each element in the 
# object filelist
# lapply( ) syntax: name of list, name of function
images <- lapply(filelist, downloadImages)

# it takes a while to download 48 images - be patient

# we now have all 48 images in a list 
# to show us that R knows this is a list we can use the function mode( )
mode(images) # output to console says "list"

# to access a single image in the list we use double square brackets
display(images[[2]])  # display the second image in the list

# we can display the image brighter
display(images[[2]]*4, method = "raster")

# it might be good to save this locally so that it doesn't need 
# to be downloaded again
# CHECK where you want to save it first
# setwd("SomeWhereYouWantToSaveTheData")
# setwd("/Users/paulbrennan/Desktop")  # example
save(images, file = "Dros_Images")

# this can then be reloaded using the load( ) function
load("Dros_Images")

# Next step to alter the list of images a a group
# make all the images greyscale
# http://stackoverflow.com/questions/20347025/convert-rgb-image-to-one-channel-gray-image-in-r-ebimage
# create another function...
makeGrey <- function(img){
  print("One image grey")
  channel(img, mode = "grey")
}

# then apply the function to the list of images using lapply
# lapply( ) syntax: (name of list, name of function)
images.grey <-lapply(images, makeGrey)

# we have now created a list of grey scale images. 
display(images.grey[[2]]*4, method = "raster)  # allows to check. 

# write a function to create image masks
maskCells <- function(img1){   
  # blur the image
  img1 <- gblur(img1*2, sigma = 5)
  
  # apply a Otsu threshold 
  img1.mask <- img1 > otsu(img1)  
  
  # generate an image with different values for each connected region
  # key here is the bwlabel( ) function
  img1.mask <- bwlabel(img1.mask)
  
  # show this as a coloured blobs 
  display(colorLabels(img1.mask), method = "raster")
  
  
  # the bwlabel() function 'counts' the blobs
  count <- max(bwlabel(img1.mask))
  
  # this outputs the count to us
  cat('Number of cells in this image =', count,'\n')
  
  return(img1.mask)
}

# make the image masks using the function and the lapply
images.mask <- lapply(images.grey, maskCells)
# result is a list of image masks

# write function to extract number of objects masked
countRegions <- function(img1.mask){return(max(bwlabel(img1.mask)))}

# to extract numbers from the image masks use lapply and function again
numbers <- lapply(images.mask, countRegions)
# result is  a list of 48 numbers in this case

# use the numbers to make some plots
# in order to graph the numbers we put them into a data frame
numb <- unlist(numbers)
df <- data.frame(numb)
df$type <- c("memb", "st1", "nuc", "st2")
df$layer <- c(1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4,5,5,5,5,6,6,6,6,7,7,7,7,8,8,8,8,9,9,9,9,
              10,10,10,10,11,11,11,11,12,12,12,12)

df.st1 <- subset(df, df$type == "st1")
df.st2 <- subset(df, df$type == "st2")
df.nuc <- subset(df, df$type == "nuc")

df.2 <- as.data.frame(cbind(df.nuc$layer, df.nuc$numb, df.st1$numb, df.st2$numb))
colnames(df.2) <- c("layer", "nuc", "st1", "st2")


# using the dataframe, I am making plots of the number of stained cells across the various layers
# make the plot for the numbers of nuclei
nucHist <- ggplot(df.2, 
                   aes(x=as.factor(layer), 
                       y=nuc)) +
            geom_bar(stat="identity", fill = "blue") +
            ggtitle("Nuclei per layer") +
            ylab("Number of regions") +
            xlab("Layers from top to bottom") +
            theme_bw()
# show the plot
nucHist

# make the plot for the numbers of regions stained with stain number 1
stain1Hist <- ggplot(df.2, 
                     aes(x=as.factor(layer), 
                         y=st1)) +
              geom_bar(stat="identity", fill = "green") +
              ggtitle("Stain number 1") +
              ylab("Number of regions") +
              xlab("Layers from top to bottom") +
              theme_bw()

# make the plot for the numbers of regions stained with stain number 2
stain2Hist <- ggplot(df.2, 
                     aes(x=as.factor(layer),
                         y=st2)) +
              geom_bar(stat="identity", fill = "grey") + 
              ggtitle("Stain number 2") +
              ylab("Number of regions") +
              xlab("Layers from top to bottom") + 
              theme_bw()



# from: http://www.cookbook-r.com/Graphs/Multiple_graphs_on_one_page_(ggplot2)/
# Multiple plot function
#
# ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects)
# - cols:   Number of columns in layout
# - layout: A matrix specifying the layout. If present, 'cols' is ignored.
#
# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),
# then plot 1 will go in the upper left, 2 will go in the upper right, and
# 3 will go all the way across the bottom.
#
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)
  
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  
  numPlots = length(plots)
  
  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  
  if (numPlots==1) {
    print(plots[[1]])
    
  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    
    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}


# run code above to make multiplot function first
multiplot(stain1Hist, nucHist, stain2Hist, cols = 3)
SCRIPT END

Extracting fluorescence of objects from an Image with R...

In the last blog post, I showed how to identify regions in an image of cells.
This is the image:

It shows about 40 or so cells stained for a protein of interest.

Using the R package, EBImage, it's possible to identify these bright images, count the numbers and then extract and calculate the fluorescence signal in each region.

Here is a graph of the fluorescence signal for each of the 41 regions:

The script that does this is here:

SCRIPT START
# showing how to extract features from an image with EBImage

library(EBImage)
library(ggplot2)

# load up image
# the readImage( ) function will work on URLs to import image into R
c2 <- readImage("https://raw.githubusercontent.com/brennanpincardiff/RforBiochemists/master/data/images/Dros_c2.tif") 
# creates an object called c2

# show a brighter version of the image in the R graphics window
display(c2*4, method = "raster")  #"raster" method means within R

# create the "mask" of this image by blur, threshold and counting
c2.b <- gblur(c2*2, sigma = 5)  #blur
c2.t <- c2.b > otsu(c2.b) # apply otsu threshold 
c2.t.cnt <- bwlabel(c2.t) # count the 'regions'

# show this as a coloured blobs 
display(colorLabels(c2.t.cnt), method = "raster") 
display(colorLabels(c2.t.cnt))  # in a browser

# next step is to extract some of the features about these objects. 
# first extract the basic features (as defined by EBImage). 
# first of these is mean intensity (fluorescence per pixel for the object. )
# this function creates a matrix with each object in a row
# these are calculated using the thresholded image and the original image
ftb <- computeFeatures.basic(c2.t.cnt, c2)

# to get a vector of the mean intensities
m.intent <- ftb[,1]
m.intent # show in Console
# there is one value for each shape...

# next category is shapes
# these a calculated from the thresholded image
fts <- computeFeatures.shape(c2.t.cnt)
area <- fts[,1]  # select first column from object called fts
perimeter <- fts[,2] # select second colum from object called fts

# put the data into a dataframe (a different kind of R object)
df <- as.data.frame(m.intent)
df$area <- area
df$signal <- df$m.intent * df$area
df$regions <- seq(1:nrow(df))

# lots of ways to plot but one of the best is ggplot2
# plot the area of each region on the mask
ggplot(df, 
       aes(x = regions, y=area)) +
  geom_bar(stat="identity")

# plot the signal of each region on the mask
ggplot(df, 
       aes(x = regions, y=signal)) +
  geom_bar(stat="identity")

# plot mean intensity
ggplot(df, 
       aes(x = regions, y = m.intent)) +
  geom_bar(stat="identity")

# dot plot of area v signal
ggplot(df, 
       aes(x=area, y=signal)) +
  geom_point()

# plot the signal of each region on the mask 
# in a slightly nicer way
ggplot(df, 
       aes(x = regions, y=signal)) +
  geom_bar(stat="identity") +
  ggtitle("Fluorescence within each region of the image") + 
  xlab("Regions from top right to bottom left of the image") +
  ylab("Total fluorescence in the region") +
  theme_bw()

# for some help on ggplot2, this R-Cookbook site is good. 
# http://www.cookbook-r.com/Graphs/Bar_and_line_graphs_(ggplot2)/
SCRIPT END

There is lots of extra information about the features that can be extracted here:

• EBImage Reference Manual

Counting and identifying stained cells step by step

On Friday, Dr Joaquín de Navascués and I held the first of two Image Analysis Workshops at Cardiff University. We had a good audience with PhD students and researchers. Joaquin led the first session which explored the fundamentals of images and using Image J or Fiji to open and explore those images.

Key points included:

• good images start at the microscope - try not to over saturate your images during collection.
• protect your original files as it can be easy to alter them permanently and lose data.
• digital pictures are arrays of numbers. 
• these numbers can be visualised and transformed in lots of different ways. 

In just under two weeks time, it's my turn. I am going to talk about using R to analyse images. I like R because it allows the development of reproducible, sharable and scalable workflows. The aim of the second workshop is to show an automated workflow that attendees can adapt to their own work if they want.

The key package that is useful for analysing images in R is EBImage developed by members of the EMBL and the EBI. It's a powerful and useful package.

This script below is one that will probably be used during the workshop. It is about using R to count the numbers of cells in an image. I'm going to go through the process step by step. The steps are as follows:

1. Download the image from Github
2. Make the image a little brighter
3. Blur the image using a Gaussian filter
4. Apply an Otsu threshold the image to convert to a binary image
5. Count the connected regions - i.e. the cells.
6. Display the regions in a multi-coloured output. 

Here is a confocal microscope picture and the image with the cell regions identified:

Here is the script that identifies and counts the regions:

SCRIPT START
# aim of this script is to identify and count cells....
library(EBImage)
library(RCurl)

# this is the link to the image
link <- "https://raw.githubusercontent.com/brennanpincardiff/RforBiochemists/master/data/images/Dros_c2.tif"
# the download.file() function downloads and saves the file with the name given
download.file(url=link, destfile="file.tif", mode="wb")

# import image of stained cells into R
c2 <- readImage("file.tif")
# creates an object called c2

# show the image in the R graphics window
display(c2, method = "raster")  #"raster" method means within R

# display a brighter image
display(c2*4, method = "raster") 
# because an image is numbers we just multiply to make it brighter

# check details by writing the name of the object
c2 # shows some information
# it's a greyscale image 

# make a brighter image by multiplying all the values by 2
c2.b <- c2*2

# gaussian blur
c2.b.blur <- gblur(c2.b, sigma = 5)
display(c2.b.blur, method = "raster")

Can you spot the difference?

## threshold using Otsu'smethod 
otsu(c2.b.blur) # gives a threshold value using Otsu algorithm
# value = 0.1777344
c2.b.blur.thres <- c2.b.blur > otsu(c2.b.blur) # apply this value 
display(c2.b.blur.thres, method = "raster")

# we see that the image is starkly black and white
# all pixels have been turned into either 0 or 1 - a binary image.

# generate an image with different values for each connected region
# key here is the bwlabel( ) function
c2.b.blur.thres.cnt <- bwlabel(c2.b.blur.thres)

# show this as a coloured blobs 
display(colorLabels(c2.b.blur.thres.cnt), method = "raster")

# count by giving us the max value in the bwlabel() function
nucNo <- max(bwlabel(c2.b.blur.thres))
# output this number to the Console
nucNo  # count = 41. 

# do we need to brighten the image (probably not in this case)?
c2.blur.thres <- gblur(c2, sigma = 5) > otsu(gblur(c2, sigma = 5))
display(c2.blur.thres, method = "raster")
max(bwlabel(c2.blur.thres))
# answer is 42 - so not a big difference - good staining. 

# what happens if we don't blur?
# we can calculate the Otsu threshold to the original image
otsu(c2)
c2.thres <- c2 > otsu(c2)
display(colorLabels(bwlabel(c2.thres)), method = "raster")

# doesn't look that different to the eye but...
# if we try counting....
max(bwlabel(c2.thres))
# answer is 1398 - lots of dots causing problems... 

# there are other ways to apply thresholds but that's for later
SCRIPT END


Resources:

• Bioconductor link about EBImage include Introduction & Reference Manual
• Supplementary information from a Nature Methods paper (2013) which using EBImage and R to generate lots of nice data - it's very detailed but interesting. 

Overview of gender breakdown across types of staff

I have spent a large part of the last month grappling with Athena SWAN data for the School of Medicine at Cardiff University. To apply for this equality award, the Athena SWAN Self Assessment Team needed to look at the gender breakdown of students and staff in the School.

Here is a bar chart I generated with the data to give an overview of the School:

Here's is the script I used to generate this data and some other graphs that explore using ggplot to make bar charts.

START
# overview of the School in terms of gender for our Athena SWAN application
library(readxl)
library(ggplot2)
library(reshape2)

# the data is here
link <- "https://raw.githubusercontent.com/brennanpincardiff/AthenaSWANBenchmarkData/master/asStaffOverview_comparedtype_20160505.xlsx"

# the download.file() function downloads and saves the file with the name given
download.file(url=link, destfile="file.xlsx", mode="wb")

# read in the data...
data <- read_excel("file.xlsx")

# calculate Female and Male as a percentage
data$tot <- data$f + data$m
data$Female <- (data$f / data$tot) *100 
data$Male <- (data$m / data$tot) *100

# first a simple plot of the Cardiff School of Medicine Data
# plot the number of females of each type of member of the School
# pull out the Cardiff data
data.c <- data[data$uni == "Cardiff",2:5]

# make the object with the Cardiff subset  
g <- ggplot(data.c,
       aes(x = type,
           y = f)) +
     geom_bar(stat="identity")+
     theme(axis.text.x = element_text(size = 12))

# show the graph
g

# The order of the bars is not what I want.
# take the order from the data file
types <- data.c$type

# apply this order to the graph
# http://stackoverflow.com/questions/5208679/order-bars-in-ggplot2-bar-graph
g + scale_x_discrete(limits = types)

# I want to put female and male on the same graph so...
# pull out and melt the Cardiff data so that have f and m for each type in long format
data.melt <- melt(data[data$uni == "Cardiff",2:5], 
                  id.var = c("type", "year"))

# put the data in an object
g <- ggplot(data.melt, 
       aes(x = type, 
           y = value,
           fill = variable)) + 
      scale_x_discrete(limits = types)

# show with geom_bar
# this is a stacked bar chart showing both female and male numbers
g + geom_bar(stat="identity")

# show female and male in separate bars using position arguement
g + geom_bar(stat="identity", position="dodge")

# have as stacked bar chart but normalised to one. 
g +  geom_bar(stat="identity", position="fill")

# I really want a percentage of each stacked 
# to allow different parts of the School to be compared
# remove the raw numbers
data.s <- data[, -c(3:6)]
data.melt.s <- melt(data.s[data.s$uni == "Cardiff",2:4], 
                  id.var = c("type"))  # id = or id.var are same
colnames(data.melt.s) <- c("type", "Gender", "Percent")

# make a new graph
p <- ggplot(data.melt.s, 
           aes(x = type, 
               y = Percent,
               fill = Gender)) +
      geom_bar(stat="identity") + 
      scale_x_discrete(limits = types) +
      ylab("Percentage of staff") +
      xlab("") +
      theme_bw() +
      theme(axis.text.x = element_text(size = 14))
p # show the graph

# adjust the colours using scale_fill_brewer()
p + scale_fill_brewer(palette = 16, direction=-1) +
  # and add a title
    ggtitle("Gender breakdown of different members of the School of Medicine") +
    theme(axis.title.y = element_text(size = 14 ))

# use data from Imperial & Leeds as comparison
# use the facet_wrap function to make all looks the same. 
# remove Prof & Support as we only have for Cardiff
data.s <- data.s[-5,]
types <- types[-5]

# melt the data into the format for the bar chart. 
data.melt.s <- melt(data.s, 
                    id.var = c("uni", "type"))
colnames(data.melt.s) <- c("uni", "type", "Gender", "Percent")

# this is a final plot with three Medical Schools (it's a bit big). 
p <- ggplot(data.melt.s, 
            aes(x = type, 
                y = Percent,
                group = uni,
                fill = Gender)) +
  geom_bar(stat="identity") +
  facet_wrap(~uni) +    # this function gives three plots
  scale_x_discrete(limits = types) +
  ylab("Percentage of staff") +
  xlab("") +
  theme_bw() +
  scale_fill_brewer(palette = 16, direction=-1) +
  ggtitle("Gender breakdown of three Schools of Medicine") +
  theme(strip.text.x = element_text(size = 14, colour = "black")) +
  theme(axis.text.x = element_text(size = 11))

p # show the plot

END of SCRIPT

I looked up these pages for help:

• http://stackoverflow.com/questions/3290330/facet-label-font-size
• http://www.sthda.com/english/wiki/ggplot2-axis-ticks-a-guide-to-customize-tick-marks-and-labels

As a biochemist, where would you like to work from an Athena SWAN point of view...

As part of my role in the School of Medicine at Cardiff University, I am helping to prepare an application for an Athena SWAN Award. As part of this, I have prepared a list of Athena SWAN Awards across the UK over the last few years and looking at other applications to learn more about the process. The data is on github. I have been using this data to find Medical Schools that have Athena SWAN awards to use as comparison data for our application and to use as a style guide.

Athena SWAN is a quality mark that focuses on equality. As such, departments and Institutions with Athena SWAN awards should be good places to work. Furthermore the better the award, the better the institution as a work place. I decided to analyse the Department names from the point of view of a biochemist.

Here is a graph, I produced (similar to one I produced before and published here):

To do the analysis, I made a user-defined function using function(). The function is called countAwardType(). To use the function, you put in an argument - the name of a department or institution in the brackets: countAwardType("Biochemistry").

There are a two Gold Awards for "Bio" Departments. These are:

• Queen’s University Belfast, School of Biological Sciences
• University of York, Department of Biology

There are three for Chemistry:

• University of York, Department of Chemistry
• Imperial College London, Department of Chemistry
• School of Chemistry, University of Edinburgh

Oxford University has the most awards - see below for Russell group graph.

SCRIPT START
# As a biochemist, where would you like to work for from an Athena SWAN point of view?

library(RCurl)
library(readxl)

# this is the link to the data
link <- "https://raw.githubusercontent.com/brennanpincardiff/AthenaSWANBenchmarkData/master/ASAwards/AthenaSWANAwardList2012_201520160316.xlsx"

# the download.file() function downloads and saves the file with the name given
download.file(url=link, destfile="file.xlsx", mode="wb")

# read in the file with the read_excel() function
awards<- read_excel("file.xlsx", col_names = TRUE)

award.type <- c("Gold", "Silver", "Bronze")

# I have created a function to count the different types of Award - Gold, Silver & Bronze
# and then return a data.frame with specific names
countAwardType <- function(x){
  
  # to test this function  make x <- c("Medicine") and go through it line by line. 
  
  # list of the different types of Athena SWAN awards
  award.type <- c("Gold", "Silver", "Bronze")
  
  # find the pattern supplied in the function as x e.g. "Medicine"
  # this creates a data.frame. 
  x.Awards <- awards[grep(x, awards$org),]
  
  # use split() function which separates the data.frame into a list by award level
  # there is a data.frame within each list.
  # this can be called by the dollar sign - $. 
  x.split <- split(x.Awards, x.Awards$Level)
  
  # count the number of rows in each list 
  gold <- nrow(x.split$Gold)
  silver <- nrow(x.split$Silver)
  bronze <- nrow(x.split$Bronze)
  
  # if there are no entries then a NULL value will be returned 
  # zero is better therefore recode the NULL values
  if(is.null(gold) == TRUE) {gold <- 0}
  if(is.null(silver) == TRUE) {silver <- 0}
  if(is.null(bronze) == TRUE) {bronze <- 0}
  
  # combine these variables into a vector
  aW <- c(gold, silver, bronze)
  
  # combine the vectors into a data.frame
  awardsDepts.x <- data.frame(x, aW, award.type)
  
  # rename the columns
  colnames(awardsDepts.x) <- c("Dept.Type", "Award.Count", "Award.Type")
  
  # return the answer. 
  return(awardsDepts.x)
}

# with this function, I can answer the following questions....

# How many Departments of "Biochemistry" have Athena SWAN awards?
# I can answer this using the grep() function 
# this searches for the word "Biochemistry" in the org column from the awards data.frame
# it creates a data.frame called biochemAwards
biochemAwards <- awards[grep("Biochemistry", awards$org),]
nrow(biochemAwards)  # answer 5

# or
# I can use my function countAwardType()
countAwardType("Biochemistry")

# returns this output:
# Dept.Type Award.Count Award.Type
# 1 Biochemistry           0       Gold
# 2 Biochemistry           1     Silver
# 3 Biochemistry           4     Bronze


# As a biochemist you might be willing to work in any "Bio" department
# these include "Biochemistry" of course. 
# or a "Chem" Department 
# or a "Medic" Department like I do...
# or perhaps a Pharmacology department "Pharm"

# lets make a list. 

biochemistry.options <- c("Biochemistry", "Bio", "Chem", "Medic", "Pharm", "Immun")

# loop through the list - make a data.frame for plotting
awardsDepts <- NULL
for(i in 1:length(biochemistry.options)){
  awardsDepts <- rbind(awardsDepts, countAwardType(biochemistry.options[i])) 
}

# re-order the data.frame
awardsDepts$Award.Type <- factor(awardsDepts$Award.Type, levels = c("Bronze", "Silver", "Gold"))

# let's draw a graph
# make the first version of the plot
p <- ggplot(awardsDepts, 
            aes(x=Dept.Type, 
                y=Award.Count, 
                fill=Award.Type)) + 
            geom_bar(stat="identity") +  # makes the barplot
     scale_fill_manual(values=c("#956C3E", 
                                "#757576", 
                                "#A28D30")) + 
     # colours as per AS website - see below
     xlab("") +  # no need for the "Depts" x-axis
     theme_few() # nice clean theme  

# I want to add a y-axis label
# and increase the size of the text
p <- p + ylab("Number of Awards") +
  theme(axis.title.y = element_text(size = 14 )) + 
  theme(axis.text = element_text(size = 14))

# make the legend a bit bigger and move it
p <- p + theme(legend.position=c(0.9,0.85), # top right
               legend.text=element_text(size = 12), # inc title size
               legend.title=element_text(size = 12)) # labels

# have a look at the bar chart
p + ggtitle("Departments Types with Athena SWAN Awards")

# Don't forget some of the quality Research Institutes: 
# Babraham and Institue of Cancer Research both have Silver Awards. 


# by way of contrast let's look at some other department types:
departments <- c("Bio", "Chem", "Medic", "Engin", "Comp", "Math", "Psych", "Physic")
awardsDepts <- NULL
for(i in 1:length(departments)){
  awardsDepts <- rbind(awardsDepts, countAwardType(departments[i])) 
}
awardsDepts$Award.Type <- factor(awardsDepts$Award.Type, levels = c("Bronze", "Silver", "Gold"))

# push in the new data
d <- p %+% awardsDepts
d + ggtitle("Departments Types with Athena SWAN Awards")

# we can also compare Institutions
russell <- c("University of Birmingham", "University of Bristol", "University of Cambridge", 
             "Cardiff University", "Durham University", "University of Edinburgh",
             "University of Exeter", "University of Glasgow", "Imperial College London",
             "King’s College London", "University of Leeds", "University of Liverpool",
             "London School of Economics and Political Science", "University of Manchester",
             "Newcastle University", "University of Nottingham", "University of Oxford",
             "Queen Mary, University of London", "Queen’s University Belfast", 
             "University of Sheffield", "University of Southampton", "University College London",
             "University of Warwick", "University of York")

awardsDepts <- NULL
for(i in 1:length(russell)){
  awardsDepts <- rbind(awardsDepts, countAwardType(russell[i])) 
}
awardsDepts$Award.Type <- factor(awardsDepts$Award.Type, levels = c("Bronze", "Silver", "Gold"))

# push in the new data
inst <- p %+% awardsDepts
inst + ylim(0,50) +  # extend the y axis
    theme(legend.position=c(0.1,0.85)) +   # move the legend
    theme(axis.text.x=element_text(angle=90, hjust=1, vjust=0.5 )) +                    # vertical names
    ggtitle("Russell group Universities with Athena SWAN Awards")

# finally Who's got Gold?
awards[grep("Gold", awards$Level),]
Type New_Renew Level Year Month                                                       org
# 70  Dept   Renewal  Gold 2015 Apr University of York, Department of Chemistry
# 360 Dept       N/A  Gold 2013 Nov Queen’s University Belfast, School of Psychology
# 361 Dept       N/A  Gold 2013 Nov University of Cambridge, Department of Physics
# 362 Dept       N/A  Gold 2013 Nov University of York, Department of Biology
# 452 Dept       N/A  Gold 2013 Apr Imperial College London, Department of Chemistry
# 539 Dept       New  Gold 2012   Nov Queen’s University Belfast, School of Biological Sciences
# 540 Dept       New  Gold 2012   Nov School of Chemistry, University of Edinburgh


# colours from AthenaSWAN website
# https://athena-swan.medschl.cam.ac.uk/wp-content/uploads/2014/03/Athena-SWAN_style-guide_October-2012.pdf
# Bronze: R149 G108 B62
# Silver: R117 G117 B118
# Gold: R162 G141 B48
# Convert RGB to hex http://www.rgbtohex.net/