Day 2: Data visualization
In-class worksheet, solutions
June 22nd, 2021
1. Plotting x-y relationships
1.1 Correlations
We will work with the iris data set available in R. This data set gives the measurements in centimeters of the variables sepal length, sepal width, petal length and petal width for 50 flowers from each of 3 species of iris. The species are Iris setosa, Iris versicolor, and Iris virginica:
# view the first several rows of the iris data set
head(iris)## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosaUsing ggplot, make a scatter plot of petal length vs. sepal length for the three species; the function you need for this is geom_point().
# plot petal length vs sepal length, colored by species
ggplot(iris, aes(x = Sepal.Length, y = Petal.Length, color = Species)) +
geom_point()
Now, do the same plot but facet by species using facet_wrap().
# plot petal length vs sepal length, faceted by species
ggplot(iris, aes(x = Sepal.Length, y = Petal.Length)) +
geom_point() +
facet_wrap(~Species)
1.2 Time series
The omp data set contains a subset of DNA microarray data measuring the differential expression of E. coli outer membrane proteins (omp) in nutrient-limited chemostatic cultures. In this particular experiment, the medium is glucose-limited. The gene column denotes one 8 genes that code for outer membrane proteins, time_min denotes the time point sampled in minutes, and au denotes the change in gene expression detected by the microarray chip (arbitrary units).
# download the `omp` data set
omp <- read_csv("https://rachaelcox.github.io/classes/datasets/ecoli_omp_expression.csv")##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## gene = col_character(),
## time_min = col_double(),
## au = col_double()
## )head(omp)## # A tibble: 6 x 3
## gene time_min au
## <chr> <dbl> <dbl>
## 1 ompA 5 -0.32
## 2 ompA 15 0.46
## 3 ompA 30 0.18
## 4 ompA 60 -0.23
## 5 ompC 5 -0.48
## 6 ompC 15 -0.56Plot the expression of each gene over time using geom_line(), coloring each line by gene. Notice anything off about this plot?
# plot the expression of each gene over time
ggplot(omp, aes(x = time_min, y = au, color = gene)) +
geom_line()
Now, make the same plot but make the size of the lines thicker by specifying the size argument inside of geom_line(). Use scale_color_colorblind() to convert the legend colors to a colorblind-friendly palette. Use xlab() and ylab to give the figure pretty labels.
# plot the expression of each gene over time & make it pretty
ggplot(omp, aes(x = time_min, y = au, color = gene)) +
geom_line(size = 2) +
scale_color_colorblind() +
xlab("Time (min)") +
ylab("Relative abundance (au)")
2. Plotting amounts
2.1 Bar plots
The bacteria data set contains data from tests of the presence of the bacterium H. influenzae in children with otitis media in the Northern Territory of Australia. We are interested in two columns of this data set: presence reports the presence (y) or absence (n) of the bacterium, treatment reports the treatment, which was placebo, drug, or drug+ (drug plus high adherence).
# download the `bacteria` data set
bacteria <- read_csv("https://rachaelcox.github.io/classes/datasets/bacteria.csv")##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## presence = col_character(),
## ap = col_character(),
## hilo = col_character(),
## week = col_double(),
## ID = col_character(),
## treatment = col_character()
## )head(bacteria)## # A tibble: 6 x 6
## presence ap hilo week ID treatment
## <chr> <chr> <chr> <dbl> <chr> <chr>
## 1 y p hi 0 X01 placebo
## 2 y p hi 2 X01 placebo
## 3 y p hi 4 X01 placebo
## 4 y p hi 11 X01 placebo
## 5 y a hi 0 X02 drug+
## 6 y a hi 2 X02 drug+Using geom_bar(), make a bar plot that shows the absolute number of cases with or without the bacterium, stacked on top of each other, for each treatment.
ggplot(bacteria, aes(x = treatment, fill = presence)) +
geom_bar()
Now modify the plot so that bars representing the absolute number of cases with or without the bacterium are shown side-by-side. Use scale_fill_brewer() to change the plot colors. Hint: This requires the argument position='dodge' in geom_bar().
ggplot(bacteria, aes(x = treatment, fill = presence)) +
geom_bar(position = 'dodge') +
scale_fill_brewer()
Now modify the plot so that bars represent the relative number of cases with or without the bacterium. What is the appropriate position option in geom_bar() to achieve this effect? Use ?geom_bar to find out. Also, try setting an alternative color palette for scale_fill_brewer() to use by specifying the type and palette arguments (see ?scale_fill_brewer for details).
ggplot(bacteria, aes(x = treatment, fill = presence)) +
geom_bar(position = 'fill') +
scale_fill_brewer(type = 'qual', palette = 'Set1')
2.2 Heat maps
The dandelion data set contains RNA-seq reads for a subset of genes differentially expressed in response to five conditions. Plot a heat map using geom_tile where each condition is on the x-axis and each gene (either transcript_dandelion or loci_arabidopsis) is on the y-axis. Fill the color of the heat map by the log2_foldchange value. Use scale_fill_distiller() to specify a continuous diverging color palette.
# download the dandelion differential expression data set
dandelion <- read_csv("https://rachaelcox.github.io/classes/datasets/dandelion_diffexp_tidy.csv")##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## transcript_dandelion = col_character(),
## baseMean = col_double(),
## z_score = col_double(),
## loci_arabidopsis = col_character(),
## protein_annotation = col_character(),
## gene_names_primary = col_character(),
## condition = col_character(),
## log2_foldchange = col_double()
## )head(dandelion)## # A tibble: 6 x 8
## transcript_dand… baseMean z_score loci_arabidopsis protein_annotat…
## <chr> <dbl> <dbl> <chr> <chr>
## 1 DN42754_c0_g1_i1 167. 1.28 AT5G67060 Transcription f…
## 2 DN42754_c0_g1_i1 167. 1.28 AT5G67060 Transcription f…
## 3 DN42754_c0_g1_i1 167. 1.28 AT5G67060 Transcription f…
## 4 DN42754_c0_g1_i1 167. 1.28 AT5G67060 Transcription f…
## 5 DN42754_c0_g1_i1 167. 1.28 AT5G67060 Transcription f…
## 6 DN42754_c0_g1_i1 167. 1.28 AT5G67060 Transcription f…
## # … with 3 more variables: gene_names_primary <chr>, condition <chr>,
## # log2_foldchange <dbl>(Note that I have told R Markdown to make a larger figure, by starting the code block with {r fig.height=6, fig.width=10} instead of {r}, because the default figure size is too narrow to show the resulting axes and map.)
ggplot(dandelion, aes(x = condition, y = transcript_dandelion, fill = log2_foldchange)) +
geom_tile() +
scale_fill_distiller(type = 'div',
palette = 'Spectral',
guide = 'colourbar',
direction = 1)
Make the same plot again, but rotate the text on the x-axis so they aren’t mashing into each other. Try a different color palette this time.
ggplot(dandelion, aes(x = condition, y = transcript_dandelion, fill = log2_foldchange)) +
geom_tile() +
scale_fill_viridis_c() +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
guides(title = 'log2(fc)')
3. Plotting distributions
3.1 Boxplots, violin plots
Using the biopsy data set, make side-by-side boxplots of clump_thickness for each outcome. The geom you need to use is geom_boxplot().
# download the biopsy data set
biopsy <- read_csv("https://rachaelcox.github.io/classes/datasets/biopsy.csv")##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## clump_thickness = col_double(),
## uniform_cell_size = col_double(),
## uniform_cell_shape = col_double(),
## marg_adhesion = col_double(),
## epithelial_cell_size = col_double(),
## bare_nuclei = col_double(),
## bland_chromatin = col_double(),
## normal_nucleoli = col_double(),
## mitoses = col_double(),
## outcome = col_character()
## )head(biopsy)## # A tibble: 6 x 10
## clump_thickness uniform_cell_si… uniform_cell_sh… marg_adhesion
## <dbl> <dbl> <dbl> <dbl>
## 1 5 1 1 1
## 2 5 4 4 5
## 3 3 1 1 1
## 4 6 8 8 1
## 5 4 1 1 3
## 6 8 10 10 8
## # … with 6 more variables: epithelial_cell_size <dbl>, bare_nuclei <dbl>,
## # bland_chromatin <dbl>, normal_nucleoli <dbl>, mitoses <dbl>, outcome <chr>ggplot(biopsy, aes(y = clump_thickness, x = outcome, fill = outcome)) +
geom_boxplot() +
scale_fill_viridis_d(direction = -1, alpha = 0.75)
Now, make the same plot again, instead using geom_violin(). What do you notice?
ggplot(biopsy, aes(y = clump_thickness, x = outcome, fill = outcome)) +
geom_violin(draw_quantiles = c(0.25, 0.5, 0.75)) +
scale_fill_viridis_d(direction = -1, alpha = 0.75, option = "E")
3.2 Histograms, density plots
Make a histogram plot of sepal lengths in the iris data set, using the default histogram settings. Then make two more such plots, with different bin widths. Use geom_histogram(). You can change the bin width by specifying the argument binwidth = or bins =. See ?geom_histogram for more information.
# default settings
ggplot(iris, aes(x = Sepal.Length)) +
geom_histogram()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# wider bins
ggplot(iris, aes(x = Sepal.Length)) +
geom_histogram(binwidth = 0.2)
# even wider bins
ggplot(iris, aes(x = Sepal.Length)) +
geom_histogram(binwidth = 0.4)
# divide all observations by a set # of bins
ggplot(iris, aes(x = Sepal.Length)) +
geom_histogram(bins = 10)
Instead of geom_histogram(), now use geom_density() and fill the area under the curves by species identity.
ggplot(iris, aes(x = Sepal.Length, fill = Species)) +
geom_density()
Now make the areas under the curve partially transparent, so the overlap of the various distributions becomes clearly visible.
ggplot(iris, aes(x = Sepal.Length, fill = Species)) +
geom_density(alpha = 0.7)
4. Layering geoms
Bonus challenge: For the iris data set, make a plot of the 2d distribution of petal length vs. sepal length, by making an x-y plot that shows the individual data points as well as contour lines indicating the density of points in a given spatial region.
ggplot(iris, aes(x = Sepal.Length, y = Petal.Length, color = Species)) +
geom_point() +
geom_density2d()
Now instead of contour lines, add a fitted straight black line (not a curve, and no confidence band!) to each group of points. You’ll need to check ?geom_smooth to see which arguments you’ll need to specify.
ggplot(iris, aes(x = Sepal.Length, y = Petal.Length, color = Species)) +
geom_point() +
geom_smooth(
aes(group = Species), # make a line for each species
method = lm, # fit a linear model
color = 'black',
se = FALSE) # no confidence interval## `geom_smooth()` using formula 'y ~ x'
In this last example, because we are manually overriding the color of the lines, we need to set the group aesthetic to tell ggplot2 to draw a separate line for each species.