Using drawProteins with BioMart to draw schematic of human MAP kinases...

Update (5th April, 2018): 'unigene' is not a valid attribute(s) so I have changed the code here.


Last week, I submitted my drawProteins R package to Bioconductor. One of the recommended things to do as a Bioconductor developer is to ensure that the package can work with other Bioconductor packages.

I've created this blog post to show how it's possible to use bioMart to pull out the UniProt accession numbers for the Gene Ontology (GO) term, "MAP kinase activity". This has a GO number of GO:0004707.

My script borrows heavily on the biomaRt users guide written by Steffen Durinck, Wolfgang Huber, Mike Smith. My thanks to them.

Using the R script below, I created this protein schematic:

START
# install bioMart if you haven't used it before
# remove the hash tags...
# source("http://www.bioconductor.org/biocLite.R ")
# biocLite()
# biocLite("biomaRt")
library(biomaRt)

# steps
# 1 choose a Mart with useMart()
# 2 choose a dataset with useMart()
# 3. make query with getBM()
# with
# A. filter  -  restriction on the query
# e.g. Interpro ID(s) [e.g. IPR000001]
# e.g. chromosome_name
# B. attributes  -  values we are interested in to retrieve.
# C. values - values you want to

listMarts()
# from https://bioconductor.org/packages/release/bioc/vignettes/biomaRt/inst/doc/biomaRt.html

# chosing a database = MART and a dataset - gets more focussed each step...
ensembl = useMart("ensembl",
                  dataset="hsapiens_gene_ensembl")

# Retrieve all entrezgene identifiers and HUGO gene symbols of genes which have
# a “MAP kinase activity” GO term associated with it.
# this is the GO:0004707
getBM(attributes = c('entrezgene','hgnc_symbol'),
      filters = 'go',
      values = 'GO:0004707',
      mart = ensembl)

# this is 14 proteins....

# create output in a dataframe and add uniprotswissprot
# which is the UniProt ID
output <- getBM(attributes = c('uniprotswissprot',
                               'hgnc_symbol'),
                filters = 'go',
                values = 'GO:0004707',
                mart = ensembl)


# returns a dataframe... pull out uniprotIDs for drawing...
uniprotIDs <- output$uniprotswissprot
# get rid of blank entries - turn into NA
uniprotIDs[uniprotIDs==""] <- NA
# remove NA
uniprotIDs <- na.omit(uniprotIDs)
# make the IDs characters
uniprotIDs <- as.character(uniprotIDs)
# just the unique ones
uniprotIDs <- unique(uniprotIDs)
# combine into one element
uniprotIDs <- paste(uniprotIDs, collapse = " ")
# this can now be used in drawProteins

# now get features from Uniprot
library(magrittr)

# install drawProteins from Github
# devtools::install_github("brennanpincardiff/drawProteins")
library(drawProteins)

uniprotIDs %>%
  drawProteins::get_features() %>%
  drawProteins::feature_to_dataframe() ->
  prot_data
# data frame with 722 observations

library(ggplot2)
# basic drawing
p <- draw_canvas(prot_data)
p <- draw_chains(p, prot_data)
p <- draw_domains(p, prot_data)
draw_repeat(p, prot_data)

p <- draw_motif(p, prot_data)
p <- draw_phospho(p, prot_data, size = 4)

# background and y-axis
p <- p + theme_bw(base_size = 20) +  # white background and change text size
  theme(panel.grid.minor=element_blank(),
    panel.grid.major=element_blank()) +
  theme(axis.ticks = element_blank(),
    axis.text.y = element_blank()) +
  theme(panel.border = element_blank())

# add titles
p <- p + labs(title = "Schematic of human MAP kinases",
  subtitle = "circles = phosphorylation sites\nsource:Uniprot")


# move legend to top
p <- p + theme(legend.position="top") + labs(fill="")

p
# AND it works!!

END





Resources and citations:

• Using drawProteins for draw NF-kappaB proteins
• More detail here: Drawing NFkappaB protein schematics with ggplot2
• bioMart home page on Bioconductor
• The biomaRt users guide 
• Durinck, Steffen , Paul T. Spellman, Ewan Birney and Wolfgang Huber "Mapping identifiers for the integration of genomic datasets with the R/Bioconductor package biomaRt." (2009) Nature Protocols 4, 1184-1191
• Durinck, Steffen, Yves Moreau, Arek Kasprzyk, Sean Davis, Bart De Moor, Alvis Brazma and Wolfgang Huber "BioMart and Bioconductor: a powerful link between biological databases and microarray data analysis." Bioinformatics 21, 3439-3440 (2005).

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. 

Opening and plotting some flow cytometry data...

Edited: 20 October 2020. Access to these files have changed and Bioconductor has changed too. 

Bioconductor is a very valuable resource for biochemists, biologists and bioinformaticians using R. It contains over 1000 software packages for the opening, analysing and visualising biomedical and biological data. 

Today, I want to show how this can be used to open a file containing some flow cytometry data and how to generate a plot with the data. 

This is the plot, I am going to generate:

It uses the Bioconductor package flowCore. The how-to manual is here and the full reference manual is here. I'm just showing how to open and plot one file which is just a starting point for flow cytometry data. There is a longer example with more plots here. 

A full flow cytometry work flow will allow opening multiple files, normalising and detailed statistics. There is lots of other packages for flow cytometry on Bioconductor.  

I found this interesting: bioinformatics.ca, the home to all things bioinformatics in Canada. They ran a course in 2013 entitled Flow Cytometry Data Analysis using R (2013) and the material is available to view for free. 

Here is the script:

### START 

## ----download_from_Bioconductor-----------------------------------

# install BiocManager

# install.packages("BiocManager")

# BiocManager::install("flowCore")

library("flowCore")

library(ggplot2)

link <- "https://github.com/brennanpincardiff/R4Biochemists201/blob/master/data/cfse_data_20111028_Bay_d7/A01.fcs?raw=true"

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

data <- flowCore::read.FCS("file.fcs", alter.names = TRUE)

#with colours indicating density

colfunc <- colorRampPalette(c("white", "lightblue", "green", "yellow", "red"))

# this colour palette can be changed to your taste 

vals <- as.data.frame(exprs(data))

ggplot(vals, aes(x=FSC.A, y=SSC.A)) +

    ylim(0, 500000) +

    xlim(0,5000000) +

    stat_density2d(geom="tile", aes(fill = ..density..), contour = FALSE) +

    scale_fill_gradientn(colours=colfunc(400)) + # gives the colour plot

    geom_density2d(colour="black", bins=5) # draws the lines inside

# there are lots of small pieces of debris. There is a blue cloud of cells too.  

# exclude debris using the filter package

vals_f <- dplyr::filter(vals, FSC.A>1000000)

# we still have 25,499 events. 

# repeat the plot

ggplot(vals_f, aes(x=FSC.A, y=SSC.A)) +

    ylim(0, 500000) +

    xlim(0,5000000) +

    stat_density2d(geom="tile", aes(fill = ..density..), contour = FALSE) +

    scale_fill_gradientn(colours=colfunc(400)) + # gives the colour plot

    geom_density2d(colour="black", bins=5) # draws the lines inside

=== END  ===