Introduction to ggplot2

Objectives

1. Learn the ggplot2 syntax.

2. Build a ggplot2 general template.

By the end of the course, students should be able to create simple, pretty, and effective figures.

Data Visualization in the tidyverse

In lesson 3, we learned how to read and save excel spreadsheet data to a R object using the tidyverse package readxl. Today we will use some example data from an excel spreadsheet to learn the basics of ggplot2, a tidyverse core package.

#data import from excel
exdata<-readxl::read_xlsx("./data/RNASeq_totalcounts_vs_totaltrans.xlsx",
                        1,.name_repair = "universal", skip=3)

## New names:
## • `Sample Name` -> `Sample.Name`
## • `Number of Transcripts` -> `Number.of.Transcripts`
## • `Total Counts` -> `Total.Counts`

exdata

## # A tibble: 8 × 4
##   Sample.Name Treatment     Number.of.Transcripts Total.Counts
##   <chr>       <chr>                         <dbl>        <dbl>
## 1 GSM1275863  Dexamethasone                 10768     18783120
## 2 GSM1275867  Dexamethasone                 10051     15144524
## 3 GSM1275871  Dexamethasone                 11658     30776089
## 4 GSM1275875  Dexamethasone                 10900     21135511
## 5 GSM1275862  None                          11177     20608402
## 6 GSM1275866  None                          11526     25311320
## 7 GSM1275870  None                          11425     24411867
## 8 GSM1275874  None                          11000     19094104

These data include total transcript read counts summed by sample and the total number of transcripts recovered by sample that had at least 100 reads. These data derive from a bulk RNAseq experiment described by Himes et al. (2014) and introduced in lesson 3. As a reminder, the authors "characterized transcriptomic changes in four primary human ASM cell lines that were treated with dexamethasone," a common therapy for asthma. Each cell line included a treated and untreated negative control resulting in a total sample size of 8.

Plotting with Excel

We could plot this data in Excel and get something like this:

While this isn't bad, it took an unnecessary amount of time to create, and there weren't a lot of options for customization.

Why ggplot2?

Outside of base R plotting, one of the most popular packages used to generate graphics in R is ggplot2, which is associated with a family of packages collectively known as the tidyverse. GGplot2 allows the user to create informative plots quickly by using a 'grammar of graphics' implementation, which is described as "a coherent system for describing and building graphs" R4DS. We will see this in action shortly. The power of this package is that plots are built in layers and few changes to the code result in very different outcomes. This makes it easy to reuse parts of the code for very different figures.

Being a part of the tidyverse collection, ggplot2 works best with data organized so that individual observations are in rows and variables are in columns.

Getting started with ggplot2

To begin plotting, we need to load our ggplot2 library. Package libraries must be loaded every time you open and use R. If you haven't yet installed the ggplot2 package on your local machine, you will need to do that using install.packages("ggplot2").

#load the ggplot2 library; you could also load library(tidyverse)
library(ggplot2)

Getting help

The R community is extensive and getting help is now easier than ever with a simple web search. If you can't figure out how to plot something, give a quick web search a try. Great resources include internet tutorials, R bookdowns, and stackoverflow. You should also use the help features within RStudio to get help on specific functions or to find vignettes. Try entering ggplot2 in the help search bar in the lower right panel under the Help tab.

The ggplot2 template

The following represents the basic ggplot2 template.

ggplot(data = <DATA>) + 
  <GEOM_FUNCTION>(mapping = aes(<MAPPINGS>))

Let's see this template in practice.

What is the relationship between total transcript sums per sample and the number of recovered transcripts per sample?

#let's plot our data
ggplot(data=exdata) + 
  geom_point(aes(x=Number.of.Transcripts, y = Total.Counts)) 

We can easily see that there is a relationship between the number of transcripts per sample and the total transcripts recovered per sample. ggplot2 default parameters are great for exploratory data analysis. But, with only a few tweaks, we can make some beautiful, publishable figures.

Let's take a closer look at the above code
The first step in creating this plot was initializing the ggplot object using the function ggplot(). Remember, we can look further for help using ?ggplot(). The function ggplot() takes data, mapping, and further arguments. However, none of this needs to actually be provided at the initialization phase, which creates the coordinate system from which we build our plot. But, typically, you should at least call the data at this point.

The data we called was from the data frame exdata, which we created above. Next, we provided a geom function (geom_point()), which created a scatter plot. This scatter plot required mapping information, which we provided for the x and y axes. More on this in a moment.

Let's break down the individual components of the code.

#What does running ggplot() do?
ggplot(data=exdata)

#What about just running a geom function?
geom_point(data=exdata,aes(x=Number.of.Transcripts, y = Total.Counts)) 

## mapping: x = ~Number.of.Transcripts, y = ~Total.Counts 
## geom_point: na.rm = FALSE
## stat_identity: na.rm = FALSE
## position_identity

#what about this
ggplot() +
geom_point(data=exdata,aes(x=Number.of.Transcripts, y = Total.Counts))

Geom functions

A geom is the geometrical object that a plot uses to represent data. People often describe plots by the type of geom that the plot uses. --- R4DS

There are multiple geom functions that change the basic plot type or the plot representation. We can create scatter plots (geom_point()), line plots (geom_line(),geom_path()), bar plots (geom_bar(), geom_col()), line modeled to fitted data (geom_smooth()), heat maps (geom_tile()), geographic maps (geom_polygon()), etc.

ggplot2 provides over 40 geoms, and extension packages provide even more (see https://exts.ggplot2.tidyverse.org/gallery/ for a sampling). The best way to get a comprehensive overview is the ggplot2 cheatsheet, which you can find at http://rstudio.com/resources/cheatsheets. --- R4DS

You can also see a number of options pop up when you type geom into the console, or you can look up the ggplot2 documentation in the help tab.

We can see how easy it is to change the way the data is plotted. Let's plot the same data using geom_line().

ggplot(data=exdata) + 
  geom_line(aes(x=Number.of.Transcripts, y = Total.Counts)) 

Mapping and aesthetics (aes())

The geom functions require a mapping argument. The mapping argument includes the aes() function, which "describes how variables in the data are mapped to visual properties (aesthetics) of geoms" (ggplot2 R Documentation). If not included it will be inherited from the ggplot() function.

An aesthetic is a visual property of the objects in your plot.---R4DS

Mapping aesthetics include some of the following:
1. the x and y data arguments
2. shapes
3. color
4. fill
5. size
6. linetype
7. alpha

This is not an all encompassing list.

Let's return to our plot above. Is there a relationship between treatment ("dex") and the number of transcripts or total counts?

#adding the color argument to our mapping aesthetic
ggplot(exdata) +
  geom_point(aes(x=Number.of.Transcripts, y = Total.Counts, 
                 color=Treatment))

There is potentially a relationship. ASM cells treated with dexamethasone in general have lower total numbers of transcripts and lower total counts.

Notice how we changed the color of our points to represent a variable, in this case. To do this, we set color equal to 'Treatment' within the aes() function. This mapped our aesthetic, color, to a variable we were interested in exploring.

#map the shape aesthetic to the variable "dex"
#use the color purple across all points (NOT mapped to a variable)
ggplot(exdata) + 
  geom_point(aes(x=Number.of.Transcripts, y = Total.Counts, 
                shape=Treatment), color="purple") 

We can also see from this that 'Treatment' could be mapped to other aesthetics. In the above example, we see it mapped to shape rather than color. By default, ggplot2 will only map six shapes at a time, and if your number of categories goes beyond 6, the remaining groups will go unmapped. This is by design because it is hard to discriminate between more than six shapes at any given moment. This is a clue from ggplot2 that you should choose a different aesthetic to map to your variable. However, if you choose to ignore this functionality, you can manually assign more than six shapes.

We could have just as easily mapped it to alpha, which adds a gradient to the point visibility by category, or we could map it to size. There are multiple options, so feel free to explore a little with your plots.

R objects can also store figures

As we have discussed, R objects are used to store things created in R to memory. This includes plots.

scatter_plot<-ggplot(exdata) + 
  geom_point(aes(x=Number.of.Transcripts, y = Total.Counts,
                 color=Treatment)) 

scatter_plot

We can add additional layers directly to our object. We will see how this works by defining some colors for our 'dex' variable.

Colors

ggplot2 will automatically assign colors to the categories in our data. Colors are assigned to the fill and color aesthetics in aes(). We can change the default colors by providing an additional layer to our figure. To change the color, we use the scale_color functions: scale_color_manual(), scale_color_brewer(), scale_color_grey(), etc. We can also change the name of the color labels in the legend using the labels argument of these functions

scatter_plot +
  scale_color_manual(values=c("red","black"),
                     labels=c('treated','untreated'))

scatter_plot +
  scale_color_grey()

scatter_plot +
  scale_color_brewer(palette = "Paired")

Similarly, if we want to change the fill, we would use the scale_fill options. To apply scale_fill to shape, we will have to alter the shapes, as only some shapes take a fill argument. Refer to the shapes in the red box in the figure below.

ggplot(data=exdata) + 
  geom_point(aes(x=Number.of.Transcripts, y = Total.Counts,fill=Treatment),
             shape=21,size=2) + #increase size and change points
  scale_fill_manual(values=c("purple", "yellow"))

There are a number of ways to specify the color argument including by name, number, and hex code. Here is a great resource from the R Graph Gallery for assigning colors in R.

There are also a number of complementary packages in R that expand our color options. One of my favorites is viridis, which provides colorblind friendly palettes. randomcoloR is a great package if you need a large number of unique colors.

library(viridis) #Remember to load installed packages before use

## Loading required package: viridisLite

scatter_plot + scale_color_viridis(discrete=TRUE, option="viridis")

Paletteer contains a comprehensive set of color palettes, if you want to load the palettes from multiple packages all at once. See the Github page for details.

Returning to our grammar of graphics

Remember, to create a plot all you you need are the data, geom_function(s), and mapping arguments.

However, there are additional components that can be added to our core components to enable us to generate even more diverse plot types.

Our grammar of graphics:

• one or more datasets,

• one or more geometric objects that serve as the visual representations of the data, – for instance, points, lines, rectangles, contours,

• descriptions of how the variables in the data are mapped to visual properties (aesthetics) of the geometric objects, and an associated scale (e. g., linear, logarithmic, rank),

• a facet specification, i.e. the use of multiple similar subplots to look at subsets of the same data,

• one or more coordinate systems,

• optional parameters that affect the layout and rendering, such text size, font and alignment, legend positions,

• statistical summarization rules

---(Holmes and Huber, 2021)

Our template can therefore be expanded to include these additional components:

ggplot(data = <DATA>) + 
  <GEOM_FUNCTION>(
     mapping = aes(<MAPPINGS>),
     stat = <STAT>
  ) +
  <FACET_FUNCTION> +
  <COORDINATE SYSTEM> +
  <THEME>

Facets

A way to add variables to a plot beyond mapping them to an aesthetic is to use facets or subplots. There are two primary functions to add facets, facet_wrap() and facet_grid(). If faceting by a single variable, use facet_wrap(). If multiple variables, use facet_grid(). The first argument of either function is a formula, with variables separated by a ~ (See below). Variables must be discrete (not continuous).

Let's return to the airway count data to see how facets are useful. Here, we are going to compare scaled and unscaled count data using a density plot.

A density plot shows the distribution of a numeric variable. --- R Graph Gallery

In our example data, density_data, the gene counts were scaled to account for technical and composition differences using the trimmed mean of M values (TMM) from EdgeR (Robinson and Oshlack 2010), but non-normalized values remained for comparison. Thus, we can compare scaled vs unscaled counts by sample using faceting.

#density plot
#let's grab the data and take a look
density_data<-read.csv("./data/density_data.csv",
                       stringsAsFactors=TRUE)

head(density_data)

##           feature sample SampleName   cell   dex albut        Run avgLength
## 1 ENSG00000000003    508 GSM1275862 N61311 untrt untrt SRR1039508       126
## 2 ENSG00000000003    508 GSM1275862 N61311 untrt untrt SRR1039508       126
## 3 ENSG00000000419    508 GSM1275862 N61311 untrt untrt SRR1039508       126
## 4 ENSG00000000419    508 GSM1275862 N61311 untrt untrt SRR1039508       126
## 5 ENSG00000000457    508 GSM1275862 N61311 untrt untrt SRR1039508       126
## 6 ENSG00000000457    508 GSM1275862 N61311 untrt untrt SRR1039508       126
##   Experiment    Sample    BioSample transcript ref_genome .abundant      TMM
## 1  SRX384345 SRS508568 SAMN02422669     TSPAN6       hg38      TRUE 1.055278
## 2  SRX384345 SRS508568 SAMN02422669     TSPAN6       hg38      TRUE 1.055278
## 3  SRX384345 SRS508568 SAMN02422669       DPM1       hg38      TRUE 1.055278
## 4  SRX384345 SRS508568 SAMN02422669       DPM1       hg38      TRUE 1.055278
## 5  SRX384345 SRS508568 SAMN02422669      SCYL3       hg38      TRUE 1.055278
## 6  SRX384345 SRS508568 SAMN02422669      SCYL3       hg38      TRUE 1.055278
##   multiplier        source abundance
## 1   1.415149        counts  679.0000
## 2   1.415149 counts_scaled  960.8864
## 3   1.415149        counts  467.0000
## 4   1.415149 counts_scaled  660.8748
## 5   1.415149        counts  260.0000
## 6   1.415149 counts_scaled  367.9388

#plot 
ggplot(data= density_data)+
  aes(x=abundance, 
             color=SampleName)+ #initialize ggplot
  geom_density() + #call density plot geom
  facet_wrap(~source) + #use facet_wrap; see ~source
  scale_x_log10()#scales the x axis using a base-10 log transformation  

## Warning: Transformation introduced infinite values in continuous x-axis

## Warning: Removed 140 rows containing non-finite values (`stat_density()`).

The distributions of sample counts did not differ greatly between samples before scaling, but regardless, we can see that the distributions are more similar after scaling.

Here, faceting allowed us to visualize multiple features of our data. We were able to see count distributions by sample as well as normalized vs non-normalized counts.

Note the help options with ?facet_wrap(). How would we make our plot facets vertical rather than horizontal?

ggplot(data= density_data)+ #initialize ggplot
  geom_density(aes(x=abundance, 
             color=SampleName)) + #call density plot geom
  facet_wrap(~source, ncol=1) + #use the ncol argument
  scale_x_log10()

## Warning: Transformation introduced infinite values in continuous x-axis

## Warning: Removed 140 rows containing non-finite values (`stat_density()`).

We could plot each sample individually using facet_grid()

ggplot(data= density_data)+ #initialize ggplot
  geom_density(aes(x=abundance, 
             color=SampleName)) + #call density plot geom
  facet_grid(as.factor(sample)~source) + # formula is sample ~ source
  scale_x_log10()

## Warning: Transformation introduced infinite values in continuous x-axis

## Warning: Removed 140 rows containing non-finite values (`stat_density()`).

Labels, legends, scales, and themes

How do we ultimately get our figures to a publishable state? The bread and butter of pretty plots really falls to the additional non-data layers of our ggplot2 code. These layers will include code to label the axes, scale the axes, and customize the legends and theme.

Let's make our original figure publishable

ggplot(exdata) + 
  geom_point(aes(x=Number.of.Transcripts, y = Total.Counts,
                 fill=Treatment),
             shape=21,size=3) + 
    #can change labels of fill levels along with colors
  scale_fill_manual(values=c("purple", "yellow"), 
                    labels=c('treated','untreated'))+ 

  labs(x="Recovered transcripts per sample", 
       y="Total sequences per sample", fill="Treatment")+
  scale_y_continuous(trans="log10") + #log transform the y axis
  theme_bw() #add a complete theme black / white

Saving plots (ggsave())

Finally, we have a quality plot ready to publish. The next step is to save our plot to a file. The easiest way to do this with ggplot2 is ggsave(). This function will save the last plot that you displayed by default. Look at the function parameters using ?ggsave().

ggsave("Plot1.png",width=5.5,height=3.5,units="in",dpi=300)

Resource list

1. ggplot2 cheatsheet
2. The R Graph Gallery
3. The R Graphics Cookbook

Acknowledgements

Material from this lesson was adapted from Chapter 3 of R for Data Science and from a 2021 workshop entitled Introduction to Tidy Transciptomics by Maria Doyle and Stefano Mangiola.