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.

Exploring some corona virus time courses...

Updated 12 Mar 2020 - updating names including "United Kingdom" and adding system date to plot.

It's hard to avoid and the instinct is to get stressed and apathetic in alternating cycles.
I found a interesting graph on twitter and decided to explore some of the data.

There is an online interactive dashboard from John Hopkins University and they have shared the data on Github.

So here is a script to pull this data into R and make a few graphs.

There also seems to be two R packages (one by Rami Krispin on CRAN) and one by GuangchuangYu on Github but I haven't explored them yet...

I explored the data to compare the number of COVID 19 cases in the UK with France and Italy. Here is the graph:

Note the y-axis is a log scale not a linear one.


START
# exploring corona virus data...
library(tidyverse)
library(lubridate)
library(ggthemes)
# import the data
url <- "https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_19-covid-Confirmed.csv"

data <- read_csv(url)

# all the cases in the world
# plot cases in the world
data %>%
    select(-`Province/State`, - Lat, -Long) %>%
    pivot_longer(-`Country/Region`, names_to = "date", values_to = "cases")%>%
    group_by(mdy(date)) %>%
    summarise(total_cases = sum(cases)) %>%
    ggplot(aes(x=`mdy(date)`, y=total_cases)) +
    geom_point() +
    scale_y_continuous(trans='log10') +
    labs(x = "Date",
    y = "Number of COVID-19 cases", 
    title = "COVID-19 Cases worldwide",
    subtitle = "source: https://github.com/CSSEGISandData/COVID-19") + 
    theme_economist_white() -> p1

p1

# plot cases in China
data %>%
    select(-`Province/State`, - Lat, -Long) %>%
    filter(`Country/Region` == "China")%>% 
    pivot_longer(-`Country/Region`, names_to = "date", values_to = "cases")%>%
    group_by(mdy(date)) %>%
    summarise(total_cases = sum(cases)) -> data2

p2 <- p1 %+% data2
p2 + ggtitle("COVID-19 cases in China")

# plot cases in Hubei - capital is Wuhan
data %>%
    select(-`Country/Region`, - Lat, -Long) %>%
    filter(`Province/State` == "Hubei") %>% 
    pivot_longer(-`Province/State`, names_to = "date", values_to = "cases")%>%
    group_by(mdy(date)) %>%
    summarise(total_cases = sum(cases)) -> data3

p2 <- p1 %+% data3
p2 + ggtitle("COVID-19 cases Hubei Province")

# plot cases outside of China
data %>%
    select(-`Province/State`, - Lat, -Long) %>%
    # exclude Mainland China
    filter(!`Country/Region` == "China") %>%
    pivot_longer(-`Country/Region`, names_to = "date", values_to = "cases")%>%
    group_by(mdy(date)) %>%
    summarise(total_cases = sum(cases)) -> world_less_china
p_wlc <- p1 %+% world_less_china
p_wlc + ggtitle("COVID-19 in world excluding China")
# fewer cases but increasing...

# let's plot some other countries
# make a function
plot_country <- function(country){
    data %>%
        select(-`Province/State`, - Lat, -Long) %>%
        filter(`Country/Region` == country) %>%
        pivot_longer(-`Country/Region`, names_to = "date", values_to = "cases")%>%
        group_by(mdy(date)) %>%
        summarise(total_cases = sum(cases)) %>%
        ggplot(aes(x=`mdy(date)`, y=total_cases)) +
        geom_point() +
        scale_y_continuous(trans='log10') +
        labs(x = "Date",
            y = "Number of COVID-19 cases", 
            title = paste("Cases in",country),
            subtitle = "source: https://github.com/CSSEGISandData/COVID-19") + 
        theme_economist_white()
    
}


plot_country("Korea, South")

# something interesting happened here...


plot_country("Italy")

plot_country("Iran")

plot_country("United Kingdom")

# compare multiple countries
three_countries <- c("United Kingdom", "France", "Italy")
data %>%
    select(-`Province/State`, - Lat, -Long) %>%
    filter(`Country/Region` %in% three_countries) %>% 
    pivot_longer(-`Country/Region`, names_to = "date", values_to = "cases") %>% 
    filter(cases>0)  %>%
    ggplot(aes(x=mdy(date), y=cases, color = `Country/Region`)) +
    geom_point() +
    scale_y_continuous(trans='log10') +
    labs(x = "Date",
        y = "Number of COVID-19 cases", 
        title = paste("Cases in UK, France & Italy"),
        subtitle = paste("source: https://github.com/CSSEGISandData/COVID-19. Updated",
            Sys.Date())) + 
    theme_economist_white()

END

Some resources:

Analyzing 2019-nCoV (Coronavirus) data — Why Quality & Sense checks are so important in Data Science?

John Hopkins Interactive Dashboard again

Treads of links to science research papers curated by Dr Muge Cevik

More dot plots exploring the Chance of dying for males and females using the UK Life Tables

Key packages: tidyverse, dplyr, lubridate