Monday, 1 July 2019

Exploring data from Nature journal about fibroblast subsets in arthritis...

Recently, I attended a conference entitled Data-driven Systems Medicine at Cardiff University. An interesting meeting organised by Dr Barbara Szomolay & Dr Tom Connor who work in the Systems Immunity Institute.

Prof. Mark Coles, Professor at Kennedy Institute of Rheumatology at Oxford University talked about a recent paper published in Nature. The title of the paper was "Distinct fibroblast subsets drive inflammation and damage in arthritis" by Croft et al . A Biorxiv of the manuscript is also available. Nature requires authors to make a "Data availability statement" which encourages authors to share their data.

This gave me the opportunity to explore some of the data in R. It was relatively straight forward workflow: import, tidy and visualise. To make a figure exactly like the paper takes a bit more effort. 

To quote from the paper:
"Rheumatoid arthritis (RA) is a prototypic IMID6 in which synovial fibroblasts contribute to both joint damage and inflammation. We found that expression of FAPα, a cell-membrane dipeptidyl peptidase10, was significantly higher in both synovial tissue and cultured synovial fibroblasts isolated from patients who fulfilled classification criteria for RA compared to patients in whom joint inflammation resolved (Fig. 1a–c), suggesting that FAPα expression may associate with a pathogenic fibroblast phenotype."

Here is a version of Figure 1b which shows FAPalpha in patient samples with different types of disease. It a quantification of staining of cells in synovial tissue.

The code to make this is below as is the code for a version of Figure 1c and 1j. 


### START of script to explore data in Figure 1b, 1c and 1j
# the packages that we need
library(readxl)
library(dplyr)
library(ggplot2)
library(tidyr)
library(ggpubr)

# It's possible to explore the data  
# This is an Excel file for Figure 1 - each part a different worksheet
link <- "https://static-content.springer.com/esm/art%3A10.1038%2Fs41586-019-1263-7/MediaObjects/41586_2019_1263_MOESM8_ESM.xlsx"

# the download.file() function downloads and saves the file with the name given
download.file(url=link, destfile="file.xlsx", mode="wb")

# check worksheet names 
sheet_names <- excel_sheets("file.xlsx")
# import the data for Figure 1b
data <- read_excel("file.xlsx", sheet_names[1])

# data in wide format...
# put into long tidy format... with gather() from dplyr
data %>% gather(na.rm = TRUE) %>%
# make patient samples factors - important for graphing and stats
 mutate(key = factor(key, c("Control", "Resolving", "RA"))) -> data_3

# First plot - a type of dot plot
ggplot(data_3, aes(key, value)) +
    geom_point()

# geom_boxplot shows the statistics - mean and quartiles
# geom_jitter() shows points off centre across the x-axis...
plot <- ggplot(data_3, aes(key, value)) +
    geom_boxplot(outlier.size=0) +
    geom_jitter(width = 0.1)
plot
# we want to put a p-value on the graph
# using the stat_compare_means() function from ggpubr package 
my_comparisons <- list( c("Control", "RA"))
plot + ylim(0,0.2) + 
    stat_compare_means(data = data_3, comparisons = my_comparisons)
# We want to replace the number with four stars as per the paper...
plot <- plot + ylim(0,0.2) + 
    stat_compare_means(data = data_3, comparisons = my_comparisons,
        symnum.args = list(cutpoints = c(0,0.0005, 1), 
            symbols = c("****","ns")))
# and add labels
plot + labs(x = "",
    y = "FAPalpha expression (pixels per unit area)", 
    title = "Figure 1b",
    subtitle = "Croft et al, Nature, 2019") +
    theme_classic() + theme(legend.position="none")
# not the exact same but very similar to the paper

# Data 1c is very similar in format.
# Here is a full script to show it...
data_1c <- read_excel("file.xlsx", sheet_names[2])
data_1c %>% gather(na.rm = TRUE) %>%
    mutate(key = factor(key, c("Control", "Resolving", "RA"))) %>%
    ggplot(aes(key, value, shape = key, colour=key)) +
    geom_boxplot(outlier.size=0) +
    geom_jitter(width = 0.1) +
    ylim(0,2) + 
    labs(x = "",
        y = "Fap mRNA expression", 
        title = "Figure 1c",
        subtitle = "Croft et al, Nature, 2019") +
    theme_classic() + theme(legend.position="none") +
    stat_compare_means(comparisons = my_comparisons,
        symnum.args = list(cutpoints = c(0,0.0005, 1), 
            symbols = c("****","ns")))


# Data 1j is another dot plot...
sheet_names[5]
data_1j <- read_excel("file.xlsx", sheet_names[5])
data_1j %>%
    gather(na.rm = TRUE) %>%
    mutate(key = factor(key, c("Control", "STIA"))) %>%
    ggplot(aes(key, value, shape = key, colour=key)) +
    geom_boxplot(outlier.size=0) +
    geom_jitter(width = 0.1) +
    ylim(0,0.25) + 
    labs(x = "",
        y = "Bioluminescence", 
        title = "Figure 1j",
        subtitle = "Croft et al, Nature, 2019") +
    theme_classic() + theme(legend.position="none") +
    stat_compare_means()


## END











Tuesday, 16 April 2019

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 NumeracyHealth Psychology  2009, Vol. 28, No. 2, 210 –216  
  • "Infographic-style charts using the R waffle package" by N Saunders 


  • Friday, 15 February 2019

    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:

    Monday, 14 January 2019

    Making a box and whisker plot with some published proteomic data...

    Updated: 1st July 2019 - the source file has changed so some of the script had to be changed.
    I'm preparing some teaching materials for another Biochemical Society R training event with the draft title of R for Biochemists 201. Some more advanced material based on feedback for participants of R for Biochemists 101.
    In preparation, I've been looking at published proteomics data. I've come across a nice paper by a group in the Barts Cancer Institute in London. The paper is entitled "Proteomic and genomic integration identifies kinase and differentiation determinants of kinase inhibitor sensitivity in leukemia cells". It was published in the journal Leukaemia.
    I visited their lab once many years ago and I have heard the senior author, Pedro Cutillas, talk.  It is very interesting work, I think.
    They have made their data available so I've spend some time writing a script that makes one part of Figure 1a - a nice box and whisker plot.

    Here is the plot:



    and here is the script:
    START
    ## data import
    library(readxl)
    library(ggplot2)
    link <- "https://static-content.springer.com/esm/art%3A10.1038%2Fs41375-018-0032-1/MediaObjects/41375_2018_32_MOESM2_ESM.xlsx"
    download.file(link, "temp_data")
    data <- read_excel("temp_data", skip=2) 
    # skip = 2 stops the first two rows being part of the file
    # the next row is used as titles of the columns
    # skip needs to be determined by looking at the data

    # remove bottom two rows as only 36 patients
    data <- data[1:36,]


    # using the geom_boxplot() function, we can draw our graph
    ggplot(data = data,
        aes(FAB, log10(`MEKi (trametinib)...13`), colour = FAB)) +
        geom_boxplot(na.rm = TRUE)

    # we can make it look a bit more like the plot in the paper using geom_jitter()
    plot <- ggplot(data = data,
        aes(FAB, log10(`MEKi (trametinib)...13`), colour = FAB)) +
        geom_boxplot(na.rm = TRUE) +
        geom_jitter(width=0.15, na.rm = TRUE) +
        theme_bw() +
        labs( y = "Log10(EC50)nM",
            title = "Sensitivity of AML patients samples to MEK inhibitor", 
            subtitle = "Casado et al (2018) Leukaemia 32:1818–1822
    doi:10.1038/s41375-018-0032-1") +
        theme(legend.position="none")
    plot

    # to print a high resolution of this the tiff() function can be used. 
    tiff("plot.tiff", height = 12, width = 17, units = 'cm', compression = "lzw", res = 300)
    plot
    dev.off()
    END


    Some resources:

    Monday, 20 August 2018

    Exploring more immunization data...

    Last week, inspired by Factfulness, I made a graph showing the BCG immunization coverage for children at 1 year. The Factfulness graph didn't mention any specific immunization programme and World Health Organisation data monitors immunization coverage for other vaccines. For that reason, I thought it would be good to download and graph more of the available global immunization data.

    This allowed me to generate this graph which shows immunization coverage across the world for nine different vaccines. The various dates that monitoring starts shows that new immunization programs are being rolled out on a regular basis - good to see.




    Below is the code for downloading the data and making the graph...
    If you would rather just make the graph with some cleaner data, the data from Aug 20, 2018 is available on github and can be downloaded using the read_csv() code shown about half way down the script.

    ## START
    ##  download the data  
    library(tidyverse)
    # install.packages("WHO")
    library(WHO)

    # check out the codes of the WHO data...
    codes <- get_codes()
    # get codes for immunizations
    immun_codes <- codes[grepl("[Ii]mmuniz", codes$display), ]
    immun_codes$label

    # go through each of the 18 to find global data...
    # Number 1 has global data

    # download number 1
    # requires internet access
    immun_data <- as.tibble(get_data(immun_codes$label[1]))

    # filter for global data
    immun_data <- filter(immun_data, region == "(WHO) Global")

    # repeat download for next 2 to 18 WHO codes  
    # had to do this as a loop as couldn't get it to work using lapply...
    # start off with second value as first is above..
    # requires internet access and patience...
    for(i in 2:length(immun_codes$label)){
        # download the data
        data <- as.tibble(get_data(immun_codes$label[i])) 
        
        #tell you that it has downloaded...
        print(paste("Dataset",immun_codes$label[i], "downloaded."))
        
        # filter the data for Global values
        data <- filter(data, region == "(WHO) Global")
        
        # if there is some data bind_rows()
        if(nrow(data)>1){
            # bind_rows() function from dplyr
            immun_data <- bind_rows(immun_data, data)
        }else{   
            # if not just tell us....
            print("No global data in this set")
        }
    }





    # reduce columns using select() function  
    immun_data <- select(immun_data, gho, region, year, value )

    # to avoid having to download every time... save a local copy
    file_name <- paste0("global_immun_data", Sys.Date())
    write_csv(immun_data, file_name)

    ## ----read_back if you have saved to continue from here
    # immun_data <- read_csv(file_name)

    # read in data from github using read_csv() function

    # immun_data <- read_csv("https://raw.githubusercontent.com/brennanpincardiff/RforBiochemists/master/data/global_immun_data2018-08-20")


    # Let's make our plot...
    plot <- ggplot(immun_data, aes(x = year, y = value, 
        colour = gho)) +
      geom_line(size = 1)+ 
      theme(legend.position="none")

    plot

    # separate the plots with facet wrap
    plotf <- plot + facet_wrap(~gho)
    plotf


    ## The individual graph titles are difficult to read
    # Shorten them by removing text using gsub() = global substitution
    immun_data$gho_s <- gsub("immunization coverage among 1-year-olds",
                                "", immun_data$gho)
    immun_data$gho_s <- gsub("immunization coverage by the nationally recommended age",
                            "", immun_data$gho_s)

    # make the plot again
    plot <- ggplot(immun_data, aes(x = year, y = value, 
        colour = gho_s)) +
      geom_line(size = 1) + 
      theme(legend.position="none")

    # separate plots with facet_wrap
    plotf <- plot + facet_wrap(~gho_s)
    plotf <- plotf + theme_bw() + theme(legend.position="none")
    plotf

    # improve plot with y limits, titles & source
    source <- paste("Source: World Health Organisation, accessed:", Sys.Date())
    plotf <- plotf + ylim(0,100)
    plotf <- plotf + labs(x = NULL, y = "Immunization Rate",
          title = "Global immunization rates", 
        subtitle = source)
    plotf



    ## END

    Some Resources: