Seems we'll all need more COVID vaccinations...

On Wednesday, July 20t0 2021, I saw this extended twitter feed by Dr Céline Gounder that reviewed the data and concepts underpinning whether we will need more COVID vaccines. It's a super long thread that is available here through Threadreader too. The thread felt really important and showed lots of data. One of the figures that interested me most is this one from a Lancet paper:

It is a lovely violin plot that shows a decrease in virus neutralisation by variant of the current corona virus (SARS‑CoV‑2). My understanding is that blood samples were taken from various people who had received various vaccines. 

I was surprised when I looked at the paper because there was no figures in it. However, the figures are all in the Supplementary Material which you have to download as a PDF to read. The Supplementary Material also had a heart warming phrase: "All data (anonymised) and full R code to produce all figures and statistical analysis presented in this manuscript are freely-available online on Github: https://github.com/davidlvb/CrickUCLH-Legacy-VOCs-2021-05"

I was very happy to see the code and data available on Github and I forked the data to see if I could reproduce this violin plot. The short answer was yes. There were some minor challenges. I'm missing the text is at the top of the figure. However, I was able to reproduce the main body of the figure within about 20 minutes of engaging with the process. I was pleased.

Here is the graph I made which is a good reproduction of the one in the Twitter feed and Figure 1B in the supplementary data. Below is some simplified code extracted from the original code that is required to make the Figure 1B.

START

# making just Figure 1B - a violin plot

library(sp)

library(tidyverse)

library(khroma)

# download the data from my fork on Github

github_link <- "https://github.com/brennanpincardiff/Crick-UCLH-Legacy-VOCs-2021-05/blob/main/Crick_Legacy_2021-24-05_B-1-617-2_PUBLIC.Rda?raw=true"

load(url(github_link))

### from Lines 21 to 261 of original script

### Subset data for further analysis

studyData <- dtHashed %>% 

    filter(COVID_vaccStatus %in% c(1,2)) # Ignore individuals who are in seroconversion window following dose 1 and dose 2 of vaccine

### Set constants for various functions / plots

strainOrder <- c("Wuhan1", "D614G", "Kent", "SAfrica", "India2")

referenceIC50 <- 2^9.802672

dose2cohort <- studyData %>% filter(COVID_vaccStatus == 2,  sampleOrderInVaccStatus == 1)

########################################################################

#   Panel 1. Vaccine responses per strain following 2nd dose Pfizer    #

########################################################################

relevantData <- dose2cohort %>%

    pivot_longer(cols = ends_with("ic50"), names_to = "strain", values_to = "ic50")

relevantData$strain <- str_replace_all(relevantData$strain, pattern = "_ic50", replacement = "")

relevantData$strain <- fct_relevel(relevantData$strain, strainOrder)

outplot <- ggplot(relevantData, aes(x=strain, y=ic50, color = strain, label = sample_barcode)) + 

    scale_colour_muted() +

    geom_hline(yintercept = referenceIC50,linetype = 3 ) +

    geom_violin(trim=TRUE) + 

    scale_y_continuous(trans='log2', 

                       breaks=c(5, 10, 64, 256, 1024, 5120), 

                       labels=c("[0]", "[<40]", "64", "256", "1024", "[>2560]"),

                       minor_breaks=c(32, 128, 512, 2048)) +

    ylab(bquote('Virus Neutralisation, '~IC[50]~ '')) +

    geom_jitter(shape=20, position=position_jitter(0.2), alpha=0.3) + 

    stat_summary(fun=median, geom = "point", color="black",  shape=5, size=1, stroke=1) + 

    theme_bw(base_family = "Helvetica Neue Thin") +

    theme(legend.position="none")+

    theme(axis.text.x = element_text(angle = 45, hjust = 1, size=12)) + 

    theme(axis.text.y = element_text(size=12))  + 

    theme(

        panel.grid.major.x = element_blank(),

        panel.grid.minor.x = element_blank(),

        axis.title.y=element_text(size=15),

        axis.title.x=element_blank()

    )

outplot

## END of this SCRIPT

Useful resources

• Making a box and whisker plot with some published proteomic data...
• Original paper: "Neutralising antibody activity against SARS-CoV-2 VOCs B.1.617.2 and B.1.351 by BNT162b2 vaccination" by Wall et al https://doi.org/10.1016/S0140-6736(21)01290-3
• David LV Bauer on Github: https://github.com/davidlvb
• David LV Bauer lab home page
• My fork of the data from the Lancet paper on Github which includes this script

An icon plot inspired by bacon and a new book....

Last Saturday, "The Art of Statistics - Learning from Data" by Professor Sir David Spiegelhalter arrived. I ordered it after watching this video entitled "Why statistics should make you suspicious". It is a very interesting book which discusses statistics, data visualisation, data science and the challenges of teaching statistics.

I like to see if I can reproduce figures from papers and books, so I spent most of yesterday trying to reproduce Figure 1.4 with R. I have been partially successful. I am pretty sure I could have made it more quickly with Powerpoint, but then I wouldn't have learned anything about R. Here is the image from the book (I hope it is OK to reproduce this... I will contact to ask permission):

It's an icon plot or an icon array. They can be used to communicate risk.

Here is my attempt to reproduce the icon array using the ggimage and  ggwaffle packages:

It's not perfect but it's the best I can do at the moment. I'm not happy with the separation of the rows of icons. The final images are very large and cause R-Studio some problems in terms of speed of rendering. However, I've learned a lot.

A much quicker and easier method uses the personograph package. The types of icon are very restricted - only male icons :-( Also the random distribution of black icons of the original image is not possible. This is supposed to show the random nature of disease which I quite like. The images are very quick to render and the code is easy to understand.

Here are the images made with personograph:

Figure 1.4

Bacon sandwich example using a pair of icon arrays. Of 100 people who do not eat 

bacon, 6 (black icons) develop bowel cancer in the normal run of events (top panel).

Of 100 people who eat bacon every day of their lives, there is 1 additional (red) case.

I've also worked with the waffle package as shown in the code below.

Here is all the R code:

## START
# installing the packages, first remove the hash tag to run:

# install.packages("personograph") # easiest way to start

# https://github.com/liamgilbey/ggwaffle
# devtools::install_github("liamgilbey/ggwaffle")

# install.packages("waffle", "readr", "ggpubr")

# https://github.com/GuangchuangYu/ggimage
# setRepositories(ind=1:2)
# install.packages("ggimage")


# first choice of package - nice and easy but limited customisation
library(personograph)
# https://github.com/joelkuiper/personograph
# https://cran.r-project.org/web/packages/personograph/index.html

# the data is supplied in a list and is plotted in order. 
data <- list(first=0.06, second=0.94)
personograph(data,  colors=list(first="black", second="#efefef"),
    fig.title = "100 people who do not eat bacon",
    draw.legend = FALSE, dimensions=c(5,20))

data_2 <- list(first=0.06, second=0.01, third=0.93)
personograph(data_2, colors=list(first="black", second="red", third="#efefef"),
    fig.title = "100 people who eat bacon every day",
    draw.legend = FALSE,  dimensions=c(5,20))

# icons all male... no random distribution...

# second package: waffle
library(waffle)
dont_eat_bacon <- c('Cancer' = 6, 'No cancer' = 94)
waffle(dont_eat_bacon, rows = 5, colors = c("#000000", "#efefef"),
    legend_pos = "bottom", title = "100 people who do not eat bacon")

eat_bacon <- c('Cancer' = 6, 'Extra case' = 1, 'No cancer' = 93)
waffle(eat_bacon, rows = 5, colors = c("#000000", "#f90000","#efefef"),
    legend_pos = "bottom", title = "100 people who eat bacon every day")

# nice clear colours but difficult to change symbols
# without installing fonts into system...
# I would rather not have to install fonts...

# Try a third package - ggwaffle

library(readr)
library(ggwaffle)
library(ggpubr)
library(ggimage)
theme_set(theme_pubr())

# in this case, we basically encode every point as a graph. 
# download the data
link <- ("https://raw.githubusercontent.com/brennanpincardiff/RforBiochemists/master/data/ggwaffledata_mf.csv")
bacon_waf <- read_csv(link)

# have a look at data
View(bacon_waf)

# basic plot with geom_waffle() from ggwaffle package
ggplot(bacon_waf, aes(x, y, fill = bacon)) + 
    geom_waffle()

# add icons with geom_icon() from ggimage package
p1 <- ggplot(bacon_waf, aes(x, y, colour = no_bacon)) + 
    geom_icon(aes(image=icon), size = 0.1) +
    scale_color_manual(values=c("black", "grey")) +
    theme_waffle()  +
    theme(legend.position = "none") +
    labs(x = "", y = "",
        title = "100 people who do not eat bacon")

# show the plot
# this is VERY SLOW to draw...
# because it contains 100 icons each of which gets
# downloaded from the internet. 
# I feel sure there is a better way but I don't know it at the moment...
p1


p2 <- ggplot(bacon_waf, aes(x, y, colour = bacon)) + 
    geom_icon(aes(image=icon), size = 0.1) +
    scale_color_manual(values=c("black","red", "grey")) +
    theme_waffle()  +
    theme(legend.position = "none") +
    labs(x = "", y = "",
        title = "100 people who eat bacon every day")

# this is the text at the bottom of the page
text <- paste("Figure 1.4\n",
"Bacon sandwich example using a pair of icon arrays, with randomly",
"scattered icons showing the incremental risk of eating bacon every",
"day. Of 100 people who do not eat bacon, 6 (solid icons) develop",
"bowel cancer in the normal run of events. Of 100 people who eat",
"bacon every day of their lives, there is 1 additional (red) case.",
sep = " ")

# format the text as ggplot object
# with ggparagraph() from the ggpubr package
text_p <- ggparagraph(text = text, size = 12, color = "black")

# arrange the images and text with ggarrange() from the ggpubr package
together <- ggarrange(p1, p2, text_p, 
    ncol = 1, nrow = 3,
    heights = c(1, 1, 0.3))

together
# AGAIN VERY SLOW...!!

# ggsave("together_3.pdf", together)
# works but these are very large PDF file and
# take a lot of time to render

# at the end it seems good to clear memory...

gc()

## END

Some resources:

Reproducing immunization data from Factfulness - another good book from Hans Rosling.  

Spiegelhalter, DJ (2008) "Understanding Uncertainty" Ann Fam Med 6:196-197 doi: 10.1370/afm.848

Galesic, M & Garcia-Retamero, R (2009) "Using Icon Arrays to Communicate Medical Risks: Overcoming Low Numeracy" Health Psychology  2009, Vol. 28, No. 2, 210 –216  

"Infographic-style charts using the R waffle package" by N Saunders