Exploring Phantom ticket prices with Tidy Tuesday...

Not biological data but an interesting Tidy Tuesday data set.

It is about Broadway Shows. Here is a plot that shows that the Phantom of the Opera seems to have become relative cheaper in a more diverse market.

START
# a script for Tidy Tuesday and Tidy Thursday at CaRdiff UseR group 
# 30 April 2020

library(tidyverse)
# pull in the data
grosses <- readr::read_csv('https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2020/2020-04-28/grosses.csv', guess_max = 40000)

# How many unique shows...
length(unique(grosses$show))
# 1122

# inspration from here: https://github.com/teunbrand/tidytuesdayscripts/blob/master/scripts/2020_04_28_Broadway_Grosses.R
# and image from here: https://twitter.com/TeunvandenBrand/status/1255253561535074306/photo/1

# first make basic graph = average ticket prices over time...

ggplot(grosses, aes(week_ending, avg_ticket_price)) +
    geom_point(alpha = 0.25)
# $500 shows depresses the axis a bit...

grosses %>%
    filter(avg_ticket_price > 400) -> expensive
# this show is Springsteen on Broadway at $500!!!
# can exclude by limiting the y-axis

# filter for the Phantom
grosses %>%
    filter(show == "The Phantom of the Opera") %>%
    ggplot(aes(week_ending, avg_ticket_price, colour=pct_capacity)) +
    geom_point() +
    geom_smooth() + 
    labs(x = "",
        y = "Average Ticket Price ($)", 
        title = "Phantom of the Opera",
        subtitle = "Data from Playbill via Tidy Tuesday") 

# why does the show seem very expensive for some weeks rather than others

grosses %>%
    filter(show == "Mamma Mia!") %>%
    ggplot(aes(week_ending, avg_ticket_price, colour=pct_capacity)) +
    geom_point() +
    geom_smooth() + 
    labs(x = "",
        y = "Average Ticket Price ($)", 
        title = "Phantom of the Opera",
        subtitle = "Data from Playbill via Tidy Tuesday") 

# how does Phantom of the Opera compare to the others?
# make one plot first
plot <- ggplot(grosses, aes(week_ending, avg_ticket_price)) +
    geom_point(alpha = 0.25, size = 0.5) + ylim(0,300)

# filter Phantom
grosses %>%
    filter(show == "The Phantom of the Opera") -> p_of_op

# add to the first plot...
plot2 <- plot +
    geom_point(data = p_of_op, aes(week_ending, avg_ticket_price),
        colour = "red") 

# show with titles
plot2 + theme_bw() +
    labs(x = "",
        y = "Average Ticket Price ($)", 
        title = "Phantom of the Opera (red): relatively cleaper in a more diverse market",
        subtitle = "More expensive some weeks!\nData from Playbill via Tidy Tuesday")

END

Resources:

• Official Repo for Tidy Tuesday
• Tidy Tuesday hash tag on twitter for inspiration
• Another Tidy Tuesday about Passwords

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

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. 

Chance of dying for males and females using the UK Life Tables

I am currently answering comments to R for Biochemists 101 - online training that I've created for the Biochemical Society. One of the participants asked about overlaying the different chances of dying for male and females inspired by the blog piece Exploring the chance of dying using the UK Life Tables.

As is often the case in R, there is a quick and simple answer which works ok if the data is very self explanatory. However, it's difficult to draw a legend.

There is also a better way that builds on the strengths of ggplot2 but it requires reformatting the data. I've used the melt() function from the reshape package.

Here is the graph made with the better answer:

The code below also shows the quick way...

START
library(readxl)
library(ggplot2)

# The UK National Life Tables
# these are available through the Office of National Statistics
# as an excel file from this page:
# https://www.ons.gov.uk/peoplepopulationandcommunity/birthsdeathsandmarriages/lifeexpectancies/datasets/nationallifetablesunitedkingdomreferencetables
# the can be downloaded from here:
ons_url <- c("https://www.ons.gov.uk/file?uri=/peoplepopulationandcommunity/birthsdeathsandmarriages/lifeexpectancies/datasets/nationallifetablesunitedkingdomreferencetables/current/nltuk1315reg.xls")
# to guarantee undistrubed access, I've put the file on github
github_url <- "https://raw.githubusercontent.com/brennanpincardiff/RforBiochemists/master/data/nltuk1315reg.xls"

# the download.file() function downloads and saves the file with the name given
download.file(url=ons_url,destfile="file.xls", mode="wb")
# if this doesn't work try replacing the url with the github_url
# then we can open the file and extract the data using the read_excel() function.
data<- read_excel("file.xls", sheet = 5, skip = 6)

# need to remove four bottom rows - these are all blank.
data <- data[1:101,]

colnames(data)

# x = age in years
# first set is males and second set is females
# mx is the central rate of mortality, defined as the number of deaths at age x last birthday in the three year period to which the National Life Table
# relates divided by the average population at that age over the same period.

# qx is the mortality rate between age x and (x +1), that is the probability that a person aged x exact will die before reaching age (x +1).

# to allow us to separate the male and female data, I am relabelling the column names

colnames(data) <- c("age",
                    "m.mx", "m.qx", "m.lx", "m.dx", "m.ex",
                    "f.mx", "f.qx", "f.lx", "f.dx", "f.ex")


# the life table allow us to draw some interesting curves
# the first one I want to explore is my chance of dying this year.

# qx is the key value.
# the chance of my dying before my next birthday.

# the ending point from the previous graph
# from http://rforbiochemists.blogspot.co.uk/2017/03/exploring-chance-of-dying-using-uk-life.html

p <- ggplot(data = data,
       aes(x = age, y = m.qx)) +
  geom_line(colour = "red") +
  geom_point(colour = "red", size = 0.75) +
  ylab("Chance of dying") +
  xlab("Current Age") +
  ggtitle("Chance of Men Dying before next Birthday") +
  scale_y_log10(
    breaks = c(0.0001, 0.001, 0.01, 0.1, 0.5),  # where to put labels
    labels = c("1/10,000", "1/1000", "1/100", "1/10", "1/2")) +  # the labels
  theme_bw()

# show the object
p

# quick way
# add geom_line() and geom_point() with more data
p +
  geom_line(data = data,
                aes(x = age, y = f.qx),
                colour = "blue") +
  geom_point(data = data,
             aes(x = age, y = f.qx),
             colour = "blue", size = 0.75) +
  # and change title
ggtitle("Chance of Dying before next Birthday")

# However, it's really difficult to add a proper legend

# Longer way - reformat the data into 'long' format
# using melt() function from reshape() package

library(reshape2)
# subset the columns we want - just three
data_s <- data.frame(cbind(data$age, data$m.qx, data$f.qx))
# rename columns
colnames(data_s) <- c("age", "male", "female")

# create new object with melt in long format
data_long <- melt(data_s, id.vars = "age")
# rename columns...
colnames(data_long) <- c("age", "gender", "chance_die")

# now draw the object
# using the colour aesthetics
# you get a legend automatically.
ggplot(data = data_long,
       aes(x = age,
           y = chance_die,
           colour = gender)) +
  geom_line() +
  geom_point() +
  ylab("Chance of dying") +
  xlab("Current Age") +
  ggtitle("Chance of Dying before next Birthday") +
  scale_y_log10(
    breaks = c(0.0001, 0.001, 0.01, 0.1, 0.5),
    labels = c("1/10,000", "1/1000", "1/100", "1/10", "1/2")) + 
  theme_bw()

END of SCRIPT

Resources

• For more about melting and overlay, here is an example of some enzymatic data http://rforbiochemists.blogspot.co.uk/2015/06/plotting-enzyme-data-with-ggplot-part-i.html
• An example of oven time courses http://rforbiochemists.blogspot.co.uk/2015/08/temperature-time-course.html
• Plotting numbers of publications http://rforbiochemists.blogspot.co.uk/2015/07/how-many-publications-in-your-area-of.html