Introduction to Solving Biological Problems Using R - Week 3

Lucy Liu

The contents of this workshop are adapted from the University of Cambridge Bioinformatics training workshop “Data Manipulation and Visualisation using R”

Introduction

This section is on the R packages dplyr and ggplot.

Base R: We call R by itself (without adding any packages) ‘base R’. It comes with a lot of functions but you can also add to these functions by installing packages. You can think of this as like installing apps on your phone - it adds extra functions.

Packages: You can think of packages as a bundle of functions. The dplyr package is a bundle of functions for manipulating your data. The ggplot2 package is a bundle of functions for making plots and it is one of the most popular packages! They were both written by the same person and they work well together.

Package installation

You need to install packages first before you use them. To install a package that you have not used before requires 2 steps:

1. install.packages(‘dplyr’) - This command downloads the package onto your computer (if you are using the server, you don’t need to do this as it would already have been done). You only need to do this ONCE.
2. library(‘dplyr’) - This command ‘opens’ this package up and you need to do this every time you open up R. People using the server you will need to use this command.

Data import

The data that we are using today can be found here. It’s important to have a look at the data before you import it. Ahat symbol is separating each data value and which command should you use to import the data?

To import the data:

1. Save the data as a text file into your working directory.
2. Import your data into R.

patients <- read.csv("patient-data-cleaned.txt", sep = "\t")
library(dplyr)
library(ggplot2)

dplyr

A vital function in dplyr is called the ‘pipe’ %>% You can think about this as a pipe, that sends data through to the next function. We’ll use the function select() to demonstrate this.

select

The select() function let’s you select columns in your data.

patients %>%
  select(Name)

What’s happening here, is that we start with our data, the variable called ‘patients’, we are piping the data to the first function called select() which is selecting just the column named ‘Name’.

To select all columns apart from Name:

patients %>%
  select(-Name) %>%
  head()

You can also get columns within a particular range. This code below selects all columns between Name and Sex.

patients %>%
  select(Name:Sex) %>%
  head()

We can also use ends_with() to select all columns whose name ends with “t”.

patients %>%
  select(ends_with("t")) %>%
  head()

Exercise: select()

1. Select all the columns between Height and Grade_Level
2. Select all the columns between Height and Grade_Level, but NOT Pet
3. Select the columns Height and Weight

Filter

filter() function let’s you choose rows from the data frame. We can use the comparison operators ==, >, < and != to help use choose the rows we want.

To choose all the rows where the patient is male we say that we want the Sex column to be equal (==) to ‘Male’, like this:

patients %>%
  filter(Sex == 'Male') %>%
  head()

To choose all the rows where the patient’s name is NOT ‘Michael’, we say that we want the Name column to not be equal to (!=) ‘Michael’:

patients %>%
  filter(Name != 'Michael') %>%
  head()

To choose all the rows where the patient is Male and their Pet is a Dog we do this:

patients %>%
  filter(Sex == 'Male', Pet == 'Dog') %>%
  head()

You can also use & here. It is exactly the same as a comma.

To choose all the rows where the patient’s height is above 170 or (symbol for or is ‘|’) their pet is a dog, we use:

patients %>%
  filter(Height > '170' | Pet == 'Dog') %>%
  head()

To choose all the rows where the patient’s pet is EITHER a dog or a bird we use:

patients %>%
  filter(Pet %in% c("Dog", "Bird")) %>%
  head()

The %in% function finds matches. It looks for all the rows where the column Pets, matches Dog or bird.

Note that there is a difference between %in% and == This:

patients %>%
  filter(Pet %in% c("Dog", "Bird"))

will not return the same result as this:

patients %>%
  filter(Pet == c("Dog", "Bird"))

The tip is, when you are looking for more than one thing (e.g. pet is either dog or bird), use %in%. When you are looking for only one thing (e.g. pet is a dog), use ==.

You can combine select and filter together like this:

patients %>%
  filter(Pet %in% c("Dog", "Bird")) %>%
  select(Name:Sex) %>%
  head()

This first chooses the rows where the patient’s pet is a dog or a bird and then selects the columns from Name to Sex.

mutate

mutate() allows you to create a new column and add it to the end of your data frame.

We can make a new column where the height is in metres instead of centimetres. We have given our new column the name “height_in_m” (we can name this anything that we choose) and we have specified that we want the number in the height column divided by 100.

patients %>%
  mutate(height_in_m = Height/100) %>%
  head()

Challenge: mutate and select

1. Create a new column called BMI which gives the BMIs of all the patients
2. Choose the rows where the overweight (BMI>25) smokers that are still alive
3. Choose the rows where the Female patients from New York or New Jersey
4. Choose the columns Name and Sex and the rows where patients are overweight and smokers

Remember that the order you get R to perform the functions matters.

group_by & summarise

The group_by() and summarise() functions are very useful when you use them together. They allow you to compute summary statistics for selected columns in our dataset.

Let’s first look at what summarise does by itself.

patients %>%
  summarise(mean_height = mean(Height))

It’s taken all the heights and given us the mean. Note that it returns just one number.

Now let’s have a look at what it does when you use it with group_by:

patients %>%
  group_by(Sex) %>%
  summarise(mean_height = mean(Height))

It groups the rows into males and females, and then it has given us the mean height for each group. Note that it returns two numbers as there were two groups.

You can calculate more than one summary statistic as well. Below we are grouping our patients by sex and race (there are 2 sexes and 4 difference races, thus we have 8 groups in total) and then calculating the mean height and weight for each group. Note that we will get 8 numbers because there are 8 groups.

patients %>%
  group_by(Sex, Race) %>%
  summarise(mean_height = mean(Height), mean_weight = mean(Weight))

Challenge: group_by & summarise

1. Use summarise to find the mean height and weight of patients in the patients dataset.
2. Find the mean weight of male and female patients from each state
3. Find the proportion of patients of each sex AND smoking status that are deceased

Hint: TRUE’s are stored as the number 1 and FALSE’s are stored as the number 0 in R. Thus if we create a vector of TRUE’s and FALSE’s and then take the mean, we get the proportion of TRUE’s in the vector.

a <- c(TRUE, FALSE, FALSE, TRUE, TRUE)
mean(a)

Prepare data for ggplot

For ggplot, we will use the original data we started with and add a BMI column to the end. We will use this dataframe to make plots with ggplot.

patients <- patients %>%
  mutate(BMI = Weight/(Height/100)^2)

ggplot2

ggplot2 is one of the most popular packages in R and allows you to create good looking graphs easily. It builds on the concept of the “Grammar of Graphics” (Wilkinson 2005, Bertin 1983) which describes a consistent syntax for the construction of a wide range of complex graphics by a concise description of their components.

The structured syntax and high level of abstraction used by ggplot2 should allow for the user to concentrate on the visualisations instead of creating the underlying code.

On top of this central philosophy ggplot2 has:

• Increased flexible over many plotting systems.
• An advanced theme system for professional/publication level graphics.
• Large developer base – Many libraries extending its flexibility.
• Large user base – Great documentation and active mailing list.

Let’s start with the basic anatomy of ggplot code - or how to give R instructions to make your plot. First you must tell ggplot what data you want to plot.

Here we are telling ggplot that we want to plot the dataframe patients:

ggplot(patients)

Now we must tell ggplot how to make the graph using aes (“aesthetics”). Here you map variables to aesthetics - i.e. what do you want on the x and y axis.

ggplot(patients, aes(x=Height, y=Weight))

Depending on the type of graph, there are other aesthetics you can map variables to. Let’s start with a scatter graph. Apart from x and y axes, you can also map variables to size of points, colour of points and shape of points. We can map each of these aesthetics to a variable:

Finally, we must tell ggplot what kind of graph we want to make. To do this we simply ‘add’ a geom like so:

ggplot(patients, aes(x=Height, y=Weight)) +
  geom_point()

geom_point is a scatter graph. Note the plus at the end of the first line - this is very important!

Challenge: ggplot 1

Recreate the above plot

Geoms - Plot types

The geom describes the type of plot used. Several geoms are available in ggplot2:

• geom_point() - Scatter plots
• geom_line() - Line plots
• geom_smooth() - Fitted line plots
• geom_bar() - Bar plots
• geom_boxplot() - Boxplots
• geom_jitter() - Jitter to plots
• geom_histogram() - Histogram plots
• geom_density() - Density plots
• geom_text() - Text to plots
• geom_errorbar() - Errorbars to plots
• geom_violin() - Violin plots

Line plots

pcPlot <- ggplot(data=patients,
        mapping=aes(x=Height,y=Weight))

pcPlot_line <- pcPlot+geom_line() 

pcPlot_line

pcPlot <- ggplot(data=patients,
        mapping=aes(x=Height,y=Weight))

pcPlot_smooth <- pcPlot+geom_smooth() 

pcPlot_smooth

Bar and frequency plots

pcPlot <- ggplot(data=patients,
        mapping=aes(x=Sex))

pcPlot_bar <- pcPlot+geom_bar() 

pcPlot_bar

pcPlot <- ggplot(data=patients,
        mapping=aes(x=Height))

pcPlot_hist <- pcPlot+geom_histogram() 

pcPlot_hist

pcPlot <- ggplot(data=patients,
        mapping=aes(x=Height))

pcPlot_density <- pcPlot+geom_density() 

pcPlot_density

Box and violin plots

pcPlot <- ggplot(data=patients,
        mapping=aes(x=Sex,y=Height))

pcPlot_boxplot <- pcPlot+geom_boxplot() 

pcPlot_boxplot

pcPlot <- ggplot(data=patients,
        mapping=aes(x=Sex,y=Height))

pcPlot_violin <- pcPlot+geom_violin() 

pcPlot_violin

An overview of geoms and thier arguments can be found at ggplot2 documentation or within the ggplot2 cheatsheet.

-ggplot2 documentation

-ggplot2 Cheatsheet

Facets

One very useful feature of ggplot is faceting. This allows you to produce plots subset by variables in your data.

To facet our data into multiple plots we can use the facet_wrap or facet_grid function specifying the variable we split by.

The facet_grid function is well suited to splitting the data by two factors.

Here we can plot the data with the Smokes variable as rows and Sex variable as columns.

facet_grid(Rows~Columns)

pcPlot <- ggplot(data=patients,aes(x=Height,y=Weight,colour=Sex))+geom_point()
pcPlot + facet_grid(Smokes~Sex)

To split by one factor we can apply the facet_grid() function ommiting the variable before the “~”" to facet along columns in plot.

facet_grid(~Columns)

pcPlot <- ggplot(data=patients,aes(x=Height,y=Weight,colour=Sex))+geom_point()
pcPlot + facet_grid(~Sex)

To split along rows in plot, the variable is placed before the “~.”.

facet_grid(Rows~.)

pcPlot <- ggplot(data=patients,aes(x=Height,y=Weight,colour=Sex))+geom_point()
pcPlot + facet_grid(Sex~.)

The facet_wrap() function offers a less grid based structure but is well suited to faceting data by one variable.

For facet_wrap() we follow as similar syntax to facet_grid()

pcPlot <- ggplot(data=patients,aes(x=Height,y=Weight,colour=Sex))+geom_point()
pcPlot + facet_wrap(~Smokes)

For more complex faceting both facet_grid and facet_wrap can accept combinations of variables.

Using facet_wrap

pcPlot <- ggplot(data=patients,aes(x=Height,y=Weight,colour=Sex))+geom_point()
pcPlot + facet_wrap(~Pet+Smokes+Sex)

Or in a nice grid format using facet_grid() and the Smokes variable against a combination of Gender and Pet.

pcPlot + facet_grid(Smokes~Sex+Pet)

Adding titles for plot and labels.

So far no plot titles have been specified. Plot titles can be specified using the labs functions.

pcPlot <- ggplot(data=patients,
        mapping=aes(x=Weight,y=Height))+geom_point()
pcPlot+labs(title="Weight vs Height",y="Height (cm)")

or specified using the ggtitle and xlab/ylab functions.

pcPlot <- ggplot(data=patients,
        mapping=aes(x=Height,y=Weight))+geom_point()
pcPlot+ggtitle("Weight vs Height")+ylab("Height (cm)")

Change the plotting order in a boxplot

We will shortly discuss how to change various aspects of the plot layout and appearance. However, a common-asked question is how to change the order in which R plots a categorical variable. Consider the boxplot to compare weights of males and females:-

ggplot(patients, aes(x=Sex, y=Weight)) + geom_boxplot()

Here, R decides the order to arrange the boxes according to the levels of the categorical variable. By default this is the alphabetical order. i.e. Female before Male.

patients$Sex

Depending on the message we want the plot to convey, we might want control over the order of boxes. The factor functions allows us to explictly change the order of the levels.

factor(patients$Sex,levels=c("Male","Female"))

With the pipe syntax we just learnt about, we can change the variable on-the-fly and generate the plot

patients %>% 
  mutate(Sex = factor(Sex,levels=c("Male","Female"))) %>% 
  ggplot(aes(x=Sex, y=Weight)) + geom_boxplot()

Scales

Scales and their legends have so far been handled using ggplot2 defaults. ggplot2 offers functionality to have finer control over scales and legends using the scale methods.

Scale methods are divided into functions by combinations of

• the aesthetics they control.

• the type of data mapped to scale.

scale_aesthetic_type

Try typing in scale_ then tab to autocomplete. This will provide some examples of the scale functions available in ggplot2.

Although different scale functions accept some variety in their arguments, common arguments to scale functions include -

• name - The axis or legend title

• limits - Minimum and maximum of the scale

• breaks - Label/tick positions along an axis

• labels - Label names at each break

Controlling the X and Y scale.

Both continous and discrete X/Y scales can be controlled in ggplot2 using the

scale_(x/y)_(continous/discrete)

In this example we control the continuous sale on the x-axis by providing a name, X-axis limits, the positions of breaks (ticks/labels) and the labels to place at breaks.

pcPlot +
  geom_point() +
  facet_grid(Smokes~Sex)+
  scale_x_continuous(name="height ('cm')",
                     limits = c(100,200),
                     breaks=c(125,150,175),
                     labels=c("small","justright","tall"))

Similary control over discrete scales is shown below.

pcPlot <- ggplot(data=patients,aes(x=Sex,y=Height))
pcPlot +
  geom_violin(aes(x=Sex,y=Height)) +
  scale_x_discrete(labels=c("Women", "Men"))

Multiple X/Y scales can be combined to give full control of axis marks.

pcPlot <- ggplot(data=patients,aes(x=Sex,y=Height,fill=Smokes))
pcPlot +
  geom_violin(aes(x=Sex,y=Height)) +
  scale_x_discrete(labels=c("Women", "Men"))+
  scale_y_continuous(breaks=c(160,180),labels=c("Short", "Tall"))

Controlling other scales.

When using fill,colour,linetype, shape, size or alpha aesthetic mappings the scales are automatically selected for you and the appropriate legends created.

pcPlot <- ggplot(data=patients,
                 aes(x=Height,y=Weight,colour=Sex))
pcPlot + geom_point(size=4)
                  

In the above example the discrete colours for the Sex variable was selected by default.

Manual discrete colour scale

Manual control of discrete variables can be performed using scale_aes_Of_Interest_manual with the values parameter. Additionally in this example an updated name for the legend is provided.

pcPlot <- ggplot(data=patients,
                 aes(x=Height,y=Weight,colour=Sex))
pcPlot + geom_point(size=4) + 
  scale_color_manual(values = c("Green","Purple"),
                     name="Gender")
                  

Colorbrewer for discrete colour scale

Here we have specified the colours to be used (hence the manual) but when the number of levels to a variable are high this may be impractical and often we would like ggplot2 to choose colours from a scale of our choice.

The brewer set of scale functions allow the user to make use of a range of palettes available from colorbrewer.

• Diverging

BrBG, PiYG, PRGn, PuOr, RdBu, RdGy, RdYlBu, RdYlGn, Spectral

• Qualitative

Accent, Dark2, Paired, Pastel1, Pastel2, Set1, Set2, Set3

• Sequential

Blues, BuGn, BuPu, GnBu, Greens, Greys, Oranges, OrRd, PuBu, PuBuGn, PuRd, Purples, RdPu, Reds, YlGn, YlGnBu, YlOrBr, YlOrRd

pcPlot <- ggplot(data=patients,
                 aes(x=Height,y=Weight,colour=Pet))
pcPlot + geom_point(size=4) + 
  scale_color_brewer(palette = "Set2")
                  

For more details on palette sizes and styles visit the colorbrewer website and ggplot2 reference page.

• colorbrewer

• ggplot2 colour scales

Continuous colour scales

So far we have looked a qualitative scales but ggplot2 offers much functionality for continuous scales such as for size, alpha (transparancy), colour and fill.

• scale_alpha_continuous() - For Transparancy

• scale_size_continuous() - For control of size.

Both these functions accept the range of alpha/size to be used in plotting.

Below the range of alpha to be used in plot is limited to between 0.5 and 1

pcPlot <- ggplot(data=patients,
                 aes(x=Height,y=Weight,alpha=BMI))
pcPlot + geom_point(size=4) + 
  scale_alpha_continuous(range = c(0.5,1))
                  

Below the range of sizes to be used in plot is limited to between 3 and 6

pcPlot <- ggplot(data=patients,
                 aes(x=Height,y=Weight,size=BMI))
pcPlot + geom_point(alpha=0.8) + 
  scale_size_continuous(range = c(3,6))
                  

The limits of the scale can also be controlled but it is important to note data outside of scale is removed from plot.

pcPlot <- ggplot(data=patients,
                 aes(x=Height,y=Weight,size=BMI))
pcPlot + geom_point() + scale_size_continuous(range = c(3,6), limits = c(25,40))
                  

What points of scale to be labeled and labels text can also be controlled.

pcPlot <- ggplot(data=patients,
                 aes(x=Height,y=Weight,size=BMI))
pcPlot + geom_point() + 
  scale_size_continuous(range = c(3,6), 
                        breaks=c(25,30), 
                        labels=c("Good","Good but not 25"))

Control of colour/fill scales can be best achieved through the gradient subfunctions of scale.

• scale_(colour/fill)_gradient - 2 colour gradient (eg. low to high BMI)

• scale_(colour/fill)_gradient2 - Diverging colour scale with a midpoint colour (e.g. Down, No Change, Up)

Both functions take a common set of arguments:-

• low - Colour for low end of gradient scale
• high - Colour for high end of gradient scale.
• na.value - Colour for any NA values.

An example using scale_colour_gradient below sets the low and high end colours to White and Red respectively

pcPlot <- ggplot(data=patients,
                 aes(x=Height,y=Weight,colour=BMI))
pcPlot + geom_point(size=4,alpha=0.8) + 
  scale_colour_gradient(low = "White",high="Red")
                  

Similarly we can use the scale_colour_gradient2 function which allows for the specification of a midpoint value and its associated colour.

pcPlot <- ggplot(data=patients,
                 aes(x=Height,y=Weight,colour=BMI))
pcPlot + geom_point(size=4,alpha=0.8) + 
  scale_colour_gradient2(low = "Blue",mid="Black", high="Red", midpoint = median(patients$BMI))

As with previous continous scales, limits and custom labels in scale legend can be added.

pcPlot <- ggplot(data=patients,
                 aes(x=Height,y=Weight,colour=BMI))
pcPlot + geom_point(size=4,alpha=0.8) + 
  scale_colour_gradient2(low = "Blue",
                         mid="Black",
                         high="Red",
                         midpoint = median(patients$BMI),
                         breaks=c(25,30),labels=c("Low","High"),
                         name="Body Mass Index")

Multiple scales may be combined to create high customisable plots and scales

pcPlot <- ggplot(data=patients,
                 aes(x=Height,y=Weight,colour=BMI,shape=Sex))
pcPlot + geom_point(size=4,alpha=0.8)+ 
  scale_shape_discrete(name="Gender") +
  scale_colour_gradient2(low = "Blue",mid="Black",high="Red",
                         midpoint = median(patients$BMI),
                         breaks=c(25,30),labels=c("Low","High"),
                         name="Body Mass Index")

Statistical transformations.

In ggplot2 many of the statistical transformations are performed without any direct specification e.g. geom_histogram() will use stat_bin() function to generate bin counts to be used in plot.

An example of statistical methods in ggplot2 which are very useful include the stat_smooth() and stat_summary() functions.

The stat_smooth() function can be used to fit a line to the data being displayed.

pcPlot <- ggplot(data=patients,
        mapping=aes(x=Weight,y=Height))
pcPlot+geom_point()+stat_smooth()

By default a “loess” smooth line is plotted by stat_smooth. Other methods available include lm, glm,gam,rlm.

pcPlot <- ggplot(data=patients,
        mapping=aes(x=Weight,y=Height))
pcPlot+geom_point()+stat_smooth(method="lm")

A useful feature of ggplot2 is that it uses previously defined grouping when performing smoothing.

If colour by Sex is an aesthetic mapping then two smooth lines are drawn, one for each sex.

pcPlot <- ggplot(data=patients,
        mapping=aes(x=Weight,y=Height,colour=Sex))
pcPlot+geom_point()+stat_smooth(method="lm")

This behaviour can be overridden by specifying an aes within the stat_smooth() function and setting inherit.aes to FALSE.

pcPlot <- ggplot(data=patients,
        mapping=aes(x=Weight,y=Height,colour=Sex))
pcPlot+geom_point()+stat_smooth(aes(x=Weight,y=Height),method="lm",inherit.aes = F)

Another useful method is stat_summary() which allows for a custom statistical function to be performed and then visualised.

The fun.y parameter specifies a function to apply to the y variables for every value of x.

In this example we use it to plot the quantiles of the Female and Male Height data

pcPlot <- ggplot(data=patients,
        mapping=aes(x=Sex,y=Height))+geom_jitter()
pcPlot+stat_summary(fun.y=quantile,geom="point",colour="purple",size=8)