R for Biologists course

R takes time to learn, like a spoken language. No one can expect to be an R expert after learning R for a few hours. This course has been designed to introduce biologists to R, showing some basics, and also some powerful things R can do (things that would be more difficult to do with Excel). The aim is to give beginners the confidence to continue learning R, so the focus here is on tidyverse and visualisation of biological data, as we believe this is a productive and engaging way to start learning R. After this short introduction you could use this book to dive a bit deeper.

Intro to R and RStudio

RStudio is an interface that makes it easier to use R. There are four windows in RStudio. The screenshot below shows an analogy linking the different RStudio windows to cooking.

R script vs console

There are two ways to work in RStudio in the console or in a script. We can type a command in the console and press Enter to run it. Try running the command below in the console.

1 + 1

## [1] 2

Or we can use an R script. To create a script, from the top menu in RStudio: File > New File > R Script. Now type the command below in the script. This time, to run the command, you use Ctrl + Enter for Windows/Linux or Cmd + Enter for MacOS. This sends the command where the cursor is from the script to the console. You can highlight multiple commands and then press Cmd/Ctrl + Enter to run them one after the other.

2 + 2

## [1] 4

As the RStudio screenshot above explains, if we work in the console we don’t have a good record (recipe) of what we’ve done. We can see commands we’ve run in the History panel (top right window), and we can go backwards and forwards through our history in the console using the up arrow and down arrow. But the history includes everything we’ve tried to run, including our mistakes so it is good practice to use an R script.

We can also add comments to a script. These are notes to ourself or others about the commands in the script. Comments start with a # which tells R not to run them as commands.

# testing R
2 + 2

## [1] 4

Keeping an accurate record of how you’ve manipulated your data is important for reproducible research. Writing detailed comments and documenting your work are useful reminders to your future self (and anyone else reading your scripts) on what your code does.

Working directory

Opening an RStudio session launches it from a specific location. This is the ‘working directory’. R looks in the working directory by default to read in data and save files. You can find out what the working directory is by using the command getwd(). This shows you the path to your working directory in the console. In Mac this is in the format /path/to/working/directory and in Windows C:\path\to\working\directory. It is often useful to have your data and R scripts in the same directory and set this as your working directory. We will do this now.

Make a folder for this course somewhere on your computer that you will be able to easily find. Name the folder for example, Intro_R_course. Then, to set this folder as your working directory:

In RStudio click on the ‘Files’ tab and then click on the three dots, as shown below.

In the window that appears, find the folder you created (e.g. Intro_R_course), click on it, then click ‘Open’. The files tab will now show the contents of your new folder. Click on More > Set As Working Directory, as shown below.

Save the script you created in the previous section as intro.R in this directory. You can do this by clicking on File > Save and the default location should be the current working directory (e.g. Intro_R_course).

Note: You can use an RStudio project as described here to automatically keep track of and set the working directory.

Packages

If it’s not already installed on your computer, you can use the install.packages function to install a package. A package is a collection of functions along with documentation, code, tests and example data.

install.packages("tidyverse")

We will see many functions in this tutorial. Functions are “canned scripts” that automate more complicated sets of commands. Many functions are predefined, or can be made available by importing R packages. A function usually takes one or more inputs called arguments. Here tidyverse is the argument to the install.packages() function.

Note: functions require parentheses after the function name.

Getting help

To see what any function in R does, type a ? before the name and help information will appear in the Help panel on the right in RStudio. Or you can search the function name in the Help panel search box. Google and Stack Overflow are also useful resources for getting help.

?install.packages

ⓘ Tab completion
A very useful feature is Tab completion. You can start typing and use Tab to autocomplete code, for example, a function name.

Common R errors

R error messages are common and can sometimes be cryptic. You most likely will encounter at least one error message during this tutorial. Some common reasons for errors are:

Case sensitivity. In R, as in other programming languages, case sensitivity is important. ?install.packages is different to ?Install.packages.
Missing commas
Mismatched parentheses or brackets
Not quoting file paths
Not finishing a command so seeing “+” in the console. If you need to, you can press ESC to cancel the command.

To see examples of some R error messages with explanations see here

Getting started with data

Data files

The data files required for this workshop are available on GitHub. To download the data.zip file, you can click here. Unzip the file and store this data folder in your working directory.

GREIN (GEO RNA-seq Experiments Interactive Navigator)

In this tutorial, we will learn some R through creating plots to visualise data from an RNA-seq experiment. RNA-seq counts file can be obtained from the GREIN platform. GREIN provides >6,500 published datasets from GEO that have been uniformly processed. It is available at http://www.ilincs.org/apps/grein/. You can search for a dataset of interest using the GEO code. We obtained the dataset used here using the code GSE60450. GREIN provide QC metrics for the RNA-seq datasets and both raw and normalized counts. We will use the normalized counts here. These are the counts of reads for each gene for each sample normalized for differences in sequencing depth and composition bias. Generally, the higher the number of counts the more the gene is expressed.

RNA-seq dataset

Here we will create some plots using RNA-seq data from the paper by Fu et al. 2015, GEO code GSE60450. This study examined expression in basal and luminal cells from mice at different stages (virgin, pregnant and lactating). There are 2 samples per group and 6 groups, 12 samples in total.

Tidyverse

In this course we will use the tidyverse. The tidyverse is a collection of R packages that includes the extremely widely used ggplot2 package.

The tidyverse makes data science faster, easier and more fun.

Why tidyverse? Why tidy data? Why is it such a game-changer?

Loading the data

We use library() to load in the packages that we need. As described in the cooking analogy in the first screenshot, install.packages() is like buying a saucepan, library() is taking it out of the cupboard to use it.

library(tidyverse)

The files we will use are csv comma-separated, so we will use the read_csv() function from the tidyverse. There is also a read_tsv() function for tab-separated values.

We will use the counts file called GSE60450_GeneLevel_Normalized(CPM.and.TMM)_data.csv that’s in a folder called data i.e. the path to the file should be data/GSE60450_GeneLevel_Normalized(CPM.and.TMM)_data.csv.

We can read the counts file into R with the command below. We’ll store the contents of the counts file in an object called counts. This stores the file contents in R’s memory making it easier to use.

# read in counts file
counts <- read_csv("data/GSE60450_GeneLevel_Normalized(CPM.and.TMM)_data.csv")

## New names:
## Rows: 23735 Columns: 14
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (2): ...1, gene_symbol dbl (12): GSM1480291, GSM1480292, GSM1480293,
## GSM1480294, GSM1480295, GSM148...
## ℹ 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.
## • `` -> `...1`

# read in metadata
sampleinfo <- read_csv("data/GSE60450_filtered_metadata.csv")

## New names:
## Rows: 12 Columns: 4
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (4): ...1, characteristics, immunophenotype, developmental stage
## ℹ 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.
## • `` -> `...1`

There is some information output by read_csv on “column specification”. It tells us that there is a missing column name in the header and it has been filled with the name “…1”, which is how read_csv handles missing column names by default. It also tells us what data types read_csv is detecting in each column. Columns with text characters have been detected (col_character) and also columns with numbers (col_double). We won’t get into the details of R data types in this tutorial but they are important to know and you can read more about them in the R for Data Science book.

In R we use <- to assign values to objects. <- is the assignment operator. It assigns values on the right to objects on the left. So to create an object, we need to give it a name (e.g. counts), followed by the assignment operator <-, and the value we want to give it. We can give an object almost any name we want but there are some rules and conventions as described in the tidyverse R style guide

We can read in a file from a path on our computer on on the web and use this as the value. Note that we need to put quotes (““) around file paths.

ⓘ Assignment operator shortcut
In RStudio, typing Alt + - (holding down Alt at the same time as the - key) will write <- in a single keystroke in Windows, while typing > Option + - (holding down Option at the same time as the - key) does the same in a Mac.

Exercise

Test what happens if you type Library(tidyverse)
What is wrong and how would you fix it?
Test what happens if you type libary(tidyverse)
What is wrong and how would you fix it?
Test what happens if you type library(tidyverse
What is wrong and how would you fix it?
Test what happens if you type
read_tsv("data/GSE60450_filtered_metadata.csv")
What is wrong and how would you fix it?
Test what happens if you type
read_csv("data/GSE60450_filtered_metadata.csv)
What is wrong and how would you fix it?
Test what happens if you type
read_csv("GSE60450_filtered_metadata.csv)
What is wrong and how would you fix it?
What is the name of the first column you get with each of these 2 commands?
read.csv("data/GSE60450_filtered_metadata.csv")
and
read_csv("data/GSE60450_filtered_metadata.csv")
If you run
read_csv("data/GSE60450_filtered_metadata.csv")
what is the difference between the column header you see developmental stage and ‘developmental stage’?

Getting to know the data

When assigning a value to an object, R does not print the value. For example, here we don’t see what’s in the counts or sampleinfo files. But there are ways we can look at the data. We will demonstrate using the sampleinfo object.

We can type the name of the object and this will print the first few lines and some information, such as number of rows.

sampleinfo

We can also use dim() to see the dimensions of an object, the number of rows and columns.

dim(sampleinfo)

## [1] 12  4

This show us there are 12 rows and 4 columns.

In the Environment Tab in the top right panel in RStudio we can also see the number of rows and columns in the objects we have in our session.

We can also take a look the first few lines with head(). This shows us the first 6 lines.

head(sampleinfo)

We can look at the last few lines with tail(). This shows us the last 6 lines. This can be useful to check the bottom of the file, that it looks ok.

tail(sampleinfo)

Or we can see the whole file with View().

View(sampleinfo)

In the Environment tab we can see how many rows and columns the object contains and we can click on the icon to view all the contents in a tab. This runs the command View() for us.

We can see all the column names with colnames().

colnames(sampleinfo)

## [1] "...1"                "characteristics"     "immunophenotype"    
## [4] "developmental stage"

We can access individual columns by name using the $ symbol. For example we can see what’s contained in the characteristics column.

sampleinfo$characteristics

##  [1] "mammary gland, luminal cells, virgin"            
##  [2] "mammary gland, luminal cells, virgin"            
##  [3] "mammary gland, luminal cells, 18.5 day pregnancy"
##  [4] "mammary gland, luminal cells, 18.5 day pregnancy"
##  [5] "mammary gland, luminal cells, 2 day lactation"   
##  [6] "mammary gland, luminal cells, 2 day lactation"   
##  [7] "mammary gland, basal cells, virgin"              
##  [8] "mammary gland, basal cells, virgin"              
##  [9] "mammary gland, basal cells, 18.5 day pregnancy"  
## [10] "mammary gland, basal cells, 18.5 day pregnancy"  
## [11] "mammary gland, basal cells, 2 day lactation"     
## [12] "mammary gland, basal cells, 2 day lactation"

If we just wanted to see the first 3 values in the column we can specify this using square brackets. Obtaining a selection of values this way is called ‘subsetting’.

sampleinfo$characteristics[1:3]

## [1] "mammary gland, luminal cells, virgin"            
## [2] "mammary gland, luminal cells, virgin"            
## [3] "mammary gland, luminal cells, 18.5 day pregnancy"

In the previous section, when we loaded in the data from the csv file, we noticed that the first column had a missing column name and by default, read_csv function assigned a name of “...1” to it. Let’s change this column to something more descriptive now. We can do this by combining a few things we’ve just learnt.

First, we use the colnames() function to obtain the column names of sampleinfo. Then we use square brackets to subset the first value of the column names ([1]). Last, we use the assignment operator (<-) to set the new value of the first column name to “sample_id”.

colnames(sampleinfo)[1] <- "sample_id"

Let’s check if this has been changed correctly.

sampleinfo

The first column is now named “sample_id”.

We can also do the same to the counts data. This time, we rename the first column name from “...1” to “gene_id”.

colnames(counts)[1] <- "gene_id"

ⓘ Multiple methods
There are multiple ways to rename columns. We’ve covered one way here, but another way is using the rename() function. When programming, you’ll often find many ways to do the same thing. Often there is one obvious method depending on the context you’re in.

Other useful commands for checking data are str() and summary().

str() shows us the structure of our data. It shows us what columns there are, the first few entries, and what data type they are e.g. character or numbers (double or integer).

str(sampleinfo)

## spc_tbl_ [12 × 4] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##  $ sample_id          : chr [1:12] "GSM1480291" "GSM1480292" "GSM1480293" "GSM1480294" ...
##  $ characteristics    : chr [1:12] "mammary gland, luminal cells, virgin" "mammary gland, luminal cells, virgin" "mammary gland, luminal cells, 18.5 day pregnancy" "mammary gland, luminal cells, 18.5 day pregnancy" ...
##  $ immunophenotype    : chr [1:12] "luminal cell population" "luminal cell population" "luminal cell population" "luminal cell population" ...
##  $ developmental stage: chr [1:12] "virgin" "virgin" "18.5 day pregnancy" "18.5 day pregnancy" ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   ...1 = col_character(),
##   ..   characteristics = col_character(),
##   ..   immunophenotype = col_character(),
##   ..   `developmental stage` = col_character()
##   .. )
##  - attr(*, "problems")=<externalptr>

summary() generates summary statistics of our data. For numeric columns (columns of type double or integer) it outputs statistics such as the min, max, mean and median. We will demonstrate this with the counts file as it contains numeric data. For character columns it shows us the length (how many rows).

summary(counts)

##    gene_id          gene_symbol          GSM1480291          GSM1480292       
##  Length:23735       Length:23735       Min.   :    0.000   Min.   :    0.000  
##  Class :character   Class :character   1st Qu.:    0.000   1st Qu.:    0.000  
##  Mode  :character   Mode  :character   Median :    1.745   Median :    1.891  
##                                        Mean   :   42.132   Mean   :   42.132  
##                                        3rd Qu.:   29.840   3rd Qu.:   29.604  
##                                        Max.   :12525.066   Max.   :12416.211  
##    GSM1480293         GSM1480294         GSM1480295          GSM1480296       
##  Min.   :    0.00   Min.   :    0.00   Min.   :     0.00   Min.   :     0.00  
##  1st Qu.:    0.00   1st Qu.:    0.00   1st Qu.:     0.00   1st Qu.:     0.00  
##  Median :    0.92   Median :    0.89   Median :     0.58   Median :     0.54  
##  Mean   :   42.13   Mean   :   42.13   Mean   :    42.13   Mean   :    42.13  
##  3rd Qu.:   21.91   3rd Qu.:   19.92   3rd Qu.:    12.27   3rd Qu.:    12.28  
##  Max.   :49191.15   Max.   :55692.09   Max.   :111850.87   Max.   :108726.08  
##    GSM1480297          GSM1480298          GSM1480299       
##  Min.   :    0.000   Min.   :    0.000   Min.   :    0.000  
##  1st Qu.:    0.000   1st Qu.:    0.000   1st Qu.:    0.000  
##  Median :    2.158   Median :    2.254   Median :    1.854  
##  Mean   :   42.132   Mean   :   42.132   Mean   :   42.132  
##  3rd Qu.:   27.414   3rd Qu.:   26.450   3rd Qu.:   24.860  
##  Max.   :10489.311   Max.   :10662.486   Max.   :15194.048  
##    GSM1480300          GSM1480301          GSM1480302       
##  Min.   :    0.000   Min.   :    0.000   Min.   :    0.000  
##  1st Qu.:    0.000   1st Qu.:    0.000   1st Qu.:    0.000  
##  Median :    1.816   Median :    1.629   Median :    1.749  
##  Mean   :   42.132   Mean   :   42.132   Mean   :   42.132  
##  3rd Qu.:   23.443   3rd Qu.:   23.443   3rd Qu.:   24.818  
##  Max.   :17434.935   Max.   :19152.728   Max.   :15997.193

Formatting the data

Converting from wide to long format

We will first convert the data from wide format into long format to make it easier to work with and plot with ggplot. We want just one column containing all the expression values instead of multiple columns with counts for each sample, as shown in the image below.

We can use pivot_longer() to easily change the format into long format.

seqdata <- pivot_longer(counts, cols = starts_with("GSM"), names_to = "Sample", 
                        values_to = "Count")

We use cols = starts_with("GSM") to tell the function we want to reformat the columns whose names start with “GSM”. pivot_longer() will then reformat the specified columns into two new columns, which we’re naming “Sample” and “Count”. The names_to = "Sample" specifies that we want the new column containing the columns we sepcified with cols to be named “Sample”, and the values_to = "Count" specifies that we want the new column contining the values to be named “Count”.

We could also specify a column range to reformat. The command below would give us the same result as the previous command.

seqdata <- pivot_longer(counts, cols = GSM1480291:GSM1480302, 
                        names_to = "Sample", values_to = "Count")

Alternatively, we could specify the columns we don’t want to reformat and pivot_longer() will reformat all the other columns. To do that we put a minus sign “-” in front of the column names that we don’t want to reformat. This is a pretty common way to use pivot_longer() as sometimes it is easier to exclude columns we don’t want than include columns we do. The command below would give us the same result as the previous command.

seqdata <- pivot_longer(counts, cols = -c("gene_id", "gene_symbol"), 
                        names_to = "Sample", values_to = "Count")

Here we see the function c() for the first time. We use this function extremely often in R when we have multiple items that we are combining. We will see it again in this tutorial.

Let’s have a look at the data.

seqdata

Joining two tables

Now that we’ve got just one column containing sample ids in both our counts and metadata objects we can join them together using the sample ids. This will make it easier to identify the categories for each sample (e.g. if it’s basal cell type) and to use that information in our plots.

We will use the function full_join() and give it the two tables we want to join. We add by = join_by(Sample == sample_id) to say we want to join on the column called “Sample” in the first table (seqdata) and the column called “sample_id” in the second table (sampleinfo) when the values match.

allinfo <- full_join(seqdata, sampleinfo, by = join_by(Sample == sample_id))