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

A simple web scrape of protein interaction data from Uniprot...

There is a vast amount of biochemical, biological and molecular data available on the internet. Lots of it is free and some of it is well curated. It's useful to be able to 'scrape' data from good quality databases and web pages.

One of my favourite sources of molecular information is the Uniprot database. According to itself, it is "a comprehensive, high-quality and freely accessible resource of protein sequence and functional information".

A little while ago, I showed how to use R to draw a schematic of the receptor of the pro-inflammatory cytokine, TNFR1. This post shows how to scrape the list of interacting proteins from the Uniprot page in order to draw a graph of the number of experiments showing each interaction. 

For scraping the webpage, I use the xml version of the page and the base R function readLines(). There are more complex ways to download the page including packages and functions that will read in xml data. I wanted to start simply in the first instance...

Here is the graph: 

Here is the script:

# START

library(ggplot2)

# extract interaction data from the TNFR1 Uniprot page

# this is the URL for the TNFR1 Uniprot page

url <- c("http://www.uniprot.org/uniprot/P19438.xml")

# scrape a Uniprot page into R

# this readLines() function is a base R function which allows reading webpages

data <- readLines(url)

print(paste(url, "has just been scraped"))

# this readLines() function is a base R function which allows reading webpages

data <- readLines(url)

print(paste(url, "has just been scraped"))

# this function will clean up the data by removing the xml code

# https://stackoverflow.com/questions/17227294/removing-html-tags-from-a-string-in-r/17227415#17227415

cleanXml <- function(x) {

  return(gsub("<.*?>", "", x))

}

# identify where the proteins names are on the scraped data

name.list <- grep("<name>", data)

name <- cleanXml(data[name.list[1]])

# identify where the interaction data lies on the scraped data

int.list <- grep("<comment type=\"interaction", data)

# 5 line numbers supplied in answer

# corresponds to five interactions on the Uniprot page...

# these are the details we want to extract:

features <- c("uniprot_id",    # 3 lines after returned number

              "uniprot_label", # 4 lines after returned number

              "uniprot_exp")   # 7 lines after returned number

# this function will extract the information from the scraped data

extInteractDataFromPage <- function(x){

  uniprot_id <- cleanXml(data[x + 3])

  uniprot_label <- cleanXml(data[x + 4])

  uniprot_exp <- cleanXml(data[x + 7])

  three.val <- c(uniprot_id, uniprot_label, uniprot_exp)

  return(three.val)

}

# use lapply() to apply the function to each interaction

output <- lapply(int.list, extInteractDataFromPage)

# generates a list of 5 proteins that interact with TNFR1

# to graph this using ggplot, we need to turn the list into a data.frame

output.1 <- as.data.frame(t(as.data.frame(output)))

str(output.1)

# but all the values are Factors so need to turn them into characters...

output.1$V1 <- as.character(output.1$V1)

output.1$V2 <- as.character(output.1$V2)

# ...and into numbers. 

output.1$V3 <- as.numeric(as.character(output.1$V3))

str(output.1) 

# looks better with characters (chr) and numbers (num)

# let's draw the graph with ggplot

g <- ggplot(output.1,         # this is the dataframe

            aes(x = V2,       # names of the proteins

                y = V3)) +    # number of experiments showing the interaction

  geom_bar(stat="identity") + # make a bar chart

  theme_bw()

g   # output the graph

# improve the labels

g <- g + xlab("") +

         ylab("Number of experiments") +  

         theme(axis.title.x = element_text(size = 16)) +

         theme(axis.title.y = element_text(size = 16)) +

         theme(axis.text.x = element_text(size = 16)) +

         theme(axis.text.y = element_text(size = 16))

g

# add a graph title and a source for the data

g <- g + ggtitle(paste("Interacting proteins for", name,"\n Source:", url))

g

g + coord_flip()  # turn the bar chart on it's side.