Using the Uniprot API to illustrate protein domains...

Update (2nd Nov 2017): I've updated the drawProteins package and it's  totally using ggplot so this script doesn't work. For more details see this blog post, please.

I'm a member of the Cardiff R User group and we're using R to explore APIs - application programming interfaces which allow R to access data. I chose the explore the API that allows me to download data from Uniprot, "a comprehensive, high-quality and freely accessible resource of protein sequence and functional information".

Tomorrow, we have our show and tell about what we have learned about using APIs with R. So in preparation and inspired by Stef Locke, the lead organiser for the Cardiff R User group, I have now created my first R package. The R package, called drawProteins, uses the API output to create a schematic protein domains.

Here are two examples of the output. The first is a schematic of a single protein using Base R plotting (script below). The second example is a schematic of three proteins generated using ggplot2 (script in a little while).

Here is the script to plot a schematic of the transcription factor p65 (top panel):

SCRIPT START
# OK so I have to write some kind of a script to use the new 
# Uniprot API
# This is R-UserGroup homework for Thursday....

# references: https://www.ebi.ac.uk/proteins/api/doc/#!/proteins/search
# references: https://academic.oup.com/nar/article-lookup/doi/10.1093/nar/gkx237

# needs the package httr
# install.packages("httr")
library(httr)

# my package with functions for extracting and graphing
# might require: install.packages("devtools")
devtools::install_github("brennanpincardiff/drawProteins")
library(drawProteins)

# so I want to use the Protein API 
uniprot_acc <- c("Q04206")  # change this for your fav protein

# Get UniProt entry by accession 
acc_uniprot_url <- c("https://www.ebi.ac.uk/proteins/api/proteins?accession=")

comb_acc_api <- paste0(acc_uniprot_url, uniprot_acc)

# basic function is GET() which accesses the API
# requires internet access
protein <- GET(comb_acc_api,
           accept_json())

status_code(protein)  # returns a 200 means it worked

# use content() function from httr to give us a list 
protein_json <- content(protein) # gives a Large list
# with 14 primary parts and lots of bits inside

# function from my package to extract names of protein
names <- drawProteins::extract_names(protein_json)

# I like a visualistion so I want to draw the features of this molecule
# https://www.ebi.ac.uk/proteins/api/features?offset=0&size=100&accession=Q04206

# Get features 
feat_api_url <- c("https://www.ebi.ac.uk/proteins/api/features?offset=0&size=100&accession=")

comb_acc_api <- paste0(feat_api_url, uniprot_acc)

# basic function is GET() which accesses the API
prot_feat <- GET(comb_acc_api,
           accept_json())

status_code(prot_feat)  # returns a 200. 

# so to process this into a data.frame
prot_feat %>%
  content() %>%
  flatten() %>%
  drawProteins::extract_feat_acc() -> 
  features

# clean up for plotting
# focus on those that we want to plot...
features_plot <- features[features$length > 0, ]  
features_plot <- features_plot[complete.cases(features_plot),]
features_plot <- features_plot[features_plot$type!= "CONFLICT",]

# add colours 
library(randomcoloR)   # might need install
colours <- distinctColorPalette(nrow(features_plot)-1)
features_plot$col <- c("white", colours)

# now draw this...  with BaseR
# here is a function that does that....
drawProteins::draw_mol_horiz(names, features_plot)

# add phosphorylation sites
# get sites first. 
features %>%
  drawProteins::phospho_site_info() %>%
  drawProteins::draw_phosphosites(5, "yellow")  # 5 circle radius
text(250,0, "Yellow circles = phosphorylation sites", cex=0.8)

SCRIPT END

Resources

• Stef's 3 lines of code to make an R package: https://itsalocke.com/the-making-of-datasaurus/
• Another useful blog about writing a package: 
  http://tinyheero.github.io/jekyll/update/2015/07/26/making-your-first-R-package.html
• Key reference: Hadley Wickham's book on Packages
• Another R script exploring proteins: Drawing a simple phylogenetic tree of the human rel homology domain family

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

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/

Counting cell nuclei in an image

I have been working with a colleague, Dr Joaquin de Navascues from the European Cancer Stem Cell Research Institute, to develop a workshop entitled an "Introduction to Biological Image Analysis". We plan to discuss the fundamentals of digital images and how to work with them in FIJI (Image J) and R, two open source data analysis tools. We aim to deliver this workshop during March or April next year at Cardiff University.

Joaquin is taking the lead on working with FIJI and I am developing the R material. For biological image analysis there is a useful package called EBImage. This has a nice introduction available and a detailed handbook - (version: 4.13.5).

I have been using this package to count cell nuclei in an image. The image is a microscope picture of a Drosophila gut. Counting nuclei involves mathematical transformations of the digital data. A digital image is a matrix of numbers.

In non-technical language the key steps are:

• blur the image 
• apply a threshold to turn nuclei into 'blobs'
• count the 'blobs'

The output from this script is:

Number of nuclei in this image = 92

The script below downloads an image from Github, opens the image, displays it, transforms it and then counts the nuclei. Because I plan to count nuclei from more than one image, I have made a function and then applied it to the downloaded file. Using user defined functions to automate your workflow is a very good use of R. 

SCRIPT:

# to install use this:
# source("http://bioconductor.org/biocLite.R")
# biocLite("EBImage")

library(EBImage)  # you might need to install - see above

# the image is on Github 
# it is from a set of cells that are stained to detect the nuclei

# this is the link to the data

link <- "https://raw.githubusercontent.com/brennanpincardiff/RforBiochemists/master/data/seq/seq_z015_c003.tif"

# the download.file() function downloads and saves the file

download.file(url=link, destfile="file.tif", mode="wb")

# EBImage uses the readImage() function to load the file. 

img1 <- readImage("file.tif")

display(img1, method = "raster")  # shows the image within R. 

display(img1*4, method = "raster") # multiply the image to make brighter 

# I have written a function to count nuclei

# includes blurring the image, applying a threshold and counting....

# it displays the image as it is changed. 

# it's not perfect and overestimates the number of nuclei. 

# it's an example that can be done. 

# improving and customizing the various options is very feasible. 

countNuclei <- function(img1){   

  # blur the image

  w = makeBrush(size = 11, shape = 'gaussian', sigma = 5)  # makes the blurring brush

  img_flo = filter2(img1*2, w) # apply the blurring filter

  display(img_flo * 4, method = "raster") # display the blurred image - brighter for display only. 

  

  # apply a threshold 

  nmaskt = thresh(img_flo *2, w=10, h=10, offset=0.05) 

  display(nmaskt, method = "raster")

  

  # the bwlabel() function 'counts' the blobs

  nucNo <- max(bwlabel(nmaskt))

  # this outputs the count to us

  cat('Number of nuclei in this image =', max(bwlabel(nmaskt)),'\n')

  return(nucNo)

}

# this applies the function to the image

nucNo <- countNuclei(img1)

END OF SCRIPT