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.

Bar chart of common mental disorders...

My day job in the School of Medicine at Cardiff University involves facilitating learning around various medical conditions including mental health. I like a few statistics so I have been exploring the prevalence of mental health disorders. I found a report about mental health from Our World in Data which shares all the data it uses on Github - making it open source. There is lots of interesting data.
As well as mental health, there is data and reports about cancer and the burden of disease.

Inspired by the mental health report from Our World in Data, I downloaded some data and generated a graph which shows the prevalence of Mental Health Disorders in the UK.

Here is the graph:

Here is the R script that generated the graph and a few other graph along the way.

===  START ===
# looking at some mental health data...
# source: https://ourworldindata.org/mental-health

library(readr)
library(dplyr)
library(tidyr)
library(ggplot2)

# download the data from Github
data <- read_csv("https://raw.githubusercontent.com/owid/owid-datasets/master/datasets/Mental%20health%20prevalence%20(IHME)/Mental%20health%20prevalence%20(IHME).csv")

# pull out data for UK and wrangle using pipes and dplyr
data %>% 
    # filter() by country and year
    filter(Entity == "United Kingdom", Year == 2016) %>%
    # select() prevalence - percentage 3rd to 13th column
    select(3:13) %>%
    # turn from wide format to long for better plotting using gather()
    gather(key = "CMHD", value = "prevalence") -> data1

# now have new object data1

# first bar chart...
ggplot(data1, aes(x = CMHD, y = prevalence)) +
    geom_bar(stat = "identity")

# plot horizontally with coord_flip()
ggplot(data1, aes(x = CMHD, y = prevalence)) +
    geom_bar(stat = "identity") +
    coord_flip()

# remove the text "- both sexes (percent)" gsub() function
data1$CMHD <- gsub(" \\- both sexes \\(percent\\)", "", data1$CMHD)
# the \\ are escape characters for minus and brackets 

# AND

# reorder the categories as factors by size of prevalence
# https://www.reed.edu/data-at-reed/resources/R/reordering_geom_bar.html
data1$CMHD <- factor(data1$CMHD, levels = data1$CMHD[order(data1$prevalence)])

p <- ggplot(data1, aes(x = CMHD, y = prevalence)) +
    geom_bar(stat = "identity") +
    coord_flip()
p

# add some labels and source....
p <- p +
    theme_bw() +
    labs(x = "",
        y = "Prevalence (%)",
        title = "Prevalence of Common Mental Health Disorders in UK (2016)", 
        subtitle = "https://ourworldindata.org/mental-health")
p

# Our World in Data website has the numbers on the plot...
p <- p +
    geom_text(aes(label=round(prevalence, 2)))
p

# Our World website has different coloured bars on the plot...
# by altering fill in the aes() of ggplot
p <- ggplot(data1, aes(x = CMHD, y = prevalence, fill = CMHD)) +
    geom_bar(stat = "identity") +
    coord_flip() +
    theme_bw() +
    labs(x = "",
        y = "Prevalence (%)",
        title = "UK Prevalence of Common Mental Health Disorders (2016)", 
        subtitle = "https://ourworldindata.org/mental-health") +
    geom_text(aes(label=round(prevalence, 1)))
p

#  Which adds a legend... so remove the legend...
p + theme(legend.position="none")

=== END ===

Some resources:

• Exploring diseases in Wales for SQL Relay...
• A bar chart for Athena SWAN awards
• ggplot documentation for geom_bar()
• http://www.sthda.com/english/wiki/ggplot2-barplots-quick-start-guide-r-software-and-data-visualization
• Another source of data: https://www.mentalhealth.org.uk/statistics

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.