In rnabioco/eda: Practical Biological Data Analysis in R/RStudio

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.align = "center"
)

Goals for the course

• Master the basics of data import, processing, plotting, and report generation in RStudio.

• Get hooked on RStudio instead of MS Excel, Prism, etc as your go-to analysis software.

• Learn productive analysis habits.

• Understand the analysis ecosystem so that new types of data are easily approachable

RStudio Overview

The class has its own R package, so you should begin by loading it into your session. You will need to reinstall the package when we make updates.

# devtools::install_github('rnabioco/practical-data-analysis')
library(pbda)

Environment

The Environment panel enables data import and environment browsing.

There are two main ways to import data within RStudio. One is interactive, using the "Import Dataset" button in the Environment panel.

read functions

Alternatively you can use functions in readr and readxl to load data on the console or within an RMarkdown report. All of the functions have read_* prefixes.

library(readr)
read_csv(readr_example("mtcars.csv"))

library(readxl)
read_excel(readxl_example("datasets.xlsx"))

You can also explore your environment with the "Global Environment" button.

Console

The Console panel enables interactive analysis.

library(tidyverse)

This one invocation loads several libraries including dplyr and ggplot2.

Viewer

The tidyverse contains ggplot2, which you covered in your pre-requisites. Let's make a simple plot in RStudio. Note that the plot appears in the "Plots" panel. From there, you can Zoom, Export to a file, or Publish on RPubs.

library(ggplot2)
ggplot(mtcars, aes(x = hp, y = mpg, color = factor(cyl))) + geom_point()

Editor

Finally we'll create an RMarkdown document in the Editor panel. Select File > New File > RMarkdown... and fill the Title and Author fields. Then click OK to create an RMarkdown document. When you have completed editing the document, click the "Knit" button to render the output (e.g., as an HTML page).

Rmarkdown is customizable with tables of contents, tabbed views and more. See the Rmarkdown documentation for examples.

Projects

RStudio provides a "Project" concept to keep related analyses together. I suggest starting each class by making your own project (File > New Project).

Tidyverse cheatsheets

See Help > Cheatsheets for helpful references. Print them out and pin them up next to your screen. The dplyr and ggplot2 cheatsheets are especially useful for beginners.

Git and Github

Github provides several advantages:

• Revision control. You have a history of all changes to your documents and can go back to any point.
• Collaboration. It is easy to share analyses with your colleagues and is built to enable multiple people to work an a common project.
• Publicity. Your github page reflects your interests and productivity. Others might visit to evaluate your code samples.

Jenny Bryan's "Happy Git with R" tutorial gets you up and running by integrating Git and R Studio.

We'll be using github to submit problem sets, so make sure this is set up by class time on Thursday.

Formatting data for efficient analysis

Reformatting your data into "tidy" format is a useful practice that allows you efficiently manipulate and analyze data. We will start with some toy examples and then will reformat a commonly used data type from a quantitative PCR experiment.

Tidy data adopts three simple conventions:

1. Columns are variables.
2. Rows are observations.
3. Store related data in separate tables.

The tibble

A tibble is similar to a data.frame with some important differences. A tibble provides a quick overview of your data including its structure (rows and columns) and column types. Use as_tibble() to convert a data.frame:

head(iris)

library(tidyverse)
as_tibble(iris)

The tidyr package

The tidyr package has useful functions to format data for efficient analysis. See the cheat sheet for reference.

gather and spread

These are the workhorse functions in tidyr. Use gather when variables are in columns and you want them in tidy format.

library(pbda)
# column names in "wide" format are variables
mtcars_tbl

# gather variables and observations into `key` and `value` columns - so-called "long" format
xx <- mtcars_tbl %>% gather(key, value, -name)
xx

# and convert back to "wide" format
xx %>% spread(key, value)

Long format

Each form has its advantages. For example, long format makes it easy to calculate summaries:

xx %>% group_by(key) %>% summarize(median(value))

Wide format

And wide format is useful for plotting.

library(ggrepel)

set.seed(42)
ggplot(mtcars_tbl) +
  geom_point(aes(hp, mpg, size = 5)) +
  geom_label_repel(aes(hp, mpg, fill = factor(cyl), label = name))

Exercise

Examine the beavers data set (beaver1 and beaver2) to identify its variables. Is it in tidy format? If not, how would you rearrange it to make it tidy?

separate and unite

separate takes a column containing multiple variables on input and returns multiple columns, each with a new variable. For example, a column with year/month/day information can be separated into invidual columns.

ys <- 1999:2002
ms <- c('Jan', 'Feb', 'Mar')
ds <- 1:10

dates <- tidyr::crossing(ys, ms, ds) %>% unite(date, ys:ds, sep = '-')

separate

dates

# separate is the inverse of unite
dates %>% separate(date, into = c('year', 'month', 'day'), sep = '-')

The sep argument can take:

• a character (split rep_value using sep = '_' into rep and value)
• a position (split a1 using sep = 1 into a and 1)

Finally the extra and fill arguments to separate control what happens when there are too many and not enough variables.

crossing and complete

crossing is useful for generating combinations of variables in tibble format. For example, use crossing to generate combinations of experimental varaibles including sample names, gene names, reaction conditions, and replicates.

genotype <- c('wt', 'mut')
gene <- c('IFN', 'ACTIN')
time <- c(0, 12, 24, 48)
rt <- c('+', '-') # reverse transcriptase added?
rep <- 1:3

samples <- tidyr::crossing(genotype, gene, time, rep, rt)

samples

Data in the 96-well plate format.

Now we'll use tidy data principles to analyze some qPCR data.

Many biological assays make use of the 96 (or 384) well plate. Note the similarity between the plate and a tibble: there are rows and columns, and each well contains a reaction that will generate one or more data points.

Sample names

All variables should be systematically listed in your sample names, i.e. name_rep_time_RT. Systematic naming makes it easy to extract relevant information.

Take this example, where the sample names are a combination of a genotype (WT and MT), a time point (0,4,8,24 hour), and a replicate (1,2,3), separated by a hyphen.

library(tidyverse)

# for reproducible `sample`
set.seed(47681)

samples <-
  tidyr::crossing(
    name = c('WT', 'MT'),
    hours = c('t0', 't4', 't8', 't24'),
    reps = 1:3
    ) %>%
  mutate(
    value = sample(1:100, n(), replace = TRUE),
    .id = row_number()
    ) %>%
  unite('sample.name', name, hours, reps, sep = '-') %>%
  select(-.id)

samples

Extracting sample names

Because the samples have systematic names, it is easy to separate this information into individual columns.

sample_info <- samples %>%
  tidyr::separate(
    sample.name,
    into = c('sample', 'hour', 'rep'),
    sep = "-"
  )

sample_info

Data manipulation

Now we can use dplyr and tidyr functions to manipulate the data.

# calculate summary statistics
sample_info %>% group_by(sample, hour) %>% summarize(mean(value))

# subtract a background value. N.B.: rearranging the table makes this calculation easy.
sample_info %>% spread(hour, value) %>% mutate(t24_norm = t24 - t0)

qPCR data

The class library provides two related tibbles that describe a simulated qPCR experiment called qpcr_names and qpcr_data.

library(pbda)

qpcr_names

qpcr_data

We will use tidying concepts to prepare this data for efficient analysis and visualization.

qPCR data tidying

• Tidy qpcr_data and qpcr_names into a structure like:

qpcr_names_tidy <- qpcr_names %>% gather(col, value, -row)
# `exp` is the relative expression level
qpcr_data_tidy <- qpcr_data %>% gather(col, exp, -row)

qpcr_data_tidy

Sample names

• Separate variables into new columns in qpcr_names_tidy.

qpcr_names_tidy <- separate(
  qpcr_names_tidy,
  value,
  into = c('sample', 'time', 'gene', 'rt', 'rep'),
  sep = '_'
)

qpcr_names_tidy

Data joining

• Join the tidied data together.

qpcr_tidy <- left_join(qpcr_names_tidy, qpcr_data_tidy)

qpcr_tidy

Statistical summary

• Calculate summary statistics for each gene, cell and time point across replicates.

qpcr_tidy %>%
  filter(rt == "+") %>%
  group_by(sample, gene, time) %>%
  summarize(mean_exp = mean(exp), var_exp = var(exp))

Plots

• Plot the expression for each gene over time.

• Calculate a fold-change for IFN over ACTIN and re-plot.

ggplot(qpcr_tidy, aes(x = time, y = exp, color = rt)) + geom_point(size=3) + facet_wrap(~gene)

Exercises

1. Tidy the mtcars data set. Note that car names are rownames in the built-in data, so they need to be moved to their own column prior to tibble conversion. You can use tibble:rownames_to_column() for this purpose.

The tidied data should look something like:

mtcars %>%
  mutate(name = rownames(.)) %>%
  as_tibble() %>%
  gather(var, value, -name)

1. For each car in the tidy mtcars data set, calculate the mean (mean()) and variance (var()) for each variable. Try using summarize(), and summarize_at() or summarize_all().

2. Plot mpg vs cyl for the mtcars data set. Which format should you use? The original data set, or the tidied one? Why?

3. Using the provided qpcr data, plot the changes in gene expression over time. Use colors to represent genotypes and facets for the different genes. If that's too easy, add error bars (geom_errorbar()) and connect each point with a line (geom_line()).

Too Easy?

If you found those exercises too easy, try these:

1. Tidy the anscombe data set. Your tidy data should look like:

# A tibble: 44 x 3
   set       x     y
   <chr> <dbl> <dbl>
 1 1         4  4.26
 2 1         5  5.68
 3 1         6  7.24
 4 1         7  4.82
 5 1         8  6.95
 6 1         9  8.81
 7 1        10  8.04
 8 1        11  8.33
 9 1        12 10.8 
10 1        13  7.58
# ... with 34 more rows

1. Calculate summary statistics in your tidy anscombe data. Calculate n, mean, and var for each x and y variable in the sets. Try using summarize_each for this calculation.

Calculate the cor between x and y for each set. What do you infer from these summary statistics?

1. Plot the tidy anscombe data. Plot each x/y point on separate facets for each set, and plot a smoothing line. Examine the plots. How do they compare to the summary statistics?

ggplot(anscombe_tidy, aes(x, y)) +
  geom_point() +
  facet_wrap(~set) +
  geom_smooth(method = 'lm', se = FALSE)

rnabioco/eda documentation built on July 12, 2022, 2:17 p.m.