Welcome to the tidyverse

≈ 15 min

The tidyverse is a collection of R packages for reading, tidying, transforming and plotting data, all designed to work together. This course builds on R basics, so if you haven’t done that course yet, do it first.

What the tidyverse is

R comes with everything you need to work with data, but over the years its functions have grown up in different styles, with different names and orders of arguments for similar jobs. The tidyverse is a collection of packages that share a single design. Data are kept in tibbles, laid out in a consistent (“tidy”) way, and almost every function takes a tibble as its first argument and gives you a tibble back, so that the functions fit together. The core packages are:

Package	What it is for
tibble	Tibbles, the tidyverse’s tables of data
readr	Reading data from text files into tibbles, and writing tibbles to files
tidyr	Tidying data, so that it is laid out in the right shape
dplyr	Choosing, filtering, changing and summarising data
ggplot2	Drawing plots and graphs

There are a few more, such as stringr for working with text, forcats for factors and lubridate for dates, which you won’t need in this course.

Installing and loading packages

A package is a bundle of functions (and sometimes data) that adds to what R can do. There are two steps to using one:

install.packages("tidyverse") downloads the package and installs it. You only need to do this once: on your own computer it stays installed, and here it is kept in your browser for the next time you come back.
library(tidyverse) loads the package into R’s memory, so that you can use its functions. You need to do this in every new R session: each time you start RStudio, restart R, or come back to this course.

When it loads, the tidyverse tells you which packages it has attached, and their versions. It also lists some conflicts, such as dplyr::filter() masks stats::filter(). This means that two packages have a function with the same name, and that, from now on, filter means the one from dplyr. That is what you want, so you can ignore the message. If you ever need the other one, you can name its package, as in stats::filter().

How this course works

This course works in the same way as R basics. The steps of each activity tick themselves when R sees that you have done them (or you can tick them yourself), questions with a Check button are marked straight away, and Show answer appears once you have had a go at a question. Parts marked If you have time don’t count towards your certificate. When you have done everything else, the last chapter has your certificate, with the name you give below.

Most of the steps in this course ask you to write some R yourself, and the text before each activity shows you examples of what you need.

Keep your commands in a script

Whenever R starts again (when you come back to the course, or restart R), it starts afresh: the tidyverse isn’t loaded, the working directory is your home folder, and none of your data have been read in. If you keep your commands in a script as you go, starting with library(tidyverse), you can get back to where you were by clicking Source at the top of the script, which runs the whole of it.

AI tools and your data

The AI tab has an AI assistant that you can ask about R, if you have a key for an AI service. Anything you type into it is sent to an AI company. Never paste personal data or patient data into it, or into any other AI tool.

Get ready

Add your name for your certificate (if you did R basics in this browser, it is already there):
Start a new script (File ▸ New File ▸ R Script) and save it as tidyverse.R. Write the commands from this course in it as you go.
Install the tidyverse: install.packages("tidyverse") . This takes about half a minute the first time.
Load it with library(tidyverse) , and read the messages it gives you. The first time, this takes about another half a minute.

Reading data from files

≈ 25 min

Most data reach R in a file. readr, one of the packages of the tidyverse, reads files of data into tibbles, and it tells you how it has interpreted each column.

readr

readr has a function for each of the common kinds of text file of data, and functions that write tibbles back out to files:

Function	What it does
`read_csv("file.csv")`	reads a file of comma-separated values
`read_tsv("file.txt")`	reads a file of tab-separated values (often ending in `.txt` or `.tsv`)
`write_csv(data, "file.csv")`	writes a tibble to a file of comma-separated values
`write_tsv(data, "file.txt")`	writes a tibble to a file of tab-separated values

Like any other function, read_csv() only gives you the data: it doesn’t keep them for you. To work with the data afterwards, store them in an object, as in child <- read_csv("Child_Variants.csv").

File paths and the working directory

A file path tells R where a file is. ~/data/small_file.txt is the file small_file.txt in the data folder of your home folder, which is what ~ stands for. You can give R the full path of a file every time, but it soon becomes tedious. Instead, you can set R’s working directory: the folder in which R looks for files, unless you tell it otherwise. Once it is set, the name of a file in that folder is enough.

getwd() tells you what the working directory is, and setwd("~/data") changes it to your data folder.
In RStudio, you can also choose Session ▸ Set Working Directory ▸ Choose Directory…, or go to the folder in the Files pane and choose More ▸ Set As Working Directory. RStudio shows you the setwd() command it runs, which you can copy into your script.

RStudio can also fill in file names for you. Type read_tsv(""), put the cursor between the quotes, and press Tab: RStudio lists the files in the working directory, and you can pick one rather than typing its name (and perhaps mistyping it).

What readr tells you

When it reads a file, readr guesses what type of data each column holds, and it tells you what it decided. Try it with a small file of seven people, giving its full path, because you haven’t set the working directory yet: trumpton <- read_tsv("~/data/trumpton.txt")

> trumpton <- read_tsv("~/data/trumpton.txt")
Rows: 7 Columns: 5
── Column specification ────────────────────────────────────────────────────────
Delimiter: "\t"
chr (2): LastName, FirstName
dbl (3): Age, Weight, Height

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
> trumpton
# A tibble: 7 × 5
  LastName FirstName   Age Weight Height
  <chr>    <chr>     <dbl>  <dbl>  <dbl>
1 Hugh     Chris        26     90    175
2 Pew      Adam         32    102    183
3 Barney   Daniel       18     88    168
4 McGrew   Chris        48     97    155
5 Cuthbert Carl         28     91    188
6 Dibble   Liam         35     94    145
7 Grub     Doug         31     89    164

chr means character (text), and dbl means a number (it is short for double, which is how computers store numbers with decimal places). Always check readr’s guesses: a column of numbers that has been read as chr usually means that something in the column isn’t a number, such as a word, or a comma used as a decimal point. When you print a tibble, you see the type of each column under its name, and only the first ten rows; View(), or a click on the tibble in the Environment pane, shows you all of it.

As in R basics, $ takes a single column out of a tibble, as a vector, so that you can give it to functions such as mean():

> mean(trumpton$Height)
[1] 168.2857

Read some files

Set the working directory to your data folder.
Read small_file.txt, which is a tab-separated file, into a new tibble called small. Use Tab to fill in the name of the file.
View small, and check that the types that readr guessed for its columns make sense.
Find the median of the log2 of the lengths in small. Use $ to take out the Length column, log2() to transform the values, and median() to find the median of the result.
Read Child_Variants.csv, which is a file of comma-separated values, into a tibble called child.
Print child in the Console to see what it looks like, and then open it in the viewer by clicking on it in the Environment pane.
Calculate the mean and the standard deviation of the MutantReadPercent column of child.

What is the median of the log2 of the lengths?

6.227664, from median(log2(small$Length)).

What is the mean of the MutantReadPercent column?

59.88219, from mean(child$MutantReadPercent).

And its standard deviation?

24.11675, from sd(child$MutantReadPercent).

readr tells you chr (6): CHR, dbSNP, REF, ALT, GENE, ENST. Why has it read the column CHR, the chromosome, as text rather than as numbers?

Because some chromosomes are called X, Y or MT, which aren’t numbers Because readr reads every column as text Because CSV files can only hold text

Tidy data

≈ 10 min

A table of data can often be laid out in more than one way. The tidyverse has a clear opinion about which way is best, and its functions are designed around it.

Think of a simple experiment, in which two genes (ABC1 and DEF1) are measured in two genotypes (WT and KO), with three replicates of each. Here are the results in wide format, with a row for each gene:

Gene	WT_1	WT_2	WT_3	KO_1	KO_2	KO_3
ABC1	8.86	4.18	8.90	4.00	14.52	13.39
DEF1	29.60	41.22	36.15	11.18	16.68	1.64

This table is compact and easy to read, and you can see the measurements of each gene at a glance. But the genotype and the replicate of each value are hidden in the names of the columns, and the values are spread over six columns. If you wanted to record something else about each sample, such as the batch it was processed in, there would be nowhere obvious to put it.

Here are the same results in long format:

Gene	Genotype	Replicate	Value
ABC1	WT	1	8.86
ABC1	WT	2	4.18
ABC1	WT	3	8.90
ABC1	KO	1	4.00
ABC1	KO	2	14.52
ABC1	KO	3	13.39
DEF1	WT	1	29.60
DEF1	WT	2	41.22
DEF1	WT	3	36.15
DEF1	KO	1	11.18
DEF1	KO	2	16.68
DEF1	KO	3	1.64

This is more verbose, because values such as the gene names are repeated. But the genotype and the replicate are now columns of their own, all the measurements are in a single column, and recording the batch of each sample would simply mean adding a Batch column.

What makes data tidy

The long format is what the tidyverse calls tidy data:

each variable has its own column,
each observation has its own row, and
each value has its own cell.

The functions of the tidyverse expect data to be tidy. For a plot of the values coloured by genotype, for example, ggplot2 needs a column that says what the genotype of each value is, and the wide table doesn’t have one. Data often arrive in wide format, because that is how people lay out tables for reading, and the tidyr package has functions for converting between the two formats: pivot_longer() and pivot_wider().

In the long format, how would you record the batch in which each sample was processed?

Add a column called Batch Add a new row for each batch Put the batch into the names of the columns, as in WT_1_batch2

A table has one row for each patient, with the columns Patient, Monday, Tuesday and Wednesday, which hold the patient’s blood pressure on each day. Is it tidy?

No: the day is a variable, so a tidy version would have the columns Patient, Day and BloodPressure Yes: each patient has a row of their own Yes, as long as there are no missing values

If you have time Making wide data long

This part is optional: it doesn’t count towards your certificate.

Type in the wide table as a tibble, and then make it long with pivot_longer(). Its arguments say which columns to gather up (all except Gene), what to call the new columns that the old column names go into, how to split the old names (at the _), and what to call the column of values:

wide <- tibble(
  Gene = c("ABC1", "DEF1"),
  WT_1 = c(8.86, 29.60), WT_2 = c(4.18, 41.22), WT_3 = c(8.90, 36.15),
  KO_1 = c(4.00, 11.18), KO_2 = c(14.52, 16.68), KO_3 = c(13.39, 1.64)
)

> pivot_longer(wide, -Gene, names_to = c("Genotype", "Replicate"), names_sep = "_", values_to = "Value")
# A tibble: 12 × 4
   Gene  Genotype Replicate Value
   <chr> <chr>    <chr>     <dbl>
 1 ABC1  WT       1          8.86
 2 ABC1  WT       2          4.18
 3 ABC1  WT       3          8.9
 4 ABC1  KO       1          4
 5 ABC1  KO       2         14.5
 6 ABC1  KO       3         13.4
 7 DEF1  WT       1         29.6
 8 DEF1  WT       2         41.2
 9 DEF1  WT       3         36.2
10 DEF1  KO       1         11.2
11 DEF1  KO       2         16.7
12 DEF1  KO       3          1.64

Tidy it

Make the wide table long, as above.

Choosing columns and rows

≈ 30 min

dplyr, the tidyverse’s package for manipulating data, has a function for choosing columns, select(), and another for choosing rows, filter().

All of dplyr’s functions take a tibble as their first argument, and they all give you back a new tibble. The tibble you started with isn’t changed, so if you want to keep a result, store it in an object. Inside the functions, you refer to the columns of the tibble just by their names, without quotes and without $.

These examples use trumpton, the small tibble from the last chapter. If it isn’t in your Environment pane, read it in again: trumpton <- read_tsv("~/data/trumpton.txt")

> trumpton
# A tibble: 7 × 5
  LastName FirstName   Age Weight Height
  <chr>    <chr>     <dbl>  <dbl>  <dbl>
1 Hugh     Chris        26     90    175
2 Pew      Adam         32    102    183
3 Barney   Daniel       18     88    168
4 McGrew   Chris        48     97    155
5 Cuthbert Carl         28     91    188
6 Dibble   Liam         35     94    145
7 Grub     Doug         31     89    164

Choosing columns with select()

Give select() the names of the columns you want, in the order you want them. You can also give their positions, or put a minus sign in front of a column that you want to leave out:

> select(trumpton, FirstName, LastName, Weight)
# A tibble: 7 × 3
  FirstName LastName Weight
  <chr>     <chr>     <dbl>
1 Chris     Hugh         90
2 Adam      Pew         102
3 Daniel    Barney       88
4 Chris     McGrew       97
5 Carl      Cuthbert     91
6 Liam      Dibble       94
7 Doug      Grub         89
> select(trumpton, 2, 4)
# A tibble: 7 × 2
  FirstName Weight
  <chr>      <dbl>
1 Chris         90
2 Adam         102
3 Daniel        88
4 Chris         97
5 Carl          91
6 Liam          94
7 Doug          89
> select(trumpton, -LastName)
# A tibble: 7 × 4
  FirstName   Age Weight Height
  <chr>     <dbl>  <dbl>  <dbl>
1 Chris        26     90    175
2 Adam         32    102    183
3 Daniel       18     88    168
4 Chris        48     97    155
5 Carl         28     91    188
6 Liam         35     94    145
7 Doug         31     89    164

Choosing rows with filter()

filter() keeps the rows for which a test is TRUE. The tests are the comparisons you used with vectors in R basics; text goes in quotes, as always.

Test	Or	Meaning
`Height > 170`	`Height < 170`	greater than, less than
`Height >= 170`	`Height <= 180`	greater than or equal to, less than or equal to
`Age == 26`	`FirstName == "Chris"`	equal to
`Age != 26`	`FirstName != "Chris"`	not equal to

> filter(trumpton, Height >= 170)
# A tibble: 3 × 5
  LastName FirstName   Age Weight Height
  <chr>    <chr>     <dbl>  <dbl>  <dbl>
1 Hugh     Chris        26     90    175
2 Pew      Adam         32    102    183
3 Cuthbert Carl         28     91    188
> filter(trumpton, FirstName == "Chris")
# A tibble: 2 × 5
  LastName FirstName   Age Weight Height
  <chr>    <chr>     <dbl>  <dbl>  <dbl>
1 Hugh     Chris        26     90    175
2 McGrew   Chris        48     97    155

A test can also use a function to transform a column first. In this table, difference can be large in either direction, and abs(), which gives the absolute value of a number (its size, ignoring its sign), finds the rows where the difference is more than 5 either way:

> transform.data <- read_tsv("transform_data.txt", show_col_types = FALSE)
> transform.data
# A tibble: 10 × 3
       WT      KO difference
    <dbl>   <dbl>      <dbl>
 1 -5.11   -3.29       1.81
 2  1.12   -1.85      -2.97
 3 -3.99   -3.77       0.222
 4 -4.18   -2.46       1.72
 5 -1.93  -10.0       -8.10
 6 -8.69   -2.38       6.31
 7 -0.670   2.73       3.40
 8 -1.15   -2.59      -1.43
 9 -1.98    1.83       3.80
10 -1.06    0.372      1.43
> filter(transform.data, abs(difference) > 5)
# A tibble: 2 × 3
     WT     KO difference
  <dbl>  <dbl>      <dbl>
1 -1.93 -10.0       -8.10
2 -8.69  -2.38       6.31

(show_col_types = FALSE tells read_tsv() not to print the column types.) You don’t need to save the results of the steps below, but do look at each one, to check that it has worked.

Filter the small file

From small, extract all the rows where the Category is A.
Extract the samples with a Length of more than 80.
Remove the Sample column, to leave just Length and Category.

Filter the child’s variants

From child, extract the variants on the mitochondrial chromosome, which is called MT.
Extract the variants with a MutantReadPercent of 70 or more.
Extract the variants with a quality of exactly 200. (Look at the data to find the name of the column.)
Extract all the variants in the gene IGFN1.
Remove the ENST and dbSNP columns.

How many samples in small have a Length of more than 80?

16, from filter(small, Length > 80): the tibble it prints has 16 rows.

How many of the child’s variants are in IGFN1?

17, from filter(child, GENE == "IGFN1").

You want only the variants whose quality is 200. Which function do you need?

filter() select()

Why does filter(small, Category == A) give an error?

A is text, so it needs quotes: without them, R looks for a column or an object called A Because a test inside filter() needs = rather than == Because the column is called category, with a small c

If you have time Fix a broken script

This part is optional: it doesn’t count towards your certificate.

The file data/broken_selects.R is a short script that tries to do some of the things you have just done, but almost every line has a mistake in it. Open it in RStudio , and work through it a line at a time, fixing each line until it runs and gives the right result. Check the result of every line before you move on to the next.

Fix it

Fix every line of broken_selects.R, and save the script.

What was wrong with each line?

The path in setwd() is for someone else’s computer (here it should be setwd("~/data")); Child_Variants.csv is a CSV file, so it needs read_csv(), not read_tsv(); Category==A needs quotes around A; select(small, length>80) should be filter(), because it chooses rows, and the column is Length, with a capital L; select() needs the tibble small, not the name of the file, as its first argument; and the tibble is called child, not Child.

Combining steps with pipes

≈ 25 min

Real questions usually take several steps: filter on one thing, then on another, then choose some columns. Pipes let you write those steps one after another, in the order in which they happen.

Suppose you want to find the people in trumpton who are at least 170 cm tall and called Chris, and see just their age and weight. That takes three steps. One way to do it is to store the result of each step in an object, and use that object in the next step:

> filter(trumpton, Height >= 170) -> answer1
> filter(answer1, FirstName == "Chris") -> answer2
> select(answer2, Age, Weight)
# A tibble: 1 × 2
    Age Weight
  <dbl>  <dbl>
1    26     90

This works, but it leaves objects you don’t need (answer1 and answer2) lying around in R’s memory, and it is hard to see what is going on.

The pipe: %>%

Every dplyr function takes a tibble as its first argument and gives you a tibble back. That means you can pass the result of one function straight into the next, without storing it. The pipe, %>%, does exactly that: it takes whatever is on its left and gives it to the function on its right as that function’s first argument. These two commands are the same:

select(trumpton, -LastName)
trumpton %>% select(-LastName)

So the three steps above become a single command, which you can read from top to bottom, saying “and then” for each %>%: take trumpton, and then keep the people who are at least 170 cm tall, and then keep those called Chris, and then select their age and weight.

> trumpton %>%
+   filter(Height >= 170) %>%
+   filter(FirstName == "Chris") %>%
+   select(Age, Weight)
# A tibble: 1 × 2
    Age Weight
  <dbl>  <dbl>
1    26     90

Writing each step on a line of its own makes a long command much easier to read. Put the %>% at the end of each line rather than the start of the next, so that R knows the command carries on. In RStudio, Ctrl+Shift+M (⌘+Shift+M on a Mac) types a pipe for you. Depending on RStudio’s settings (Tools ▸ Global Options ▸ Code), it may type |> instead: that is R’s own pipe, which works in the same way for everything in this course.

Several tests at once

You can also give filter() several tests, separated by commas, and it keeps the rows for which all of them are TRUE. These two commands are the same as the two filters above:

trumpton %>% filter(Height >= 170, FirstName == "Chris")
trumpton %>% filter(Height >= 170 & FirstName == "Chris")

& means “and”, | means “or”, and ! means “not”. %in% tests whether a value is one of a set of values, so filter(trumpton, FirstName %in% c("Chris", "Liam")) keeps the people called Chris or Liam, and filter(trumpton, !FirstName %in% c("Chris", "Liam")) keeps everyone else.

Multi-stage filtering

From child, extract the variants with a quality of 200, a coverage of more than 50, and a MutantReadPercent of more than 70.
Remove all the variants on the X, Y and MT chromosomes.
Show just the chromosome and position (CHR and POS) of the variants that have an ID in dbSNP. (Variants without an ID have a full stop, ., in the dbSNP column.)

How many variants pass all three quality tests?

3961, from child %>% filter(QUAL == 200, COVERAGE > 50, MutantReadPercent > 70).

How many of the variants have an ID in dbSNP?

23690, from child %>% filter(dbSNP != ".") %>% select(CHR, POS).

What does small %>% filter(Length > 80) mean?

The same as filter(small, Length > 80) Keep the rows of small whose Length is more than 80 per cent of the largest Store the filtered rows in small

If you have time Transforming data in a filter

This part is optional: it doesn’t count towards your certificate.

In R basics you used nchar(), which counts the characters in each piece of text, and substr(), which extracts part of a piece of text. Both work inside a filter.

Transform and filter

Use nchar() on the REF column to select the variants that are deletions: those where the reference sequence is more than one base long.
Use substr() to find all the variants in genes whose names start with Q.

How many deletions are there, and how many variants in genes whose names start with Q?

602 deletions, from child %>% filter(nchar(REF) > 1), and 30 variants in genes starting with Q, from child %>% filter(substr(GENE, 1, 1) == "Q").

Your first plots with ggplot2

≈ 25 min

ggplot2 is the tidyverse’s package for plots. It is powerful and flexible, and it works in the same way as the rest of the tidyverse: the data come in a tibble, in tidy format, as the first argument.

How a ggplot is put together

Every plot in ggplot2 is built in the same way:

Start with ggplot(), and give it the tibble that holds your data.
Say which columns of the data should control which parts of the plot, with aes(). These are the aesthetics: the position along the x axis, the position along the y axis, the colour, the size, and so on.
Add a geometry with +: the kind of plot you want, such as points, lines or bars.
Add anything else you need, such as titles or labels, with more +.

The geometries all have names that start with geom_:

Geometry	What it draws
`geom_point()`	points: scatter plots and strip charts
`geom_line()`	line graphs
`geom_boxplot()`	box plots
`geom_col()` and `geom_bar()`	bar plots
`geom_histogram()` and `geom_density()`	distributions

These examples use a tibble of gene expression levels, from expr_data.txt. To try them yourself, read it in first: expression <- read_tsv("~/data/expr_data.txt")

> expression
# A tibble: 12 × 4
   Gene       WT     KO pValue
   <chr>   <dbl>  <dbl>  <dbl>
 1 Mia1     5.83  3.24  0.1
 2 Snrpa    8.59  5.02  0.001
 3 Itpkc    8.49  6.16  0.04
 4 Adck4    7.69  6.41  0.2
 5 Numbl    8.37  6.81  0.1
 6 Ltbp4    6.96 10.4   0.001
 7 Shkbp1   7.57  5.83  0.1
 8 Spnb4   10.7   9.38  0.2
 9 Blvrb    7.32  5.29  0.05
10 Pgam1    0     0.285 0.5
11 Sertad3  8.13  3.02  0.0001
12 Sertad1  7.69  4.34  0.01

To plot the expression in the KO samples against the expression in the WT samples, the tibble is expression, the aesthetics are x = WT and y = KO, and the geometry is points:

ggplot(expression, aes(x = WT, y = KO)) + geom_point()

Change geom_point() to geom_line(), and the same data are drawn as a line instead. Because the tibble is the first argument of ggplot(), you can also pipe the data into it.

Fixed and mapped aesthetics

An aesthetic can be set in two ways. A fixed value is the same for every point: all the points red, say, or all of them size 4. You give fixed values to the geometry, outside aes():

expression %>%
  ggplot(aes(x = WT, y = KO)) +
  geom_point(colour = "red2", size = 5)

A mapped aesthetic comes from a column of your data: points coloured by genotype, or sized by the number of observations. Mappings go inside aes(), as in aes(x = WT, y = KO, colour = genotype), and ggplot2 adds a legend to explain them. A mapping can be categorical (a column of text, which gives a separate colour for each value) or quantitative (a column of numbers, which gives a scale of colours from low to high).

The + goes at the end of the line

Like the pipe, the + that adds a layer must go at the end of a line, not at the start of the next one. If a line ends without one, R thinks the command has finished there: it draws the plot without the layer, and then complains about the + on the next line.

Your first ggplot

Read brain_bodyweight.txt, a tab-separated file, into a tibble called brain. It holds the weights of the brains and bodies of 27 species of animal.
Draw a scatter plot (with geom_point()) of the weights of the brains and bodies, with brain on the x axis and body on the y axis.
Make all the points blue2 (colour = "blue2") and size 3.
The points are all squashed into the bottom left-hand corner of the plot. Use the log.brain and log.body columns instead, and see whether the plot looks better.
Change the colour of the points from a fixed blue to a mapping, so that they are coloured by Category.

What do you notice about the extinct species?

Their bodies are far heavier, for the size of their brains, than those of the other animals, so they sit well above the rest They follow the same trend as the other animals They had the largest brains of all

What is the difference between geom_point(colour = "red") and ggplot(brain, aes(x = log.brain, y = log.body, colour = Category))?

The first makes every point red; the second colours each point according to the value in its Category column There is no difference: both set the colour of the points The first is for points and the second is for lines

If you have time Leave out the extinct species

This part is optional: it doesn’t count towards your certificate.

Draw your last plot again, but without the extinct species. Use what you learnt about filter() to remove them from the data, and pipe the filtered data straight into ggplot().

Filter, then plot

Plot the data without the extinct species.

Bar plots and distributions

≈ 25 min

Points are only the start. ggplot2 has geometries for bar plots, histograms, density plots and many more, and they all work in the same way.

Bar plots: geom_col() and geom_bar()

There are two geometries for bar plots, and which one you need depends on your data. Use geom_col() when your data already contain the heights of the bars. Then you need two aesthetics: x for the bars, and y for their heights. To plot the expression of each gene in the WT samples:

expression %>%
  ggplot(aes(x = Gene, y = WT)) +
  geom_col(fill = "red2")

Use geom_bar() when you want ggplot2 to count things for you: how many rows there are for each value of x. It only needs the x aesthetic. For example, ggplot(small, aes(x = Category)) + geom_bar() draws a bar for each category, as tall as the number of samples in it.

Distributions: histograms and density plots

A histogram divides the range of the values into bins, and draws a bar for each bin, as tall as the number of values in it. The argument binwidth sets how wide the bins are. These examples use many_values.txt, which has 100,000 values from two genotypes:

many.values <- read_tsv("many_values.txt", show_col_types = FALSE)
many.values %>%
  ggplot(aes(x = values)) +
  geom_histogram(binwidth = 0.1, fill = "yellow", colour = "black")

A density plot (geom_density()) draws a smooth curve instead, which is often easier to read when you are comparing groups. Here the fill is mapped to the genotype, and alpha, the opacity, is set to 0.5, so that each curve shows through the other:

many.values %>%
  ggplot(aes(x = values, fill = genotype)) +
  geom_density(colour = "black", alpha = 0.5)

fill and colour

For a shape with an inside, such as a bar or the area under a density curve, fill sets the colour of the inside and colour sets the colour of its outline. For points and lines, colour is the colour of the point or line itself. The help page of each geometry (for example ?geom_histogram) lists all the aesthetics it understands.

More plots

Draw a bar plot of the lengths of the samples in category A of small. Filter the data to keep only the category A samples, pass the filtered tibble to ggplot(), use Sample for x and Length for y, and use geom_col(), because the data give the heights of the bars.
Plot the distribution of the MutantReadPercent values of child as a histogram, and change its fill and colour to make it look better.
Try different values of binwidth in geom_histogram() to change how finely the distribution is divided up.
Plot the same distribution as a density plot.
Draw a bar plot of how often each base is the reference base of a variant. Filter child to keep only the variants whose REF is a single letter (use nchar()), pass them to ggplot() with REF as the x aesthetic, and use geom_bar(): there is no y aesthetic, so ggplot2 counts the variants for each base.

Which base is most often the reference base of a single-base variant?

G, with 6982 variants (C has 6813, T 5718 and A 5707). child %>% filter(nchar(REF) == 1) %>% count(REF) gives you the numbers.

You have a tibble with one row for each sample, and you want a bar for each category, as tall as the number of samples in it. Which geometry do you need?

geom_bar(), with Category as x geom_col(), with Category as x geom_histogram(), with Category as x

What happens to a histogram when you make binwidth smaller?

There are more, narrower bars, which show more detail but are more ragged There are fewer, wider bars The bars get shorter, but there are just as many

Titles, labels and layers

≈ 15 min

A plot that someone else will read needs a title and clear labels. ggplot2 builds a plot up in layers, so you add these, and other annotations, with +.

ggtitle() adds a title to a plot, and xlab() and ylab() change the labels of its axes. A plot can also have more than one geometry, each of them a layer drawn on top of the ones before. geom_text() writes text at each point, taken from the label aesthetic. Here the names of the genes are written just below the points, because vjust = 1.2 moves each label down a little:

expression %>%
  ggplot(aes(x = WT, y = KO, label = Gene)) +
  geom_point() +
  ggtitle("Expression level comparison") +
  xlab("WT expression level (log2 RPM)") +
  ylab("KO expression level (log2 RPM)") +
  geom_text(vjust = 1.2)

The arguments hjust and vjust move the labels sideways and up or down. Try a few values to see what they do.

Label your plot

Go back to your plot of log.body against log.brain, coloured by Category, and give it a title and proper labels for its axes.
Add a geom_text() layer that writes the name of each species next to its point. It needs a mapping of label to Species, which you can add to the aes() in ggplot() or give to geom_text() in an aes() of its own. Then use hjust to move the names off the points.

Why does this code draw a plot with no points, and then give an error?

ggplot(brain, aes(x = log.brain, y = log.body))
  + geom_point()

The + is at the start of the second line, so R thinks the command ended with the first line geom_point() needs its own aes() The data should be piped into ggplot()

If you have time A line of best fit

This part is optional: it doesn’t count towards your certificate.

The plot shows that animals with heavier brains tend to have heavier bodies. You can test whether this relationship is significant, and measure how strong it is, with two functions that ask similar questions of the same data, but which are used in quite different ways. That is a good reminder to read the help page of a function before you use it.

cor.test() takes two vectors of the same length and tests the correlation between them. It gives you the correlation coefficient and a p-value. Use $ to take the log.brain and log.body columns out of brain, and note the p-value, so that you can add it to the title of your plot.

lm() fits a linear model. Its two main arguments are a formula, which describes the relationship you want to model, and the data. The formula log.body ~ log.brain means “log.body as predicted by log.brain”, so the command is lm(log.body ~ log.brain, data = brain). It gives you the intercept and the slope of the line of best fit, which you can draw on your plot with geom_abline(intercept = …, slope = …).

Test, and draw the line

Test the correlation between log.brain and log.body with cor.test().
Fit a linear model of log.body against log.brain with lm().
Add the line of best fit to your plot with geom_abline(), and put the p-value in the title.

What p-value does cor.test() give for the correlation?

2.44e-06, which is 0.00000244: a highly significant correlation. The correlation coefficient is 0.77, and lm() gives an intercept of −2.283 and a slope of 1.215.