session7.Rmd

---
title: "Session 7 -- Restructuring data for analysis"
---

> #### Learning objectives
>
> * Understand what makes a data set 'tidy' and why you'd want your data to be structured this way
> * Use `pivot_longer()` and `pivot_wider()` operations to restructure data frames
> * Tease apart columns containing multiple variables using `separate()`
> * Modify character variables using string manipulation functions from the stringr package
> * Customize non-data components of plots created using ggplot2 by changing the theme

---

# Restructuring data

The data you collect or obtain from a third party is not always in a form that
is suitable for exploratory analysis and visualization and may need to be
restructured before you can fully make use of it.

This is particularly true of the plotting and summarizing tools we've been
looking at in this course, which are designed specifically to work on data in a
format referred to as 'tidy'. This is where the tidyverse gets its name.

In this session, we will look at what it means for data to be 'tidy' and how
you can transform your data, if necessary, into this form. We'll also look at
useful functions for handling compound variables, that is columns that contain
more than one type of measurement or attribute (you'd be surprised how common
this is) and some of the string manipulation functions from the stringr
package that can help with cleaning up values within a data frame.

Finally, we'll take another look at customizing plots created with ggplot2 by
changing various non-data components that are largely cosmetic.

The functions we're mostly focusing on in this session are from the **tidyr**
and **stringr** packages both of which get loaded as part of the tidyverse.

```{r}
library(tidyverse)
```

---

# Tidy data

So what is 'tidy data' and why should we care about it?

To answer these questions, we'll look at different ways in which a simple data
set can be represented and consider the challenges associated with each. The
data set in question is a subset of data from the
[World Health Organization Global Tuberculosis Report](https://www.who.int/health-topics/tuberculosis).

The tidyr package contains a series of tibbles that represent the same set of
information on the number of new cases of tuberculosis (TB) recorded in each of
3 countries in the years 1999 and 2000 as well as the populations of those
countries.

Here is the first table, `table1`, that contains a row for every combination of
country and year and separate columns containing the numbers of TB cases and the
population.

```{r}
table1
```

Two alternative ways of representing the same information are given in `table2`
and `table3`. We'll consider each of these in turn, shortly.

```{r}
table2
table3
```

The final representation has the data split across two tables, a scenario that
is actually quite likely given that population data will almost certainly have
been collected separately from the recording of TB cases.

```{r}
table4a
table4b
```

Time series data like this are very commonly represented as they are in
`table4a` and `table4b` with a series of dates or years as columns extending
across a spreadsheet. You will find numerous examples of this if you seek out
various data sets made available by the UK Office for National Statistics and
various other national and international organizations.

Tables 1 to 4 are all different representations of the same underlying data but
one of these tables is structured in such a way as to be most readily used in
the tidyverse.

## Rules for tidy data

```{block type = "rmdblock"}
**Tidy data**

A tidy data set is a data frame (or table) for which the following are true:

1. Each **variable** has its own column
2. Each **observation** has its own row
3. Each **value** has its own cell

![](images/tidy_data.png)

A **variable** contains all values that measure the same underlying attribute
(like height, temperature, duration) across units.

An **observation** contains all values measured on the same unit (like a person
or a day) across attributes.

_From ['Tidy Data'](https://www.jstatsoft.org/index.php/jss/article/view/v059i10/v59i10.pdf) by Hadley Wickham._
```

**Question:** _Which of the above representations of TB cases is tidy?_

Another way of framing the question is to consider what are the variables in the
TB data set, i.e. what are the things that vary and for which we can attach a
value for each observation?

Take another look at tables 4a and 4b. Do each of the columns correspond to a
variable? The country is certainly a variable. In this data set it takes one of
three values: Afghanistan, Brazil or China.

But what about the other two columns, '1999' and '2000'? These are values, not
variables. The variable in this case would be 'year' and could take a value of
1999 or 2000. Tables 4a and 4b are not tidy.

There is also another rather obvious problem with tables 4a and 4b -- the data
are contained in two separate data frames. The data would almost certainly have
been collected separately, so it's hardly surprising, but whenever numbers of
people affected by a disease, or engaging in an activity, are compared between
countries we almost always want to be comparing the rate (the percentage within
the population) and not the absolute number, so that the comparison is fair. We
need to combine the data from these two tables in order to calculate the rate.

The only table that is truly tidy is `table1`. It contains one column for each
of the variables in the data set, namely country, year, the number of new TB
cases and the population. We'll look at tables 2 and 3 shortly and why these
aren't tidy and what we can do about it, but first we'll see how we can convert
tables 4a and 4b into the tidy format.

----

# Pivoting operations

## `pivot_longer()`

Tables 4a and 4b are examples of what is often referred to as 'wide format'.
While neither table looks especially wide, you can imagine that the more complete
WHO data set contains data for very many years and that if each had its own
column, the table would be very much wider.

What we need is a column for 'year' so that we have a count, whether it is the
number of TB cases or the population, for each unique combination of country and
year. Transforming the table in this way is known as 'pivoting' and the tidyr
package provides the `pivot_longer()` function for just such an operation.

```{r}
table4a
table4a_long <- pivot_longer(table4a,
                             c(`1999`, `2000`),
                             names_to = "year",
                             values_to = "cases")
table4a_long
```

As with almost all tidyverse operations the first argument is the data frame
we're working on. The second specifies which columns we are operating on. Here
we've used a vector with `c()` but it is also quite customary and normally more
convenient to use a range of columns, e.g. `` `1999`:`2000` ``. Remember that we
have to use backticks because R doesn't like variable names starting with
numbers.

```{block type = "rmdblock"}
**Pivoting operations**

**`pivot_longer(data, cols, names_to = "name", values_to = "value")`**

Pivot data from wide to long, increasing the number of rows and decreasing the
number of columns.

`table4a_long <- pivot_longer(table4a, c("1999", "2000"), names_to = "year", values_to = "cases")`

![](images/pivot_longer.png)

**`pivot_wider(data, names_from = name, values_from = value)`**

Pivot data from long to wide, increasing the number of columns and decreasing
the number of rows.

`table4a_wide <- pivot_wider(table4a_long, names_from = year, values_from = cases)`
```

The `names_to` and `values_to` arguments are so called because we are taking
the names of the columns we specified (1999 and 2000) and putting these in a new
column with a name given by `names_to`. Likewise, we are taking values in each
of our specified columns and putting these in a new column whose name is given
by `values_to`.

We can do the same with `table4b` and join the two resulting tables together to
recreate `table1`.

```{r}
table4b_long <- pivot_longer(table4b,
                             -country,
                             names_to = "year",
                             values_to = "population")
table4b_long
# join the two pivoted tables to recreate table1
left_join(table4a_long, table4b_long, by = c("country", "year"))
```

## `pivot_wider()`

`pivot_wider()` has the opposite effect of `pivot_longer()` and can be used
to reverse the pivoting operation we just performed on `table4a`.

```{r}
table4a_long
pivot_wider(table4a_long, names_from = "year", values_from = "cases")
```

`pivot_wider()` has looked for all distinct values in the `year` column and
created a new column for each. The values are taken from the `cases` column.
An `NA` would be inserted into the new columns for any instances where there
is no observation for a given year. That isn't the case in this example but we
can show this by removing the last row from `table4a_long`.

```{r}
table4a_long <- head(table4a_long, 5)
table4a_long
pivot_wider(table4a_long, names_from = "year", values_from = "cases")
```

In some cases the wide format is the more human-readable form and usually it is
a more compact way of representing the data. In this example, there is no
duplication of the year value in the wide format of `table4a` and `table4b`. We
will see later that this is much more apparent with larger data tables, e.g.
gene expression matrices.

Let's look again at `table2`.

```{r}
table2
```

Are the type and count columns true variables? And what is the observational
unit in this case?

If we consider the observational unit to be a country in a specific year then
`table2` is not tidy because observations are split across two rows. Also the
count variable contains counts of what are essentially different things, the
number of cases of TB or the total population.

In tricky situations like this, a tell-tale sign that the data is not in a tidy
format is when we attempt to perform some exploratory data analysis and
visualization and find we're having to do quite a bit of work to calculate or
plot what we want.

For example, if we wanted to create a bar plot of the numbers of TB cases in
each country, we would have to first remove the rows corresponding to the
populations using a filter operation. Surely that wouldn't be necessary if the
data were tidy.

Similarly, the rate of new TB cases, i.e. the proportion of the population
infected with TB, is something we should be able to calculate easily in a
simple operation. However, this is actually quite difficult to do with the data
structured as they are in `table2`.

We can use `pivot_wider()` to sort this out. The `type` column contains the
variable names so we'd need to set `names_from = "type"`, while the values will
be taken from the `count` column.

```{r}
table2_fixed <- pivot_wider(table2, names_from = "type", values_from = "count")
table2_fixed
```

The resulting table is exactly the same as `table1` and now the rate of
infection can be calculated rather straightforwardly.

```{r}
mutate(table2_fixed, rate = cases / population)
```

----

# Splitting columns

Table 3 contains an example of a column that contains multiple values. It is a
somewhat convoluted example but occasionally you may come across data like this.

```{r}
table3
```

The `rate` column contains both the number of TB cases and the population
separated by a '/' character. The rate column is a character type so not
terribly useful for doing anything of a mathematical nature in its current
guise.

## `separate_wider_delim()`

The `separate_wider_delim()` function allows us to split a character column into
multiple columns based on a delimiter or separator.

```{r}
table3_separated <- separate_wider_delim(table3,
                                         rate,
                                         names = c("cases", "population"),
                                         delim = "/")
table3_separated
```

The `delim` argument takes a vector of characters to split the string at.  In
this case our separator is just the '/' character.

The resulting data frame is still not quite what we want though. This becomes
apparent as soon as we try to do anything with the new `cases` and `population`
columns.

```{r error = TRUE}
mutate(table3_separated, rate = cases / population)
```

By default the separated values are character types. We could convert these
using `mutate()`.

```{r}
table3_separated <- separate_wider_delim(table3,
                                         rate,
                                         names = c("cases", "population"),
                                         delim = "/") %>%
    mutate(across(c(cases, population), as.integer)) %>%
    mutate(rate = cases / population)
table3_separated
```

In addition to `separate_wider_delim()`, the tidyr package provides a two 
other functions for splitting columns:

* `separate_wider_regex()` - for splitting columns based on a regular expression
* `separate_wider_fixed()` - for splitting columns based on a fixed number of
                             characters

---

# Example 1: METABRIC gene expression

Although tables 1 to 4 contain real data they are, of course, 'toy' data frames
created for demonstration and teaching purposes. We'll now turn our attention to
the METABRIC expression data and see how this needs to be transformed into a
tidier format to open up different avenues for exploring the data.

We'll first load the table and then select just the columns we're going to need.

```{r message = FALSE}
metabric <- read_csv("data/metabric_clinical_and_expression_data.csv") %>%
    select(Patient_ID, ER_status, ESR1:MLPH)
metabric
```

When we first looked at visualization using ggplot2 we created the following
box plot.

```{r box_plot_1, fig.width = 5, fig.height = 4}
ggplot(metabric) +
    geom_boxplot(mapping = aes(x = ER_status, y = GATA3))
```

But what if we would like to create a series of plots using the faceting
functions in ggplot2 with one plot for each gene?

Faceting requires a categorical variable, which is used to divide the data into
subsets to be used for each plot. In this case we'd need a gene column. Clearly
our data are not structured in this way.

We have gene names for column headings. Are these variables? Well, maybe,
although a more correct name for each of these variables or column headings
would be 'Expression of ESR1', 'Expression of ERBB2', etc.

But we could consider that these gene column headings are actually values of a
variable called 'gene' or 'gene symbol'. In this regard, what we have is a 'wide
format' table.

Most gene expression matrices have a similar form, although usually there have
rows for each gene and columns for each sample. It should be said that the gene
expression matrix format is a very compact way of representing the data which
could be a consideration when dealing with tens of thousands of genes and
anywhere between a few tens of samples to a few thousand, such is the case for
METABRIC.

Furthermore, there are lots of tools for working with gene expression data in
the form of matrices, including many packages in the Bioconductor project.
Fortunately, as we've seen, `pivot_longer()` and `pivot_wider()` provide a very
convenient means of converting between tidy and matrix-like formats.

We'll convert our table of ER status and gene expression data to the tidy
format.

```{r}
metabric <- pivot_longer(metabric,
                         ESR1:MLPH,
                         names_to = "Gene",
                         values_to = "Expression")
metabric
```

Note how we specified a range of columns between `ESR1` and `MLPH`, which is a
lot easier than naming each column individually.

We're now in a position to create our faceted box plot chart.

```{r box_plot_2, fig.height = 6}
ggplot(data = metabric,
       mapping = aes(x = ER_status, y = Expression)) +
    geom_boxplot() +
    facet_wrap(vars(Gene))
```

In carrying out this transformation, the observational unit has changed. The
tidy format has _one-row-per-gene-per-sample_, while wide format was
_one-row-per-sample_. The tidy format is much less compact and involves
considerable duplication of values in the first three columns (`Patient_ID`,
`ER_status` and `Gene`).

One of the other plot types we've used in exploring these data was a scatter
plot comparing the expression of two genes across all the samples. For this,
the _one-row-per-sample_ representation is the more appropriate and being able
to convert back to this format allows us to create the plot.

```{r scatter_plot_1}
metabric %>%
    pivot_wider(names_from = "Gene", values_from = "Expression") %>%
    ggplot() +
    geom_point(mapping = aes(x = GATA3, y = ESR1, colour = ER_status))
```

---

# Example 2: Protein levels in MCF-7 after treatment with tamoxifen

Our second real example features another data set generated by CRUK CI
scientists
([Papachristou _et al._, Nature Communications 9:2311, 2018](https://pubmed.ncbi.nlm.nih.gov/29899353))
in which the dynamics of endogenous chromatin-associated protein complexes were
investigated using quantitative mass spectrometry.

We'll look at just one of several tabular data sets made available as
supplementary data, which contains the total level of protein in MCF-7 cells at
various time points after treatment with the drug tamoxifen. MCF-7 is a breast
cancer cell line isolated in 1970 from a 69-year old woman and established for
use in breast cancer research by the Michigan Cancer Foundation-7 institute in
Detroit. Tamoxifen is a hormone therapy used in the treatment of estrogen
receptor-positive breast cancer.

The table in question is
[supplementary data 9](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5998130/bin/41467_2018_4619_MOESM11_ESM.xlsx).

```{r message = FALSE}
library(readxl)
protein_levels <- read_excel("data/41467_2018_4619_MOESM11_ESM.xlsx", skip = 1)
select(protein_levels, `Uniprot Accession`, `Gene Name`, ends_with("rep01"))
```

This has a very similar structure to a gene expression matrix having one row for
each protein (or gene) and a column for each sample. We've only shown columns
for the first of four replicates for each group defined by a treatment and a
time point. The control group is, in effect, the 'untreated' group in which the
cells are treated with the vehicle (ethanol) alone.

If we wanted to calculate the mean protein levels within each group, i.e. the
average level for the protein measured in the 4 replicates, or we wanted to
show the spread of values for the replicates as a box plot, then the data as
currently structured is not in the most suitable form. In what follows, we'll
transform the table to allow us to do both these analyses.

To simplify matters, we're going to focus on just a few proteins, those whose
levels are markedly reduced 24 hours after treatment with tamoxifen compared
with the vehicle.

```{r}
protein_levels <- protein_levels %>%
    filter(`log2FC(24h/veh)` < -0.75) %>%
    select(accession = `Uniprot Accession`,
           gene = `Gene Name`,
           vehicle.rep01:tam.24h.rep04)
protein_levels
```

This is a fairly typical example of a table with columns for each sample where
the sample names contain quite a lot of information, in this case:

* the treatment
* the time at which the protein levels are measured after treatment
* the number of the replicate sample for that treatment and time point

To make use of this information we need to pivot the table such that the sample
name is in a column and protein levels are in another column, and then to split
the sample name column into its constituent parts.

For the first part of this transformation we'll use `pivot_longer()`.

```{r}
protein_levels <- pivot_longer(protein_levels,
                               vehicle.rep01:tam.24h.rep04,
                               names_to = "sample",
                               values_to = "protein_level")
protein_levels
```

Now we can use `separate()` to disentangle the components of the sample name.

It looks like we need to split on the '.' character but that has a special
meaning in a regular expression, i.e. match _any_ character. Fortunately,
the default regular expression used by `separate()` splits on any character
that isn't a letter or a number, so will do just fine.

But there is another problem. The vehicle sample names don't follow the
pattern of "treatment.time.replicate". In actual fact the vehicle measurements
were taken 24 hours after treatment with ethanol alone. What we should do to
correct matters is add the '24h' component. For that, we're going to use one
of several really useful string manipulation functions provided by the
stringr package.

```{block type = "rmdblock"}
**Some useful stringr functions**

**`str_c(..., sep = "")`**

Join multiple strings into a single string.

`str_c("tidyverse", "fab", sep = " is ")`

` [1] "tidyverse is fab"`

`str_c(c("Ashley", "Matt"), c("hiking", "swimming"), sep = " loves ")`

` [1] "Ashley loves hiking" "Matt loves swimming"`

---

**`str_replace(string, pattern, replacement)`**

**`str_replace_all(string, pattern, replacement)`**

Substitute a matched pattern in a string (or character vector).

`str_replace("Oscar is the best cat in the world", "best", "loveliest")`

` [1] "Oscar is the loveliest cat in the world"`

`str_replace_all("the cat sat on the mat", "at", "AT")`

` [1] "the cAT sAT on the mAT"`

---

**`str_remove(string, pattern)`**

**`str_remove_all(string, pattern)`**

Remove matched patterns in a string.

Equivalent to `str_replace(string, pattern, "")` and
`str_replace_all(string, pattern, "")`.

---

**`str_sub(string, start, end)`**

Extract substrings from a character vector at the given start and end positions.

`str_sub("Ashley", start = 1, end = 3)`

` [1] "Ash"`

`str_sub(c("tamoxifen", "vehicle"), 1, 3)`

` [1] "tam" "veh"`

```

`str_replace()` looks like the function we need so let's see how it works on our
sample names.

```{r}
str_replace("vehicle.rep01", "vehicle", "vehicle.24h")
```

The first argument is the character vector we're working on, in this case a
single character string. The second argument is the pattern or substring we want
to replace and the third is the string we want to replace it with.

Looking at the help we can see that, like very many R functions, it works on a
vector of character values (or strings), so let's try this on a few of our
sample names, say the first 10.

```{r}
str_replace(protein_levels$sample[1:10], "vehicle", "vehicle.24h")
```

Finally, we can modify the values in our data frame using `mutate()`.

```{r}
protein_levels <- mutate(protein_levels,
                         sample = str_replace(sample, "vehicle", "vehicle.24h"))
protein_levels
```

Now we're ready to separate the `sample` column into its component parts.

```{r}
protein_levels <- separate_wider_delim(protein_levels,
                                       sample,
                                       names = c("treatment", "time", "replicate"),
                                       delim = ".")
protein_levels
```

The groups we want to compare are in fact the combination of the treatment and
the time point. We'll create a new column called 'group' in which we
concatenate the treatment and time values using another stringr function,
`str_c()`. But we'll also shorten 'vehicle' to 'veh' so we have similar length
treatment labels ('veh' and 'tam') in the plot we'll eventually create. For that
we'll use `str_sub()`. Our group variable is categorical so we'll convert this
to a factor while we're at it.

```{r}
protein_levels <- protein_levels %>%
    mutate(group = str_sub(treatment, 1, 3)) %>%
    mutate(group = str_c(group, time)) %>%
    mutate(group = as_factor(group))
protein_levels
```

Computing the mean protein levels within each group is now very straightforward.

```{r}
protein_levels %>%
    group_by(accession, gene, group) %>%
    summarize(protein_level = mean(protein_level))
```

Plotting the protein levels in each group as box plots with each data point
overlaid is only slightly more involved. We can use faceting to create separate
plots for each of the proteins (we'll use gene symbols as they're more
recognizable).

```{r box_plot_3, fig.height = 7}
ggplot(data = protein_levels, mapping = aes(x = group, y = protein_level)) +
    geom_boxplot(colour = "grey60") +
    geom_jitter(width = 0.15) +
    facet_wrap(vars(gene), scales = "free_y", ncol = 3)
```

---

# Further customization of plots with ggplot2

Finally, we'll turn our attention back to visualization using ggplot2 and how we
can customize our plots by changing the theme and consider various other odds
and ends.

## Themes

Themes can be used to customize non-data components of a plot. Let's create a
plot showing the expression of estrogen receptor alpha (ESR1) for each of the
Integrative cluster breast cancer subtypes.

```{r message = FALSE}
# read in the METABRIC data, convert the Integrative_cluster variable into a
# categorical variable with the levels in the correct order, and select just
# the columns and rows we're going to use
cluster_order <- c("1", "2", "3", "4ER-", "4ER+", "5", "6", "7", "8", "9", "10")
metabric <- read_csv("data/metabric_clinical_and_expression_data.csv") %>%
    mutate(Integrative_cluster = factor(Integrative_cluster,
                                        levels = cluster_order)) %>%
    mutate(`3-gene_classifier` = replace_na(`3-gene_classifier`,
                                            "Unclassified")) %>%
    select(Patient_ID,
           ER_status,
           PR_status,
           `3-gene_classifier`,
           Integrative_cluster,
           ESR1) %>%
    filter(!is.na(Integrative_cluster), !is.na(ESR1))
metabric
```

```{r box_plot_4}
# plot the ESR1 expression for each integrative cluster
plot <- ggplot(data = metabric) +
    geom_boxplot(mapping = aes(x = Integrative_cluster, y = ESR1, fill = Integrative_cluster)) +
    labs(x = "Integrative cluster", y = "ESR1 expression")
plot
```

The default theme has the characteristic grey background which isn't
particularly suitable for printing on paper. We can change to one of a number
of alternative themes available in the ggplot2 package, e.g. the black and
white theme.

```{r box_plot_5}
plot + theme_bw()
```

The help page for the available themes, which can be accessed using `?theme_bw`,
lists each and when you might want to use them. It states that `theme_bw` may
work better for presentations displayed on a projector.

Each of these themes is really just a collection of attributes relating to how
various non-data elements of the plot will be displayed. We can override any of
these individual settings using the `theme()` function. A look at the help page
(`?theme`) shows that there are a very large number of settings that you can
change. The following example demonstrates a few of these.

```{r box_plot_6}
plot +
  theme_bw() +
  theme(
    panel.grid.major.x = element_blank(),
    axis.ticks.x = element_blank(),
    legend.position = "none"
  )
```

Here's another example that also involves customizing the labels, scales and
colours.

```{r bar_plot_1, fig.width = 5}
ggplot(data = metabric) +
    geom_bar(mapping = aes(x = `3-gene_classifier`, fill = ER_status)) +
    scale_y_continuous(limits = c(0, 700),
                       breaks = seq(0, 700, 100),
                       expand = expand_scale(mult = 0)) +
    scale_fill_manual(values = c("firebrick2", "dodgerblue2")) +
    labs(x = NULL, y = "samples", fill = "ER status") +
    theme_bw() +
    theme(panel.border = element_blank(),
          panel.grid = element_blank(),
          axis.ticks.x = element_blank(),
          axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
          axis.line.y = element_line(),
          axis.ticks.length.y = unit(0.2, "cm"),
          legend.position = "bottom")
```

The ggthemes package contains some extra themes and might be fun to check out.
Here's an example of a plot that uses the `theme_gdocs` theme that resembles the
default look of charts in Google Docs.

```{r box_plot_7, fig.width = 5, fig.height = 4}
library(ggthemes)
metabric %>%
    filter(`3-gene_classifier` == "HER2+") %>%
    ggplot(mapping = aes(x = PR_status, y = ESR1)) +
    geom_boxplot() +
    geom_jitter(mapping = aes(colour = PR_status),
                width = 0.25,
                alpha = 0.4,
                show.legend = FALSE) +
    scale_colour_brewer(palette = "Set1") +
    labs(x = "PR status", y = "ESR1 expression") +
    theme_gdocs()
```

## Position adjustments

All geoms in ggplot2 have a position adjustment that can be set using the
`position` argument. This has different effects for different types of plot but
essentially this resolves how overlapping geoms are displayed.

For example, let's consider the stacked bar plot we created earlier showing the
numbers of patients in each of the 3-gene classifications subdivided by ER
status. The default position value for `geom_bar()` is "stack" which is why the
plot is shown as a stacked bar chart. An alternative way of representing these
data would be to show separate bars for each ER status side-by-side by setting
`position = "dodge"`.

```{r bar_plot_2}
ggplot(data = metabric) +
    geom_bar(mapping = aes(x = `3-gene_classifier`, fill = ER_status),
             position = "dodge") +
    scale_y_continuous(limits = c(0, 700),
                       breaks = seq(0, 700, 100),
                       expand = expand_scale(mult = 0)) +
    scale_fill_manual(values = c("firebrick2", "dodgerblue2")) +
    labs(x = NULL, y = "samples", fill = "ER status") +
    theme_bw() +
    theme(panel.border = element_blank(),
          panel.grid = element_blank(),
          axis.ticks.x = element_blank(),
          axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
          axis.line.y = element_line(),
          axis.ticks.length.y = unit(0.2, "cm"))
```

Another position adjustment we've come across is `geom_jitter()`, which is just
a convenient shortcut for `geom_point(position = "jitter")`. A variation on this,
`position_jitterdodge()`, comes in handy when we are overlaying points on top of
a box plot. We show an example of just such a plot in which first use
`postion = "jitter"`.

```{r box_plot_8}
ggplot(data = metabric,
       mapping = aes(x = `3-gene_classifier`, y = ESR1, colour = PR_status)) +
    geom_boxplot() +
    geom_point(position = "jitter", size = 0.5, alpha = 0.3) +
    labs(x = "3-gene classification",
         y = "ESR1 expression",
         colour = "PR status") +
    scale_color_brewer(palette = "Set1") +
    theme_minimal() +
    theme(panel.grid.major.x = element_blank(),
          axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
          axis.ticks.x = element_blank())
```

The PR-negative and PR-positive points have distinct colours but are overlapping
in a way that is aesthetically displeasing. What we want is for the points to
have both jitter and to be dodged in the same way as the boxes. With
`position_jitterdodge()` we get a better looking plot.

```{r box_plot_9}
ggplot(data = metabric,
       mapping = aes(x = `3-gene_classifier`, y = ESR1, colour = PR_status)) +
    geom_boxplot() +
    geom_point(position = position_jitterdodge(), size = 0.5, alpha = 0.3) +
    labs(x = "3-gene classification",
         y = "ESR1 expression",
         colour = "PR status") +
    scale_color_brewer(palette = "Set1") +
    theme_minimal() +
    theme(panel.grid.major.x = element_blank(),
          axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
          axis.ticks.x = element_blank())
```

## Saving plot images

Use `ggsave()` to save the last plot you displayed.

```{r eval = FALSE}
ggsave("esr1_expression.png")
```

You can alter the width and height of the plot and can change the image file
type.

```{r eval = FALSE}
ggsave("esr1_expression.pdf", width = 20, height = 12, units = "cm")
```

You can also pass in a plot object you have created instead of using the last
plot displayed. See the help page (`?ggsave`) for more details.

---

# Summary

In this session we have covered the following:

* Restructuring data frames using pivoting operations
* Splitting compound variables into separate columns
* Modifying character columns using various string manipulation functions
* Customizing plots created with ggplot2 by changing the theme and position adjustments
* Saving plots created with ggplot2 as image files

---

# Exercises

This week's involves transforming two data sets in order to ask a specific
question of each. The first is a more complete set of data on the number of new
tuberculosis cases each year from 1980 to 2013 in 219 countries, taken from the
[World Health Organization Global Tuberculosis Report](https://www.who.int/health-topics/tuberculosis).

For the second "real world" example we return to the METABRIC data set and
introduce the other set of assay data obtained as part of the original 2012
study that concerns genomic copy number. We'll transform the copy number data
downloaded from [cBioPortal](https://www.cbioportal.org/study/summary?id=brca_metabric)
so that we can look at the copy number states for some of our favourite genes
and then combine that with the expression data to see how copy number influences
gene expression.

The clinical data, gene expression values and mutation data were all downloaded
from [cBioPortal](https://www.cbioportal.org/study/summary?id=brca_metabric).

```{r message = FALSE}
library(tidyverse)
```

---

**1. Tidy the WHO tuberculosis data set**

In this week's materials we demonstrated the restructuring of different
representations of a very small subset of the WHO tuberculosis data set (tables
1 to 4a and b). The more complete data set is also available in its original
format in another data frame available as part of the **tidyr** package, named
`who`.

```{r}
who
```

What are the variables in this data set?

<details><summary>Answer</summary>
* country
* iso2 and iso3 country codes
* year
* gender
* age group
* diagnosis method
* number of cases of tuberculosis
</details>


_Hint: see the help page for `who` to understand the meaning of the column headings._

_Hint: the Markdown Quick Reference available from the Help menu in RStudio is a very useful guide to formatting text in R markdown documents._

Why isn't this data set tidy and what kinds of analysis are not straightforward
with the data in this format?

Transform the data into a tidy format.

The column names `new_sp_m014` to `newrel_f65` are values for a compound
variable, i.e.  multiple variables joined together: 

* new/new_ - all are new, this is redundant and can be removed
* diagnosis - "ep", "rel", "sn", or "sp"
* gender - "m" or "f"
* age group - "014", "1524", "2534", "3544", "4554", "5564", "65" (0 - 14 years,
  15 - 24 years ... 65+ years)


The restructured data frame should have columns for each of the separated values
in these column headings.

Also remove all rows with missing values for the number of tuberculosis cases.

<details><summary>Answer</summary>

```{r}
# the hints below lead to the following solution
who_tidy <- who %>%
    pivot_longer(new_sp_m014:newrel_f65,
                 names_to = "group",
                 values_to = "cases") %>%
    mutate(group = str_remove(group, "new_")) %>%
    mutate(group = str_remove(group, "new")) %>%
    separate_wider_delim(group,
                         names = c("diagnosis", "group"),
                         delim = "_") %>%
    mutate(group = str_replace_all(group, c("m" = "m_", "f" = "f_"))) %>%
    separate_wider_delim(group,
                         names = c("gender", "age_group"),
                         delim = "_") %>%
    filter(!is.na(cases))

# this is how it could be done removing missing values as part of the pivot and
# using a more sophisticated regular expression to remove the new_? prefix in a
# single step
who_tidy <- who %>%
    pivot_longer(starts_with("new"),
                 names_to = "group",
                 values_to = "cases",
                 values_drop_na = TRUE) %>%
    mutate(group = str_remove(group, "^new_?")) %>%
    separate_wider_delim(group,
                         names = c("diagnosis", "group"),
                         delim = "_") %>%
    mutate(group = str_replace_all(group, c("m" = "m_", "f" = "f_"))) %>%
    separate_wider_delim(group,
                         names = c("gender", "age_group"),
                         delim = "_")

# alternative using str_replace_all to get groups into suitable state for single
# separate operation
who_tidy <- who %>%
    pivot_longer(starts_with("new"),
                 names_to = "group",
                 values_to = "cases",
                 values_drop_na = TRUE) %>%
    mutate(group = str_replace_all(group, 
                                   c("newrel" = "rel",
                                     "new_" = "",
                                     "_m" = "_m_",
                                     "_f" = "_f_"))) %>%
    separate_wider_delim(group,
                         names = c("diagnosis", "gender", "age_group"),
                         delim = "_")

# alternative using separate_wider_regex to split the group column in one step
who_tidy <- who %>%
    pivot_longer(starts_with("new"),
                 names_to = "group",
                 values_to = "cases",
                 values_drop_na = TRUE) %>%
    separate_wider_regex(group,
                         patterns = c("new_?",
                                      diagnosis = "[spnerl]+",
                                      "_",
                                      gender = "[mf]",
                                      age_group = "[0-9]+$"))

# alternative using names_pattern in pivot_longer (need to know regular
# expressions for this) avoiding need for separate operation
who_tidy <- who %>%
    pivot_longer(starts_with("new"),
                 names_pattern = "new_?(.*)_(.)(.*)",
                 names_to = c("diagnosis", "gender", "age_group"),
                 values_to = "cases",
                 values_drop_na = TRUE)

# convert the age group into a factor and change the levels to a more
# human-readable form
# (could also use factor or recode_factor for this)
who_tidy <- mutate(who_tidy,
                   age_group = fct_recode(age_group,
                                          "0-14" = "014",
                                          "15-24" = "1524",
                                          "25-34" = "2534",
                                          "35-44" = "3544",
                                          "45-54" = "4554",
                                          "55-64" = "5564",
                                          "65+" = "65"))

who_tidy
```

</details>

*Hint 1: the column headers from new_sp_m014 onwards contain a superfluous
prefix "new" or "new_" that you will want to remove; this can be done using one
of the stringr functions in a single step if you know about regular expressions
but is also achievable in two steps if you don't use any regular expression
magic.*

_Hint 2: separating two of the variables is a bit tricky because there is no
separator character; you could modify the variable before separating._

---

**2. Recreate `table1` from your tidy version of the WHO tuberculosis data set**

Recreate the first 3 columns of `table1`, containing a small subset of the WHO
tuberculosis data set, from your tidy version of `who`.

```{r}
select(table1, 1:3)
```

<details><summary>Answer</summary>
```{r}
who_tidy %>%
    filter(country %in% c("Afghanistan", "Brazil", "China")) %>%
    filter(year %in% c(1999, 2000)) %>%
    group_by(country, year) %>%
    summarise(cases = sum(cases))
```
</details>

_Hint: this will involve some filtering, grouping and summarization._

---

**3. Create time series plots for cases of tuberculosis**

Create a time series plot showing the number of tuberculosis cases in the United
Kingdom with separate lines for each age group. The year should be on the x axis
and the number of cases on the y axis.

<details><summary>Answer</summary>
```{r}
who_tidy <- who_tidy %>%
    group_by(country, year, age_group) %>%
    summarise(cases = sum(cases))

uk <- "United Kingdom of Great Britain and Northern Ireland"
who_tidy %>%
    filter(country == uk) %>%
    ggplot() +
    geom_line(mapping = aes(x = year, y = cases, colour = age_group))
```
</details>

_Hint: we're not distinguishing between diagnoses and genders in this plot so you'll need to sum the numbers of cases for each country, year and age group._

Use faceting in ggplot2 to create a series of plots showing the same information
for France, Brazil, Japan and Namibia.

The numbers of cases in these countries are on quite different scales. Look at
help page for `facet_wrap()` and try changing the `scales` argument so we can
more easily compare the overall patterns between countries.

<details><summary>Answer</summary>
```{r}
who_tidy %>%
    filter(country %in% c("France", "Brazil", "Japan", "Namibia")) %>%
    ggplot() +
    geom_line(mapping = aes(x = year, y = cases, colour = age_group)) +
    facet_wrap(vars(country), scales = "free_y")
```
</details>

---