-
Notifications
You must be signed in to change notification settings - Fork 10
/
week6.Rmd
842 lines (636 loc) · 30.8 KB
/
week6.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
---
title: "Week 6 -- 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, c(`1999`, `2000`), 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()`
The `separate()` function allows us to split a character column into multiple
columns based on a delimiter or separator.
```{r}
table3_separated <- separate(table3, rate, into = c("cases", "population"), sep = "/")
table3_separated
```
The `sep` argument takes a regular expression that defines how to split the
values. We've mentioned regular expressions before -- they are a language for
specifying search patterns used to find sequences of characters within text and
well worth learning. 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_at()`.
```{r}
mutate_at(table3_separated, vars(cases, population), as.integer)
```
But another option is to specify `convert = TRUE` when carrying out the
`separate()` operation, in which case it will deduce the type of the values and
convert the column accordingly.
```{r}
table3_separated <- separate(table3, rate, into = c("cases", "population"), sep = "/", convert = TRUE)
mutate(table3_separated, rate = cases / population)
```
---
# 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) +
geom_boxplot(mapping = aes(x = ER_status, y = Expression)) +
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", "beer"), sep = " loves ")`
` [1] "Ashley loves hiking" "Matt loves beer"`
---
**`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.
Alias for `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("Matthew", start = 1, end = 4)`
` [1] "Matt"`
`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(protein_levels, sample, into = c("treatment", "time", "replicate"))
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
metabric <- read_csv("data/metabric_clinical_and_expression_data.csv") %>%
mutate(Integrative_cluster = factor(Integrative_cluster, levels = c("1", "2", "3", "4ER-", "4ER+", "5", "6", "7", "8", "9", "10"))) %>%
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 `black_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}
metabric %>%
ggplot(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}
metabric %>%
ggplot(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
---
# Assignment
Assignment: [assignment6.Rmd](assignments/assignment6.Rmd)
Solutions: [assignment6_solutions.Rmd](assignments/assignment6_solutions.Rmd)
and [assignment6_solutions.html](assignments/assignment6_solutions.html)