R Crash Course

Introduction to R - Part III

Recap

library(tidyverse)
gapminder <- read_csv("raw_data/gapminder.csv")

Summarising and grouping with dplyr

The summarise function can take any R function that takes a vector of values (i.e. a column from a data frame) and returns a single value. Some of the more useful functions include:

• min minimum value
• max maximum value
• sum sum of values
• mean mean value
• sd standard deviation
• median median value
• IQR the interquartile range
• n_distinct the number of distinct values
• n the number of observations (Note: this is a special function that doesn’t take a vector argument, i.e. column)

summarise(gapminder, min(lifeExp), max(gdpPercap), mean(pop))

It is also possible to summarise using a function that takes more than one value, i.e. from multiple columns. For example, we could compute the correlation between year and life expectancy. Here we also assign names to the table that is produced.

gapminder %>% 
summarise(MinLifeExpectancy = min(lifeExp), 
          MaximumGDP = max(gdpPercap), 
          AveragePop = mean(pop), 
          Correlation = cor(year, lifeExp))

However, it is not particularly useful to calculate such values from the entire table as we have different continents and years. The group_by function allows us to split the table into different categories, and compute summary statistics for each year (for example).

gapminder %>% 
    group_by(year) %>% 
    summarise(MinLifeExpectancy = min(lifeExp), 
              MaximumGDP = max(gdpPercap), 
              AveragePop = mean(pop))

Other summary statistics that can be useful include first and last which are used to report the first and last values for a particular group. For instance, we might want to look at the increase in wealth over time for each country by extracting their gdpPercap in 1952 and 2007.

gapminder %>% 
  group_by(country) %>% 
  arrange(year) %>% 
summarise(StartGDP = first(gdpPercap), EndGDP = last(gdpPercap), GDPIncrease = EndGDP - StartGDP)

The nice thing about summarise is that it can followed up by any of the other dplyr verbs that we have met so far (select, filter, arrange..etc).

Returning to the correlation between life expectancy and year, we can summarise as follows:-

gapminder %>%     
    group_by(country) %>% 
    summarise(Correlation = cor(year , lifeExp))

We can then arrange the table by the correlation to see which countries have the lowest correlation

gapminder %>%      
    group_by(country) %>% 
    summarise(Correlation = cor(year , lifeExp)) %>% 
    arrange(Correlation)

We can filter the results to find obsevations of interest

gapminder %>%      
    group_by(country) %>% 
    summarise(Correlation = cor(year , lifeExp)) %>% 
    filter(Correlation < 0)

The countries we identify could then be used as the basis for a plot.

filter(gapminder, country %in% c("Rwanda","Zambia","Zimbabwe")) %>% 
  ggplot(aes(x=year, y=lifeExp,col=country)) + geom_line()

---

---

---

Exercise

• Produce a plot to show the change in average gdpPercap for each continent over time.
• see below for a suggestion
  • HINT: you will need to specifiy the stat=identity option when creating the bar plot

......+ geom_bar(stat="identity")

---

---

---

Joining

In many real life situations, data are spread across multiple tables or spreadsheets. Usually this occurs because different types of information about a subject, e.g. a patient, are collected from different sources. It may be desirable for some analyses to combine data from two or more tables into a single data frame based on a common column, for example, an attribute that uniquely identifies the subject.

dplyr provides a set of join functions for combining two data frames based on matches within specified columns. For those familiar with such SQL, these operations are very similar to carrying out join operations between tables in a relational database.

As a toy example, lets consider two data frames that contain the names of various bands, and the instruments that they play:-

band_instruments

band_members

There are various ways in which we can join these two tables together. We will just consider the case of a “left join”.

Animated gif by Garrick Aden-Buie

left_join returns all rows from the first data frame regardless of whether there is a match in the second data frame. Rows with no match are included in the resulting data frame but have NA values in the additional columns coming from the second data frame.

Animations to illustrate other types of join are available at https://github.com/gadenbuie/tidy-animated-verbs

left_join(band_members, band_instruments)

Joining, by = "name"

right_join is similar but returns all rows from the second data frame that have a match with rows in the first data frame based on the specified column.

right_join(band_members, band_instruments)

Joining, by = "name"

inner_join only returns those rows where matches could be made

inner_join(band_members, band_instruments)

Joining, by = "name"

---

---

---

Exercise (open-ended)

• The file medal_table.csv in the raw_data/ project sub-directory contains data about how many medals how been won by various countries at the Beijing summer olympics of 2008.
• Read this csv file into R and join with the gapminder data from 2007
• What interesting summaries / plots can you make from the data? For example…
• what countries have the greatest proportion of gold medals (ignore countries with too few medals)
• calculate the number of medals won per million people and re-arrange by this new measure. What countries perform best?
• how similar is the distribution of total medals between continents?
• do countries with a larger population tend to win more medals?
• do countries with larger GDP tend to win more medals?
• are these trends consistent among different continents?

---

---

---

Using R to analyse RNAi data

You should have access to an example RNAi screen dataset; rnai.csv. In RStudio, create a new project in the directory that this csv is saved.

You can create a new markdown file to document the analysis through the menus.

File -> New File -> R markdown >

The first code chunk should make sure we have the tidyverse package loaded.

library(tidyverse)

As a quality control measure, we can look at the signal from each plate and make sure there are no systematic differences. This could be visualised as a boxplot with geom_boxplot as we have seen before.

ggplot(rnai, aes(x=as.factor(Plate), y=`FL/RL_rep1`)) + geom_boxplot()

We also have some positional information at our disposal that could be of interest. This is encoded in the well name, but we can split the well name into a row and column ID using a mutate command with substr to extract the letter and number of the well.

rnai <- mutate(rnai,Row = substr(Well,1,1), Col=substr(Well,2,3))
rnai

This can be visualised using a type of geom that we haven’t seen before called geom_tile that produces a “heatmap” style of plot.

ggplot(rnai, aes(x=Row,y=Col,fill=`FL/RL_rep1`)) + geom_tile() + facet_wrap(~Plate)

---

---

---

Exercise

• Try and re-produce the geom_boxplot and geom_tile above with the FL/RL_Av column instead of FL/RL_Rep1
  • R will struggle to produce the plot. Look at the contents of this column and think about why this might happen?
• Apply appropriate filtering to the data in order to produce a boxplot of the FL/RL_Av values for each Plate
  • HINT: the R code !is.na can be used to check identify what rows do not contain an NA value in a particular column

---

---

---

Before proceeding to further analysis we can clean the dataset so that we only have the gene identifiers and average FL and RL values as columns.

rnai_cleaned <- filter(rnai, !is.na(BKN)) %>% 
  mutate(FL=as.numeric(`FL_Av`), RL=as.numeric(`RL_Av`), Ratio = as.numeric(`FL/RL_Av`)) %>% 
  select(GeneSymbol, GeneID,FL,RL,Ratio)
rnai_cleaned

The Gene identifiers are not particularly useful by themselves but we can use resources such as Biomart to map these to more recognisable names, and also obtain homologous human genes. The biomaRt package is able to perform the mapping, and without going into too much detail here is the code to do so.

library(biomaRt)
genes_of_interest <-rnai_cleaned$GeneID
ensembl <- useEnsembl(biomart = "ensembl",
                      dataset = "dmelanogaster_gene_ensembl")
anno <- getBM(c("ensembl_gene_id","external_gene_name", "hsapiens_homolog_ensembl_gene","hsapiens_homolog_associated_gene_name","hsapiens_homolog_orthology_confidence"), "external_gene_name", genes_of_interest, ensembl) %>% 
  rename(GeneID = external_gene_name)

This can now be joined to the cleaned dataset using the left_join command from before.

rnai_cleaned <- left_join(rnai_cleaned, anno)

Joining, by = "GeneID"

rnai_cleaned

A basic scatter plot of FL and RL reveals some extreme values in the dataset.

ggplot(rnai_cleaned, aes(x = FL, y = RL,col=Ratio)) + geom_point()

---

---

---

Exercise

• Identify the names of the genes that have extreme RL or FL values
• Produce a version of the plot that has the name of these genes labelled
  • check the identity of these genes. Do they have any biological significance?

• Decide on suitable cut-offs on the ratio of RL to FL to identify genes of interest
• Write the rows corresponding to these genes to a csv file for further processing