Introduction to R

Tidyverse

Dr. Yu Cheng Hsu

Learning objectives

By the end of this session, you will

Wrangle (dplyr)

  • Read tabular data into R
  • Use core verbs: filter, mutate, summarise (with the pipe %>%)
  • Reflect on good coding practice in data wrangling

Plot (ggplot2)

  • Understand the basic grammar and layered structure of ggplot2
  • Create simple plots using data, aesthetics, and geometric layers

We DO NOT expect you to Memorize every function in ggplot2 or dplyr

Let’s collect some data

https://forms.cloud.microsoft/r/N29JZEAxFS

https://forms.cloud.microsoft/r/N29JZEAxFS

Package

  • A package is a bundle of functions, data, and documentation. You can import them through
library(package_name)
  • while using functions under the package you can use
package_name::function_name()

or

library(package_name)
function_name()
  • You can check what packages have been installed, imported to your environment using the tab Packages

Checking packages

Checking packages

Installing package

Golden rule: refer to the documentation of the specific package

  • CRAN (the Comprehensive R Archive Network) (Major package channel)
install.packages("ggplot2")
install.packages("tidyverse")
  • Bioconductor (Biomedical packages)
BiocManager::install("DESeq2")
  • GitHub (code at the developmental stage)
devtools::install_github("tidyverse/ggplot2")

Why wrangling matters

Raw experiment tables are almost never analysis-ready.

  • Columns may be wide, messy, or split across files
  • We need a clean table: one row per observation, clear columns for predictors/outcomes

dplyr is how we reshape and filter those tables so plots and models can run.

Keep a before/after picture in mind: messy spreadsheet → tidy analysis table.

Transforming data with dplyr

Using tidyverse an “umbrella-package” that installs several packages useful for data analysis which work together well such as tidyr, dplyr, ggplot2, tibble, etc.

# Title: Load tidyverse packages
# tidyverse loads multiple data science packages like ggplot2, dplyr, and tidyr.
library(tidyverse)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.2.0     ✔ readr     2.2.0
✔ forcats   1.0.1     ✔ stringr   1.6.0
✔ ggplot2   4.0.2     ✔ tibble    3.3.1
✔ lubridate 1.9.5     ✔ tidyr     1.3.2
✔ purrr     1.2.1     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

In this class, we actually use dplyr is a package for making tabular data manipulation easier. You can choose using

# Title: Load dplyr package
# dplyr provides functions for easy data frame manipulation.
library(dplyr)

Example data - Yeast gene expression

  • Rows = one gene measured in one nutrient-limitation condition
  • Columns = gene name, limiting nutrient, growth rate, expression
  • Numbers = log₂ expression relative to a reference RNA (how strongly that gene was “on”)

We will practice dplyr verbs on this table.

  • Source: Brauer et al. (2008), Coordination of Growth Rate, Cell Cycle, Stress Response, and Metabolic Activity in Yeast.
  • Goal: How does yeast gene expression change with growth rate when a single nutrient is limited?
  • Organism: budding yeast (Saccharomyces cerevisiae)

  • Exposure: rich media except one limiting nutrient — Glucose, Leucine, Ammonia, Sulfate, Phosphate, or Uracil

  • Observations: microarray expression at six growth rates (0.05 to 0.3 h⁻¹)

  • Running example: LEU1 (leucine biosynthesis) is turned up when leucine is starved, then falls as the restriction is relaxed.

  • brauer2007_tidy.csv: teaching subset of the expression table (one row per gene × nutrient × growth rate)

Reading data

# Load the yeast expression table and gene annotation from CSV files for later use.

ydat <- read_csv("https://bdsr.stephenturner.us/data/brauer2007_tidy.csv")
Rows: 198430 Columns: 7
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (5): symbol, systematic_name, nutrient, bp, mf
dbl (2): rate, expression

ℹ 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.

A reminder: Variable assignment

  • Only I use <- to assign or update the value of a variable, then the variable will be updated to the new value.
  • Othersies you just do a calculation and the variable will not be updated.
# Assign the value '50' to the weight object.
weight <- 50
# Print out weight to the screen (50)
weight
# What's the value of weight plus 10?
weight + 10
# Weight is still 50
weight
# Weight is only modified if we *reassign* weight to the modified value
weight <- weight+10
# Weight is now 60
weight
  • For rest of the lecture, I won’t update the ydat object, but instead will show the output of each operation. You can choose to update it if you want to keep the changes.

Inspecting the data

Getting the first 3 rows of the data frame by head()

# Show the first 3 rows of the expression data.
head(ydat, 3)
# A tibble: 3 × 7
  symbol systematic_name nutrient  rate expression bp                      mf   
  <chr>  <chr>           <chr>    <dbl>      <dbl> <chr>                   <chr>
1 SFB2   YNL049C         Glucose   0.05      -0.24 ER to Golgi transport   mole…
2 <NA>   YNL095C         Glucose   0.05       0.28 biological process unk… mole…
3 QRI7   YDL104C         Glucose   0.05      -0.02 proteolysis and peptid… meta…

Getting the summary of the data frame by summary()

# Show a summary of the expression data.
summary(ydat)
    symbol          systematic_name      nutrient              rate       
 Length:198430      Length:198430      Length:198430      Min.   :0.0500  
 Class :character   Class :character   Class :character   1st Qu.:0.1000  
 Mode  :character   Mode  :character   Mode  :character   Median :0.2000  
                                                          Mean   :0.1752  
                                                          3rd Qu.:0.2500  
                                                          Max.   :0.3000  
   expression             bp                 mf           
 Min.   :-6.500000   Length:198430      Length:198430     
 1st Qu.:-0.290000   Class :character   Class :character  
 Median : 0.000000   Mode  :character   Mode  :character  
 Mean   : 0.003367                                        
 3rd Qu.: 0.290000                                        
 Max.   : 6.640000                                        

Inspecting the data frame structure by str()

# Show the structure of the expression data.
str(ydat)
spc_tbl_ [198,430 × 7] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
 $ symbol         : chr [1:198430] "SFB2" NA "QRI7" "CFT2" ...
 $ systematic_name: chr [1:198430] "YNL049C" "YNL095C" "YDL104C" "YLR115W" ...
 $ nutrient       : chr [1:198430] "Glucose" "Glucose" "Glucose" "Glucose" ...
 $ rate           : num [1:198430] 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 ...
 $ expression     : num [1:198430] -0.24 0.28 -0.02 -0.33 0.05 -0.69 -0.55 -0.75 -0.24 -0.16 ...
 $ bp             : chr [1:198430] "ER to Golgi transport" "biological process unknown" "proteolysis and peptidolysis" "mRNA polyadenylylation*" ...
 $ mf             : chr [1:198430] "molecular function unknown" "molecular function unknown" "metalloendopeptidase activity" "RNA binding" ...
 - attr(*, "spec")=
  .. cols(
  ..   symbol = col_character(),
  ..   systematic_name = col_character(),
  ..   nutrient = col_character(),
  ..   rate = col_double(),
  ..   expression = col_double(),
  ..   bp = col_character(),
  ..   mf = col_character()
  .. )
 - attr(*, "problems")=<externalptr> 

Pipeline mindset

Without dplyr, if we want to manipulate data, we can still use the following vanilla way to perform the task :

# Title: Base R data manipulation example
# Show how row selection, counting, and column updates work without dplyr.

# Select rows 2 to 5 from the expression table.
ydat <- ydat[2:5,]

# Count how many Leucine-limited and Glucose-limited rows are in the selection.
sum(ydat$nutrient == 'Leucine')
sum(ydat$nutrient == 'Glucose')

# Overwrite the nutrient column with a logical value indicating leucine limitation.
ydat$nutrient = ydat$nutrient == 'Leucine'

Discussion question

  • What are the common features for these operations
  • What are the pitfalls in using this approach ?

Pipe operator |>

To prevent us from keep doing like a matryoshka doll

# Title: Nested function call example
# This shows the same work as a pipeline, but in nested form.
h(g(f(data_frame)))

R by default supports Pipe operator |> We can rewrite the above code as

# Pipe operator example
# Use the pipe operator to write the same 
# steps in a clearer, left-to-right order.
data_frame |> 
  f() |> 
  g() |> 
  h()
# This is the equivalent code of the left one.

data_frame <-f(data_frame)
data_frame <-g(data_frame)
data_frame <-h(data_frame)

‘|>’ takes the output from the left hand side and then pipes it to be the first argument of the right function.

Discussion question

  • What are the benefits of using pipe?

Select columns or rows

  • Select the designated columns by select(COLUMN_NAMES)
# Title: Select expression columns
ydat |>
    select(symbol, nutrient, rate, expression) |> # Keep gene, nutrient, growth rate, and expression.
    head(3) # This line is for display control, you can omit it and the previous |>.
# A tibble: 3 × 4
  symbol nutrient  rate expression
  <chr>  <chr>    <dbl>      <dbl>
1 SFB2   Glucose   0.05      -0.24
2 <NA>   Glucose   0.05       0.28
3 QRI7   Glucose   0.05      -0.02

Select by conditions

  • Select the designated rows by filter(CONDITION)
# Title: Filter rows by gene
ydat|>
    filter(symbol == 'LEU1') |> # Keep only the leucine biosynthesis gene LEU1.
    head(3) # This line is for display control, you can omit it and the previous |>.
# A tibble: 3 × 7
  symbol systematic_name nutrient  rate expression bp                   mf      
  <chr>  <chr>           <chr>    <dbl>      <dbl> <chr>                <chr>   
1 LEU1   YGL009C         Glucose   0.05      -1.12 leucine biosynthesis 3-isopr…
2 LEU1   YGL009C         Glucose   0.1       -0.77 leucine biosynthesis 3-isopr…
3 LEU1   YGL009C         Glucose   0.15      -0.67 leucine biosynthesis 3-isopr…

Common error

filter() function has different meaning in other R default library. If the error message occurs, you can try replaced it by dplyr::filter() instead

Manipulating columns

  • You can create new column by mutate()
# Title: Create or modify columns
ydat |>
    mutate( # expression is log2(sample/reference); 2^expression recovers the raw ratio.
              signal = 2^expression,
              restriction  = "Nutrient-limited"
              ) |>
    head(3)
# A tibble: 3 × 9
  symbol systematic_name nutrient  rate expression bp               mf    signal
  <chr>  <chr>           <chr>    <dbl>      <dbl> <chr>            <chr>  <dbl>
1 SFB2   YNL049C         Glucose   0.05      -0.24 ER to Golgi tra… mole…  0.847
2 <NA>   YNL095C         Glucose   0.05       0.28 biological proc… mole…  1.21 
3 QRI7   YDL104C         Glucose   0.05      -0.02 proteolysis and… meta…  0.986
# ℹ 1 more variable: restriction <chr>

Perform summary statistics

Frequently we will need to provide summary statistics for your data, we can use summarize() to do it

# Title: Compute summary statistics

ydat |>
    summarize( # Summarize expression with count, mean, standard deviation, and LEU1 row count.
              n=n(),
              mean_expr=mean(expression),
              sd_expr=sd(expression),
              n_leu1=sum(symbol=='LEU1',na.rm = TRUE)
    )
# A tibble: 1 × 4
       n mean_expr sd_expr n_leu1
   <int>     <dbl>   <dbl>  <int>
1 198430   0.00337   0.668     36

Some common summary function could be

  • mean()
  • sd()
  • max(), min()
  • sum()
  • n()

Summary statistics by group

  • group_by() and summarize() provide summary statistics by group.
# Title: Grouped summary by nutrient
ydat |>
# Group the expression table by limiting nutrient,
    group_by(nutrient) |> 
#  and compute summary statistics for each group.
    summarize(      
              n=n(),
              mean_expr=mean(expression)
                  )
# A tibble: 6 × 3
  nutrient      n mean_expr
  <chr>     <int>     <dbl>
1 Ammonia   33141 -0.000671
2 Glucose   33138  0.0223  
3 Leucine   33178  0.00577 
4 Phosphate 33068 -0.00207 
5 Sulfate   32897 -0.00381 
6 Uracil    33008 -0.00141 

Condition matching

if_else(CONDITION, 'TRUE_VALUE', 'FALSE_VALUE') for dichotomous processing

# Title: Create binary restriction status
# Use if_else() to add a restriction label based on growth rate.
ydat |>
    select(symbol, nutrient, rate, expression) |>
    mutate( restriction = if_else(
        rate == 0.05, "Severe restriction", "Relaxed"
    )) |>
    head(3)
# A tibble: 3 × 5
  symbol nutrient  rate expression restriction       
  <chr>  <chr>    <dbl>      <dbl> <chr>             
1 SFB2   Glucose   0.05      -0.24 Severe restriction
2 <NA>   Glucose   0.05       0.28 Severe restriction
3 QRI7   Glucose   0.05      -0.02 Severe restriction

Condition matching

  • case_when(Condition ~~ "Label") to handle sequential logical control
# Title: Create categorical restriction stage
# Use case_when() to assign descriptive labels for selected growth rates.
ydat |>
    select(symbol, nutrient, rate, expression) |>
    mutate( restriction = case_when(
        rate <= 0.05 ~ "Severe restriction",
        rate <= 0.15 ~ "Moderate restriction",
        rate <= 0.3 ~ "Relaxed restriction"
    ))|>
    head(3)
# A tibble: 3 × 5
  symbol nutrient  rate expression restriction       
  <chr>  <chr>    <dbl>      <dbl> <chr>             
1 SFB2   Glucose   0.05      -0.24 Severe restriction
2 <NA>   Glucose   0.05       0.28 Severe restriction
3 QRI7   Glucose   0.05      -0.02 Severe restriction

Best practice of data wrangling

  1. The “Read-Write-Verify” Loop

    • Write a step
    • Run it immediately
    • Verify the output

    DO NOT propagate error and raise your hand

  2. Look up documentations

  3. Make your code readable

  4. Trial and error

Advance topics

Renaming dataframe columns

  • You can rename column by rename()
# Title: Rename columns
ydat |> # Rename selected expression columns to make them more readable.
    rename(
              Gene = symbol,
              Nutrient = nutrient,
              Rate = rate
              ) |>
    head(3)
# A tibble: 3 × 7
  Gene  systematic_name Nutrient  Rate expression bp                       mf   
  <chr> <chr>           <chr>    <dbl>      <dbl> <chr>                    <chr>
1 SFB2  YNL049C         Glucose   0.05      -0.24 ER to Golgi transport    mole…
2 <NA>  YNL095C         Glucose   0.05       0.28 biological process unkn… mole…
3 QRI7  YDL104C         Glucose   0.05      -0.02 proteolysis and peptido… meta…

Groupby multiple conditions

You can group by multiple groups together, but note that they need to be stated in the same group_by() brackets

# Title: Multi-group summary
ydat |>
    group_by(nutrient, rate) |> # Group by both limiting nutrient and growth rate,
    summarize(                   #  then summarize the combined groups.
              n=n(),
              mean_expr=mean(expression),
              sd_expr=sd(expression),
    )

Checkpoint — three verbs for today

If you remember only three dplyr moves after the wrangling block, make them:

  1. filter() — keep rows that match a condition
  2. mutate() — add / change a column
  3. summarise() (+ group_by()) — collapse rows into a summary

Then we switch to ggplot2 and build plots layer by layer.

Introduction to ggplot2

About ggplot2

  • Visualization library in R for creating elegant and complex visualizations.
  • ggplot2 is part of the tidyverse collection of R packages.
  • Developed by Hadley Wickham [@wickham2010layered]. He received the COPSS Presidents’ Award for his contributions to the tidyverse collection.
  • Flexible and intuitive once you understand its core concepts.
library(tidyverse)
# or you can load ggplot2 directly
library(ggplot2)

About the data

Gapminder data is a dataset for the survey on each country with its life expectancy, GDP per capita, and population size. The data is collected by Gapminder Foundation.

  • country: a categorical variable 142 levels
  • continent :a categorical variable with 5 levels
  • year: going from 1952 to 2007 in increments of 5 years
  • pop : population
  • gdpPercap: GDP per capita
  • lifeExp: life expectancy

Read the data

gm <- read_csv("https://bdsr.stephenturner.us/data/gapminder.csv")
Rows: 1704 Columns: 6
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): country, continent
dbl (4): year, lifeExp, pop, gdpPercap

ℹ 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.

ggplot2 example

  • What are the elements in a graph ?

Components in the graph

Common components of a figure, figure from @ggplot2

Common components of a figure, figure from @ggplot2

Choosing the Visualization

Thinking process and decision flowchart will help you determine which type of graph to use (at least within the scope of this course):

  1. What is the purpose of displaying the graph?
  2. What types of data are you going to present?

flowchart LR
  A{Data type} -->|Continuous| B{Purpose}
  A{Data type}  -->|Discrete| C{Purpose}

  B{Purpose}-->|Exploration| D((Histogram/Boxplot))
  B{Purpose} -->|Association| E((Scatter plot))
  B{Purpose} -->|Association+time| T((Line plot))

  C{Purpose}-->|Exploration| F((Bar chart))
  C{Purpose} -->|Association| G((Tree map))

Overview of ggplot2 logics

ggplot2 is designed as if you are plotting by hand, so the process looks like

  1. Defining coordinate systems
  2. Plotting dots, lines, bars, and etc.
  3. Adding legend and axis information

Layer features of ggplot2

Layer features of ggplot2

Layer 1: Coordinate system

First, we define the coordinate system for our gm data in R

The function is ggplot(), which takes two arguments 1

  • data: The data you are going to use
  • mapping: Should be wrapped in aes() specifying how each column maps to the coordinate.
# Initialize a ggplot object with 
# data and aesthetic mappings
p <- ggplot(data=gm, 
  mapping = aes(x = gdpPercap))
p  # Display the base plot (no layers added yet)

# Alternative 
# p <- ggplot(gm, aes(gdpPercap))
# p
Figure 1: Showing coordinate systems in ggplot2

Layer 2: Plotting data

The general format for plotting items will be

geom_XXX(mapping=aes())

In aes(), you can map anything to the figure, common terms could be

  • x
  • y
  • color, or colour
  • fill, shape, etc.

These functions are named in the format geom_XXXX:

Function name
Histogram geom_histogram()
Box chart geom_boxplot()
Bar chart geom_bar()
Scatter chart geom_point()
Line chart geom_line()

Layer 2: Example

p+geom_point(aes(y = lifeExp))

Scatter plot for GDP per capita and life expectancy

Scatter plot for GDP per capita and life expectancy
p+geom_boxplot()

Box plot for life expectancy

Box plot for life expectancy
  • + operator to stack on the original ggplot()
  • Adding additional mapping in the geom_XXX()
  • Overwriting mapping in the ggplot()

More examples

ggplot(data=gm)+
  geom_histogram(aes(x=gdpPercap, color=continent))
`stat_bin()` using `bins = 30`. Pick better value `binwidth`.

Histogram for GDP per capita

Histogram for GDP per capita
ggplot(data=gm)+
  geom_bar(aes(x=continent,fill=continent))

Number of records in each continent

Number of records in each continent

Some notes 1: color affects the border color (dots/lines), fill affects the rectangle color

Layer 3: Legends and axis

Scales are functions indicate how each variable behave and presents on the graph. The function usually in the format scale_(AES)_(datatype):

Discrete Continuous
X scale_x_discrete() scale_x_continuous()
Y scale_y_discrete() scale_y_continuous()
color scale_color_discrete() scale_color_continuous()

Elements in the scale

Elements in the scale
Argument name Axis Legend
name Label Title
breaks Ticks Key
labels Tick label Key Label

Example

# Create a scatter plot with customized axis and color legend labels
   p +
  geom_point(aes(y = lifeExp, color=continent)) +  # Add points layer for data representation
  scale_x_continuous(name="GDP per capita") +  # Customize x-axis label
  scale_y_continuous(name="Life Expectancy",breaks=c(25,50,75)) + # Customize y-axis label
  scale_colour_discrete(name="Continent",labels=c("Afr", "Ameri", "Asia", "Eur", "Ocean"))  # Customize color legend title
Figure 2: Scatter plot of GDP per capita vs life expectancy with customized axis and legend labels.

Save figure

Add the function ggsave("FILE_PATH+FILE_NAME.png") after your ggplot instructions:

# Saving plots as figure_theme_grey.png at the folder figure
ggsave("figure_1.png",p)
  • It supports most of the common image formats like pdf, jpeg, png

Saved figure location

Saved figure location

Final notes

  • There are far more things you can do in ggplot2
  • Try to look up on the reference (or even use AI) when you need it.