Data visualization with the programming language R - published.

Today, I published a piece in The Biochemist about data visualisation with R. I enjoyed writing the piece and it includes some of the code and blog posts that I have written here. 

The code I used to create the figures has all been uploaded on Github. To reproduce the figures you can cut and paste the code from there into a script window on R-Studio. 

Here is a list of pages with code for the figures. 

Figure 1 - R allows reproducible data visualisations. 

• Figure 1A - My first published R data visualization – a cluster dendrogram from a proteomic data set. First published in 2014 and reproduced in 2021 (Alsagaby et al, 2014, https://doi.org/10.1021/pr5002803).
• Figure 1B - A violin plot showing neutralizing antibody titres against five SARS-CoV-2 strains reproduced from shared data and code from Wall et al (2021, https://doi.org/10.1016/S0140-6736(21)01290-3).
• Figure 1C - Phylogenetic tree of human kinase domains inspired by visualization by Manning et al (2002, https://doi.org/10.1126/science.1075762).

Figure 2 - Three most viewed data viz from this blog

• Figure 2A - A volcano plot. 
• Figure 2B - A drug dose–response curve and LD50 calculation.
• Figure 2C - A simple flow cytometry plot.

Figure 3 - Tidy Tuesday data visualisations

• Figure 3A - A volcanic activity time line inspired by @ijeamaka_a
• Figure 3B - Illustrating the importance of numbers in password strength. Across a range of password types, inclusion of numbers increases password strength.
• Figure 3C (Github only) - Showing the proportion of female culprits in Scooby Doo shows from 1960s to 2020s.

Figure 4 - Showcasing drawProteins (below)

• Figure 4A - The Hexsticker for the Bioconductor package drawProteins. (b) 
• Figure 4B (Github only) - Schematic of SARS-CoV S1 protein variants. Protein and variant information from Uniprot (https://www.uniprot.org/uniprot/P0DTC2).

There are lots more scripts on this blog. I hope you find it useful for learning R. 

Do reach out if you need any help or if some of the code doesn't work. Either comment here, contact me through Github or twitter: @brennanpcardiff

Programmatic visualization of UK SARS-CoV-2 spike protein variant

Building on my visualization of SARS-CoV-2 spike proteins,  this script provides a R script to allow you to draw a schematic of the corona virus S1 spike protein and the UK variant that has changes within the S1 protein. 

Here is the visualisation and below is the code to make it. 

START

# viz the changes of the UK variant in S1 spike protein....

library(drawProteins)

library(ggplot2)

library(tidyverse)

# download protein data from

# Uniprot link: https://www.uniprot.org/uniprot/P0DTC2

drawProteins::get_features("P0DTC2") -> spike_sars

drawProteins::feature_to_dataframe(spike_sars) -> spike_data

# pull out S1 chain... begins 13 ends: 685

spike_data %>%

    filter(begin > 12 & end < 686) -> s1_bot

# duplicate this and put order = 2

s1_top <- s1_bot

s1_top$order <- 2

# combine these two 

s1_both <- rbind(s1_top, s1_bot)

# draw canvas, chains & regions

draw_canvas(s1_both) -> p

p <- draw_chains(p, s1_both, labels = c("S1 protein", "B.1.1.7 variant"))

p <- draw_regions(p, s1_both)

# here are the details of the changes...

# https://virological.org/t/preliminary-genomic-characterisation-of-an-emergent-sars-cov-2-lineage-in-the-uk-defined-by-a-novel-set-of-spike-mutations/563

uk_variant <- tribble(

    ~type, ~description, ~begin, ~end, ~length, ~accession, ~entryName, ~taxid,

    ~order,

    "B.1.1.7", "deletion", 69, 70, 1, "P0DTC2","SPIKE_SARS2", 2697049, 1,

    "B.1.1.7", "deletion", 144, 144, 1, "P0DTC2","SPIKE_SARS2", 2697049, 1,

    "B.1.1.7", "substitution", 501, 501, 1, "P0DTC2","SPIKE_SARS2", 2697049, 1,

    "B.1.1.7", "substitution", 570, 570, 1, "P0DTC2","SPIKE_SARS2", 2697049, 1,   

    "B.1.1.7", "substitution", 681, 681, 1, "P0DTC2","SPIKE_SARS2", 2697049, 1,

    "B.1.1.7", "substitution", 716, 716, 1, "P0DTC2","SPIKE_SARS2", 2697049, 1,

    "B.1.1.7", "substitution", 982, 982, 1, "P0DTC2","SPIKE_SARS2", 2697049, 1,

    "B.1.1.7", "substitution", 1118, 1118, 1, "P0DTC2","SPIKE_SARS2", 2697049, 1,

)

# overlay information about the variants

p <- p + geom_point(data = filter(uk_variant, begin < 686),

                          aes(x = begin,

                              y = order+0.2, 

                              shape = description), size = 5)

# style the plot a bit...

p <- p + theme_bw(base_size = 14) + # white background

    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()) +

    theme(legend.position = "bottom")

p <- p + labs(title = "Schematic of SARS-CoV-2 S1 Protein and UK variant",

              subtitle = "Source: Uniprot (https://www.uniprot.org/uniprot/P0DTC2)")

p

END

I feel this could, and maybe will, be better but I'm stopping for now :-)

Some Resources

• Plotting disease causing myosin variants using drawProteins
• Using drawProteins to draw MAP kinase proteins
• Using drawProteins step by step guide on Bioconductor
• Reference page for geom_point()

For more help, bug reports or to suggest features

• drawProteins on Github
• Bioconductor forum for questions
• If you use drawProteins in a publication, please cite my paper.

Programmatic visualization of SARS-CoV-2 Spike Protein

The SARS-CoV-2 Spike protein and its variants are very important at the moment. I thought it would be interesting to showcase using my Bioconductor package drawProteins to programmatically draw a visualization of the Spike protein. This helped me understand a little more about the protein too. 

The data for making the visualization is from Uniprot: https://www.uniprot.org/uniprot/P0DTC2

Here is the visualisation and below is the code to make it. 

START

library(drawProteins)

library(ggplot2)

library(tidyverse)

# Uniprot link: https://www.uniprot.org/uniprot/P0DTC2

drawProteins::get_features("P0DTC2") -> spike_sars

drawProteins::feature_to_dataframe(spike_sars) -> spike_data

# From the Uniprot entry, it say that the Spike protein

# is made as a single protein and then processed into 

# S1 and S2 protein. 

# thus the Uniprot entry has multiple chains

# Processing Uniprot data to create different proteins

# pull out full length chain

spike_data %>% 

    filter(begin < 685 & end == 1273) -> spike_data_1

# want this at the top... 

spike_data_1$order = 3

# pull out S1 chain... begins 13 ends: 685

spike_data %>%

    filter(begin < 685 & end < 686) -> spike_data_2

# want this next... 

spike_data_2$order = 2

# pull out S2 chain...  begins 686; ends: 1273

spike_data %>%

    filter(begin > 685 & end < 1274) -> spike_data_3

# want this at the bottom

spike_data_3$order = 1

# combine all back for plotting

spike_data_o <- rbind(spike_data_1, spike_data_2, spike_data_3)

# pull out names for chains

spike_data_o %>%

    filter(type == "CHAIN") -> chains

chain_names <- c(chains$description[1:3], "")

# draw canvas and chains...

draw_canvas(spike_data_o) -> p

p <- draw_chains(p, spike_data_o,

                 labels = chain_names) 

# add regions to S1 and S2

p <- draw_regions(p, spike_data_o)

p <- p + theme_bw(base_size = 14) + # white background

    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()) +

    theme(legend.position = "bottom")

p <- p + labs(title = "Schematic of SARS-CoV-2 Spike Protein",

    subtitle = "Source: Uniprot (https://www.uniprot.org/uniprot/P0DTC2")

p

END

Some Resources

• Plotting disease causing myosin variants using drawProteins
• Using drawProteins to draw MAP kinase proteins
• Using drawProteins step by step guide on Bioconductor
• Reference page for geom_point()

For more help, bug reports or to suggest features

• drawProteins on Github
• Bioconductor forum for questions
• If you use drawProteins in a publication, please cite my paper.

Creating a drawProteins hex sticker...

Lots of Bioconductor packages have created hex stickers - those fun stickers that go on the top of your laptop. They are a bit nerdy but, hey, who doesn't like nerdy 😀

I thought I should create a hex sticker for drawProteins since it's now available in Bioconductor.

It's been created totally in R using the hexSticker package from the talented GuangchuangYu.

Here is the sticker:

Below is the R code that makes the sticker. I've chosen to show three of NF-kappaB proteins but the code could be adapted to illustrate your favourite protein. You can also customise the sticker colour and lots of other things.

START

# script to create drawProteins sticker

# setwd("/Users/paulbrennan/Documents/drawProteins_hexsticker")

library(ggplot2)

# install.packages("hexSticker")

library(hexSticker)

library(drawProteins)

data("five_rel_data")

p <- draw_canvas(five_rel_data)

p <- draw_chains(p, five_rel_data, label_chains = FALSE, size =0.2)

p <- draw_domains(p, five_rel_data,

  show.legend = FALSE,

  label_size = 2)

p

p <- draw_phospho(p, five_rel_data)

p <- p  + theme_void()

sticker(p, package="drawProteins", p_size=6, s_x=.98,

  s_y=.88, s_width=1.4, s_height=1.1,h_fill = "#114466", h_color="#001030")

# this doesn't quite work :-(

# size of text is wrong... no labels... cramped!

# I think five proteins is probably too many...

# reduce to three....

# P19838

# Q04206

# Q04864

# get from UniProt rather than using internal data

prot_data <- drawProteins::get_features("Q04206 Q04864 P19838")

# turn the data into a dataframe

prot_data <- drawProteins::feature_to_dataframe(prot_data)

# create canvas

p <- draw_canvas(prot_data)

# customise the labels

p <- draw_chains(p, prot_data, label_size = 1.4,

  labels = c("c-Rel",

    "p50/p105",

    "p50/p105",

    "p65/Rel A"),

  size=0.2)

# draw domains without a legend or labels

p <- draw_domains(p, prot_data, show.legend = FALSE, label_domains = FALSE)

# add text and control size with geom_text rather than labels from drawProteins

p <- p + ggplot2::geom_text(data = prot_data[prot_data$type == "DOMAIN", ],

  ggplot2::aes(x = begin + (end-begin)/2,

               y = order,

               label = description),

           size = 1.3)

# keep phosphorylation sites small

p <- draw_phospho(p, prot_data, size = 1)

# a completely blank theme using theme_void()

p <- p  + theme_void()

p

# make the sticker using the ggplot object entitled p

# creates a png file

sticker(p,

  package="drawProteins",

  p_size=6,  # font size for package name

  s_x=.95,   # x position for ggplot object

  s_y=.86,   # y position for ggplot object

  s_width=1.4,

  s_height=1.1,

  h_fill = "#e8a935",   # background colour for hexagon

  h_color="#001030",    # border colour

  url = "www.bioconductor.org")

# automatically creates the PNG file with the name drawProteins.

# make the PDF

sticker(p,

  package="drawProteins",

  p_size=6,  # font size for package name

  s_x=.95,   # x position for ggplot object

  s_y=.86,   # y position for ggplot object

  s_width=1.4,

  s_height=1.1,

  h_fill = "#e8a935",   # background colour for hexagon

  h_color="#001030",    # border colour

  url = "www.bioconductor.org",

  filename="drawProteins.pdf")

END

Some Resources:

• Using drawProteins to draw MAP kinase proteins
• Using drawProteins step by step guide on Bioconductor

For more help, bug reports or to suggest features

• drawProteins on Github

Plotting disease causing myosin variants using drawProteins

My drawProteins package has been accepted by Bioconductor. It's currently available as a development version -  a nice conclusion of lots of effort.

Last week, I had my first query about drawProteins from Kevin T. Booth, a graduate student at the University of Iowa. Kevin wanted to plot disease causing protein variants.

Kevin's working on the MYO7A gene (Uniprot ID “Q13402”). Kevin told me that "alterations in MYO7A can result in 3 different phenotypes-Autosomal Dominant non-syndromic hearing loss (DFNA11),  Autosomal Recessive non-syndromic hearing loss (DFNB2) or Usher Syndrome Type 1 (USH1B), another type of inherited deafness.

Looking at the Uniprot page for MYO7A gene, it contains a long list of disease associated variants. We created a visualisation of these variants using the geom_point() function to apply a ggplot2 layer to a schematic of the protein.

Here is the script that generated this and a few other images on the way. Feedback welcome.

SCRIPT START
# exploring plotting variants
# using data from Q13402

# download the protein data first...

library(drawProteins)
library(ggplot2)
# pulling down data for Kevin's protein - MYO7A
prot_data <- drawProteins::get_features("Q13402")
# produce data frame
prot_data <- drawProteins::feature_to_dataframe(prot_data)

# make protein schematic
p <- draw_canvas(prot_data)
p <- draw_chains(p, prot_data)
p <- draw_domains(p, prot_data)
p

# nice domain structure
# note that it is a long protein > 2000 amino acids. 

p <- p + theme_bw(base_size = 14) + # white background
    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())
p

# EXPLORE VARIANTS
# many VARIANTS
# count them...
nrow(prot_data[prot_data$type == "VARIANT",])
# ans = 74

# these can be plotted using geom_point() function
p + geom_point(data = prot_data[prot_data$type == "VARIANT",],
    aes(x = begin,
        y = order+0.25))

# wanted focus on disease variants and colour the by disease...
# so need to change the info...
# colour or shape by info within it...

# duplicate prot_data
prot_data2 <- prot_data
# duplicate description..
prot_data2$disc2 <- prot_data$description
# remove word "in"
prot_data2$disc2 <- gsub("in ", "", prot_data2$disc2)

# pull out the disease name...
library(stringr)
temptable <- str_split_fixed(prot_data2$disc2, ";|:", 2)
prot_data2 <- cbind(prot_data2, temptable)
# name columns
tempnames <- colnames(prot_data2)
tempnames[11] <- "disease"
colnames(prot_data2) <- tempnames

# add colour argument within aesthetics
p + geom_point(data = prot_data2[prot_data2$type == "VARIANT",],
    aes(x = begin,
        y = order+0.25,
        colour = disease))

# exclude some of the VARIANTS as not disease causing..
# dbSNPs, rare polymorphism and the blank one 
prot_data2$disease[prot_data2$disease=="dbSNP"] <- NA
prot_data2$disease[prot_data2$disease==""] <- NA
prot_data2$disease[prot_data2$disease=="rare polymorphism"] <- NA

# clean up variant list for easier legend
disease_var_data <- prot_data2[prot_data2$type == "VARIANT",]
disease_var_data <- na.omit(disease_var_data)
disease_var_data$disease <- gsub("found patients with ", 
    "", disease_var_data$disease)
disease_var_data$disease <- gsub("found a patient with ", 
    "", disease_var_data$disease)

# add colour argument within aesthetics
p + geom_point(data = disease_var_data,
    aes(x = begin,
        y = order+0.25,
        colour = disease))

# over plotting is an problem 
# enlarge canvas and explore geom_jitter()
p <- p + ylim(0.75, 3)
p <- p + geom_jitter(data = disease_var_data,
            aes(x = begin,
                y = 2,
                colour = disease),
            width = 0, height = 0.75, 
            size = 2)

# add titles
myo7A_subtitle <- paste0("circles = Disease Variants\n",
    "source: Uniprot accessed 8 Feb 2018")

p <- p + labs(title = "Schematic of human Unconventional myosin-VIIa",
    subtitle = myo7A_subtitle)
p

SCRIPT END

Some Resources

• Using drawProteins to draw MAP kinase proteins
• Using drawProteins step by step guide on Bioconductor
• Reference page for geom_point()

For more help, bug reports or to suggest features

• drawProteins on Github
• Bioconductor forum for questions 
• You can also email me like Kevin did - thanks Kevin. 

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).

Using drawProteins for draw NF-kappaB proteins...

Update, 4 Jan 2018: Following advice from Bioconductor, the package magrittr is not integral to drawProteins so needs to be activated separately to use this script. 

Update,  3 Nov 2017: I've renamed the functions e.g. geom_chains is now draw_chains and added a new function draw_canvas which allows one to draw without the chain, if you want.

Building on what I did previously, more work done on the drawProteins package allows the simplier script below to draw this protein schematic:

SCRIPT START
# devtools::install_github("brennanpincardiff/drawProteins")
library(drawProteins)
library(ggplot2)
library(magrittr)

# accession numbers of five NF-kappaB proteins
"Q04206 Q01201 Q04864 P19838 Q00653" %>%
  drawProteins::get_features() %>%
  drawProteins::feature_to_dataframe() ->
  prot_data

p <- draw_canvas(prot_data)
p <- draw_chains(p, prot_data)
p <- draw_domains(p, prot_data)
p <- draw_repeat(p, prot_data)
p <- draw_motif(p, prot_data)
p <- draw_phospho(p, prot_data, size = 8)

# 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(x = "Amino acid number",         # label x-axis
              y = "",  # label y-axis
              title = "Schematic of human NF-kappaB proteins",
              subtitle = "circles = phosphorylation sites\nRHD = Rel Homology Domain\nsource:Uniprot")


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

p
SCRIPT END

Resources

• Github repository for the drawProteins package