Data Visualization Workshop
+Han Yu, Meghana Munipalle
+ +Introduction
+Data visualization is a crucial skill in data science and analytics. +In this tutorial, we will explore how to create various types of +visualizations in R using the ggplot2 package, a powerful tool for +creating informative and aesthetically pleasing graphics.
+Installing and Loading Required Packages
+Before we begin, ensure you have the necessary packages installed. To +avoid any errors, install these at minimum to ensure their most recent +versions:
+install.packages(c("xfun", "knitr", "rmarkdown"))Then check that these packages are here:
+library(dplyr)
+library(tidyverse)
+library(ggplot2)
+library(ggridges)
+library(datasauRus)
+library(colorspace)
+library(pheatmap)
+library(survival)
+library(survminer)
+library(qqman)
+library(plotly)
+library(gganimate)
+library(gifski)
+library(palmerpenguins)
+library(gapminder)Module 1: Brief Review of Basic Visualization in R
+Let’s start with the simple built-in penguins data set to get +familiar with the structure of ggplot2.
+data(penguins)
+head(penguins)penguins <- penguins %>% drop_na()
+
+# Basic scatter plot
+ggplot(penguins) +
+ geom_point(aes(x = bill_length_mm,
+ y = bill_depth_mm))
For this workshop, we will mostly use a “simple” RNA-seq +dataset on the transcriptional response to SARS-CoV-2 infection in +eye tissue from human donors.
+Run this cell first to do some quick data cleanup.
+counts <- read_csv("data/GSE164073_Eye_count_matrix.csv.gz") # reads .gz directly## Rows: 27946 Columns: 19
+## ── Column specification ────────────────────────────────────────────────────────
+## Delimiter: ","
+## chr (1): Gene
+## dbl (18): MW1_cornea_mock_1, MW2_cornea_mock_2, MW3_cornea_mock_3, MW4_corne...
+##
+## ℹ Use `spec()` to retrieve the full column specification for this data.
+## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.# Some genes have zero reads in every sample - let's drop them
+counts <- counts %>% filter(rowSums(across(-Gene)) > 0)
+
+expr <- counts %>%
+ pivot_longer(-Gene, names_to = "sample", values_to = "count") %>%
+ separate_wider_delim(sample, delim = "_",
+ names = c("id", "tissue", "condition", "replicate"),
+ cols_remove = FALSE) %>%
+ mutate(condition = factor(condition, levels = c("mock", "CoV2"))) %>%
+ group_by(sample) %>%
+ mutate(cpm = count / sum(count) * 1e6, # counts per million
+ log_cpm = log2(cpm + 1)) %>%
+ ungroup()Now let’s explore the dataset.
+head(expr, 5) # check the first 5 rowsWe can also create a per-sample summary table. The whole dataset is +one row per gene per sample, which is a lot (tens of thousands) of rows +for each sample. Here, we take the 18 distinct samples and summarize +their info into one row for each.
+samples <- expr %>%
+ group_by(sample, tissue, condition, replicate) %>%
+ summarise(library_size = sum(count), genes_detected = sum(count > 0), .groups = "drop")We want to group by sample, but we also want to keep tissue,
+condition, and replicate in the output. So, we use
+group_by(sample, tissue, condition, replicate). In
+summarize, we condense the info for each sample into
+library_size - how many reads across every gene in the
+sample?
We’ll use both samples and expr for our
+visualizations depending on what we want to compare going forward. Now,
+let’s use the dataset to build different useful plots!
1.1 bar plot
+Bar plot is used to visually compare categorical data by representing +values with rectangular bars.
+# bar plot: discrete variable
+ggplot(samples, aes(x = sample, y = library_size / 1e6, fill = condition)) +
+ geom_col() +
+ theme(axis.text.x = element_text(angle = 45, hjust = 1),
+ plot.margin = margin(t = 5, r = 5, b = 5, l = 40))
1.2 histogram/ frequency/ density plot
+A histogram, frequency plot, and density plot all visualize the +distribution of a dataset, but they differ in how they represent data: a +histogram groups data into bins and shows counts, a frequency plot +connects points representing counts for each bin, and a density plot +estimates the probability density function, providing a smooth +representation of distribution shape.
+# histogram plot: continuous variable
+ggplot(expr, aes(x = log_cpm)) + geom_histogram(binwidth = 0.5)
# frequency plot
+ggplot(expr, aes(x = log_cpm, colour = tissue)) + geom_freqpoly(binwidth = 0.5)
# density plot (smooth histogram = density plot)
+ggplot(expr, aes(x = log_cpm, colour = sample)) + geom_density()
1.3 box plot
+A box plot is used to summarize the distribution of a dataset by +displaying its median, quartiles, and potential outliers, helping to +identify variability and skewness.
+# box plot: one discrete variable + one continuous variable
+ggplot(expr, aes(x = sample, y = log_cpm, fill = tissue)) +
+ geom_boxplot() +
+ theme(axis.text.x = element_text(angle = 45, hjust = 1))
1.4 violin plot
+A violin plot combines aspects of a box plot and a density plot to +show the distribution, probability density, and variability of a +dataset, making it easier to identify patterns and multimodal +distributions.
+# violin plot
+ggplot(expr, aes(x = tissue, y = log_cpm, fill = tissue)) + geom_violin(trim = FALSE, alpha = 0.5)
stat_summary(
+ fun = mean,
+ fun.max = function(x) {
+ mean(x) + sd(x)
+ },
+ fun.min = function(x) {
+ mean(x) - sd(x)
+ },
+ geom = "pointrange"
+ )## geom_pointrange: na.rm = FALSE, orientation = NA
+## stat_summary: fun.data = NULL, fun = function (x, ...)
+## UseMethod("mean"), fun.max = function (x)
+## {
+## mean(x) + sd(x)
+## }, fun.min = function (x)
+## {
+## mean(x) - sd(x)
+## }, fun.args = list(), na.rm = FALSE, orientation = NA
+## position_identity1.5 Example 1
+Let’s see an example of how to improve a data visualization plot step +by step. We’ll look at our RNA-seq dataset to investigate how SARS-CoV-2 +infection changes the expression of one gene, SOD2, in +three tissues of the human eye.
+SOD2 (superoxide dismutase 2) is an enzyme that protects cells from +oxidative stress and is commonly switched on during inflammation. It is +one of the genes that changes most clearly after infection in this +dataset.
+1.5.1 Create the dataset we will use in this example
+The count matrix has one row per gene and one column per sample. To +compare samples fairly, we convert counts to counts per million (CPM) +(correcting for each sample’s sequencing depth) and then log-transform +them. We need every gene to calculate each sample’s library size, but we +only keep SOD2 for this example.
+# total reads per sample (library size), named by sample
+lib_size <- colSums(counts[, -1])
+
+sod2_df <- counts %>%
+ filter(Gene == "SOD2") %>%
+ pivot_longer(-Gene, names_to = "sample", values_to = "count") %>%
+ # sample names look like "MW4_cornea_CoV2_1"
+ separate_wider_delim(
+ sample,
+ delim = "_",
+ names = c("id", "tissue", "condition", "replicate"),
+ cols_remove = FALSE
+ ) %>%
+ mutate(
+ log_cpm = log2(count / lib_size[sample] * 1e6 + 1),
+ tissue = str_to_title(tissue),
+ condition = fct_recode(
+ factor(condition, levels = c("mock", "CoV2")),
+ Mock = "mock",
+ `SARS-CoV-2` = "CoV2"
+ ),
+ # one x-axis group per tissue x condition, in a sensible order
+ group = factor(
+ paste(tissue, condition, sep = "\n"),
+ levels = paste(rep(c("Cornea", "Limbus", "Sclera"), each = 2),
+ c("Mock", "SARS-CoV-2"), sep = "\n")
+ )
+ ) %>%
+ select(sample, tissue, condition, group, log_cpm)
+
+head(sod2_df)Let’s also define a theme once here so that we don’t have to keep +repeating it for following code blocks:
+theme_example <- theme_classic(base_size = 14) +
+ theme(
+ axis.text = element_text(color = "black", size = 12),
+ plot.margin = margin(3, 7, 3, 1.5)
+ )
+
+y_label <- "SOD2 expression (log2 CPM + 1)"1.5.2 What’s the problem of the following plot?
+# points-errorbars
+sod2_errbar <- sod2_df %>%
+ ggplot(aes(x = group, y = log_cpm)) +
+ stat_summary(
+ fun = mean,
+ fun.max = function(x) {
+ mean(x) + 2 * sd(x)
+ },
+ fun.min = function(x) {
+ mean(x) - 2 * sd(x)
+ },
+ geom = "pointrange"
+ ) +
+ xlab(NULL) +
+ ylab(y_label) +
+ theme_example
+
+sod2_errbar
1.5.3 Use box/violin plot to show the distribution
+# box plot
+sod2_box <- sod2_df %>%
+ ggplot(aes(x = group, y = log_cpm)) +
+ geom_boxplot(fill = "grey90") +
+ xlab(NULL) +
+ ylab(y_label) +
+ theme_example
+
+sod2_box
# violin plot
+sod2_violin <- sod2_df %>%
+ ggplot(aes(x = group, y = log_cpm)) +
+ geom_violin(fill = "grey90") +
+ xlab(NULL) +
+ ylab(y_label) +
+ theme_example
+
+sod2_violin
1.5.4 Use points to show the original data
+# points
+sod2_points <- sod2_df %>%
+ ggplot(aes(x = group, y = log_cpm)) +
+ geom_point(size = 2) +
+ xlab(NULL) +
+ ylab(y_label) +
+ theme_example
+
+sod2_points
# jitter points
+sod2_jitter <- sod2_df %>%
+ ggplot(aes(x = group, y = log_cpm)) +
+ geom_point(
+ position = position_jitter(width = .15, height = 0, seed = 320),
+ size = 2
+ ) +
+ xlab(NULL) +
+ ylab(y_label) +
+ theme_example
+
+sod2_jitter
1.5.5 Use sina plot to combine violin plot and jitter points?
+What if we combine two things, and add in the data points to a violin +plot to see if the points show the distribution’s shape?
+# violin + jitter
+sod2_violin_jitter <- sod2_df %>%
+ ggplot(aes(x = group, y = log_cpm)) +
+ geom_violin(color = "transparent", fill = "gray90") +
+ geom_jitter(width = 0.15, height = 0, size = 2) +
+ xlab(NULL) +
+ ylab(y_label) +
+ theme_example
+
+sod2_violin_jitter
With only 3 points per group, the violin plot is shapes that aren’t +necessarily reflecting the data. With small sample sizes, using a +technique that shows the raw points is the best. Let’s try something +else.
+1.5.6 A better plot for small groups: points with a mean
+Instead of a violin, we can add a short line at each group’s mean and +use color to separate the conditions.
+sod2_final <- sod2_df %>%
+ ggplot(aes(x = group, y = log_cpm, color = condition)) +
+ stat_summary(
+ fun = mean,
+ geom = "crossbar",
+ width = 0.5,
+ linewidth = 0.4,
+ color = "grey40"
+ ) +
+ geom_point(
+ position = position_jitter(width = .1, height = 0, seed = 320),
+ size = 2.5
+ ) +
+ scale_color_manual(values = c(Mock = "grey50", `SARS-CoV-2` = "firebrick")) +
+ xlab(NULL) +
+ ylab(y_label) +
+ labs(
+ color = NULL,
+ caption = "Points are individual samples (n = 3 per group). Lines are group means."
+ ) +
+ theme_example +
+ theme(legend.position = "top")
+
+sod2_final
Compared with the earlier plot, the reader can now see every +measurement, how many samples there are, and how consistent the +replicates are.
+Module 2: Advanced Customization and Data Transformation
+2.0 Prepare the dataset
+Let’s continue using the same RNA-seq dataset, but we’ll create one +more small table to use in our plotting. This one is genes by sample and +each gene’s variance across samples.
+library(tidyverse)
+library(colorspace)# long format: one row per gene per sample
+expr <- counts %>%
+ pivot_longer(-Gene, names_to = "sample", values_to = "count") %>%
+ separate_wider_delim(
+ sample,
+ delim = "_",
+ names = c("id", "tissue", "condition", "replicate"),
+ cols_remove = FALSE
+ ) %>%
+ mutate(condition = factor(condition, levels = c("mock", "CoV2"))) %>%
+ group_by(sample) %>%
+ mutate(
+ cpm = count / sum(count) * 1e6, # counts per million
+ log_cpm = log2(cpm + 1)
+ ) %>%
+ ungroup() %>%
+ mutate(group = factor(
+ paste(tissue, condition),
+ levels = c("cornea mock", "cornea CoV2",
+ "limbus mock", "limbus CoV2",
+ "sclera mock", "sclera CoV2")
+ ))
+
+# one row per sample
+samples <- expr %>%
+ group_by(sample, tissue, condition, replicate) %>%
+ summarise(
+ library_size = sum(count) / 1e6, # millions of reads
+ genes_detected = sum(count > 0),
+ .groups = "drop"
+ )
+
+# gene x sample matrix of log_cpm, and each gene's variance across samples
+expr_mat <- expr %>%
+ select(Gene, sample, log_cpm) %>%
+ pivot_wider(names_from = sample, values_from = log_cpm) %>%
+ column_to_rownames("Gene") %>%
+ as.matrix()
+
+gene_var <- apply(expr_mat, 1, var)2.1 Data Transformations and Faceting
+In data visualization, data transformation is a crucial step to +ensure that the data is in the right format and structure for effective +visual representation.
+2.1.1 Log Transformation
+There are many different data transformation methods, and the +following is an example of log transformation.
+We would like to investigate the relationship between how highly a +gene is expressed and how much it changes after SARS-CoV-2 infection, +across different tissues. (This is also called an MA plot.)
+# Step 1: calculate the mean CPM of each gene in mock and infected samples, per tissue
+gene_summary <- expr %>%
+ group_by(tissue, Gene, condition) %>%
+ summarise(mean_cpm = mean(cpm), .groups = "drop") %>%
+ pivot_wider(names_from = condition, values_from = mean_cpm) %>%
+ mutate(
+ avg_cpm = (mock + CoV2) / 2, # average expression
+ log2FC = log2((CoV2 + 1) / (mock + 1)) # log2 fold change, infected vs mock
+ ) %>%
+ filter(avg_cpm >= 1) # drop barely expressed genes
+gene_summary# Step 2: draw the plot
+gene_summary %>%
+ ggplot(aes(x = avg_cpm, y = log2FC)) +
+ geom_point(aes(color = tissue), size = 0.5, alpha = 0.3)
# Step 3: log transformation of x-axis (avg_cpm)
+gene_summary %>%
+ ggplot(aes(x = avg_cpm, y = log2FC)) +
+ geom_point(aes(color = tissue), size = 0.5, alpha = 0.3) +
+ scale_x_log10()
# Step 4 (Optional): add comma format or unit format to label the axis for better interpretation
+gene_summary %>%
+ ggplot(aes(x = avg_cpm, y = log2FC)) +
+ geom_point(aes(color = tissue), size = 0.5, alpha = 0.3) +
+ scale_x_log10(breaks = c(1, 10, 100, 1000, 10000),
+ labels = scales::label_comma())
gene_summary %>%
+ ggplot(aes(x = avg_cpm, y = log2FC)) +
+ geom_point(aes(color = tissue), size = 0.5, alpha = 0.3) +
+ scale_x_log10(
+ name = "Average expression",
+ breaks = c(1, 10, 100, 1000, 10000),
+ labels = scales::unit_format(unit = "CPM"))
2.1.2 Faceting
+Data faceting is used in data visualization to break down data into +subsets and display them across multiple smaller plots. This technique +helps in revealing patterns, trends, and relationships that might be +obscured in a single, aggregated plot. The following is an example of +data faceting.
+# facet_grid(): create a grid of graphs, by rows and columns
+expr %>%
+ ggplot(aes(x = log_cpm)) +
+ geom_density() +
+ facet_grid(. ~ tissue)
expr %>%
+ ggplot(aes(x = log_cpm)) +
+ geom_histogram() +
+ facet_grid(tissue ~ .)## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.
# facet_wrap(): create small multiples by "wrapping" a series of plots
+gene_summary %>%
+ ggplot(aes(x = avg_cpm, y = log2FC)) +
+ geom_point(size = 0.5, alpha = 0.3, show.legend = FALSE) +
+ scale_x_log10() +
+ facet_wrap(~tissue)
2.2 Advanced Aesthetic Customization - color
+Color plays a crucial role in data visualization. We will explore how +to use ColorSpace to create plots tailored to different data types.
+Show all available color in colorspace.
+hcl_palettes(plot = TRUE)
hcl_palettes("qualitative", plot = TRUE)
hcl_palettes("sequential (single-hue)", n = 7, plot = TRUE)
hcl_palettes("sequential (multi-hue)", n = 7, plot = TRUE)
hcl_palettes("diverging", n = 7, plot = TRUE)
2.2.1 qualitative colors
+Qualitative (PCA of samples): principal component analysis (PCA) +summarizes the expression of many genes into a few new variables (PC1, +PC2, …) that capture the biggest differences between samples. Create a +scatter plot of PC1 vs. PC2, coloring points by tissue.
+Qualitative (iris dataset) : create a scatter plot of Sepal.Length +vs. Sepal.Width, coloring points by Species.
+# Example - RNA-seq data
+# PCA on the 500 most variable genes
+top500 <- names(sort(gene_var, decreasing = TRUE))[1:500]
+pca <- prcomp(t(expr_mat[top500, ]))
+
+pca_df <- as_tibble(pca$x[, 1:2], rownames = "sample") %>%
+ left_join(samples, by = "sample")
+
+# qualitative colors (Categorical Data)
+ggplot(pca_df, aes(x = PC1, y = PC2, color = tissue)) +
+ geom_point(size = 3) +
+ scale_color_discrete_qualitative(palette = "Dark 3") + # Using discrete qualitative color scale
+ ggtitle("Qualitative Colors Example: Samples by Tissue")
# Example: built-in iris data set
+ggplot(iris, aes(x = Sepal.Length, y = Sepal.Width, color = Species)) +
+ geom_point(size = 3) +
+ scale_color_discrete_qualitative(palette = "Dark 3") + # Using discrete qualitative color scale
+ ggtitle("Qualitative Colors Example: Iris Species")
2.2.2 Sequential(single-hue) colors
+Sequential single-hue (samples table): create a bar plot of samples +based on their library size (total reads, in millions), using a +single-hue sequential color gradient.
+Sequential single-hue (mtcars dataset): create a bar plot of car +models (rownames(mtcars)) based on their miles per gallon (mpg), using a +single-hue sequential color gradient.
+# sequential (single-hue) colors (Ordered Data - RNA-seq example)
+ggplot(samples, aes(x = reorder(sample, library_size), y = library_size, fill = library_size)) +
+ geom_bar(stat = "identity") +
+ scale_fill_continuous_sequential(palette = "Blues 3") + # Sequential color for fill
+ coord_flip() + # Flip for better readability
+ ylab("Library size (millions of reads)") +
+ xlab("Sample") +
+ ggtitle("Sequential (Single-Hue) Colors Example: Library Size per Sample") +
+ theme_minimal()
# sequential (single-hue) colors (Ordered Data)
+# Order data by mpg
+mtcars$car <- rownames(mtcars)
+mtcars <- mtcars[order(mtcars$mpg), ]
+# Create a bar plot
+ggplot(mtcars, aes(x = reorder(car, mpg), y = mpg, fill = mpg)) +
+ geom_bar(stat = "identity") +
+ scale_fill_continuous_sequential(palette = "Blues 3") + # Sequential color for fill
+ coord_flip() + # Flip for better readability
+ ggtitle("Sequential (Single-Hue) Colors Example: MPG in mtcars") +
+ theme_minimal()
2.2.3 Sequential (multi-hue) colors
+Sequential multi-hue (expression heatmap): the heatmap shows the +expression of the 30 most variable genes in every sample using color +intensity. It highlights which genes are highly expressed in which +samples.
+Sequential multi-hue (volcano dataset): The volcano heatmap provides +a 2D topographic view of elevation data using color intensity. It +highlights how terrain elevation changes over space.
+# Example: RNA-seq dataset
+# sequential (multi-hue) colors (Ordered Data)
+top30 <- names(sort(gene_var, decreasing = TRUE))[1:30]
+
+# Plot the heatmap
+expr %>%
+ filter(Gene %in% top30) %>%
+ ggplot(aes(x = sample, y = Gene, fill = log_cpm)) +
+ geom_tile() +
+ scale_fill_continuous_sequential(palette = "YlOrRd") + # Sequential multi-hue color scale
+ ggtitle("Sequential (Multi-Hue) Colors Example: Top 30 Variable Genes") +
+ xlab("Sample") +
+ ylab("Gene") +
+ theme_minimal() +
+ theme(axis.text.x = element_text(angle = 45, hjust = 1))
# sequential (multi-hue) colors (Ordered Data)
+# Convert volcano dataset to long format for ggplot
+volcano_df <- as.data.frame(as.table(volcano))
+colnames(volcano_df) <- c("x", "y", "elevation")
+# Plot the heatmap
+ggplot(volcano_df, aes(x = x, y = y, fill = elevation)) +
+ geom_tile() +
+ scale_fill_continuous_sequential(palette = "YlOrRd") + # Sequential multi-hue color scale
+ ggtitle("Sequential (Multi-Hue) Colors Example: Volcano Elevation") +
+ xlab("X Coordinate") + # Add x-axis label
+ ylab("Y Coordinate") + # Add y-axis label
+ theme_minimal() +
+ theme(
+ axis.text = element_blank(), # Remove axis text
+ axis.ticks = element_blank() # Remove axis ticks
+ )
2.2.4 Diverging colors
+Diverging (cornea fold changes): map the log2 fold change of the +genes that changed most after infection in cornea, using a diverging +color scale. The goal is to understand which genes go up and which go +down after infection. A positive value means the gene is higher in +infected samples, while a negative value means it is lower. Zero (no +change) is a meaningful middle point.
+Diverging (USArrests dataset): map UrbanPop deviation from the median +across states using a diverging color scale. The goal is to understand +how each state’s urban population percentage compares to the median +urban population percentage. A positive value means the state’s urban +population is above the median, while a negative value means it is below +the median.
+# diverging colors (cornea fold changes)
+# 30 genes with the largest change; avg_cpm >= 10 avoids noisy, barely expressed genes
+cornea_fc <- gene_summary %>%
+ filter(tissue == "cornea", avg_cpm >= 10) %>%
+ slice_max(abs(log2FC), n = 30)
+
+# Create a bar plot
+ggplot(cornea_fc, aes(x = reorder(Gene, log2FC), y = log2FC, fill = log2FC)) +
+ geom_bar(stat = "identity") +
+ scale_fill_continuous_divergingx(palette = "RdBu", rev = TRUE) + # rev: red = up, blue = down
+ coord_flip() +
+ xlab("Gene") +
+ ggtitle("Diverging Colors Example: Fold Change after Infection in Cornea") +
+ theme_minimal()
# diverging colors (USArrest Dataset)
+# Compute deviation from median
+USArrests$State <- rownames(USArrests)
+USArrests$UrbanPop_dev <- USArrests$UrbanPop - median(USArrests$UrbanPop)
+# Create a bar plot
+ggplot(USArrests, aes(x = reorder(State, UrbanPop_dev), y = UrbanPop_dev, fill = UrbanPop_dev)) +
+ geom_bar(stat = "identity") +
+ scale_fill_continuous_divergingx(palette = "RdBu") +
+ coord_flip() +
+ ggtitle("Diverging Colors Example: UrbanPop Deviation in USArrests") +
+ theme_minimal()
2.3 Example 2
+Compare the before-and-after heatmaps below to see how color choice +enhance readability and perceptual accuracy in data visualization. The +heatmap shows the mean expression of several inflammation-related genes +in each tissue and condition.
+2.3.1 Before
+Use the default color to create the heatmap.
+# before
+inflam_genes <- c("SOD2", "C3", "CXCL1", "CXCL6", "TNFAIP6", "CFB", "TNFAIP3")
+
+expr %>%
+ filter(Gene %in% inflam_genes) %>%
+ group_by(group, Gene) %>%
+ summarise(mean_log_cpm = mean(log_cpm), .groups = "drop") %>%
+ ggplot(aes(x = group, y = Gene, fill = mean_log_cpm)) +
+ geom_tile() +
+ theme_classic()
2.3.2 After
+Use the Sequential color to show the heatmap.
+# color choice 1 (Viridis)
+expr %>%
+ filter(Gene %in% inflam_genes) %>%
+ group_by(group, Gene) %>%
+ summarise(mean_log_cpm = mean(log_cpm), .groups = "drop") %>%
+ ggplot(aes(x = group, y = Gene, fill = mean_log_cpm)) +
+ geom_tile() +
+ theme_classic() +
+ scale_fill_continuous_sequential(palette = "Viridis", rev = FALSE)
# color choice 2 (Inferno)
+expr %>%
+ filter(Gene %in% inflam_genes) %>%
+ group_by(group, Gene) %>%
+ summarise(mean_log_cpm = mean(log_cpm), .groups = "drop") %>%
+ ggplot(aes(x = group, y = Gene, fill = mean_log_cpm)) +
+ geom_tile() +
+ theme_classic() +
+ scale_fill_continuous_sequential(palette = "Inferno", begin = 0.15, rev = FALSE)
2.4 Size & Transparency
+Explore how to adjust the size and transparency of the points to +enhance readability. Each point is one gene in one tissue, comparing its +mean expression in mock (x) and infected (y) samples. When you have this +many points (like about 40,000), size and transparency make a big +difference.
+ggplot(gene_summary) +
+ geom_point(aes(x = mock + 1,
+ y = CoV2 + 1,
+ color = tissue),
+ size = 2,
+ alpha = 0.4) +
+ scale_x_log10() +
+ scale_y_log10()
ggplot(gene_summary) +
+ geom_point(aes(x = mock + 1,
+ y = CoV2 + 1,
+ color = tissue),
+ size = 4,
+ alpha = 0.8) +
+ scale_x_log10() +
+ scale_y_log10()
2.5 Legend
+The legend plays a crucial role in interpreting the plot. Below, we +will explore how to adjust its appearance, create different legends for +different layers, and learn how to delete or merge legends +effectively.
+ggplot(gene_summary, aes(x = mock + 1, y = CoV2 + 1, color = tissue)) +
+ geom_point(size = 1, alpha = 0.4) +
+ scale_x_log10() +
+ scale_y_log10()
2.5.1 Adjust the appearance of the legend
+# Adjust the appearance of the legend
+ggplot(gene_summary, aes(x = mock + 1, y = CoV2 + 1, color = tissue)) +
+ geom_point(size = 1, alpha = 0.4) +
+ scale_x_log10() +
+ scale_y_log10() +
+ guides(color = guide_legend(override.aes = list(size = 3, alpha = 1)))
2.5.2 Generate different legends for different layers
+# Generate different legends for different layers
+ggplot(gene_summary, aes(x = mock + 1, y = CoV2 + 1, color = tissue)) +
+ geom_point(size = 1, alpha = 0.4) +
+ geom_smooth(method = "lm", se = FALSE) +
+ scale_x_log10() +
+ scale_y_log10()## `geom_smooth()` using formula = 'y ~ x'
ggplot(gene_summary, aes(x = mock + 1, y = CoV2 + 1, color = tissue)) +
+ geom_point(aes(alpha = "Observed"), size = 1) +
+ geom_smooth(method = "lm", se = FALSE, aes(alpha = "Fitted")) +
+ scale_x_log10() +
+ scale_y_log10() +
+ scale_alpha_manual(
+ name = NULL,
+ values = c(Observed = 0.4, Fitted = 1),
+ breaks = c("Observed", "Fitted")
+ ) +
+ guides(alpha = guide_legend(override.aes = list(
+ linetype = c(0, 1),
+ shape = c(16, NA),
+ color = "black"
+ )))## `geom_smooth()` using formula = 'y ~ x'
2.5.3 Delete legend
+# delete legend
+pca_df %>%
+ ggplot(aes(x = PC1, y = PC2, color = tissue, size = library_size)) +
+ geom_point()
pca_df %>%
+ ggplot(aes(x = PC1, y = PC2, color = tissue, size = library_size)) +
+ geom_point() +
+ guides(
+ color = guide_legend("tissue"), # keep
+ size = "none" # remove
+ )
2.5.4 Merge legend
+# merge legend
+pca_df %>%
+ ggplot(
+ aes(x = PC1, y = PC2,
+ color = library_size, size = library_size)
+ ) +
+ geom_point(alpha = 0.6) +
+ scale_color_viridis_c()
# method 1
+pca_df %>%
+ ggplot(
+ aes(x = PC1, y = PC2,
+ color = library_size, size = library_size)
+ ) +
+ geom_point(alpha = 0.6) +
+ scale_color_viridis_c() +
+ guides(
+ color = guide_legend(),
+ size = guide_legend()
+ )
# method 2
+pca_df %>%
+ ggplot(
+ aes(x = PC1, y = PC2,
+ color = library_size, size = library_size)
+ ) +
+ geom_point(alpha = 0.6) +
+ scale_color_viridis_c(guide = "legend")
2.6 Theme
+Explore how to change the theme of the plots.
+df <- data.frame(x = 1:3, y = 1:3)
+base <- ggplot(df, aes(x, y)) + geom_point()
+# default theme
+base + theme_grey() + ggtitle("theme_grey()")
# other built-in themes in ggplot2
+base + theme_bw() + ggtitle("theme_bw()")
base + theme_linedraw() + ggtitle("theme_linedraw()")
base + theme_light() + ggtitle("theme_light()")
base + theme_dark() + ggtitle("theme_dark()")
base + theme_minimal() + ggtitle("theme_minimal()")
base + theme_classic() + ggtitle("theme_classic()")
base + theme_void() + ggtitle("theme_void()")
2.7 Exporting plots
+Once a plot looks right, we usually want to save it as a file, for +example to finish a figure in a graphic editor like Inkscape. There are +two kinds of image formats:
+-
+
- Raster formats (PNG, JPEG, TIFF) store the image as +a grid of pixels. File size doesn’t depend on the number of points, but +text and lines can’t be edited and the image gets blurry when +enlarged. +
- Vector formats (SVG, PDF) store every point as a +separate shape. The text is sharp & editable but a plot with tens of +thousands of points becomes a very large file that is slow to open and +edit. +
Our mock vs. infected scatter plot has about 40,000 points, so it is +a good example of this problem.
+library(ggrastr)p_scatter <- ggplot(gene_summary, aes(x = mock + 1, y = CoV2 + 1, color = tissue)) +
+ geom_point(size = 0.5, alpha = 0.4) +
+ scale_x_log10() +
+ scale_y_log10() +
+ labs(x = "Mock (CPM + 1)", y = "SARS-CoV-2 (CPM + 1)") +
+ theme_classic()
+
+# raster: a fixed grid of pixels (dpi = dots per inch)
+ggsave("scatter_raster.png", p_scatter, width = 6, height = 5, dpi = 300)
+
+# vector: every one of the ~40,000 points is a separate shape
+ggsave("scatter_vector.pdf", p_scatter, width = 6, height = 5)2.7.1 Rasterize only the points with ggrastr
+The ggrastr package offers a middle ground: it rasterizes only the
+layers with many data points and keeps everything else (axes, text,
+legend) in vector format. Wrap the layer you want to rasterize in
+rasterise():
p_mixed <- ggplot(gene_summary, aes(x = mock + 1, y = CoV2 + 1, color = tissue)) +
+ rasterise(geom_point(size = 0.5, alpha = 0.4), dpi = 300) +
+ scale_x_log10() +
+ scale_y_log10() +
+ labs(x = "Mock (CPM + 1)", y = "SARS-CoV-2 (CPM + 1)") +
+ theme_classic()
+
+ggsave("scatter_mixed.pdf", p_mixed, width = 6, height = 5)
+
+# compare file sizes in MB
+file.size(c("scatter_raster.png", "scatter_vector.pdf", "scatter_mixed.pdf")) / 1e6## [1] 0.282269 2.092735 0.300479Alternatively, rasterize() can be applied to a finished
+plot, choosing layers by type:
p_mixed2 <- rasterize(p_scatter, layers = "Point", dpi = 300)When the mixed PDF is opened in Inkscape, the points appear as a +single image that is quick to move and resize, and the axis labels, tick +labels, and legend are still editable.
+Some tips:
+-
+
dpi = 300is a common minimum for print figures. Higher +resolution -> sharper image but larger file size
+- To save as SVG instead of PDF, use a file name ending in
+
.svg(this requires the svglite package); Inkscape opens +both formats
+ - Set
widthandheightin +ggsave()to the size the figure will have in the final +document (in inches by default), so that text sizes stay consistent and +you don’t need to rescale the figure later
+
2.8 Exercise 1
+Improve the plot based on the following requirements.
+-
+
- Improving the axes and legend labels +
- Adding a title for the plot +
- Tweaking the color scale to your favorite palette +
- The background should be white, not pale grey +
- Major gridlines should be a pale grey and minor gridlines should be +removed +
- The plot title should be centered in 12pt bold text +
ggplot(penguins, aes(x = bill_length_mm, y = bill_depth_mm, color = species)) +
+ geom_point()
The answer is as follows. Do not open the answer code block before +completing this exercise.
+# answer
+ggplot(penguins, aes(x = bill_length_mm, y = bill_depth_mm, color = species)) +
+ geom_point() +
+ labs(
+ title = "Penguin Bill Length vs. Depth", # Add title
+ x = "Bill Length (mm)", # Improve x-axis label
+ y = "Bill Depth (mm)", # Improve y-axis label
+ color = "Species" # Improve legend label
+ ) +
+ scale_color_discrete_qualitative(palette = "Dark 3") + # Choose a favorite color palette
+ theme_minimal(base_size = 12) + # Set white background and base size
+ theme(
+ legend.justification = c(0.5, 0.5), # Adjust legend positioning
+ legend.title = element_text(size = 10), # Set legend title font size
+ legend.text = element_text(size = 9), # Set legend text font size
+ plot.title = element_text(size = 12, face = "bold", hjust = 0.5), # Title centered in 12pt bold
+ panel.grid.major = element_line(color = "grey90"), # Pale grey major gridlines
+ panel.grid.minor = element_blank(), # Remove minor gridlines
+ panel.background = element_blank() # Set background to white
+ )
Module 3. Advanced Visualization Techniques
+Next, we will explore some specialized plots and learn which types of +data are best suited for each. To avoid version errors, if you’re using +a new version of R and want to make sure, re-install these:
+install.packages("plotly")
+install.packages("ggplot2")Check for these libraries:
+library(tidyverse)
+library(colorspace)
+library(pheatmap)
+library(plotly)
+library(ggplot2)
+library(gganimate)
+library(gifski)
+library(survival)
+library(survminer)
+library(qqman)3.1 Specialized plot - heatmaps
+# base R: heatmap() function
+data <- matrix(rnorm(100), nrow = 10)
+rownames(data) <- paste("Row", 1:10)
+colnames(data) <- paste("Col", 1:10)
+heatmap(data)
# ggplot2 with geom_tile()
+data <- expand.grid(Row = LETTERS[1:5], Column = LETTERS[6:10])
+data$Value <- runif(nrow(data), 1, 100)
+ggplot(data, aes(x = Column, y = Row, fill = Value)) +
+ geom_tile() +
+ scale_fill_continuous_sequential(palette = "Inferno", rev = FALSE) +
+ labs(title = "Heatmap using ggplot2", x = "Columns", y = "Rows", fill = "Value") +
+ theme_minimal()
# pheatmap for annotated heatmaps
+data <- matrix(rnorm(100), nrow = 10)
+rownames(data) <- paste("Gene", 1:10)
+colnames(data) <- paste("Sample", 1:10)
+pheatmap(data, color = colorRampPalette(c("blue", "white", "red"))(50))
3.2 Specialized plot - KM Curve
+# The 'lung' dataset is already loaded from the 'survival' library!
+
+# Create a survival object
+surv_object <- Surv(time = lung$time, event = lung$status)
+
+# Fit a Kaplan-Meier survival curve by group
+km_fit_group <- survfit(surv_object ~ sex, data = lung)
+
+# Plot the Kaplan-Meier curves for groups
+ggsurvplot(
+ km_fit_group,
+ data = lung,
+ conf.int = TRUE,
+ pval = TRUE, # Add p-value from log-rank test
+ xlab = "Time (days)",
+ ylab = "Survival Probability",
+ title = "Kaplan-Meier Curve by Sex",
+ palette = c("blue", "red"), # Different colors for groups
+ legend.labs = c("Male", "Female"), # Custom legend labels
+ ggtheme = theme_minimal()
+)
3.3 Specialized plot - Manhattan plot
+# Load the geuvadis dataset
+data("gwasResults")
+
+# Inspect the first few rows
+head(gwasResults)# Structure of the data
+# SNP: SNP identifier
+# CHR: Chromosome
+# BP: Base-pair position
+# P: P-value
+
+manhattan(
+ gwasResults,
+ chr = "CHR", # Chromosome column
+ bp = "BP", # Base-pair position column
+ p = "P", # P-value column
+ snp = "SNP", # SNP column
+ genomewideline = -log10(5e-8), # Genome-wide significance threshold
+ suggestiveline = -log10(1e-5), # Suggestive significance threshold
+ main = "Manhattan Plot", # Title
+ annotatePval = 0.01
+)
3.4 Specialized plot - 2d density plot
+# Scatter plot
+ggplot(diamonds, aes(x = carat, y = price)) +
+ geom_point() +
+ theme_minimal() +
+ labs(title = "Simple 2D Density Plot",
+ x = "Carat",
+ y = "Price")
# 2D density plot
+ggplot(diamonds, aes(x = carat, y = price)) +
+ geom_bin2d(bins = 100) +
+ scale_fill_continuous(type = "viridis") +
+ theme_minimal() +
+ labs(title = "Simple 2D Density Plot",
+ x = "Carat",
+ y = "Price")
3.5 Interactive plot - plotly
+Interactive plots allow users to explore data more deeply through +actions like clicking, zooming, or hovering, providing richer insights +and making the data more understandable and engaging. Next, we will +explore how to make plots interactive.
+3.5.1 ggplotly
+p <- ggplot(penguins) +
+ geom_point(aes(x = bill_length_mm,
+ y = bill_depth_mm,
+ color = species))
+
+ggplotly(p)3.5.2 Customize tooltip
+What if we are interested in which island the penguin lives or any +other information?
+penguins <- penguins %>%
+ mutate(text = paste(
+ "Island: ", island,
+ "\nSpecies: ", species,
+ "\nBill Length (mm): ", bill_length_mm,
+ "\nBill Depth (mm): ", bill_depth_mm,
+ "\nFlipper Length (mm): ",flipper_length_mm,
+ "\nBody Mass (g): ", body_mass_g, sep=""
+ ))
+p <- ggplot(penguins) +
+ geom_point(aes(x = bill_length_mm,
+ y = bill_depth_mm,
+ color = species,
+ text = text)) +
+ theme(legend.position = "none")
+ggplotly(p, tooltip="text")3.5.3 More advanced customization
+Adjust the tooltip based on the following requests.
+-
+
Use bold text to make the tooltip more readable
+Add a title to the tooltip
+Change the background color of the tooltip to white
+Change the font color of the tooltip to black
+
penguins <- penguins %>%
+ mutate(text = paste(
+ "<span style='font-size:16px;'><b>", species, "</b></span><br>",
+ "<br><b>Island:</b> ", island,
+ "<br><b>Bill Length (mm):</b> ", bill_length_mm,
+ "<br><b>Bill Depth (mm):</b> ", bill_depth_mm,
+ "<br><b>Flipper Length (mm):</b> ",flipper_length_mm,
+ "<br><b>Body Mass (g):</b> ", body_mass_g,
+ sep=""
+ ))
+p <- ggplot(penguins) +
+ geom_point(aes(x = bill_length_mm,
+ y = bill_depth_mm,
+ color = species,
+ text = text)) +
+ theme(legend.position = "none")
+pp <- ggplotly(p, tooltip="text")
+pp <- pp %>%
+ layout(
+ hoverlabel = list(bgcolor = "white",
+ font = list(size = 12, color = "black"))
+ )
+pp3.6 Animated plot - gganimate
+Using gganimate to create animated plots helps visualize changes over +time, making it easier to understand dynamic patterns and trends in +data. It adds an extra layer of clarity by showing how variables evolve, +which can be especially useful for time-series data.
+# Create the plot
+plot <- ggplot(gapminder, aes(gdpPercap, lifeExp, size = pop, colour = country)) +
+ geom_point(alpha = 0.7, show.legend = FALSE) +
+ scale_colour_manual(values = country_colors) +
+ scale_size(range = c(2, 12)) +
+ scale_x_log10() +
+ facet_wrap(~continent) +
+ # Here comes the gganimate specific bits
+ labs(title = 'Year: {frame_time}', x = 'GDP per capita', y = 'life expectancy') +
+ transition_time(year) +
+ ease_aes('linear')
+
+# Save at gif:
+anim <- animate(plot, width = 500, height = 357, units = "px", res = 96, renderer = gifski_renderer())
+anim_save("output.gif")
+anim
3.7 Shiny app
+Please open the “shiny app” folder.
+