Prerequisites

Required R packages: ggpubr to easily create ggplot2-based publication ready plots.

Install from CRAN:

install.packages("ggpubr")

Or, install the latest developmental version from GitHub as follow:

if(!require(devtools)) install.packages("devtools")
devtools::install_github("kassambara/ggpubr")

Load ggpubr:

library(ggpubr)

Basic plots

Demo data set:

df <- ToothGrowth
df$dose <- as.factor(df$dose)
head(df)

##    len supp dose
## 1  4.2   VC  0.5
## 2 11.5   VC  0.5
## 3  7.3   VC  0.5
## 4  5.8   VC  0.5
## 5  6.4   VC  0.5
## 6 10.0   VC  0.5

Plot:

p <- ggdensity(df, x = "len", fill = "dose", 
               palette = "jco", 
               ggtheme = theme_light(), legend = "top")
p

Facet by one grouping variables

Divide by the levels of the supp variable in the horizontal direction:

facet(p, facet.by = "supp")

Divide by the levels of supp in the vertical direction –> use ncol = 1:

facet(p, facet.by = "supp", ncol = 1)

Facet by two grouping variables

The data can be split up by one or two variables that vary on the horizontal and/or vertical direction.

For example in facet.by = c(“supp”, “dose”):

• “supp”, the first variable, will be displayed in vertical direction
• “dose”, the second variable, will be displayed in horizontal direction.

# Divide with "supp" vertical, "dose" horizontal
facet(p, facet.by = c("supp", "dose"),
      short.panel.labs = FALSE)

Modifying panel label appearance

Additional arguments are available to customize the appearance of panel labels (see ?facet). These include:

# Divide with "supp" vertical, "dose" horizontal
facet(p, facet.by = c("supp", "dose"),
       panel.labs = list(
         supp = c("Orange Juice", "Vitamin C"),
         dose = c("D0.5", "D1", "D2")
         ),
       panel.labs.background = list(color = "steelblue", fill = "steelblue", size = 0.5),
       panel.labs.font = list(color = "white"),
       panel.labs.font.x = list(angle = 45, color = "white")
      )

Prerequisites

Install ggpubr (version >= 0.1.4), for easily creating ggplot2-based publication ready plots.

Install from CRAN:

install.packages("ggpubr")

Or, install the latest developmental version from GitHub as follow:

if(!require(devtools)) install.packages("devtools")
devtools::install_github("kassambara/ggpubr")

Load ggpubr:

library(ggpubr)

Basic plots

We start by creating a basic box plot colored by groups. To add the panel border line, we’ll use the helper function border() [in ggpubr].

# Basic plot
p <- ggboxplot(ToothGrowth, x = "dose", y = "len",
               color = "dose")
p

# Add grids
p + grids(linetype = "dashed")

# Add panel border line
p + border("black")
  
# Change background color
p + bgcolor("#BFD5E3") +
  border("#BFD5E3") 

Change titles and axis labels

Change plot titles and labels as follow:

# Change titles and axis labels
p2 <- ggpar(p, 
            title = "Box Plot created with ggpubr",
            subtitle = "Length by dose",
            caption = "Source: ggpubr",
            xlab ="Dose (mg)", 
            ylab = "Teeth length",
            legend.title = "Dose (mg)")
p2

Change the font/appearance of plot titles and labels. Use ggpar():

ggpar(p2, 
      font.title = c(14, "bold.italic", "red"),
      font.subtitle = c(10,  "orange"),
      font.caption = c(10,  "orange"),
      font.x = c(14,  "blue"),
      font.y = c(14,  "#993333")
      )

or, equivalently, use font():

p2 +
 font("title", size = 14, color = "red", face = "bold.italic")+
 font("subtitle", size = 10, color = "orange")+
 font("caption", size = 10, color = "orange")+
 font("xlab", size = 12, color = "blue")+
 font("ylab", size = 12, color = "#993333")

Note that, you can change simultaneously titles/labels text and appearance at once using the function ggpar() as follow:

# Change title texts and fonts
# line break: \n
ggpar(p, title = "Plot of length \n by dose",
      xlab ="Dose (mg)", ylab = "Teeth length",
      legend.title = "Dose (mg)",
      font.title = c(14,"bold.italic", "red"),
      font.x = c(14, "bold", "#2E9FDF"),
      font.y = c(14, "bold", "#E7B800"))

Change legend position & appearance

ggpar(p,
      legend = "right", legend.title = "Dose (mg)") + 
  font("legend.title", color = "blue", face = "bold")+ 
  font("legend.text", color = "red")

Change color palettes

# Use custom color palette
ggpar(p, palette = c("#00AFBB", "#E7B800", "#FC4E07"))

# Use brewer palette. 
# Type RColorBrewer::display.brewer.all(), to see possible palettes
ggpar(p, palette = "Dark2" )

# Use grey palette
ggpar(p, palette = "grey")
   
# Use scientific journal palette from ggsci package
# Allowed values: "npg", "aaas", "lancet", "jco", 
#   "ucscgb", "uchicago", "simpsons" and "rickandmorty".
ggpar(p, palette = "npg") # nature

# jco color palette
p + color_palette("jco")

# Custom color
p + color_palette(c("#00AFBB", "#E7B800", "#FC4E07"))

p3 <- ggscatter(mtcars, x = "wt", y = "mpg", color = "mpg",
                size = 2)

# Use one custom color
p3 + gradient_color("red")

# Two colors
p3 + gradient_color(c("blue",  "red"))

# Three colors
p3 + gradient_color(c("blue", "white", "red"))

# Use RColorBrewer palette
p3 + gradient_color("RdYlBu")

p4 <- ggscatter(mtcars, x = "wt", y = "mpg", fill = "mpg",
                size = 4, shape = 21)
p4 + gradient_fill(c("blue", "white", "red"))

Change axis limits and scales

The following arguments can be used in ggpar():

# Change y axis limits
ggpar(p, ylim = c(0, 50))

# Change y axis scale to log2
ggpar(p, yscale = "log2")

# Format axis scale
ggpar(p, yscale = "log2", format.scale = TRUE)

Alternatively, you can also use the function xscale() and yscale() [in ggpubr], as follow:

p + yscale("log2", .format = TRUE)

Customize axis text and ticks

# Change the font of x and y axis texts.
# Rotate x and y texts, angle = 45
p + 
  font("xy.text", size = 12, color = "blue", face = "bold") +
  rotate_x_text(45)+       
  rotate_y_text(45)

# remove ticks and axis texts
p + rremove("ticks")+
  rremove("axis.text")

Rotate a plot

# Horizontal box plot
p + rotate()

Change themes

The default theme in ggpubr is theme_pubr(), for publication ready theme.

The argument ggtheme can be used in any ggpubr plotting functions to change the plot theme.

Allowed values include ggplot2 official themes: theme_gray(), theme_bw(), theme_minimal(), theme_classic(), theme_void(), etc. It’s also possible to use the function “+” to add a theme.

# Gray theme
p + theme_gray()

# Black and white theme
p + theme_bw()

# Theme light
p + theme_light()

# Minimal theme
p + theme_minimal()

# Empty theme
p + theme_void()

Remove ggplot components

The function rremove() [in ggpubr] can be used to remove a specific component from a ggplot.

Usage:

rremove(object)

Object: character string specifying the plot components. Allowed values include:

• “grid” for both x and y grids
• “x.grid” for x axis grids
• “y.grid” for y axis grids
• “axis” for both x and y axes
• “x.axis” for x axis
• “y.axis” for y axis
• “xlab”, or “x.title” for x axis label
• “ylab”, or “y.title” for y axis label
• “xylab”, “xy.title” or “axis.title” for both x and y axis labels
• “x.text” for x axis texts (x axis tick labels)
• “y.text” for y axis texts (y axis tick labels)
• “xy.text” or “axis.text” for both x and y axis texts
• “ticks” for both x and y ticks
• “x.ticks” for x ticks
• “y.ticks” for y ticks
• “legend.title” for the legend title
• “legend” for the legend

Examples:

# Basic plot
p <- ggboxplot(ToothGrowth, x = "dose", y = "len",
  ggtheme = theme_gray())
p

# Remove all grids
p + rremove("grid")

# Remove only x grids
p + rremove("x.grid")

Prerequisites

if(!require(devtools)) install.packages("devtools")
devtools::install_github("kassambara/ggpubr")

install.packages("ggpubr")

library(ggpubr)

# ToothGrowth
data("ToothGrowth")
head(ToothGrowth)

##    len supp dose
## 1  4.2   VC  0.5
## 2 11.5   VC  0.5
## 3  7.3   VC  0.5
## 4  5.8   VC  0.5
## 5  6.4   VC  0.5
## 6 10.0   VC  0.5

# mtcars 
data("mtcars")
mtcars$name <- rownames(mtcars)
mtcars$cyl <- as.factor(mtcars$cyl)
head(mtcars[, c("name", "wt", "mpg", "cyl")])

##                                name   wt  mpg cyl
## Mazda RX4                 Mazda RX4 2.62 21.0   6
## Mazda RX4 Wag         Mazda RX4 Wag 2.88 21.0   6
## Datsun 710               Datsun 710 2.32 22.8   4
## Hornet 4 Drive       Hornet 4 Drive 3.21 21.4   6
## Hornet Sportabout Hornet Sportabout 3.44 18.7   8
## Valiant                     Valiant 3.46 18.1   6

Create some plots

Here, we’ll use ggplot2-based plotting functions available in ggpubr. You can use any ggplot2 functions to create the plots that you want for arranging them later.

We’ll start by creating 4 different plots:

• Box plots and dot plots using the ToothGrowth data set
• Bar plots and scatter plots using the mtcars data set

You’ll learn how to combine these plots in the next sections using specific functions.

• Create a box plot and a dot plot:

# Box plot (bp)
bxp <- ggboxplot(ToothGrowth, x = "dose", y = "len",
                 color = "dose", palette = "jco")
bxp

# Dot plot (dp)
dp <- ggdotplot(ToothGrowth, x = "dose", y = "len",
                 color = "dose", palette = "jco", binwidth = 1)
dp

• Create an ordered bar plot and a scatter plot:

Create ordered bar plots. Change the fill color by the grouping variable “cyl”. Sorting will be done globally, but not by groups.

# Bar plot (bp)
bp <- ggbarplot(mtcars, x = "name", y = "mpg",
          fill = "cyl",               # change fill color by cyl
          color = "white",            # Set bar border colors to white
          palette = "jco",            # jco journal color palett. see ?ggpar
          sort.val = "asc",           # Sort the value in ascending order
          sort.by.groups = TRUE,      # Sort inside each group
          x.text.angle = 90           # Rotate vertically x axis texts
          )
bp + font("x.text", size = 8)

# Scatter plots (sp)
sp <- ggscatter(mtcars, x = "wt", y = "mpg",
                add = "reg.line",               # Add regression line
                conf.int = TRUE,                # Add confidence interval
                color = "cyl", palette = "jco", # Color by groups "cyl"
                shape = "cyl"                   # Change point shape by groups "cyl"
                )+
  stat_cor(aes(color = cyl), label.x = 3)       # Add correlation coefficient
sp

Arrange on one page

To arrange multiple ggplots on one single page, we’ll use the function ggarrange()[in ggpubr], which is a wrapper around the function plot_grid() [in cowplot package]. Compared to the standard function plot_grid(), ggarange() can arrange multiple ggplots over multiple pages.

ggarrange(bxp, dp, bp + rremove("x.text"), 
          labels = c("A", "B", "C"),
          ncol = 2, nrow = 2)

Alternatively, you can also use the function plot_grid() [in cowplot]:

library("cowplot")
plot_grid(bxp, dp, bp + rremove("x.text"), 
          labels = c("A", "B", "C"),
          ncol = 2, nrow = 2)

or, the function grid.arrange() [in gridExtra]:

library("gridExtra")
grid.arrange(bxp, dp, bp + rremove("x.text"), 
             ncol = 2, nrow = 2)

Annotate the arranged figure

R function: annotate_figure() [in ggpubr].

figure <- ggarrange(sp, bp + font("x.text", size = 10),
                    ncol = 1, nrow = 2)

annotate_figure(figure,
                top = text_grob("Visualizing mpg", color = "red", face = "bold", size = 14),
                bottom = text_grob("Data source: \n mtcars data set", color = "blue",
                                   hjust = 1, x = 1, face = "italic", size = 10),
                left = text_grob("Figure arranged using ggpubr", color = "green", rot = 90),
                right = "I&#39;m done, thanks :-)!",
                fig.lab = "Figure 1", fig.lab.face = "bold"
                )

Align plot panels

A real use case is, for example, when plotting survival curves with the risk table placed under the main plot.

To illustrate this case, we’ll use the survminer package. First, install it using install.packages(“survminer”), then type this:

# Fit survival curves
library(survival)
fit <- survfit( Surv(time, status) ~ adhere, data = colon )

# Plot survival curves
library(survminer)
ggsurv <- ggsurvplot(fit, data = colon, 
                     palette = "jco",                              # jco palette
                     pval = TRUE, pval.coord = c(500, 0.4),        # Add p-value
                     risk.table = TRUE                            # Add risk table
                     )
names(ggsurv)

## [1] "plot"           "table"          "data.survplot"  "data.survtable"

ggsurv is a list including the following components:

• plot: survival curves
• table: the risk table plot

You can arrange the survival plot and the risk table as follow:

ggarrange(ggsurv$plot, ggsurv$table, heights = c(2, 0.7),
          ncol = 1, nrow = 2)

ggarrange(ggsurv$plot, ggsurv$table, heights = c(2, 0.7),
          ncol = 1, nrow = 2, align = "v")

Change column/row span of a plot

ggarrange(sp,                                                 # First row with scatter plot
          ggarrange(bxp, dp, ncol = 2, labels = c("B", "C")), # Second row with box and dot plots
          nrow = 2, 
          labels = "A"                                        # Labels of the scatter plot
          ) 

ggdraw()

draw_plot(plot, x = 0, y = 0, width = 1, height = 1)

draw_plot_label(label, x = 0, y = 1, size = 16, ...)

library("cowplot")
ggdraw() +
  draw_plot(bxp, x = 0, y = .5, width = .5, height = .5) +
  draw_plot(dp, x = .5, y = .5, width = .5, height = .5) +
  draw_plot(bp, x = 0, y = 0, width = 1, height = 0.5) +
  draw_plot_label(label = c("A", "B", "C"), size = 15,
                  x = c(0, 0.5, 0), y = c(1, 1, 0.5))

library("gridExtra")
grid.arrange(sp,                             # First row with one plot spaning over 2 columns
             arrangeGrob(bxp, dp, ncol = 2), # Second row with 2 plots in 2 different columns
             nrow = 2)                       # Number of rows

grid.arrange(bp,                                    # bar plot spaning two columns
             bxp, sp,                               # box plot and scatter plot
             ncol = 2, nrow = 2, 
             layout_matrix = rbind(c(1,1), c(2,3)))

library("gridExtra")
library("cowplot")

# Arrange plots using arrangeGrob
# returns a gtable (gt)
gt <- arrangeGrob(bp,                               # bar plot spaning two columns
             bxp, sp,                               # box plot and scatter plot
             ncol = 2, nrow = 2, 
             layout_matrix = rbind(c(1,1), c(2,3)))

# Add labels to the arranged plots
p <- as_ggplot(gt) +                                # transform to a ggplot
  draw_plot_label(label = c("A", "B", "C"), size = 15,
                  x = c(0, 0, 0.5), y = c(1, 0.5, 0.5)) # Add labels
p

library(grid)

# Move to a new page
grid.newpage()

# Create layout : nrow = 3, ncol = 2
pushViewport(viewport(layout = grid.layout(nrow = 3, ncol = 2)))

# A helper function to define a region on the layout
define_region <- function(row, col){
  viewport(layout.pos.row = row, layout.pos.col = col)
} 

# Arrange the plots
print(sp, vp = define_region(row = 1, col = 1:2))   # Span over two columns
print(bxp, vp = define_region(row = 2, col = 1))
print(dp, vp = define_region(row = 2, col = 2))
print(bp + rremove("x.text"), vp = define_region(row = 3, col = 1:2))

Use common legend for combined ggplots

To place a common unique legend in the margin of the arranged plots, the function ggarrange() [in ggpubr] can be used with the following arguments:

• common.legend = TRUE: place a common legend in a margin
• legend: specify the legend position. Allowed values include one of c(“top”, “bottom”, “left”, “right”)

ggarrange(bxp, dp, labels = c("A", "B"),
          common.legend = TRUE, legend = "bottom")

Scatter plot with marginal density plots

# Scatter plot colored by groups ("Species")
sp <- ggscatter(iris, x = "Sepal.Length", y = "Sepal.Width",
                color = "Species", palette = "jco",
                size = 3, alpha = 0.6)+
  border()                                         

# Marginal density plot of x (top panel) and y (right panel)
xplot <- ggdensity(iris, "Sepal.Length", fill = "Species",
                   palette = "jco")
yplot <- ggdensity(iris, "Sepal.Width", fill = "Species", 
                   palette = "jco")+
  rotate()

# Cleaning the plots
yplot <- yplot + clean_theme() 
xplot <- xplot + clean_theme()

# Arranging the plot
ggarrange(xplot, NULL, sp, yplot, 
          ncol = 2, nrow = 2,  align = "hv", 
          widths = c(2, 1), heights = c(1, 2),
          common.legend = TRUE)

Mix table, text and ggplot2 graphs

In this section, we’ll show how to plot a table and text alongside a chart. The iris data set will be used.

We start by creating the following plots:

1. a density plot of the variable “Sepal.Length”. R function: ggdensity() [in ggpubr]
2. a plot of the summary table containing the descriptive statistics (mean, sd, … ) of Sepal.Length.
   • R function for computing descriptive statistics: desc_statby() [in ggpubr].
   • R function to draw a textual table: ggtexttable() [in ggpubr].
3. a plot of a text paragraph. R function: ggparagraph() [in ggpubr].

We finish by arranging/combining the three plots using the function ggarrange() [in ggpubr]

# Density plot of "Sepal.Length"
#::::::::::::::::::::::::::::::::::::::
density.p <- ggdensity(iris, x = "Sepal.Length", 
                       fill = "Species", palette = "jco")

# Draw the summary table of Sepal.Length
#::::::::::::::::::::::::::::::::::::::
# Compute descriptive statistics by groups
stable <- desc_statby(iris, measure.var = "Sepal.Length",
                      grps = "Species")
stable <- stable[, c("Species", "length", "mean", "sd")]
# Summary table plot, medium orange theme
stable.p <- ggtexttable(stable, rows = NULL, 
                        theme = ttheme("mOrange"))

# Draw text
#::::::::::::::::::::::::::::::::::::::
text <- paste("iris data set gives the measurements in cm",
              "of the variables sepal length and width",
              "and petal length and width, respectively,",
              "for 50 flowers from each of 3 species of iris.",
             "The species are Iris setosa, versicolor, and virginica.", sep = " ")
text.p <- ggparagraph(text = text, face = "italic", size = 11, color = "black")

# Arrange the plots on the same page
ggarrange(density.p, stable.p, text.p, 
          ncol = 1, nrow = 3,
          heights = c(1, 0.5, 0.3))

Insert a graphical element inside a ggplot

The function annotation_custom() [in ggplot2] can be used for adding tables, plots or other grid-based elements within the plotting area of a ggplot. The simplified format is :

annotation_custom(grob, xmin, xmax, ymin, ymax)

density.p + annotation_custom(ggplotGrob(stable.p),
                              xmin = 5.5, ymin = 0.7,
                              xmax = 8)

# Scatter plot colored by groups ("Species")
#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
sp <- ggscatter(iris, x = "Sepal.Length", y = "Sepal.Width",
                color = "Species", palette = "jco",
                size = 3, alpha = 0.6)

# Create box plots of x/y variables
#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
# Box plot of the x variable
xbp <- ggboxplot(iris$Sepal.Length, width = 0.3, fill = "lightgray") +
  rotate() +
  theme_transparent()
# Box plot of the y variable
ybp <- ggboxplot(iris$Sepal.Width, width = 0.3, fill = "lightgray") +
  theme_transparent()

# Create the external graphical objects
# called a "grop" in Grid terminology
xbp_grob <- ggplotGrob(xbp)
ybp_grob <- ggplotGrob(ybp)


# Place box plots inside the scatter plot
#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
xmin <- min(iris$Sepal.Length); xmax <- max(iris$Sepal.Length)
ymin <- min(iris$Sepal.Width); ymax <- max(iris$Sepal.Width)
yoffset <- (1/15)*ymax; xoffset <- (1/15)*xmax

# Insert xbp_grob inside the scatter plot
sp + annotation_custom(grob = xbp_grob, xmin = xmin, xmax = xmax, 
                       ymin = ymin-yoffset, ymax = ymin+yoffset) +
  # Insert ybp_grob inside the scatter plot
  annotation_custom(grob = ybp_grob,
                       xmin = xmin-xoffset, xmax = xmin+xoffset, 
                       ymin = ymin, ymax = ymax)

# Import the image
img.file <- system.file(file.path("https://www.sthda.com/english/sthda-upload/images/ggpubr", "background-image.png"),
                        package = "ggpubr")
img <- png::readPNG(img.file)

library(ggplot2)
library(ggpubr)
ggplot(iris, aes(Species, Sepal.Length))+
  background_image(img)+
  geom_boxplot(aes(fill = Species), color = "white")+
  fill_palette("jco")

library(ggplot2)
library(ggpubr)
ggplot(iris, aes(Species, Sepal.Length))+
  background_image(img)+
  geom_boxplot(aes(fill = Species), color = "white", alpha = 0.5)+
  fill_palette("jco")

mypngfile <- download.file("https://upload.wikimedia.org/wikipedia/commons/thumb/e/e4/France_Flag_Map.svg/612px-France_Flag_Map.svg.png", 
                           destfile = "france.png", mode = &#39;wb&#39;) 

img <- png::readPNG(&#39;france.png&#39;) 

ggplot(iris, aes(x = Sepal.Length, y = Sepal.Width)) +
  background_image(img)+
  geom_point(aes(color = Species), alpha = 0.6, size = 5)+
  color_palette("jco")+
  theme(legend.position = "top")

Arrange over multiple pages

If you have a long list of ggplots, say n = 20 plots, you may want to arrange the plots and to place them on multiple pages. With 4 plots per page, you need 5 pages to hold the 20 plots.

The function ggarrange() [in ggpubr] provides a convenient solution to arrange multiple ggplots over multiple pages. After specifying the arguments nrow and ncol, the function ggarrange() computes automatically the number of pages required to hold the list of the plots. It returns a list of arranged ggplots.

For example the following R code,

multi.page <- ggarrange(bxp, dp, bp, sp,
                        nrow = 1, ncol = 2)

returns a list of two pages with two plots per page. You can visualize each page as follow:

multi.page[[1]] # Visualize page 1
multi.page[[2]] # Visualize page 2

You can also export the arranged plots to a pdf file using the function ggexport() [in ggpubr]:

ggexport(multi.page, filename = "multi.page.ggplot2.pdf")

PDF file: Multi.page.ggplot2

library(gridExtra)
res <- marrangeGrob(list(bxp, dp, bp, sp), nrow = 1, ncol = 2)

# Export to a pdf file
ggexport(res, filename = "multi.page.ggplot2.pdf")

# Visualize interactively
res

Nested layout with ggarrange()

We’ll arrange the plot created in section (@ref(mix-table-text-and-ggplot)) and (@ref(create-some-plots)).

p1 <- ggarrange(sp, bp + font("x.text", size = 9),
                ncol = 1, nrow = 2)

p2 <- ggarrange(density.p, stable.p, text.p, 
                ncol = 1, nrow = 3,
                heights = c(1, 0.5, 0.3))

ggarrange(p1, p2, ncol = 2, nrow = 1)

Export plots

R function: ggexport() [in ggpubr].

First, create a list of 4 ggplots corresponding to the variables Sepal.Length, Sepal.Width, Petal.Length and Petal.Width in the iris data set.

plots <- ggboxplot(iris, x = "Species",
                   y = c("Sepal.Length", "Sepal.Width", "Petal.Length", "Petal.Width"),
                   color = "Species", palette = "jco"
                   )
plots[[1]]  # Print the first plot
plots[[2]]  # Print the second plots and so on...

Next, you can export individual plots to a file (pdf, eps or png) (one plot per page). It’s also possible to arrange the plots (2 plot per page) when exporting them.

Export individual plots to a pdf file (one plot per page):

ggexport(plotlist = plots, filename = "test.pdf")

Arrange and export. Specify nrow and ncol to display multiple plots on the same page:

ggexport(plotlist = plots, filename = "test.pdf",
         nrow = 2, ncol = 1)

Acknoweledgment

We sincerely thank all developers for their efforts behind the packages that ggpubr depends on, namely:

• Baptiste Auguie (2016). gridExtra: Miscellaneous Functions for “Grid” Graphics. R package version 2.2.1. https://CRAN.R-project.org/package=gridExtra
• Claus O. Wilke (2016). cowplot: Streamlined Plot Theme and Plot Annotations for ‘ggplot2’. R package version 0.7.0. https://CRAN.R-project.org/package=cowplot
• H. Wickham. ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York, 2009.

Prerequisites

install.packages("ggpubr")

if(!require(devtools)) install.packages("devtools")
devtools::install_github("kassambara/ggpubr")

library(ggpubr)

Basic bar plots

Create a demo data set:

df <- data.frame(dose=c("D0.5", "D1", "D2"),
                 len=c(4.2, 10, 29.5))
print(df)

##   dose  len
## 1 D0.5  4.2
## 2   D1 10.0
## 3   D2 29.5

Basic bar plots:

# Basic bar plots with label
p <- ggbarplot(df, x = "dose", y = "len",
          color = "black", fill = "lightgray")
p

# Rotate to create horizontal bar plots
p + rotate()

Change fill and outline colors by groups:

ggbarplot(df, x = "dose", y = "len",
   fill = "dose", color = "dose", palette = "jco")

Multiple grouping variables

Create a demo data set:

df2 <- data.frame(supp=rep(c("VC", "OJ"), each=3),
                  dose=rep(c("D0.5", "D1", "D2"),2),
                  len=c(6.8, 15, 33, 4.2, 10, 29.5))
print(df2)

##   supp dose  len
## 1   VC D0.5  6.8
## 2   VC   D1 15.0
## 3   VC   D2 33.0
## 4   OJ D0.5  4.2
## 5   OJ   D1 10.0
## 6   OJ   D2 29.5

Plot y = “len” by x = “dose” and change color by a second group: “supp”

# Stacked bar plots, add labels inside bars
ggbarplot(df2, x = "dose", y = "len",
  fill = "supp", color = "supp", 
  palette = c("gray", "black"),
  label = TRUE, lab.col = "white", lab.pos = "in")

# Change position: Interleaved (dodged) bar plot
ggbarplot(df2, x = "dose", y = "len",
          fill = "supp", color = "supp", 
          palette = c("gray", "black"),
          position = position_dodge(0.9))

Ordered bar plots

Load and prepare data:

# Load data
data("mtcars")
dfm <- mtcars
# Convert the cyl variable to a factor
dfm$cyl <- as.factor(dfm$cyl)
# Add the name colums
dfm$name <- rownames(dfm)
# Inspect the data
head(dfm[, c("name", "wt", "mpg", "cyl")])

##                                name   wt  mpg cyl
## Mazda RX4                 Mazda RX4 2.62 21.0   6
## Mazda RX4 Wag         Mazda RX4 Wag 2.88 21.0   6
## Datsun 710               Datsun 710 2.32 22.8   4
## Hornet 4 Drive       Hornet 4 Drive 3.21 21.4   6
## Hornet Sportabout Hornet Sportabout 3.44 18.7   8
## Valiant                     Valiant 3.46 18.1   6

Create ordered bar plots. Change the fill color by the grouping variable “cyl”. Sorting will be done globally, but not by groups.

ggbarplot(dfm, x = "name", y = "mpg",
          fill = "cyl",               # change fill color by cyl
          color = "white",            # Set bar border colors to white
          palette = "jco",            # jco journal color palett. see ?ggpar
          sort.val = "desc",          # Sort the value in dscending order
          sort.by.groups = FALSE,     # Don't sort inside each group
          x.text.angle = 90           # Rotate vertically x axis texts
          )

Sort bars inside each group. Use the argument sort.by.groups = TRUE.

ggbarplot(dfm, x = "name", y = "mpg",
          fill = "cyl",               # change fill color by cyl
          color = "white",            # Set bar border colors to white
          palette = "jco",            # jco journal color palett. see ?ggpar
          sort.val = "asc",           # Sort the value in dscending order
          sort.by.groups = TRUE,      # Sort inside each group
          x.text.angle = 90           # Rotate vertically x axis texts
          )

Deviation graphs

The deviation graph shows the deviation of quantitative values to a reference value. In the R code below, we’ll plot the mpg z-score from the mtcars data set.

Calculate the z-score of the mpg data:

# Calculate the z-score of the mpg data
dfm$mpg_z <- (dfm$mpg -mean(dfm$mpg))/sd(dfm$mpg)
dfm$mpg_grp <- factor(ifelse(dfm$mpg_z < 0, "low", "high"), 
                     levels = c("low", "high"))
# Inspect the data
head(dfm[, c("name", "wt", "mpg", "mpg_z", "mpg_grp", "cyl")])

##                                name   wt  mpg  mpg_z mpg_grp cyl
## Mazda RX4                 Mazda RX4 2.62 21.0  0.151    high   6
## Mazda RX4 Wag         Mazda RX4 Wag 2.88 21.0  0.151    high   6
## Datsun 710               Datsun 710 2.32 22.8  0.450    high   4
## Hornet 4 Drive       Hornet 4 Drive 3.21 21.4  0.217    high   6
## Hornet Sportabout Hornet Sportabout 3.44 18.7 -0.231     low   8
## Valiant                     Valiant 3.46 18.1 -0.330     low   6

Create an ordered bar plot, colored according to the level of mpg:

ggbarplot(dfm, x = "name", y = "mpg_z",
          fill = "mpg_grp",           # change fill color by mpg_level
          color = "white",            # Set bar border colors to white
          palette = "jco",            # jco journal color palett. see ?ggpar
          sort.val = "asc",           # Sort the value in ascending order
          sort.by.groups = FALSE,     # Don't sort inside each group
          x.text.angle = 90,          # Rotate vertically x axis texts
          ylab = "MPG z-score",
          xlab = FALSE,
          legend.title = "MPG Group"
          )

Rotate the plot: use rotate = TRUE and sort.val = “desc”

ggbarplot(dfm, x = "name", y = "mpg_z",
          fill = "mpg_grp",           # change fill color by mpg_level
          color = "white",            # Set bar border colors to white
          palette = "jco",            # jco journal color palett. see ?ggpar
          sort.val = "desc",          # Sort the value in descending order
          sort.by.groups = FALSE,     # Don't sort inside each group
          x.text.angle = 90,          # Rotate vertically x axis texts
          ylab = "MPG z-score",
          legend.title = "MPG Group",
          rotate = TRUE,
          ggtheme = theme_minimal()
          )

Alternatives to bar plots

ggdotchart(dfm, x = "name", y = "mpg",
           color = "cyl",                                # Color by groups
           palette = c("#00AFBB", "#E7B800", "#FC4E07"), # Custom color palette
           sorting = "ascending",                        # Sort value in descending order
           add = "segments",                             # Add segments from y = 0 to dots
           ggtheme = theme_pubr()                        # ggplot2 theme
           )

ggdotchart(dfm, x = "name", y = "mpg",
           color = "cyl",                                # Color by groups
           palette = c("#00AFBB", "#E7B800", "#FC4E07"), # Custom color palette
           sorting = "descending",                       # Sort value in descending order
           add = "segments",                             # Add segments from y = 0 to dots
           rotate = TRUE,                                # Rotate vertically
           group = "cyl",                                # Order by groups
           dot.size = 6,                                 # Large dot size
           label = round(dfm$mpg),                        # Add mpg values as dot labels
           font.label = list(color = "white", size = 9, 
                             vjust = 0.5),               # Adjust label parameters
           ggtheme = theme_pubr()                        # ggplot2 theme
           )

ggdotchart(dfm, x = "name", y = "mpg_z",
           color = "cyl",                                # Color by groups
           palette = c("#00AFBB", "#E7B800", "#FC4E07"), # Custom color palette
           sorting = "descending",                       # Sort value in descending order
           add = "segments",                             # Add segments from y = 0 to dots
           add.params = list(color = "lightgray", size = 2), # Change segment color and size
           group = "cyl",                                # Order by groups
           dot.size = 6,                                 # Large dot size
           label = round(dfm$mpg_z,1),                        # Add mpg values as dot labels
           font.label = list(color = "white", size = 9, 
                             vjust = 0.5),               # Adjust label parameters
           ggtheme = theme_pubr()                        # ggplot2 theme
           )+
  geom_hline(yintercept = 0, linetype = 2, color = "lightgray")

ggdotchart(dfm, x = "name", y = "mpg",
           color = "cyl",                                # Color by groups
           palette = c("#00AFBB", "#E7B800", "#FC4E07"), # Custom color palette
           sorting = "descending",                       # Sort value in descending order
           rotate = TRUE,                                # Rotate vertically
           dot.size = 2,                                 # Large dot size
           y.text.col = TRUE,                            # Color y text by groups
           ggtheme = theme_pubr()                        # ggplot2 theme
           )+
  theme_cleveland()                                      # Add dashed grids

Prerequisites

if(!require(devtools)) install.packages("devtools")
devtools::install_github("kassambara/ggpubr")

install.packages("ggpubr")

library(ggpubr)

# ToothGrowth
data("ToothGrowth")
head(ToothGrowth)

##    len supp dose
## 1  4.2   VC  0.5
## 2 11.5   VC  0.5
## 3  7.3   VC  0.5
## 4  5.8   VC  0.5
## 5  6.4   VC  0.5
## 6 10.0   VC  0.5

# mtcars 
data("mtcars")
head(mtcars[, c("wt", "mpg", "cyl")])

##                     wt  mpg cyl
## Mazda RX4         2.62 21.0   6
## Mazda RX4 Wag     2.88 21.0   6
## Datsun 710        2.32 22.8   4
## Hornet 4 Drive    3.21 21.4   6
## Hornet Sportabout 3.44 18.7   8
## Valiant           3.46 18.1   6

Error plots

R function: ggerrorplot() [in ggpubr].

Simplified format:

ggerrorplot(data, x, y, desc_stat = "mean_se")

For example, the following R code uses the ToothGrowth data set and plots y = “len” by x = “dose”.

# Mean +/- standard deviation
ggerrorplot(ToothGrowth, x = "dose", y = "len", 
            desc_stat = "mean_sd")
# Change error plot type and add mean points
ggerrorplot(ToothGrowth, x = "dose", y = "len", 
            desc_stat = "mean_sd",
            error.plot = "errorbar",            # Change error plot type
            add = "mean"                        # Add mean points
            )

It’s also possible to add jitter points (representing individual points), dot plots and violin plots:

# Add jittered points
ggerrorplot(ToothGrowth, x = "dose", y = "len", 
            desc_stat = "mean_sd", color = "black",
            add = "jitter", add.params = list(color = "darkgray")
            )
# Add dot plots
ggerrorplot(ToothGrowth, x = "dose", y = "len", 
            desc_stat = "mean_sd", color = "black",
            add = "dotplot", add.params = list(color = "darkgray")
            )
# Add violin plots
ggerrorplot(ToothGrowth, x = "dose", y = "len", 
            desc_stat = "mean_sd", color = "black",
            add = "violin", add.params = list(color = "darkgray")
            )

To add p-values comparing means, use this:

# Specify the comparisons you want
my_comparisons <- list( c("0.5", "1"), c("1", "2"), c("0.5", "2") )
ggerrorplot(ToothGrowth, x = "dose", y = "len",
            desc_stat = "mean_sd", color = "black",
            add = "violin", add.params = list(color = "darkgray"))+ 
  stat_compare_means(comparisons = my_comparisons)+ # Add pairwise comparisons p-value
  stat_compare_means(label.y = 50)                  # Add global p-value

Read more at : Add P-values and Significance Levels to ggplots.

Color by a grouping variable:

# Color by "dose" (same variable used on x-axis)
ggerrorplot(ToothGrowth, x = "dose", y = "len", 
            desc_stat = "mean_sd", 
            color = "dose", palette = "jco")
# Color by another grouping variable "supp"
ggerrorplot(ToothGrowth, x = "dose", y = "len", 
            desc_stat = "mean_sd", 
            color = "supp", palette = "jco",
            position = position_dodge(0.3)     # Adjust the space between bars
            )

Line plots

You can create a line plot of mean +/- error using the function ggline()[in ggpubr]. The format is as follow:

# Basic line plots of means +/- se with jittered points
ggline(ToothGrowth, x = "dose", y = "len", 
       add = c("mean_se", "jitter"))

Color by groups:

ggline(ToothGrowth, x = "dose", y = "len", 
       add = c("mean_se", "jitter"),
       color = "supp", palette = "jco")

Bar plots

R function ggbarplot()[in ggpubr]. The format is as follow:

# Basic bar plots of means +/- se with jittered points
ggbarplot(ToothGrowth, x = "dose", y = "len", 
       add = c("mean_se", "jitter"))

Color by groups:

ggbarplot(ToothGrowth, x = "dose", y = "len", 
          add = c("mean_se", "jitter"),
          color = "supp", palette = "jco",
          position = position_dodge(0.8))

Add labels

In this section we’ll plot group means with individual information.

Data: mtcars.

# Prepare the data set
# Use row names as individual names
df <- as.data.frame(mtcars[, c("am", "hp")])
df$name <- rownames(df)
head(df)

##                   am  hp              name
## Mazda RX4          1 110         Mazda RX4
## Mazda RX4 Wag      1 110     Mazda RX4 Wag
## Datsun 710         1  93        Datsun 710
## Hornet 4 Drive     0 110    Hornet 4 Drive
## Hornet Sportabout  0 175 Hornet Sportabout
## Valiant            0 105           Valiant

Create a bar plot with individual labels.

set.seed(123)
# Bar plot of mean +/- se, add individual points
ggbarplot(df, x = "am", y = "hp",
          add = c("mean_se", "point"),
          color = "am", fill = "am", alpha = 0.5,
          palette = "jco")+
   ggrepel::geom_text_repel(aes(label = name))

Application to gene expression data

In our previous article - Facilitating Exploratory Data Visualization: Application to TCGA Genomic Data - we described how to visualize gene expression data using box plots, violin plots, dot plots and stripcharts. We also demonstrated how to combine the plot of multiples variables (genes) in the same plot.

Here we provide some R code to visualize the mean expression profile of one or multiple genes. We’ll use the gene expression data set described in our previous tutorial: Facilitating Exploratory Data Visualization: Application to TCGA Genomic Data.

expr <- read.delim("https://raw.githubusercontent.com/kassambara/data/master/expr_tcga.txt",
                   stringsAsFactors = FALSE)

The data set contains the mRNA expression for five genes of interest - GATA3, PTEN, XBP1, ESR1 and MUC1 - from 3 different data sets:

• Breast invasive carcinoma (BRCA),
• Ovarian serous cystadenocarcinoma (OV) and
• Lung squamous cell carcinoma (LUSC)

The R code below displays the mean expression of three genes - “GATA3”, “PTEN” and “XBP1”.

ggline(expr, x = "dataset",
      y = c("GATA3", "PTEN", "XBP1"),
      combine = TRUE,
      ylab = "Expression", 
      add = "mean_sd")

You can also add other geometries on the mean plot such as jitter points, dotplot or violin. To add a violin plot, type this:

ggline(expr, x = "dataset",
      y = c("GATA3", "PTEN", "XBP1"),
      combine = TRUE,
      ylab = "Expression", 
      color = "gray",                                     # Line color
      add = c("mean_sd", "violin"),                     
      add.params = list(color = "dataset"),
      palette = "jco"
      )

To add jitter points, we’ll use a small subset of data for readability:

# Subset 50 random rows
set.seed(123)
random_rows <- sample(1:nrow(expr), 50)
expr2 <- expr[random_rows, ]
# Visualize
ggline(expr2, x = "dataset",
      y = c("GATA3", "PTEN", "XBP1"),
      combine = TRUE,
      ylab = "Expression", 
      color = "gray",                           
      add = c("mean_sd", "jitter"),                     
      add.params = list(color = "dataset", size = 0.5),
      palette = "jco"
      )

As previously shown, you can merge the three plots as follow:

# Merge the three plot
ggline(expr2, x = "dataset",
      y = c("GATA3", "PTEN", "XBP1"),
      merge = TRUE,
      ylab = "Expression", 
      add = "mean_sd",
      palette = "jco")
# Add  jitter
ggline(expr2, x = "dataset",
      y = c("GATA3", "PTEN", "XBP1"),
      merge = TRUE,
      ylab = "Expression", 
      add = c("mean_sd", "jitter"),              # Add jitter points
      palette = "jco")

Show line labels:

ggline(expr2, x = "dataset",
      y = c("GATA3", "PTEN", "XBP1"),
      merge = TRUE,
      ylab = "Expression",                   
      add = "mean",   
      show.line.label = TRUE,
      repel = TRUE,
      legend = "none",
      palette = rep("black", 3)   # Black color for each line
      )

Let’s plot a complex plot with point labels:

# Add  jitter
ggline(expr2, x = "dataset",
      y = c("GATA3", "PTEN", "XBP1"),
      merge = TRUE,
      ylab = "Expression", 
      add = c("mean_se", "jitter"),              # Add mean_se and jitter points
      add.params = list(size = 0.7),             # Add point size
      label = "bcr_patient_barcode",             # Add point labels
      label.select = list(top.up = 2),           # show only labels for the top 2 points
      font.label = list(color = ".y."),          # Color labels by .y., here gene names
      repel = TRUE,                              # Use repel to avoid labels overplotting
      palette = "jco")

Plot a bar plot of means:

# Create bar plots
ggbarplot(expr2, x = "dataset",
      y = c("GATA3", "PTEN", "XBP1"),
      combine = TRUE,
      ylab = "Expression", 
      add = "mean_se")

# Merge bar plots
ggbarplot(expr2, x = "dataset",
      y = c("GATA3", "PTEN", "XBP1"),
      merge = TRUE,
      ylab = "Expression", 
      add = "mean_se", palette = "jco")

Error plots:

# Create error plots
ggerrorplot(expr2, x = "dataset",
      y = c("GATA3", "PTEN", "XBP1"),
      combine = TRUE,
      ylab = "Expression", 
      add = "mean_sd")

# Merge error plots
ggerrorplot(expr2, x = "dataset",
      y = c("GATA3", "PTEN", "XBP1"),
      merge = TRUE,
      ylab = "Expression", 
      add = "mean_sd", palette = "jco",
      position = position_dodge(0.3)
      )

Prerequisites

if(!require(devtools)) install.packages("devtools")
devtools::install_github("kassambara/ggpubr")

install.packages("ggpubr")

library(ggpubr)

# Load data
data("mtcars")
df <- mtcars
# Convert cyl as a grouping variable
df$cyl <- as.factor(df$cyl)

# Inspect the data
head(df[, c("wt", "mpg", "cyl", "qsec")])

##                     wt  mpg cyl qsec
## Mazda RX4         2.62 21.0   6 16.5
## Mazda RX4 Wag     2.88 21.0   6 17.0
## Datsun 710        2.32 22.8   4 18.6
## Hornet 4 Drive    3.21 21.4   6 19.4
## Hornet Sportabout 3.44 18.7   8 17.0
## Valiant           3.46 18.1   6 20.2

Basic plots

ggscatter(df, x = "wt", y = "mpg",
          add = "reg.line",                                 # Add regression line
          conf.int = TRUE,                                  # Add confidence interval
          add.params = list(color = "blue",
                            fill = "lightgray")
          )+
  stat_cor(method = "pearson", label.x = 3, label.y = 30)  # Add correlation coefficient

You can change the point shape, by specifying the argument shape, for example:

ggscatter(df, x = "wt", y = "mpg",
          shape = 18)

To see the different point shapes commonly used in R, type this:

show_point_shapes()

Color by groups

Grouping variable: cyl. To add a correlation coefficient per group, specify the grouping variable using the mapping function aes() as follow.

ggscatter(df, x = "wt", y = "mpg",
          add = "reg.line",                         # Add regression line
          conf.int = TRUE,                          # Add confidence interval
          color = "cyl", palette = "jco",           # Color by groups "cyl"
          shape = "cyl"                             # Change point shape by groups "cyl"
          )+
  stat_cor(aes(color = cyl), label.x = 3)           # Add correlation coefficient

# Extending the regression line --> fullrange = TRUE
# Add marginal rug (marginal density) ---> rug = TRUE
ggscatter(df, x = "wt", y = "mpg",
          add = "reg.line",                         # Add regression line
          color = "cyl", palette = "jco",           # Color by groups "cyl"
          shape = "cyl",                            # Change point shape by groups "cyl"
          fullrange = TRUE,                         # Extending the regression line
          rug = TRUE                                # Add marginal rug
          )+
  stat_cor(aes(color = cyl), label.x = 3)           # Add correlation coefficient

Add concentration ellipses

Main arguments:

ggscatter(df, x = "wt", y = "mpg",
          color = "cyl", palette = "jco",
          shape = "cyl",
          ellipse = TRUE)

# Change the ellipse type to 'convex'
ggscatter(df, x = "wt", y = "mpg",
          color = "cyl", palette = "jco",
          shape = "cyl",
          ellipse = TRUE, ellipse.type = "convex")

# Add group mean points and stars
ggscatter(df, x = "wt", y = "mpg",
          color = "cyl", palette = "jco",
          shape = "cyl",
          ellipse = TRUE, 
          mean.point = TRUE,
          star.plot = TRUE)

Add point labels

Main arguments:

# Use row names as point labels
df$name <- rownames(df)
ggscatter(df, x = "wt", y = "mpg",
   color = "cyl", palette = "jco",
   label = "name", repel = TRUE)

# Select some labels to show
ggscatter(df, x = "wt", y = "mpg",
   color = "cyl", palette = "jco",
   label = "name", repel = TRUE,
   label.select = c("Toyota Corolla", "Merc 280", "Duster 360"))

# Show labels according to some criteria: x and y values
ggscatter(df, x = "wt", y = "mpg",
   color = "cyl", palette = "jco",
   label = "name", repel = TRUE,
   label.select = list(criteria = "`x` > 4 & `y` < 15"))

Bubble chart

In a bubble chart, points size is controlled by a continuous variable, here “qsec”. In the R code below, the argument alpha is used to control color transparency. alpha should be between 0 and 1.

ggscatter(df, x = "wt", y = "mpg",
   color = "cyl", palette = "jco",
   size = "qsec", alpha = 0.5)+
  scale_size(range = c(0.5, 15))    # Adjust the range of points size

Color by a continuous variable

The R code below, will color points according to the values of a continuous variable, here “mpg”. By default, a blue gradient color is created. This can be changed using the helper function gradient_color() [in ggpubr].

# Color by continuous variable
p <- ggscatter(df, x = "wt", y = "mpg",
               color = "mpg")
p

# Change gradient color
p + gradient_color(c("blue", "white", "red"))

Add marginal plots

The function ggMarginal() [in ggExtra package], can be used to easily add a marginal histogram, density or boxplot to a scatter plot.

First, install the ggExtra package as follow: install.packages(“ggExtra”); then type the following R code:

# Add density distribution as marginal plot
library("ggExtra")
p <- ggscatter(iris, x = "Sepal.Length", y = "Sepal.Width",
               color = "Species", palette = "jco",
               size = 3, alpha = 0.6)
ggMarginal(p, type = "density")

# Change marginal plot type
ggMarginal(p, type = "boxplot")

One limitation of ggExtra is that it can’t cope with multiple groups in the scatter plot and the marginal plots. In the R code below, we provide a solution using the cowplot package.

# Scatter plot colored by groups ("Species")
sp <- ggscatter(iris, x = "Sepal.Length", y = "Sepal.Width",
                color = "Species", palette = "jco",
                size = 3, alpha = 0.6)+
  border()                                         

# Marginal density plot of x (top panel) and y (right panel)
xplot <- ggdensity(iris, "Sepal.Length", fill = "Species",
                   palette = "jco")
yplot <- ggdensity(iris, "Sepal.Width", fill = "Species", 
                   palette = "jco")+
  rotate()

# Cleaning the plots
sp <- sp + rremove("legend")
yplot <- yplot + clean_theme() + rremove("legend")
xplot <- xplot + clean_theme() + rremove("legend")

# Arranging the plot using cowplot
library(cowplot)
plot_grid(xplot, NULL, sp, yplot, ncol = 2, align = "hv", 
          rel_widths = c(2, 1), rel_heights = c(1, 2))

Add marginal boxplot:

# Scatter plot colored by groups ("Species")
sp <- ggscatter(iris, x = "Sepal.Length", y = "Sepal.Width",
                color = "Species", palette = "jco",
                size = 3, alpha = 0.6, ggtheme = theme_bw())             

# Marginal boxplot of x (top panel) and y (right panel)
xplot <- ggboxplot(iris, x = "Species", y = "Sepal.Length", 
                   color = "Species", fill = "Species", palette = "jco",
                   alpha = 0.5, ggtheme = theme_bw())+
  rotate()
yplot <- ggboxplot(iris, x = "Species", y = "Sepal.Width",
                   color = "Species", fill = "Species", palette = "jco",
                   alpha = 0.5, ggtheme = theme_bw())

# Cleaning the plots
sp <- sp + rremove("legend")
yplot <- yplot + clean_theme() + rremove("legend")
xplot <- xplot + clean_theme() + rremove("legend")

# Arranging the plot using cowplot
library(cowplot)
plot_grid(xplot, NULL, sp, yplot, ncol = 2, align = "hv", 
          rel_widths = c(2, 1), rel_heights = c(1, 2))

The problem with the above plots, is the presence of extra spaces between the main plot and the marginal density plots. Recently, in a tweet post, Claus Wilke provides the following solution for creating a perfect scatter plot with marginal density plots or histogram plots:

library(cowplot) 

# Main plot
pmain <- ggplot(iris, aes(x = Sepal.Length, y = Sepal.Width, color = Species))+
  geom_point()+
  ggpubr::color_palette("jco")

# Marginal densities along x axis
xdens <- axis_canvas(pmain, axis = "x")+
  geom_density(data = iris, aes(x = Sepal.Length, fill = Species),
              alpha = 0.7, size = 0.2)+
  ggpubr::fill_palette("jco")

# Marginal densities along y axis
# Need to set coord_flip = TRUE, if you plan to use coord_flip()
ydens <- axis_canvas(pmain, axis = "y", coord_flip = TRUE)+
  geom_density(data = iris, aes(x = Sepal.Width, fill = Species),
                alpha = 0.7, size = 0.2)+
  coord_flip()+
  ggpubr::fill_palette("jco")


p1 <- insert_xaxis_grob(pmain, xdens, grid::unit(.2, "null"), position = "top")
p2<- insert_yaxis_grob(p1, ydens, grid::unit(.2, "null"), position = "right")
ggdraw(p2)

Add 2d density estimation

# Add 2d density estimation
sp <- ggscatter(iris, x = "Sepal.Length", y = "Sepal.Width",
                color = "lightgray")
sp + geom_density_2d()

# Gradient color
sp + stat_density_2d(aes(fill = ..level..), geom = "polygon")

# Change gradient color: custom
sp + stat_density_2d(aes(fill = ..level..), geom = "polygon")+
  gradient_fill(c("white", "steelblue"))

# Change the gradient color: RColorBrewer palette
sp + stat_density_2d(aes(fill = ..level..), geom = "polygon") +
  gradient_fill("YlOrRd")

Application to gene expression data

We’ll use the gene expression data set described in our previous tutorial: Facilitating Exploratory Data Visualization: Application to TCGA Genomic Data.

expr <- read.delim("https://raw.githubusercontent.com/kassambara/data/master/expr_tcga.txt",
                   stringsAsFactors = FALSE)

The data set contains the mRNA expression for five genes of interest - GATA3, PTEN, XBP1, ESR1 and MUC1 - from 3 different data sets:

• Breast invasive carcinoma (BRCA),
• Ovarian serous cystadenocarcinoma (OV) and
• Lung squamous cell carcinoma (LUSC)

The following plots show the association between GATA3 and ESR1 genes expression.

# Association between GATA3 and ESR1
# Color points by dataset
# Add correlation coefficient by dataset
ggscatter(expr, x = "GATA3", y = "ESR1", size = 0.3, 
          rug = TRUE,                                # Add marginal rug
          color = "dataset", palette = "jco") +
  stat_cor(aes(color = dataset), method = "spearman")

Facet/split by data set, add regression line and confidence interval:

ggscatter(expr, x = "GATA3", y = "ESR1", size = 0.3,
          color = "dataset", palette = "jco",
          facet.by = "dataset", #scales = "free_x",
          add = "reg.line", conf.int = TRUE) +
  stat_cor(aes(color = dataset), method = "spearman", label.y = 6)

Combining multiple plots. Visualize the correlation of GATA3 with two other genes (ESR1 and MUC1)

ggscatter(expr, x = "GATA3", y = c("ESR1", "MUC1"), size = 0.3,
          combine = TRUE, ylab = "Expression",
          color = "dataset", palette = "jco",
          add = "reg.line", conf.int = TRUE) +
  stat_cor(aes(color = dataset), method = "spearman")

Further readings

See also the ggpmisc R package to add linear model equation to a scatter plot.

Prerequisites

install.packages("ggpubr")

if(!require(devtools)) install.packages("devtools")
devtools::install_github("kassambara/ggpubr")

library(ggpubr)

# Load the bioconductor installer. 
source("https://bioconductor.org/biocLite.R")
# Install the main RTCGA package
biocLite("RTCGA")
# Install the clinical and mRNA gene expression data packages
biocLite("RTCGA.clinical")
biocLite("RTCGA.mRNA")

library(RTCGA)
infoTCGA()

## # A tibble: 38 x 13
##   Cohort    BCR Clinical     CN   LowP Methylation   mRNA mRNASeq    miR
## *                 
## 1    ACC     92       92     90      0          80      0      79      0
## 2   BLCA    412      412    410    112         412      0     408      0
## 3   BRCA   1098     1097   1089     19        1097    526    1093      0
## 4   CESC    307      307    295     50         307      0     304      0
## 5   CHOL     51       45     36      0          36      0      36      0
## 6   COAD    460      458    451     69         457    153     457      0
## # ... with 32 more rows, and 4 more variables: miRSeq , RPPA ,
## #   MAF , rawMAF 

Gene expression data

The R function expressionsTCGA() [in RTCGA package] can be used to easily extract the expression values of genes of interest in one or multiple cancer types.

In the following R code, we start by extracting the mRNA expression for five genes of interest - GATA3, PTEN, XBP1, ESR1 and MUC1 - from 3 different data sets:

• Breast invasive carcinoma (BRCA),
• Ovarian serous cystadenocarcinoma (OV) and
• Lung squamous cell carcinoma (LUSC)

library(RTCGA)
library(RTCGA.mRNA)
expr <- expressionsTCGA(BRCA.mRNA, OV.mRNA, LUSC.mRNA,
                        extract.cols = c("GATA3", "PTEN", "XBP1","ESR1", "MUC1"))
expr

## # A tibble: 1,305 x 7
##            bcr_patient_barcode   dataset GATA3  PTEN  XBP1   ESR1  MUC1
##                                     
## 1 TCGA-A1-A0SD-01A-11R-A115-07 BRCA.mRNA  2.87 1.361  2.98  3.084  1.65
## 2 TCGA-A1-A0SE-01A-11R-A084-07 BRCA.mRNA  2.17 0.428  2.55  2.386  3.08
## 3 TCGA-A1-A0SH-01A-11R-A084-07 BRCA.mRNA  1.32 1.306  3.02  0.791  2.99
## 4 TCGA-A1-A0SJ-01A-11R-A084-07 BRCA.mRNA  1.84 0.810  3.13  2.495 -1.92
## 5 TCGA-A1-A0SK-01A-12R-A084-07 BRCA.mRNA -6.03 0.251 -1.45 -4.861 -1.17
## 6 TCGA-A1-A0SM-01A-11R-A084-07 BRCA.mRNA  1.80 1.311  4.04  2.797  3.53
## # ... with 1,299 more rows

To display the number of sample in each data set, type this:

nb_samples <- table(expr$dataset)
nb_samples

## 
## BRCA.mRNA LUSC.mRNA   OV.mRNA 
##       590       154       561

We can simplify data set names by removing the “mRNA” tag. This can be done using the R base function gsub().

expr$dataset <- gsub(pattern = ".mRNA", replacement = "",  expr$dataset)

Let’s simplify also the patients’ barcode column. The following R code will change the barcodes into BRCA1, BRCA2, …, OV1, OV2, …., etc

expr$bcr_patient_barcode <- paste0(expr$dataset, c(1:590, 1:561, 1:154))
expr

## # A tibble: 1,305 x 7
##   bcr_patient_barcode dataset GATA3  PTEN  XBP1   ESR1  MUC1
##                          
## 1               BRCA1    BRCA  2.87 1.361  2.98  3.084  1.65
## 2               BRCA2    BRCA  2.17 0.428  2.55  2.386  3.08
## 3               BRCA3    BRCA  1.32 1.306  3.02  0.791  2.99
## 4               BRCA4    BRCA  1.84 0.810  3.13  2.495 -1.92
## 5               BRCA5    BRCA -6.03 0.251 -1.45 -4.861 -1.17
## 6               BRCA6    BRCA  1.80 1.311  4.04  2.797  3.53
## # ... with 1,299 more rows

The above (expr) dataset has been saved at https://raw.githubusercontent.com/kassambara/data/master/expr_tcga.txt. This data is required to practice the R code provided in this tutotial.

If you experience some issues in installing the RTCGA packages, You can simply load the data as follow:

expr <- read.delim("https://raw.githubusercontent.com/kassambara/data/master/expr_tcga.txt",
                   stringsAsFactors = FALSE)

Box plots

Create a box plot of a gene expression profile, colored by groups (here data set/cancer type):

library(ggpubr)
# GATA3
ggboxplot(expr, x = "dataset", y = "GATA3",
          title = "GATA3", ylab = "Expression",
          color = "dataset", palette = "jco")
# PTEN
ggboxplot(expr, x = "dataset", y = "PTEN",
          title = "PTEN", ylab = "Expression",
          color = "dataset", palette = "jco")

Instead of repeating the same R code for each gene, you can create a list of plots at once, as follow:

# Create a  list of plots
p <- ggboxplot(expr, x = "dataset", 
               y = c("GATA3", "PTEN", "XBP1"),
               title = c("GATA3", "PTEN", "XBP1"),
               ylab = "Expression", 
               color = "dataset", palette = "jco")
# View GATA3
p$GATA3
# View PTEN
p$PTEN
# View XBP1
p$XBP1

To add p-values and significance levels to the boxplots, read our previous article: Add P-values and Significance Levels to ggplots. Briefly, you can do this:

my_comparisons <- list(c("BRCA", "OV"), c("OV", "LUSC"))
ggboxplot(expr, x = "dataset", y = "GATA3",
          title = "GATA3", ylab = "Expression",
          color = "dataset", palette = "jco")+
  stat_compare_means(comparisons = my_comparisons)

For each of the genes, you can compare the different groups as follow:

compare_means(c(GATA3, PTEN, XBP1) ~ dataset, data = expr)

## # A tibble: 9 x 8
##      .y. group1 group2         p     p.adj p.format p.signif   method
##                             
## 1  GATA3   BRCA     OV 1.11e-177 3.34e-177  < 2e-16     **** Wilcoxon
## 2  GATA3   BRCA   LUSC  6.68e-73  1.34e-72  < 2e-16     **** Wilcoxon
## 3  GATA3     OV   LUSC  2.97e-08  2.97e-08  3.0e-08     **** Wilcoxon
## 4   PTEN   BRCA     OV  6.79e-05  6.79e-05  6.8e-05     **** Wilcoxon
## 5   PTEN   BRCA   LUSC  1.04e-16  3.13e-16  < 2e-16     **** Wilcoxon
## 6   PTEN     OV   LUSC  1.28e-07  2.56e-07  1.3e-07     **** Wilcoxon
## # ... with 3 more rows

If you want to select items (here cancer types) to display or to remove a particular item from the plot, use the argument select or remove, as follow:

# Select BRCA and OV cancer types
ggboxplot(expr, x = "dataset", y = "GATA3",
          title = "GATA3", ylab = "Expression",
          color = "dataset", palette = "jco",
          select = c("BRCA", "OV"))
# or remove BRCA
ggboxplot(expr, x = "dataset", y = "GATA3",
          title = "GATA3", ylab = "Expression",
          color = "dataset", palette = "jco",
          remove = "BRCA")

To change the order of the data sets on x axis, use the argument order. For example order = c(“LUSC”, “OV”, “BRCA”):

# Order data sets
ggboxplot(expr, x = "dataset", y = "GATA3",
          title = "GATA3", ylab = "Expression",
          color = "dataset", palette = "jco",
          order = c("LUSC", "OV", "BRCA"))

To create horizontal plots, use the argument rotate = TRUE:

ggboxplot(expr, x = "dataset", y = "GATA3",
          title = "GATA3", ylab = "Expression",
          color = "dataset", palette = "jco",
          rotate = TRUE)

To combine the three gene expression plots into a multi-panel plot, use the argument combine = TRUE:

ggboxplot(expr, x = "dataset",
          y = c("GATA3", "PTEN", "XBP1"),
          combine = TRUE,
          ylab = "Expression",
          color = "dataset", palette = "jco")

You can also merge the 3 plots using the argument merge = TRUE or merge = “asis”:

ggboxplot(expr, x = "dataset",
          y = c("GATA3", "PTEN", "XBP1"),
          merge = TRUE,
          ylab = "Expression", 
          palette = "jco")

In the plot above, It’s easy to visually compare the expression level of the different genes in each cancer type.

But you might want to put genes (y variables) on x axis, in order to compare the expression level in the different cell subpopulations.

In this situation, the y variables (i.e.: genes) become x tick labels and the x variable (i.e.: dataset) becomes the grouping variable. To do this, use the argument merge = “flip”.

ggboxplot(expr, x = "dataset",
          y = c("GATA3", "PTEN", "XBP1"),
          merge = "flip",
          ylab = "Expression", 
          palette = "jco")

You might want to add jittered points on the boxplot. Each point correspond to individual observations. To add jittered points, use the argument add = “jitter” as follow. To customize the added elements, specify the argument add.params.

ggboxplot(expr, x = "dataset",
          y = c("GATA3", "PTEN", "XBP1"),
          combine = TRUE,
          color = "dataset", palette = "jco",
          ylab = "Expression", 
          add = "jitter",                              # Add jittered points
          add.params = list(size = 0.1, jitter = 0.2)  # Point size and the amount of jittering
          )

You can add and adjust a dotplot as follow:

ggboxplot(expr, x = "dataset",
          y = c("GATA3", "PTEN", "XBP1"),
          combine = TRUE,
          color = "dataset", palette = "jco",
          ylab = "Expression", 
          add = "dotplot",                              # Add dotplot
          add.params = list(binwidth = 0.1, dotsize = 0.3)
          )

You might want to label the boxplot by showing the names of samples with the top n highest or lowest values. In this case, you can use the following arguments:

• label: the name of the column containing point labels.
• label.select: can be of two formats:
  • a character vector specifying some labels to show.
  • a list containing one or the combination of the following components:
    • top.up and top.down: to display the labels of the top up/down points. For example, label.select = list(top.up = 10, top.down = 4).
    • criteria: to filter, for example, by x and y variables values, use this: label.select = list(criteria = “`y` > 3.9 & `y` < 5 & `x` %in% c(‘BRCA’, ‘OV’)”).

For example:

ggboxplot(expr, x = "dataset",
          y = c("GATA3", "PTEN", "XBP1"),
          combine = TRUE,
          color = "dataset", palette = "jco",
          ylab = "Expression", 
          add = "jitter",                               # Add jittered points
          add.params = list(size = 0.1, jitter = 0.2),  # Point size and the amount of jittering
          label = "bcr_patient_barcode",                # column containing point labels
          label.select = list(top.up = 2, top.down = 2),# Select some labels to display
          font.label = list(size = 9, face = "italic"), # label font
          repel = TRUE                                  # Avoid label text overplotting
          )

A complex criteria for labeling can be specified as follow:

label.select.criteria <- list(criteria = "`y` > 3.9 & `x` %in% c('BRCA', 'OV')")
ggboxplot(expr, x = "dataset",
          y = c("GATA3", "PTEN", "XBP1"),
          combine = TRUE,
          color = "dataset", palette = "jco",
          ylab = "Expression", 
          label = "bcr_patient_barcode",              # column containing point labels
          label.select = label.select.criteria,       # Select some labels to display
          font.label = list(size = 9, face = "italic"), # label font
          repel = TRUE                                # Avoid label text overplotting
          )

Violin plots

The following R code draws violin plots with box plots inside:

ggviolin(expr, x = "dataset",
          y = c("GATA3", "PTEN", "XBP1"),
          combine = TRUE, 
          color = "dataset", palette = "jco",
          ylab = "Expression", 
          add = "boxplot")

Instead of adding a box plot inside the violin plot, you can add the median + interquantile range as follow:

ggviolin(expr, x = "dataset",
          y = c("GATA3", "PTEN", "XBP1"),
          combine = TRUE, 
          color = "dataset", palette = "jco",
          ylab = "Expression", 
          add = "median_iqr")

Stripcharts and dot plots

To draw a stripchart, type this:

ggstripchart(expr, x = "dataset",
             y = c("GATA3", "PTEN", "XBP1"),
             combine = TRUE, 
             color = "dataset", palette = "jco",
             size = 0.1, jitter = 0.2,
             ylab = "Expression", 
             add = "median_iqr",
             add.params = list(color = "gray"))

For a dot plot, use this:

ggdotplot(expr, x = "dataset",
          y = c("GATA3", "PTEN", "XBP1"),
          combine = TRUE, 
          color = "dataset", palette = "jco",
          fill = "white",
          binwidth = 0.1,
          ylab = "Expression", 
          add = "median_iqr",
          add.params = list(size = 0.9))

Density plots

To visualize the distribution as a density plot, use the function ggdensity() as follow:

# Basic density plot
ggdensity(expr,
       x = c("GATA3", "PTEN",  "XBP1"),
       y = "..density..",
       combine = TRUE,                  # Combine the 3 plots
       xlab = "Expression", 
       add = "median",                  # Add median line. 
       rug = TRUE                       # Add marginal rug
)

# Change color and fill by dataset
ggdensity(expr,
       x = c("GATA3", "PTEN",  "XBP1"),
       y = "..density..",
       combine = TRUE,                  # Combine the 3 plots
       xlab = "Expression", 
       add = "median",                  # Add median line. 
       rug = TRUE,                      # Add marginal rug
       color = "dataset", 
       fill = "dataset",
       palette = "jco"
)

# Merge the 3 plots
# and use y = "..count.." instead of "..density.."
ggdensity(expr,
       x = c("GATA3", "PTEN",  "XBP1"),
       y = "..count..",
       merge = TRUE,                    # Merge the 3 plots
       xlab = "Expression", 
       add = "median",                  # Add median line. 
       rug = TRUE ,                     # Add marginal rug
       palette = "jco"                  # Change color palette
)

# color and fill by x variables
ggdensity(expr,
       x = c("GATA3", "PTEN",  "XBP1"),
       y = "..count..",
       color = ".x.", fill = ".x.",     # color and fill by x variables
       merge = TRUE,                    # Merge the 3 plots
       xlab = "Expression", 
       add = "median",                  # Add median line. 
       rug = TRUE ,                     # Add marginal rug
       palette = "jco"                  # Change color palette
)

# Facet by "dataset"
ggdensity(expr,
       x = c("GATA3", "PTEN",  "XBP1"),
       y = "..count..",
       color = ".x.", fill = ".x.", 
       facet.by = "dataset",            # Split by "dataset" into multi-panel
       merge = TRUE,                    # Merge the 3 plots
       xlab = "Expression", 
       add = "median",                  # Add median line. 
       rug = TRUE ,                     # Add marginal rug
       palette = "jco"                  # Change color palette
)

Histogram plots

To visualize the distribution as a histogram plot, use the function gghistogram() as follow:

# Basic histogram plot 
gghistogram(expr,
       x = c("GATA3", "PTEN",  "XBP1"),
       y = "..density..",
       combine = TRUE,                  # Combine the 3 plots
       xlab = "Expression", 
       add = "median",                  # Add median line. 
       rug = TRUE                       # Add marginal rug
)

# Change color and fill by dataset
gghistogram(expr,
       x = c("GATA3", "PTEN",  "XBP1"),
       y = "..density..",
       combine = TRUE,                  # Combine the 3 plots
       xlab = "Expression", 
       add = "median",                  # Add median line. 
       rug = TRUE,                      # Add marginal rug
       color = "dataset", 
       fill = "dataset",
       palette = "jco"
)

# Merge the 3 plots
# and use y = "..count.." instead of "..density.."
gghistogram(expr,
       x = c("GATA3", "PTEN",  "XBP1"),
       y = "..count..",
       merge = TRUE,                    # Merge the 3 plots
       xlab = "Expression", 
       add = "median",                  # Add median line. 
       rug = TRUE ,                     # Add marginal rug
       palette = "jco"                  # Change color palette
)

# color and fill by x variables
gghistogram(expr,
       x = c("GATA3", "PTEN",  "XBP1"),
       y = "..count..",
       color = ".x.", fill = ".x.",     # color and fill by x variables
       merge = TRUE,                    # Merge the 3 plots
       xlab = "Expression", 
       add = "median",                  # Add median line. 
       rug = TRUE ,                     # Add marginal rug
       palette = "jco"                  # Change color palette
)

# Facet by "dataset"
gghistogram(expr,
       x = c("GATA3", "PTEN",  "XBP1"),
       y = "..count..",
       color = ".x.", fill = ".x.", 
       facet.by = "dataset",            # Split by "dataset" into multi-panel
       merge = TRUE,                    # Merge the 3 plots
       xlab = "Expression", 
       add = "median",                  # Add median line. 
       rug = TRUE ,                     # Add marginal rug
       palette = "jco"                  # Change color palette
)

Empirical cumulative density function

# Basic ECDF plot 
ggecdf(expr,
       x = c("GATA3", "PTEN",  "XBP1"),
       combine = TRUE,                 
       xlab = "Expression", ylab = "F(expression)"
)

# Change color  by dataset
ggecdf(expr,
       x = c("GATA3", "PTEN",  "XBP1"),
       combine = TRUE,                 
       xlab = "Expression", ylab = "F(expression)",
       color = "dataset", palette = "jco"
)

# Merge the 3 plots and color by x variables
ggecdf(expr,
       x = c("GATA3", "PTEN",  "XBP1"),
       merge = TRUE,                 
       xlab = "Expression", ylab = "F(expression)",
       color = ".x.", palette = "jco"
)

# Merge the 3 plots and color by x variables
# facet by "dataset" into multi-panel
ggecdf(expr,
       x = c("GATA3", "PTEN",  "XBP1"),
       merge = TRUE,                 
       xlab = "Expression", ylab = "F(expression)",
       color = ".x.", palette = "jco",
       facet.by = "dataset"
)

Quantile - Quantile plot

# Basic ECDF plot 
ggqqplot(expr,
       x = c("GATA3", "PTEN",  "XBP1"),
       combine = TRUE, size = 0.5
)

# Change color  by dataset
ggqqplot(expr,
       x = c("GATA3", "PTEN",  "XBP1"),
       combine = TRUE, color = "dataset", palette = "jco",
       size = 0.5
)

# Merge the 3 plots and color by x variables
ggqqplot(expr,
       x = c("GATA3", "PTEN",  "XBP1"),
       merge = TRUE,  
       color = ".x.", palette = "jco"
)

# Merge the 3 plots and color by x variables
# facet by "dataset" into multi-panel
ggqqplot(expr,
       x = c("GATA3", "PTEN",  "XBP1"),
       merge = TRUE, size = 0.5,
       color = ".x.", palette = "jco",
       facet.by = "dataset"
)

Prerequisites

install.packages("ggpubr")

if(!require(devtools)) install.packages("devtools")
devtools::install_github("kassambara/ggpubr")

library(ggpubr)

data("ToothGrowth")
head(ToothGrowth)

##    len supp dose
## 1  4.2   VC  0.5
## 2 11.5   VC  0.5
## 3  7.3   VC  0.5
## 4  5.8   VC  0.5
## 5  6.4   VC  0.5
## 6 10.0   VC  0.5

Methods for comparing means

The standard methods to compare the means of two or more groups in R, have been largely described at: comparing means in R.

The most common methods for comparing means include:

A practical guide to compute and interpret the results of each of these methods are provided at the following links:

R functions to add p-values

Here we present two new R functions in the ggpubr package:

• compare_means(): easy to use solution to performs one and multiple mean comparisons.
• stat_compare_means(): easy to use solution to automatically add p-values and significance levels to a ggplot.

compare_means(formula, data, method = "wilcox.test", paired = FALSE,
  group.by = NULL, ref.group = NULL, ...)

stat_compare_means(mapping = NULL, comparisons = NULL hide.ns = FALSE,
                   label = NULL,  label.x = NULL, label.y = NULL,  ...)

Compare two independent groups

Perform the test:

compare_means(len ~ supp, data = ToothGrowth)

## # A tibble: 1 x 8
##     .y. group1 group2      p  p.adj p.format p.signif   method
##                       
## 1   len     OJ     VC 0.0645 0.0645    0.064       ns Wilcoxon

Returned value is a data frame with the following columns:

• .y.: the y variable used in the test.
• p: the p-value
• p.adj: the adjusted p-value. Default value for p.adjust.method = “holm”
• p.format: the formatted p-value
• p.signif: the significance level.
• method: the statistical test used to compare groups.

Create a box plot with p-values:

p <- ggboxplot(ToothGrowth, x = "supp", y = "len",
          color = "supp", palette = "jco",
          add = "jitter")
#  Add p-value
p + stat_compare_means()
# Change method
p + stat_compare_means(method = "t.test")

Note that, the p-value label position can be adjusted using the arguments: label.x, label.y, hjust and vjust.

The default p-value label displayed is obtained by concatenating the method and the p columns of the returned data frame by the function compare_means(). You can specify other combinations using the aes() function.

For example,

• aes(label = ..p.format..) or aes(label = paste0(“p =”, ..p.format..)): display only the formatted p-value (without the method name)
• aes(label = ..p.signif..): display only the significance level.
• aes(label = paste0(..method.., “\n”, “p =”, ..p.format..)): Use line break (“\n”) between the method name and the p-value.

As an illustration, type this:

p + stat_compare_means( aes(label = ..p.signif..), 
                        label.x = 1.5, label.y = 40)

If you prefer, it’s also possible to specify the argument label as a character vector:

p + stat_compare_means( label = "p.signif", label.x = 1.5, label.y = 40)

Compare two paired samples

Perform the test:

compare_means(len ~ supp, data = ToothGrowth, paired = TRUE)

## # A tibble: 1 x 8
##     .y. group1 group2       p   p.adj p.format p.signif   method
##                         
## 1   len     OJ     VC 0.00431 0.00431   0.0043       ** Wilcoxon

Visualize paired data using the ggpaired() function:

ggpaired(ToothGrowth, x = "supp", y = "len",
         color = "supp", line.color = "gray", line.size = 0.4,
         palette = "jco")+
  stat_compare_means(paired = TRUE)

Compare more than two groups

• Global test:

# Global test
compare_means(len ~ dose,  data = ToothGrowth, method = "anova")

## # A tibble: 1 x 6
##     .y.        p    p.adj p.format p.signif method
##                     
## 1   len 9.53e-16 9.53e-16  9.5e-16     ****  Anova

Plot with global p-value:

# Default method = "kruskal.test" for multiple groups
ggboxplot(ToothGrowth, x = "dose", y = "len",
          color = "dose", palette = "jco")+
  stat_compare_means()
# Change method to anova
ggboxplot(ToothGrowth, x = "dose", y = "len",
          color = "dose", palette = "jco")+
  stat_compare_means(method = "anova")

• Pairwise comparisons. If the grouping variable contains more than two levels, then pairwise tests will be performed automatically. The default method is “wilcox.test”. You can change this to “t.test”.

# Perorm pairwise comparisons
compare_means(len ~ dose,  data = ToothGrowth)

## # A tibble: 3 x 8
##     .y. group1 group2        p    p.adj p.format p.signif   method
##                           
## 1   len    0.5      1 7.02e-06 1.40e-05  7.0e-06     **** Wilcoxon
## 2   len    0.5      2 8.41e-08 2.52e-07  8.4e-08     **** Wilcoxon
## 3   len      1      2 1.77e-04 1.77e-04  0.00018      *** Wilcoxon

# Visualize: Specify the comparisons you want
my_comparisons <- list( c("0.5", "1"), c("1", "2"), c("0.5", "2") )
ggboxplot(ToothGrowth, x = "dose", y = "len",
          color = "dose", palette = "jco")+ 
  stat_compare_means(comparisons = my_comparisons)+ # Add pairwise comparisons p-value
  stat_compare_means(label.y = 50)     # Add global p-value

If you want to specify the precise y location of bars, use the argument label.y:

ggboxplot(ToothGrowth, x = "dose", y = "len",
          color = "dose", palette = "jco")+ 
  stat_compare_means(comparisons = my_comparisons, label.y = c(29, 35, 40))+
  stat_compare_means(label.y = 45)

(Adding bars, connecting compared groups, has been facilitated by the ggsignif R package )

• Multiple pairwise tests against a reference group:

# Pairwise comparison against reference
compare_means(len ~ dose,  data = ToothGrowth, ref.group = "0.5",
              method = "t.test")

## # A tibble: 2 x 8
##     .y. group1 group2        p    p.adj p.format p.signif method
##                         
## 1   len    0.5      1 6.70e-09 6.70e-09  6.7e-09     **** T-test
## 2   len    0.5      2 1.47e-16 2.94e-16  < 2e-16     **** T-test

# Visualize
ggboxplot(ToothGrowth, x = "dose", y = "len",
          color = "dose", palette = "jco")+
  stat_compare_means(method = "anova", label.y = 40)+      # Add global p-value
  stat_compare_means(label = "p.signif", method = "t.test",
                     ref.group = "0.5")                    # Pairwise comparison against reference

• Multiple pairwise tests against all (base-mean):

# Comparison of each group against base-mean
compare_means(len ~ dose,  data = ToothGrowth, ref.group = ".all.",
              method = "t.test")

## # A tibble: 3 x 8
##     .y. group1 group2        p    p.adj p.format p.signif method
##                         
## 1   len  .all.    0.5 1.24e-06 3.73e-06  1.2e-06     **** T-test
## 2   len  .all.      1 5.67e-01 5.67e-01     0.57       ns T-test
## 3   len  .all.      2 1.37e-05 2.74e-05  1.4e-05     **** T-test

# Visualize
ggboxplot(ToothGrowth, x = "dose", y = "len",
          color = "dose", palette = "jco")+
  stat_compare_means(method = "anova", label.y = 40)+      # Add global p-value
  stat_compare_means(label = "p.signif", method = "t.test",
                     ref.group = ".all.")                  # Pairwise comparison against all

A typical situation, where pairwise comparisons against “all” can be useful, is illustrated here using the myeloma data set available on Github.

We’ll plot the expression profile of the DEPDC1 gene according to the patients’ molecular groups. We want to know if there is any difference between groups. If yes, where the difference is?

To answer to this question, you can perform a pairwise comparison between all the 7 groups. This will lead to a lot of comparisons between all possible combinations. If you have many groups, as here, it might be difficult to interpret.

Another easy solution is to compare each of the seven groups against “all” (i.e. base-mean). When the test is significant, then you can conclude that DEPDC1 is significantly overexpressed or downexpressed in a group xxx compared to all.

# Load myeloma data from GitHub
myeloma <- read.delim("https://raw.githubusercontent.com/kassambara/data/master/myeloma.txt")
# Perform the test
compare_means(DEPDC1 ~ molecular_group,  data = myeloma,
              ref.group = ".all.", method = "t.test")

## # A tibble: 7 x 8
##      .y. group1           group2        p   p.adj p.format p.signif method
##                                   
## 1 DEPDC1  .all.       Cyclin D-1 0.149690 0.44907  0.14969       ns T-test
## 2 DEPDC1  .all.       Cyclin D-2 0.523143 1.00000  0.52314       ns T-test
## 3 DEPDC1  .all.     Hyperdiploid 0.000282 0.00169  0.00028      *** T-test
## 4 DEPDC1  .all. Low bone disease 0.005084 0.02542  0.00508       ** T-test
## 5 DEPDC1  .all.              MAF 0.086107 0.34443  0.08611       ns T-test
## 6 DEPDC1  .all.            MMSET 0.576291 1.00000  0.57629       ns T-test
## # ... with 1 more rows

# Visualize the expression profile
ggboxplot(myeloma, x = "molecular_group", y = "DEPDC1", color = "molecular_group", 
          add = "jitter", legend = "none") +
  rotate_x_text(angle = 45)+
  geom_hline(yintercept = mean(myeloma$DEPDC1), linetype = 2)+ # Add horizontal line at base mean
  stat_compare_means(method = "anova", label.y = 1600)+        # Add global annova p-value
  stat_compare_means(label = "p.signif", method = "t.test",
                     ref.group = ".all.")                      # Pairwise comparison against all

# Visualize the expression profile
ggboxplot(myeloma, x = "molecular_group", y = "DEPDC1", color = "molecular_group", 
          add = "jitter", legend = "none") +
  rotate_x_text(angle = 45)+
  geom_hline(yintercept = mean(myeloma$DEPDC1), linetype = 2)+ # Add horizontal line at base mean
  stat_compare_means(method = "anova", label.y = 1600)+        # Add global annova p-value
  stat_compare_means(label = "p.signif", method = "t.test",
                     ref.group = ".all.", hide.ns = TRUE)      # Pairwise comparison against all

Multiple grouping variables

• Two independent sample comparisons after grouping the data by another variable:

Perform the test:

compare_means(len ~ supp, data = ToothGrowth, 
              group.by = "dose")

## # A tibble: 3 x 9
##    dose   .y. group1 group2       p  p.adj p.format p.signif   method
##                         
## 1   0.5   len     OJ     VC 0.02319 0.0464    0.023        * Wilcoxon
## 2   1.0   len     OJ     VC 0.00403 0.0121    0.004       ** Wilcoxon
## 3   2.0   len     OJ     VC 1.00000 1.0000    1.000       ns Wilcoxon

Visualize (1/2). Create a multi-panel box plots facetted by group (here, “dose”):

# Box plot facetted by "dose"
p <- ggboxplot(ToothGrowth, x = "supp", y = "len",
          color = "supp", palette = "jco",
          add = "jitter",
          facet.by = "dose", short.panel.labs = FALSE)
# Use only p.format as label. Remove method name.
p + stat_compare_means(label = "p.format")

# Or use significance symbol as label
p + stat_compare_means(label =  "p.signif", label.x = 1.5)

Visualize (2/2). Create one single panel with all box plots. Plot y = “len” by x = “dose” and color by “supp”:

p <- ggboxplot(ToothGrowth, x = "dose", y = "len",
          color = "supp", palette = "jco",
          add = "jitter")
p + stat_compare_means(aes(group = supp))

# Show only p-value
p + stat_compare_means(aes(group = supp), label = "p.format")

# Use significance symbol as label
p + stat_compare_means(aes(group = supp), label = "p.signif")

• Paired sample comparisons after grouping the data by another variable:

Perform the test:

compare_means(len ~ supp, data = ToothGrowth, 
              group.by = "dose", paired = TRUE)

## # A tibble: 3 x 9
##    dose   .y. group1 group2      p  p.adj p.format p.signif   method
##                        
## 1   0.5   len     OJ     VC 0.0330 0.0659    0.033        * Wilcoxon
## 2   1.0   len     OJ     VC 0.0191 0.0572    0.019        * Wilcoxon
## 3   2.0   len     OJ     VC 1.0000 1.0000    1.000       ns Wilcoxon

Visualize. Create a multi-panel box plots facetted by group (here, “dose”):

# Box plot facetted by "dose"
p <- ggpaired(ToothGrowth, x = "supp", y = "len",
          color = "supp", palette = "jco", 
          line.color = "gray", line.size = 0.4,
          facet.by = "dose", short.panel.labs = FALSE)
# Use only p.format as label. Remove method name.
p + stat_compare_means(label = "p.format", paired = TRUE)

Other plot types

• Bar and line plots (one grouping variable):

# Bar plot of mean +/-se
ggbarplot(ToothGrowth, x = "dose", y = "len", add = "mean_se")+
  stat_compare_means() +                                         # Global p-value
  stat_compare_means(ref.group = "0.5", label = "p.signif",
                     label.y = c(22, 29))                   # compare to ref.group
# Line plot of mean +/-se
ggline(ToothGrowth, x = "dose", y = "len", add = "mean_se")+
  stat_compare_means() +                                         # Global p-value
  stat_compare_means(ref.group = "0.5", label = "p.signif",
                     label.y = c(22, 29))     

• Bar and line plots (two grouping variables):

ggbarplot(ToothGrowth, x = "dose", y = "len", add = "mean_se",
          color = "supp", palette = "jco", 
          position = position_dodge(0.8))+
  stat_compare_means(aes(group = supp), label = "p.signif", label.y = 29)
ggline(ToothGrowth, x = "dose", y = "len", add = "mean_se",
          color = "supp", palette = "jco")+
  stat_compare_means(aes(group = supp), label = "p.signif", 
                     label.y = c(16, 25, 29))