Data visualization in ggplot
RNA-seq edition
Kim Dill-McFarland
version September 08, 2022
3. Data visualization in ggplot
In this workshop, we introduce you to data visualization in
ggplot (docs.ggplot2.org), the tidyverse package for plots. In
it, we cover how to create:
- Volcano plot
- PCA plot
- Customization such as color, shape, theme
- Facets
During the workshop, we will build an R script together, which will be posted as ‘live_notes’ after the workshop at https://github.com/BIGslu/workshops/tree/main/2022.08.15_R.tidyverse.workshop/live_notes
We load ggplot within the tidyverse. For
more information on installation, see the setup
instructions.
library(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6 ✔ purrr 0.3.4
## ✔ tibble 3.1.8 ✔ dplyr 1.0.9
## ✔ tidyr 1.2.0 ✔ stringr 1.4.1
## ✔ readr 2.1.2 ✔ forcats 0.5.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()Intro to ggplot
Why ggplot?
ggplot2 is an implementation of Grammar of Graphics (Wilkinson 1999) for R
Benefits:
- handsome default settings
- snap-together building block approach
- automatic legends, colors, facets
- statistical overlays like regressions lines and smoothers (with confidence intervals)
Drawbacks:
- it can be hard to get it to look exactly the way you want
- requires having the input data in a certain format
ggplot layers
ggplot functions by adding layers to a plot, somewhat
like how dplyr/tidyr sequentially modify data
with the pipe %>%. Instead, ggplot uses the
+ to combine layers such as:
- data: 2D table (
data.frame) of variables - aesthetics: map variables to visual attributes (position, color, etc)
- geoms: graphical representation of data (points, lines, etc)
- stats: statistical transformations to get from data to points in the plot (binning, summarizing, smoothing)
- scales: control how to map a variable to an aesthetic
- facets: juxtapose mini-plots of data subsets, split by variables
- guides: axes, legend, etc. that reflect the variables and their values
The idea is to independently specify and combine layers to create the plot you want. There are at least three things we have to specify to create any plot:
- data
- aesthetic mappings from data variables to visual properties
- a geom describing how to draw those properties
Data
A quick reminder about our data. There are bulk RNA-seq contained in
the dat object in file “data/dat_voom.RData”. Normalized
log2 counts per million (CPM) are contained in dat$E and
library metadata are in dat$targets. You will not need the
other pieces of dat for this workshop.
load("data/dat_voom.RData")
dat$E[1:3,1:3]## Loading required package: limma## pt01_Media pt01_Mtb pt02_Media
## A1BG 1.620022 1.717603 1.086946
## A2M 7.259912 6.680758 5.939909
## A2ML1 -4.052404 -3.515057 -4.015712dat$targets[1:3,]## group lib.size norm.factors libID ptID condition age_dys sex
## pt01_Media 1 8295928 1.092872 pt01_Media pt01 Media 12410 M
## pt01_Mtb 1 5716203 0.819490 pt01_Mtb pt01 Mtb 12410 M
## pt02_Media 1 8087601 1.119352 pt02_Media pt02 Media 12775 M
## ptID_old RNAseq methylation total_seq sample.weights
## pt01_Media pt00001 TRUE FALSE 9114402 0.8943797
## pt01_Mtb pt00001 TRUE FALSE 8918699 0.9850363
## pt02_Media pt00002 TRUE FALSE 9221555 1.2582092Here is a brief description of the variables in the metadata.
libID(character): Unique library ID for each sequencing library. Formatted as ptID_conditionptID(character): Patient ID. Each patient has 2 libraries which differ in theirconditioncondition(character): Experimental condition. Media for media control. Mtb of M. tuberculosis infected.age_dys(numeric): Age in agessex(character): Biological sex, M or FptID_old(character): Old patient ID with leading zeroesRNAseq(logical): If the library has pass-filter RNA-seq datamethylation(logical): If the library has pass-filter methylation datatotal_seq(numeric): Total number of pass-filter sequences
Basic plots
Boxplot
We will start with a boxplot of the total sequences per library. First, we call a ggplot object. Since we’ve given it no data or aesthetics, this is a blank plot.
ggplot()
Then, we give the plot data and specify aesthetics for which parts of
the data we want to plot. Now, we see that ggplot knows what the scale
of the y data will be. Throughout this section, we will highlight lines
of code that change with the comment
#<--- see change.
ggplot(dat$targets) +
aes(y = total_seq) #<--- see change
Finally, we add a geom to denote how we want the data displayed, in this case as a boxplot.
ggplot(dat$targets) +
aes(y = total_seq) +
geom_boxplot() #<--- see change
Since we did not specify an x-variable, all the data are plotted in one boxplot. We add an aesthetic to specify x and split up the data.
ggplot(dat$targets) +
aes(y = total_seq, x = condition) + #<--- see change
geom_boxplot()
Data ordering with factors
Our x-variable happens to be in the correct alphabetical order with
the Media controls before Mtb infection.
However, alphabetical does not always give you want you want. You can
force a variable’s order by converting it to a factor. Let’s switch the
condition order.
dat$targets %>%
mutate(condition_fct = factor(condition,
levels = c("Mtb", "Media"))) %>% #<--- see change
ggplot() +
aes(y = total_seq, x = condition_fct) +
geom_boxplot()
Notice how we now pipe %>% the data into ggplot so
that we can modify it in dplyr first.
Dotplot
These same data could be displayed as a dotplot instead. We return to
the original condition order since it was correct.
ggplot(dat$targets) +
aes(y = total_seq, x = condition) +
geom_point() #<--- see change
And if we want to make sure we see every point, a jitter works nicely.
ggplot(dat$targets) +
aes(y = total_seq, x = condition) +
geom_jitter() #<--- see change
Multiple geoms
The building block nature of ggplot allows us to add multiple geoms on top of each other. So we can combine the two previous plots like so
ggplot(dat$targets) +
aes(y = total_seq, x = condition) +
geom_boxplot() + #<--- see change
geom_jitter() #<--- see change
But you may notice that there are some duplicate points. Looking back at the first boxplots, there are points for data outside the box. These points remain in the above plot plus appear jittered from the second geom. We can fix this by telling the boxplot to not plot the outliers.
ggplot(dat$targets) +
aes(y = total_seq, x = condition) +
geom_boxplot(outlier.shape = NA) + #<--- see change
geom_jitter()
Layer order matters
Because ggplot works in layers, the order matters. If, for example,
we put the boxplot after the jitter, the boxes cover up points! We use
alpha for transparency so you can see these points. Thus,
be careful when considering how you setup your plot!
ggplot(dat$targets) +
aes(y = total_seq, x = condition) +
geom_jitter() + #<--- see change
geom_boxplot(outlier.shape = NA, alpha = 0.8) #<--- see change
Why are my dots moving?
With jitter, you are adding random error to be able to see all the points. Because this is truly random, every time you run jitter, you will get a slightly different plot. To combat this, set a seed before making the plot.
set.seed(468) #<--- see change
ggplot(dat$targets) +
aes(y = total_seq, x = condition) +
geom_boxplot(outlier.shape = NA) +
geom_jitter()
set.seed(468) #<--- see change
ggplot(dat$targets) +
aes(y = total_seq, x = condition) +
geom_boxplot(outlier.shape = NA) +
geom_jitter()
RNA-seq plots
PCA
One dotplot of interest in RNA-seq analysis is principle component analysis (PCA). In this plot, each RNA-seq library is a single point. The axes are unitless but in general, points that are closer together represent libraries with more similar gene expression. You can see large scale changes in gene expression in PCA. However, even when you don’t see clear PCA groupings, there can be significant differentially expressed genes. So don’t despair if your data don’t cluster in this plot!
Here, we calculate PCA from the normalized counts. Note that we
transpose t() the counts data to get one PCA point per
library, instead of one per gene.
dat.pca <- prcomp(t(dat$E), scale.=TRUE, center=TRUE)
dat.pca$x[1:3,1:3]## PC1 PC2 PC3
## pt01_Media -103.91401 -33.55704 22.40058
## pt01_Mtb 52.95006 -74.68146 24.78278
## pt02_Media -60.61507 27.13460 24.58148Then we use join to add metadata to the PCA data, so we can annotate the PCA plot with variables of interest.
dat.pca.meta <- as.data.frame(dat.pca$x) %>%
rownames_to_column("libID") %>%
inner_join(dat$targets)## Joining, by = "libID"A nice trick with joins is to not specify what to join
by. In this case, the function finds all columns shared
between the data frames and uses them to match rows.
With geom_point, we have a PCA!
ggplot(dat.pca.meta) +
aes(x = PC1, y = PC2) +
geom_point()
Color and shape
But it’s not very informative, so let’s annotate with our condition groups. In general, you want to use annotations in the following order: space, color, shape, size. These are ordered by how easily the human eye can discern differences.
We already have space as libraries are separated along x and y. Now we add color.
ggplot(dat.pca.meta) +
aes(x = PC1, y = PC2, color = condition) + #<--- see change
geom_point()
I like to keep all my aesthetics in one layer, but you could keep
adding aes() functions to the plot if you prefer. This
gives the same plot.
ggplot(dat.pca.meta) +
aes(x = PC1, y = PC2) + #<--- see change
geom_point() +
aes(color = condition) #<--- see changeWe could also use shape for the same data. Which plot to you think is easier to interpret?
ggplot(dat.pca.meta) +
aes(x = PC1, y = PC2, shape = condition) + #<--- see change
geom_point()
We’re not showing size or transparency because these aesthetics are difficult to see and should be avoided unless you really need them or the differences between data points are large enough to see in these aesthetics.
Labels
We can also improve our axes and legend labels. First, we see what proportion of variation is explained by each PC.
summary(dat.pca)$importance## PC1 PC2 PC3 PC4 PC5 PC6
## Standard deviation 57.21487 38.99933 35.32881 34.21306 29.76604 27.13699
## Proportion of Variance 0.24395 0.11334 0.09301 0.08723 0.06603 0.05488
## Cumulative Proportion 0.24395 0.35729 0.45030 0.53753 0.60356 0.65844
## PC7 PC8 PC9 PC10 PC11 PC12
## Standard deviation 24.48881 23.98672 21.83296 21.23941 19.84912 18.33097
## Proportion of Variance 0.04469 0.04288 0.03552 0.03362 0.02936 0.02504
## Cumulative Proportion 0.70313 0.74601 0.78153 0.81515 0.84451 0.86955
## PC13 PC14 PC15 PC16 PC17 PC18
## Standard deviation 18.05337 17.40958 17.13667 15.67883 14.48455 14.12701
## Proportion of Variance 0.02429 0.02259 0.02188 0.01832 0.01563 0.01487
## Cumulative Proportion 0.89384 0.91642 0.93831 0.95663 0.97226 0.98713
## PC19 PC20
## Standard deviation 13.14034 1.137607e-13
## Proportion of Variance 0.01287 0.000000e+00
## Cumulative Proportion 1.00000 1.000000e+00And use that for labels as well as change the legend.
ggplot(dat.pca.meta) +
aes(x = PC1, y = PC2, color = condition) +
geom_point() +
labs(x = "PC1 (24%)", y = "PC2 (11%)", #<--- see change
color = "Infection")
Coordinates
Since PCA units are relative, it is nice to have the scales of the x
and y be the same. That is, 1 inch of space corresponding to the same
values for both. This is accomplished with coord_fixed.
Note that in this PCA, this function does not change the plot much but
it can can a big effect in other PCA! There are also several other
coordinate systems you can see with coord_.
ggplot(dat.pca.meta) +
aes(x = PC1, y = PC2, color = condition) +
geom_point() +
labs(x = "PC1 (24%)", y = "PC2 (11%)",
color = "Infection") +
coord_fixed() #<--- see change
Themes
ggplot has a number of themes that change the aesthetics
of the plot. I like theme_classic for PCA. Type
theme_ into the console to see the available themes
suggested with auto-complete. You can also checkout more themes in the
package ggthemes.
ggplot(dat.pca.meta) +
aes(x = PC1, y = PC2, color = condition) +
geom_point() +
labs(x = "PC1 (24%)", y = "PC2 (11%)",
color = "Infection") +
coord_fixed() +
theme_classic() #<--- see change
volcano plot
Let’s move on to another application of geom_point. A
volcano plot shows differentially expressed gene significance and fold
change.
Linear modeling
We’ll need to perform linear modeling to determine these results.
Below is the workflow for comparing Mtb-infected vs media samples,
corrected for sex and blocked by patient. We use our R package kimma. You
DO NOT need to run these steps. They are here in case
you want to run a similar analysis in the future.
install.packages("devtools")
devtools::install_github("BIGslu/kimma")library(kimma)fit <- kmFit(dat, model = "~ condition + sex + (1|ptID)",
run.lme = TRUE, use.weights = TRUE)lme/lmerel model: expression~condition+sex+(1|ptID)
Input: 20 libraries from 10 unique patients
Model: 20 libraries
Complete: 13419 genes
Failed: 0 genessave(fit, file="data/kimma_result.RData")Instead, load the results from the data downloaded for this workshop.
Note that this is an S3 object similar to the dat object.
In it, we have the linear mixed effects (lme) model results in
fit$lme, which includes our log2 fold change
estimate and significance FDR for each gene
and each variable.
load("data/kimma_result.RData")
fit$lme## # A tibble: 40,257 × 7
## model gene variable estimate pval FDR gene_biotype
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <chr>
## 1 lme A1BG condition -0.313 0.264 0.347 protein_coding
## 2 lme A1BG sex 0.0796 0.814 0.954 protein_coding
## 3 lme A1BG (1 | ptID) 0 1 1 protein_coding
## 4 lme A2M condition -0.0940 0.605 0.682 protein_coding
## 5 lme A2M sex 1.40 0.0140 0.158 protein_coding
## 6 lme A2M (1 | ptID) 7.79 0.00525 0.0574 protein_coding
## 7 lme A2ML1 condition 0.791 0.0191 0.0356 protein_coding
## 8 lme A2ML1 sex -2.53 0.00000156 0.000388 protein_coding
## 9 lme A2ML1 (1 | ptID) 1.39 0.239 0.498 protein_coding
## 10 lme A4GALT condition 2.66 0.0000318 0.000108 protein_coding
## # … with 40,247 more rowsColor
Let’s start with a basic dotplot like before. We first filter the fit
data to just the variable condition so that we have 1 point
per gene.
fit$lme %>%
filter(variable == "condition") %>%
ggplot() +
aes(x = estimate, y = -log10(FDR)) +
geom_point()
Now let’s color the significant points. There is no “significant” variable in the data frame, so we first make one with mutate.
fit$lme %>%
filter(variable == "condition") %>%
mutate(significance = case_when(FDR < 0.01 ~ "FDR < 0.01", #<--- see change
TRUE ~ "NS")) %>%
ggplot() +
aes(x = estimate, y = -log10(FDR), color = significance) + #<--- see change
geom_point()
Or go further and color up vs down regulated genes differently.
fit$lme %>%
filter(variable == "condition") %>%
mutate(significance = case_when(FDR < 0.01 & estimate < 0 ~ "down", #<--- see change
FDR < 0.01 & estimate > 0 ~ "up",
TRUE ~ "NS")) %>%
ggplot() +
aes(x = estimate, y = -log10(FDR), color = significance) +
geom_point()
And we can specify the colors we want with another layer.
fit$lme %>%
filter(variable == "condition") %>%
mutate(significance = case_when(FDR < 0.01 & estimate < 0 ~ "down",
FDR < 0.01 & estimate > 0 ~ "up",
TRUE ~ "NS")) %>%
ggplot() +
aes(x = estimate, y = -log10(FDR), color = significance) +
geom_point() +
scale_color_manual(values = c("blue", "grey", "red")) #<--- see change
Lines
Since we’re defining significance at a cutoff, we can plot that
cutoff with a vertical line with geom_vline. Other similar
geoms are geom_hline for a horizontal line and
geom_abline for a line with slope a and
intercept b.
fit$lme %>%
filter(variable == "condition") %>%
mutate(significance = case_when(FDR < 0.01 & estimate < 0 ~ "down",
FDR < 0.01 & estimate > 0 ~ "up",
TRUE ~ "NS")) %>%
ggplot() +
aes(x = estimate, y = -log10(FDR), color = significance) +
geom_point() +
scale_color_manual(values = c("blue", "grey", "red")) +
geom_hline(yintercept = -log10(0.01)) #<--- see change
Scales
You might also wish to center the 0 estimate or otherwise alter the plots limits.
fit$lme %>%
filter(variable == "condition") %>%
mutate(significance = case_when(FDR < 0.01 & estimate < 0 ~ "down",
FDR < 0.01 & estimate > 0 ~ "up",
TRUE ~ "NS")) %>%
ggplot() +
aes(x = estimate, y = -log10(FDR), color = significance) +
geom_point() +
scale_color_manual(values = c("blue", "grey", "red")) +
geom_hline(yintercept = -log10(0.01)) +
lims(x = c(-9, 9)) #<--- see change
If you accidentally (or purposefully) remove data point with your limits, R will warn you!
fit$lme %>%
filter(variable == "condition") %>%
mutate(significance = case_when(FDR < 0.01 & estimate < 0 ~ "down",
FDR < 0.01 & estimate > 0 ~ "up",
TRUE ~ "NS")) %>%
ggplot() +
aes(x = estimate, y = -log10(FDR), color = significance) +
geom_point() +
scale_color_manual(values = c("blue", "grey", "red")) +
geom_hline(yintercept = -log10(0.01)) +
lims(x = c(-3, 3)) #<--- see change## Warning: Removed 165 rows containing missing values (geom_point).
Label points
We’re going to use a ggplot extension package to label our most significant genes. Let’s load that package now.
library(ggrepel)First, we determine our top 2 up and down regulated genes. There are
some new dplyr functions in here! See the comments for
information.
top_deg <- fit$lme %>%
filter(variable == "condition") %>%
mutate(direction = case_when(estimate < 0 ~ "down",
estimate > 0 ~ "up")) %>%
group_by(direction) %>%
slice_min(FDR, n=2) %>% # Keeps n number of minimum FDR values
pull(gene) # Extract the variable and turns the data into a vector
top_deg## [1] "IFNAR2" "NAIP" "TERF2IP" "ZBTB21"Then we add those labels to our fit data and use that variable in
geom_text_repel. In this case, the data removal warning is
expected as we’re not labeling the majority of points.
fit$lme %>%
filter(variable == "condition") %>%
mutate(significance = case_when(FDR < 0.01 & estimate < 0 ~ "down",
FDR < 0.01 & estimate > 0 ~ "up",
TRUE ~ "NS")) %>%
mutate(label = case_when(gene %in% top_deg ~ gene)) %>% #<--- see change
ggplot() +
aes(x = estimate, y = -log10(FDR), color = significance) +
geom_point() +
scale_color_manual(values = c("blue", "grey", "red")) +
geom_hline(yintercept = -log10(0.01)) +
lims(x = c(-9, 9)) +
geom_text_repel(aes(label = label)) #<--- see change## Warning: Removed 13415 rows containing missing values (geom_text_repel).
The text legend gets overlayed with the point legend, which looks odd. Let’s remove one of them.
fit$lme %>%
filter(variable == "condition") %>%
mutate(significance = case_when(FDR < 0.01 & estimate < 0 ~ "down",
FDR < 0.01 & estimate > 0 ~ "up",
TRUE ~ "NS")) %>%
mutate(label = case_when(gene %in% top_deg ~ gene)) %>%
ggplot() +
aes(x = estimate, y = -log10(FDR), color = significance) +
geom_point() +
scale_color_manual(values = c("blue", "grey", "red")) +
geom_hline(yintercept = -log10(0.01)) +
lims(x = c(-9, 9)) +
geom_text_repel(aes(label = label), show.legend = FALSE) #<--- see change## Warning: Removed 13415 rows containing missing values (geom_text_repel).
And a small edit if you don’t want the text to be colored. When you put an aesthetic in a specific geom, it only applies that change to that one geom.
fit$lme %>%
filter(variable == "condition") %>%
mutate(significance = case_when(FDR < 0.01 & estimate < 0 ~ "down",
FDR < 0.01 & estimate > 0 ~ "up",
TRUE ~ "NS")) %>%
mutate(label = case_when(gene %in% top_deg ~ gene)) %>%
ggplot() +
aes(x = estimate, y = -log10(FDR)) + #<--- see change
geom_point(aes(color = significance)) + # <--- see change
scale_color_manual(values = c("blue", "grey", "red")) +
geom_hline(yintercept = -log10(0.01)) +
lims(x = c(-9, 9)) +
geom_text_repel(aes(label = label), show.legend = FALSE)## Warning: Removed 13415 rows containing missing values (geom_text_repel).
Theme
To round it out, let’s look at theme_bw, my favorite one
for volcano plots.
fit$lme %>%
filter(variable == "condition") %>%
mutate(significance = case_when(FDR < 0.01 & estimate < 0 ~ "down",
FDR < 0.01 & estimate > 0 ~ "up",
TRUE ~ "NS")) %>%
mutate(label = case_when(gene %in% top_deg ~ gene)) %>%
ggplot() +
aes(x = estimate, y = -log10(FDR)) +
geom_point(aes(color = significance)) +
scale_color_manual(values = c("blue", "grey", "red")) +
geom_hline(yintercept = -log10(0.01)) +
lims(x = c(-9, 9)) +
geom_text_repel(aes(label = label)) +
theme_bw() #<--- see change## Warning: Removed 13415 rows containing missing values (geom_text_repel).
There you have it! Did you think you’d reach a 11 function plot so soon?
Facets
With big data, we often want to plot more than one thing at once.
Facets allow you to make similar plots across multiple subsets of data
with only 1 additional line of code. Here, we filter for 2
variable and then facet those into 2 plots.
fit$lme %>%
filter(variable %in% c("condition", "sex")) %>% #<--- see change
mutate(significance = case_when(FDR < 0.01 & estimate < 0 ~ "down",
FDR < 0.01 & estimate > 0 ~ "up",
TRUE ~ "NS")) %>%
mutate(label = case_when(gene %in% top_deg ~ gene)) %>%
ggplot() +
aes(x = estimate, y = -log10(FDR)) +
geom_point(aes(color = significance)) +
scale_color_manual(values = c("blue", "grey", "red")) +
geom_hline(yintercept = -log10(0.01)) +
#<--- see change No lims used
geom_text_repel(aes(label = label)) +
theme_bw() +
facet_wrap(~ variable) #<--- see change## Warning: Removed 26830 rows containing missing values (geom_text_repel).
Since the a and y scales differ for the two variable, we allow them
to differ with scales. This is my we removed the
lims function earlier. scales can be set to
“free_x”, “free_y”, or “free” for both x and y.
fit$lme %>%
filter(variable %in% c("condition", "sex")) %>%
mutate(significance = case_when(FDR < 0.01 & estimate < 0 ~ "down",
FDR < 0.01 & estimate > 0 ~ "up",
TRUE ~ "NS")) %>%
mutate(label = case_when(gene %in% top_deg ~ gene)) %>%
ggplot() +
aes(x = estimate, y = -log10(FDR)) +
geom_point(aes(color = significance)) +
scale_color_manual(values = c("blue", "grey", "red")) +
geom_hline(yintercept = -log10(0.01)) +
geom_text_repel(aes(label = label)) +
theme_bw() +
facet_wrap(~ variable, scales = "free") #<--- see change## Warning: Removed 26830 rows containing missing values (geom_text_repel).
Save plots
Once you have your plot as you desire, save it to your computer for
use in presentations, publications, Twitter, and more!
ggsave allows you to save in a bunch of formats (see the
help page with ?ggsave).
Since we’ve only been printing our plots to the plot pane in RStudio,
we need to first save a plot to the environment. We’ll use the single
condtion variable volcano plot.
p1 <- fit$lme %>% #<--- see change
filter(variable == "condition") %>%
mutate(significance = case_when(FDR < 0.01 & estimate < 0 ~ "down",
FDR < 0.01 & estimate > 0 ~ "up",
TRUE ~ "NS")) %>%
mutate(label = case_when(gene %in% top_deg ~ gene)) %>%
ggplot() +
aes(x = estimate, y = -log10(FDR)) +
geom_point(aes(color = significance)) +
scale_color_manual(values = c("blue", "grey", "red")) +
geom_hline(yintercept = -log10(0.01)) +
geom_text_repel(aes(label = label)) +
theme_bw()
p1## Warning: Removed 13415 rows containing missing values (geom_text_repel).
Now we can save the object. My most commonly used plot formats are
.pdf for publication (vectorized PDF that can be edited in
Illustrator/Inkscape) and .png for all other uses (not
vectorized, smaller files). You specify the file type in the filename
and give dimensions in inches (default). Without dimensions, R used
whatever the ratio of you plot pane currently is (not
reproducible!).
ggsave(filename = "images/volcano_plot.png", plot = p1, width = 3, height = 4)## Warning: Removed 13415 rows containing missing values (geom_text_repel).ggsave(filename = "images/volcano_plot.pdf", plot = p1, width = 3, height = 4)## Warning: Removed 13415 rows containing missing values (geom_text_repel).Commented plot scripts
Here are each of the final plots fulled comments for your reference.
ggplot(dat.pca.meta) + # Make ggplot from data in dat.pca.meta
aes(x = PC1, y = PC2, # Plot PCs 1 and 2 along x and y
color = condition) + # Color by condition
geom_point() + # Plot points for each value
labs(x = "PC1 (24%)", y = "PC2 (11%)", # Relabel axes
color = "Infection") + # Relabel legend
coord_fixed() + # Force xy scales to be equal
theme_classic() # Set themefit$lme %>%
# Keep results for only the condition variable
filter(variable == "condition") %>%
# Create variable for significant up and down genes
mutate(significance = case_when(FDR < 0.01 & estimate < 0 ~ "down",
FDR < 0.01 & estimate > 0 ~ "up",
TRUE ~ "NS")) %>%
# Create variable for gene name labels
mutate(label = case_when(gene %in% top_deg ~ gene)) %>%
ggplot() + # initiate ggplot plot
aes(x = estimate, y = -log10(FDR)) + # plot estimate and FDR
geom_point(aes(color = significance)) + # plot points, color by significance
scale_color_manual(values = c("blue", "grey", "red")) + # recolor points
geom_hline(yintercept = -log10(0.01)) + # add FDR cuttof line
lims(x = c(-9, 9)) + # set axis limits
geom_text_repel(aes(label = label)) + # label points
theme_bw() # set theme ... +
facet_wrap(~ variable) # facet by variableExercises: ggplot
- Using the combined expression and metadata you made in the last
session, plot the expression of one gene in
MediavsMtb. The plot type is up to you! As a reminder, here is how we made the combined data.
full_data <- as.data.frame(dat$E) %>%
rownames_to_column("gene") %>%
pivot_longer(-gene, names_to = "libID", values_to = "log2CPM") %>%
inner_join(dat$targets)## Joining, by = "libID"full_data## # A tibble: 268,380 × 15
## gene libID log2CPM group lib.s…¹ norm.…² ptID condi…³ age_dys sex ptID_…⁴
## <chr> <chr> <dbl> <fct> <dbl> <dbl> <chr> <chr> <dbl> <chr> <chr>
## 1 A1BG pt01… 1.62 1 8.30e6 1.09 pt01 Media 12410 M pt00001
## 2 A1BG pt01… 1.72 1 5.72e6 0.819 pt01 Mtb 12410 M pt00001
## 3 A1BG pt02… 1.09 1 8.09e6 1.12 pt02 Media 12775 M pt00002
## 4 A1BG pt02… 1.10 1 5.28e6 0.912 pt02 Mtb 12775 M pt00002
## 5 A1BG pt03… 0.874 1 5.18e6 1.05 pt03 Media 11315 M pt00003
## 6 A1BG pt03… 1.71 1 5.23e6 0.958 pt03 Mtb 11315 M pt00003
## 7 A1BG pt04… 1.57 1 9.00e6 1.13 pt04 Media 8395 M pt00004
## 8 A1BG pt04… -0.264 1 5.40e6 0.855 pt04 Mtb 8395 M pt00004
## 9 A1BG pt05… 1.63 1 1.26e7 1.07 pt05 Media 7300 M pt00005
## 10 A1BG pt05… -0.0305 1 1.01e7 0.930 pt05 Mtb 7300 M pt00005
## # … with 268,370 more rows, 4 more variables: RNAseq <lgl>, methylation <lgl>,
## # total_seq <dbl>, sample.weights <dbl>, and abbreviated variable names
## # ¹lib.size, ²norm.factors, ³condition, ⁴ptID_old- Color the PCA by the numeric variable
total_seq. Explore the help pages forscale_color_continuous,scale_color_gradient, orscale_color_gradient2and recolor your PCA to more easily see differences. As a reminder, here is the base PCA plot.
ggplot(dat.pca.meta) +
aes(x = PC1, y = PC2) +
geom_point() +
labs(x = "PC1 (24%)", y = "PC2 (11%)",
color = "Infection") +
theme_classic()
- Using
geom_text_repel, label the points in the PCA by theirlibID