-
Notifications
You must be signed in to change notification settings - Fork 1
Expand file tree
/
Copy pathday5.Rmd
More file actions
697 lines (582 loc) · 22.8 KB
/
day5.Rmd
File metadata and controls
697 lines (582 loc) · 22.8 KB
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
---
title: "Intro to R 2026 day 5"
author: "Matt Cannon"
date: '2026-04-24'
output:
html_document:
code_folding: hide
toc: true
toc_depth: 2
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
Load libraries
```{r}
library(tidyverse)
```
# Day 5
## General stuff
### geoms as layers
- When you're making a plot using ggplot2, think about it as if you're making a canvas and then adding layers to change it
- The layers are added in the order you add the code
- Layers can do things like add points, lines, or more abstract things like change the font of parts of the plot
- Each layer is separated by a "+" sign, which is pretty much the only time you'll see this in R
- The documentation for ggplot2 is quite good
- https://ggplot2.tidyverse.org/reference/
- One trick: skim the documentation for a plot that looks similar to what you want, then look at the code they used to make it
- Another: You can also do a google image search for ggplot and the type of plot you want to make then go to the page where the image is from. Often times you'll find the code they used to make it.
### Anatomy of a ggplot
ggplot( # This makes your "canvas"
data_frame_here, # The data you're using - must be a long form data.frame or tibble - see day 3 code for discussion
aes( # "Aesthetics" - x, y, colors, etc - Applies to all layers unless overwritten
x = x here, # The name of the column that has your x values - note: don't put this in quotes
y = y here, # The name of the column that has your y values
other stuff # You can specify a lot of things here - check the documentation
)
) +
layer_1 + # Layers can be things like geom_point(), geom_line(), geom_boxplot(), etc
layer_2 + # You can add multiple different types of geoms, or even multiple of the same type
layer_3 +
other_stuff + # You can add other things like themes, labels, etc
more_stuff +
change_how_something_looks +
modify_the_plot_somehow
You can really customize the plot and the code can get quite long. Thinking about it as layers can help keep track of what you're doing
I'm not going over all the things in this document, just the ones marked with lots of pound signs ##############################################################
I put the rest in here for reference. If you go through it, I've tried to drop in a bunch of random bits of code that might be useful to you in the future.
## Functions for today's activity, also for reference
#### geom_point() ##############################################################
Add a scatterplot layer to the plot. This requires the x and y aesthetics to be mapped so that ggplot knows what to plot on the x and y axes.
```{r}
ggplot(
mtcars,
aes(
x = hp,
y = carb
)
) +
geom_point()
# Set the overall look of the plot - I like the white background better than gray
# This will apply to all plots you make after you run this code, so you only need to run it once per session
theme_set(theme_bw()) # Changes default colors
theme_update(plot.title = element_text(hjust = 0.5)) # Force title to be in the center
ggplot(
mtcars,
aes(
x = hp,
y = carb
)
) +
geom_point() +
theme_bw()
```
#### aes() ################################################
Change colors/shapes/etc
- fill = color of the **inside** of a shape
- color = color of the **border** of a shape, note that points use color for the color of the point
- shape = the shape of the point
- size = the size of the point or line
- alpha = transparency of the point or line
- linetype = the type of line
- linewidth = the width of the line, used in geoms that draw lines like geom_line() or geom_errorbar()
- group = the group that the data is in, used for things like boxplots
- label = Text to plot somewhere, used in geom_text(), geom_label(), etc
If you put the argument outside of aes() inside of a geom_{something}()
- applies to everything in that layer
- does not show up in the legend
If you put it inside aes() in either ggplot(aes({thing here})) or geom_{something}(aes({thing here}))
- Applies to parts of the layer/plot
- Shows up in legend
- If it is in ggplot(), it applies to all layers
- Unless another layer overwrites it
```{r}
mtcars %>%
ggplot(
aes(
x = as.factor(cyl),
y = hp
)
) +
geom_boxplot(
color = "red", # These specifications aren't in aes()
fill = "blue", # They apply to everything and we give specific values
linewidth = 4
)
ggplot(
mtcars,
aes(
x = hp, # These specifications are in aes()
y = carb, # They apply to parts of the plot
color = as.factor(vs), # They show up in the legend
shape = as.factor(gear), # We give column names that have the groups we'll use
size = disp
)
) +
geom_point()
```
#### geom_histogram()
- Better looking than hist()
- Specify how many bars to show with `bins`
- You put only x in the aes()
```{r}
hist(storms$pressure, n = 200)
ggplot(
storms,
aes(x = pressure)
) +
geom_histogram(bins = 200)
```
#### geom_density()
- Smooth your data out in a smoothed histogram-looking plot
- The adjust argument can make it more or less smooth - 0 is the raw data, 1 is very smooth
- This can actually hide some of the nuance in your data
```{r}
ggplot(
storms,
aes(
x = pressure,
fill = as.factor(category)
)
) +
geom_density(alpha = 0.1, adjust = 0.5) + # alpha sets transparency -- 0 is clear, 1 is opaque
scale_color_brewer(palette = "Set2") # The brewer color palettes are pretty, and "Set2" is color-blind friendly
```
#### geom_col() ################################################################
Add a barplot layer to the plot
- geom_col() takes both x and y aesthetics - x is the category / factor on the x axis, y is the height of the bar
- Unlike geom_bar() which counts the number of observations in each x category, geom_col() uses the y value you supply to set the height of the bar
```{r}
mtcars %>%
rownames_to_column("car_make") %>% # This makes the row names a column so we can use it in the plot
mutate(car_make = str_remove(car_make, " .+")) %>% # This removes everything after the first space in the car_make column
ggplot(
aes(
x = car_make,
y = hp
)
) +
geom_col(color = "red") +
theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
ggtitle("this is a title") +
theme_update(plot.title = element_text(hjust = 0.5))
mtcars %>%
rownames_to_column("car_make") %>%
ggplot(aes(y = car_make, x = hp)) +
geom_col()
```
#### geom_errorbar()###########################################################
Add error bars as a layer to plot
- Need to have a column of values to specify how high/low the error bars should go, such as a standard deviation or standard error
- You specify the ymin / ymax aesthetics inside aes() in geom_errorbar() to tell ggplot how far the error bars should extend
```{r}
summarized_hp <-
mtcars %>%
group_by(cyl) %>%
summarize(
mean_hp = mean(hp),
sd_hp = sd(hp)
)
summarized_hp
ggplot(
summarized_hp,
aes(
x = as.factor(cyl),
y = mean_hp
)
) +
geom_col() +
geom_errorbar(
aes(
ymin = mean_hp - sd_hp,
ymax = mean_hp + sd_hp
),
width = 0.2,
color = "red"
)
# This has error bars in the back since that layer was added first!
ggplot(
summarized_hp,
aes(
x = as.factor(cyl),
y = mean_hp
)
) +
geom_errorbar(
aes(
ymin = mean_hp - sd_hp,
ymax = mean_hp + sd_hp
),
color = "red"
) +
geom_bar(stat = "identity") +
geom_point(color = "blue")
```
#### geom_jitter()
Drop-in replacement for geom_point()
Adds noise in x and y directions by default so you can see individual points
```{r}
ggplot(
summarized_hp,
aes(
x = as.factor(cyl),
y = mean_hp
)
) +
geom_col() +
geom_errorbar(
aes(
ymin = mean_hp - sd_hp,
ymax = mean_hp + sd_hp
)
) +
geom_jitter(
data = mtcars, # Note that we can provide an individual geom it's own data
aes(
x = as.factor(cyl),
y = hp
),
color = "red", height = 0
)
```
#### ggbeeswarm package
Adds noise in x direction and pulls things to the center to make it look nicer than jitter
- Need to install the ggbeeswarm package first with `install.packages("ggbeeswarm")`
- Works similar to geom_point(), though it struggles with large datasets
```{r}
ggplot(
summarized_hp,
aes(
x = as.factor(cyl),
y = mean_hp
)
) +
geom_col() +
geom_errorbar(
aes(
ymin = mean_hp - sd_hp,
ymax = mean_hp + sd_hp
)
) +
ggbeeswarm::geom_beeswarm(
data = mtcars,
aes(
x = as.factor(cyl),
y = hp
),
color = "red"
)
```
#### geom_boxplot() ############################################################
Boxplots!
- Specify both x and y aesthetics - x is the factor / category on the x axis, y is the numeric variable that the boxplot will summarize
- Can turn the plot sideways by making y a category and x the numeric variable
```{r}
mtcars %>%
ggplot(aes(x = as.factor(cyl), y = hp)) +
geom_boxplot()
```
#### geom_violin()
- Similar to a boxplot, but shows the density of the distribution of the numeric variable along the y axis for each category on the x axis
- This can over-smooth the data if you have non-standard data distributions
- The `adjust` argument controls the bandwidth / smoothing of the density estimate - higher values smooth more, lower values less
```{r}
mtcars %>%
ggplot(aes(x = as.factor(cyl), y = hp)) +
geom_violin(adjust = 1, fill = "red", alpha = 0.5)
```
#### geom_smooth()
Add a regression line to your plot
Can specify the type of regression with `method`
Can remove gray confidence interval with `se = FALSE`
```{r}
ggplot(
mtcars,
aes(
x = mpg,
y = disp
)
) +
geom_point() +
geom_smooth(method = "lm")
```
#### Use multiple datasets in the same plot
Contents of the columns used as axes have to match
Each geom using it's own data needs to have "data = " in it
```{r}
ggplot(
summarized_hp,
aes(
x = as.factor(cyl),
y = mean_hp
)
) +
geom_col() +
geom_errorbar(
aes(
ymin = mean_hp - sd_hp,
ymax = mean_hp + sd_hp
),
width = 0.5
) +
geom_jitter(
data = mtcars,
aes(
x = as.factor(cyl),
y = hp
),
color = "red"
)
```
#### labs()
Add labels
```{r}
ggplot(mtcars) +
geom_point(aes(x = hp, y = disp)) +
labs(
x = "X label goes here",
y = "Y label goes here",
title = "This is a title",
subtitle = "Subtitle!",
caption = "Cannon et al, 2034!",
tag = "A"
)
```
#### facet_wrap() ##############################################################
Faceting lets you split a plot into multiple panels based on the values of one or more categorical variables. `facet_wrap()` will wrap the panels into a grid, and you can specify the faceting formula with `~ var1`.
- You can facet on more than one variable by using `+` in the formula - e.g. `~ var1 + var2` will make panels for each combination of the levels of var1 and var2
- You can control the scales of the facets independently with the `scales` argument - e.g. `scales = "free_y"` lets the y axis vary across the facets rather than being fixed across all panels
- You can also control the number of rows / columns in the wrapped facets with `nrow` / `ncol` - e.g. `facet_wrap(~ var1, nrow = 1)` will put all the facets in a single row rather than letting ggplot pick the layout automatically
```{r}
dplyr::storms
ggplot(
dplyr::storms,
aes(x = pressure)
) +
geom_histogram(bins = 100) +
facet_wrap(~ category + status)
ggplot(
dplyr::storms,
aes(x = pressure)
) +
geom_histogram(bins = 100) +
facet_wrap(~category, scales = "free_y")
ggplot(
dplyr::storms,
aes(x = pressure)
) +
geom_histogram(bins = 100) +
facet_wrap(
~category,
nrow = 1,
scales = "free_y"
)
```
#### theme()
A good guide to what theme elements are what:
https://henrywang.nl/ggplot2-theme-elements-demonstration/
You can change just about any aspect of how the plot looks with theme()
You'll have to look up how to change each specific bit
Generally, you call theme() and then specify the elements you want to tweak inside it - e.g. `theme(legend.position = "bottom")` etc. See the help documentation for theme() for all the elements you can modify - `?theme`.
```{r}
ggplot(
storms,
aes(
x = pressure,
fill = category
)
) +
geom_density(alpha = 0.5) + # alpha sets transparency - 0 is clear, 1 is opaque
scale_color_brewer(palette = "Set2") + # The brewer color palettes are pretty, and "Set2" is color-blind friendly
theme(
legend.position = c(0.2, 0.9),
legend.direction = "horizontal",
legend.key.width = unit(2, "cm"),
legend.title = element_text(face = "bold"),
strip.background = element_rect(color = "white", fill = "white"),
strip.text.x = element_text(size = 30, face = "bold"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank(),
axis.line.x = element_line(color = "black"),
axis.text.y = element_blank(),
axis.ticks.y = element_blank()
)
```
### ggsave() ###################################################################
Saves an image to your computer
- Specify the filename, with path to the folder if needed
- Optionally specify the width / height of the output image (in inches by default) with the `width` / `height` arguments
- You can also pass a variable with a plot output to the `plot` argument rather than relying on the last plot that was printed to the plotting pane
- If you don't specify the width / height, ggsave will use the size of the current plotting device / plot window to determine the output image dimensions, which can be problematic
```{r}
ggsave(
"examplePlot.png",
width = 10,
height = 8
)
plot_in_variable <-
ggplot(
dplyr::storms,
aes(x = pressure)
) +
geom_histogram(bins = 100) +
facet_wrap(~category)
# ggsave will get the plot size from how big your "plot" windows is in the right pane if you don't specify it!!!!!
ggsave(
"anotherExamplePlot.png",
plot = plot_in_variable
)
```
### heatmap() ##################################################################
- This function takes in a numeric matrix (or something that can be coerced to a matrix) and produces a heatmap of the values in that matrix.
- The column and row names of the matrix are used to label the axes of the heatmap.
- Can specify if you want to cluster by rows and/or columns
If you like heatmaps also check out the pheatmap package. You may need to install it first with `install.packages("pheatmap")`.
```{r}
heatmap(as.matrix(mtcars))
# Scaling the categories instead
heatmap(as.matrix(mtcars), scale = "column")
pheatmap::pheatmap(as.matrix(mtcars))
```
### PCA
- Principal Component Analysis (PCA) is a technique to reduce the dimensionality of a dataset by transforming the original variables into a new set of uncorrelated variables (the principal components) that capture the most variance in the data.
- Breaks the data down into new principal components that are linear combinations of the original variables, ordered by the amount of variance in the data that each component explains.
- The first principal component (PC1) explains the most variance in the data, the second principal component (PC2) explains the second most, and so on - so you can plot the first few PCs to see the main directions of variation in the dataset.
If you run the example code below you will need to install the ggrepel package with `install.packages("ggrepel")` first. This is a package that makes it easier to add text labels to points in a ggplot without the labels overlapping each other.
```{r}
mtcars <- mtcars
pca_raw <- prcomp(mtcars, scale = TRUE)
?prcomp
# The output of prcomp is a list
# The element named "x" contains the actual principle components
# Though, it needs to be converted from a matrix to a dataframe for ggplot
pc_values <-
pca_raw$x %>%
as.data.frame() %>%
rownames_to_column("car_make")
ggplot(
pc_values,
aes(
x = PC1,
y = PC2
)
) +
geom_point() +
ggrepel::geom_text_repel(aes(label = car_make))
```
## Activity
For the activities, lets use the txhousing dataset. You can load it like this:
```{r}
txhousing <- dplyr::txhousing
```
This has information on housing in Texas - things like sales, median price, and volume of transactions across different cities and months. It includes data from 2000 through 2015.
> Column Name Data Type Description
> city <chr> Name of the Multiple Listing Service (MLS) area or city in Texas.
> year <int> The year the data was recorded (2000-2015).
> month <int> The month the data was recorded (1-12).
> date <dbl> A continuous numeric representation of the date (e.g., 2000.000 for Jan 2000). Useful for plotting time series smoothly.
> sales <dbl> The total number of home sales that occurred during that month.
> volume <dbl> The total value of all home sales (in US Dollars) for that month.
> median <dbl> The median sale price (in US Dollars) of homes sold that month.
> listings <dbl> The total number of active home listings available.
> inventory <dbl> "Months inventory" - the estimated amount of time (in months) it would take to sell all current listings at the current pace of sales.
### Make a scatterplot using geom_point()
Lets make a scatterplot of median home price over the years. The x-axis will be the year and the y-axis will be the median home price.
```{r}
```
### Lets make a boxplot of the median housing prices by city
Make the y-axis the city and the x-axis the median housing price so that the boxplots are horizontal rather than vertical.
```{r}
```
### Now lets plot the median housing prices from 2015 over the months
I'll make a new dataset for you to use below that only contains the 2015 data.
```{r}
txhousing_2015 <-
txhousing %>%
dplyr::filter(year == 2015)
txhousing_2015
```
Lets use this to help people see the difference between adding arguments inside vs outside of `aes()`
1) Make the plot with no color specified for the points
2) Make another plot where you put the color argument inside `aes()` so that the points are colored by city
3) Now make a third plot where you keep `color = city`, but put it inside of geom_point() and not inside `aes()`. What happens? R is looking for a variable called `city` in the environment outside of the data - it can't find it so it throws an error.
4) Now replace `city` with `"city"`. R is trying to make the points the color "city" - that is, it is treating "city" as a color value rather than a variable in the data, which isn't a valid color so it errors out.
5) Now try putting `color = "blue"` inside of the aes(). Are your points blue? This is because inside of the `aes()`, color specifies what category to use. R is treating "blue" as a category for your data and so it colors them by the first default color, which is red.
6) Now try putting `color = "blue"` outside of the aes(). Are your points blue now?
```{r}
```
### Now lets make a summary of total listings each year across all cities and plot that over the years
You will need to group the data by year and summarize the total number of listings for each year.
The final plot will show the total number of listings on the y-axis and the year on the x-axis.
```{r}
txhousing_listings_by_year <-
txhousing %>%
dplyr::group_by(year) %>%
dplyr::summarise(
total_listings = sum(listings, na.rm = TRUE)
)
```
```{r}
```
### Now lets make the plot look nicer
Copy the code from the previous plot. Do each of the following changes below (one at a time) to slowly build up the plot. Generally when making a plot I don't write all of the code at once, I start with a basic plot and then add code incrementally.
- Add in a title
- Add labels for the x- and y-axes
- Use `+ theme_bw()` to give the plot a clean black and white theme
- Make the line thicker with the `linewidth` argument
- Make the points larger and red with the `size` and `color` arguments.
- Add geom_smooth() to add a smoothed trend line to the plot behind all the other geoms.
- Force the title to be centered using `theme()` with the `plot.title = element_text(hjust = 0.5))` argument.
```{r}
```
If you want to get advanced, try adding more arguments to `theme()` to do things like change the font size of the x/y labels, tick marks, legend text, etc.
## Advanced Activities
### Recreate plots from a paper
- Data is from https://doi.org/10.1371/journal.pbio.2005756
- Najafov A, Zervantonakis IK, Mookhtiar AK, Greninger P, March RJ, Egan RK, et al. (2018) BRAF and AXL oncogenes drive RIPK3 expression loss in cancer. PLoS Biol 16(8): e2005756. https://doi.org/10.1371/journal.pbio.2005756
- I downloaded their figure data and made it easier to import:
- Fig1A_partial.txt
- Fig1C.txt
- Fig2F.txt
- Fig3F.txt
- You can see what the plots should look like in the corresponding files in materials/Figure_***.png files
- You don't have to make them look exact, feel free to play around with how you want them to look
#### Figure 1A - geom_boxplot()
```{r}
```
#### Figure 1C - heatmap()
```{r}
```
#### Figure 2F - geom_point()
```{r}
```
#### Figure 3F - geom_col
use geom_errorbars()
and facet_wrap() to split the plots by cell_line
```{r}
```
## If you're super fast:
Keep going, do the rest 🙃
Download the data from the supplemental data here:
https://doi.org/10.1371/journal.pbio.2005756.s001
####
# You're done! You made it through introduction to R!
####
Suggestions for what to do next:
- Try to practice working in R as much as you can
- **Find ways** that R can be useful for you and explore packages that are useful for the type of work you do
- You're still learning! Don't get down on yourself when you struggle, because you will
- The better you get at R, the less you'll struggle, but everyone runs into problems
- Read the documentation for functions/packages you think you'll use a lot
- Make all the errors!
- Then google them and try to understand what went wrong
- Ask AI for help, but make sure to ask it to explain the answer rather than just give the code
- Go back through and re-read the materials to reinforce what you've learned and pick up little details you may have missed
- Check out other R resources/tutorials to continue learning beyond this workshop
- youtube has lots of good tutorials on R
- If you're going to the SCRGOT single cell workshop, do some reading on Seurat
- https://satijalab.org/seurat/
- I run a weekly bioinformatics meeting where we talk about bioinformatics topics and help troubleshoot problems
- See https://docs.google.com/spreadsheets/d/1yCwYNAL9MhYErl98CikoiTW6yx5gg54w66eD73WsRVc/edit?usp=sharing for example previous topics and video recordings of lectures
- Email/teams me if you want to join (OSU welcome too)