---
title: Introduction to R answers
description: Worked answers for the exercises in the Introduction to R with Tidyverse course
author: Simon Andrews, Laura Biggins
author-affiliation: Babraham Institute
title-block-banner: true
format: 
  html:
    toc: true
    df-print: paged
editor: visual
embed-resources: true
---

This document provides worked answers for all of the [exercises](https://www.bioinformatics.babraham.ac.uk/training/Introduction_to_R_tidyverse/Introduction%20to%20R%20with%20Tidyverse%20Exercises.pdf) in the [Introduction to R with Tidyverse](https://www.bioinformatics.babraham.ac.uk/training.html#rintrotidy) course.

# Exercise 1

## Numbers

We start by doing some simple calculations.

```{r}
31 * 78

697 / 41
```

We next look at how to assign data to named variables and then use those variables in calculations.

We make assignments using arrows and they can point to the right or the left depending on the ordering of our data and variable name.

```{r}
x <- 39

y <- 22
22 -> y # This is equivalent to the line above
```

We can then use these in calculations instead of re-entering the data

```{r}
x - y
```

We can also save the results directly into a new variable.

```{r}
z <- x - y

z
```

## Using functions

We want to calculate the square root of 2345, and perform a log2 transformation on the result.

This could be done in two steps, by creating an intermediate variable.

```{r}
sqrt_value <- sqrt(2345)
log2(sqrt_value)
```

Or we can nest the functions.

```{r}
log2(sqrt(2345))
```

## Text

We can also store text. The difference with text is that we need to indicate to R that this isn't something it should try to understand. We do this by putting it into quotes.

If we forget the quotes then R will look for a variable called `simon` and will produce an error when it can't find it.

```{r}
#| #error: true
#| eval: false
# do not copy this code - it will produce an error!
my_name <- simon
# Error: object 'simon' not found
```

It works nicely when we include the quotes.

```{r}
my_name <- "simon"
```

We can use the `nchar()` function to find out how many characters `my_name` has in it.

```{r}
nchar(my_name)
```

We can also use the `substr()` function to get the first letter of `my_name`.

```{r}
substr(my_name, start = 1, stop = 1)
```

# Exercise 2

## Making vectors manually

We're going to manually make some vectors to work with. For the first one there is no pattern to the numbers so we're going to make it completely manually with the `c()` function.

```{r}
some_numbers <- c(2, 5, 8, 12, 16) 
```

For the second one we're making an integer series, so we can use the colon notation to enter this more quickly.

```{r}
number_range <- 5:9
```

Now we can do some maths using the two vectors.

```{r}
some_numbers - number_range
```

Because the two vectors are the same size then the equivalent positions are matched together. Thus the final answer is:

`(2-5), (5-6), (8-7), (12-8), (16-9)`

## Vector creating functions

We're going to use some functions to create vectors.

First we're going to make a numerical sequence with the `seq()` function.

-   Use `seq()` to make a vector of 100 values starting at 2 and increasing by 3 each time and store it in a new variable called `number_series`

```{r}
number_series <- seq(from = 2, by = 3, length.out = 100)

number_series
```

-   Multiply every value in `number_series` by 1000 and overwrite the original variable.

```{r}
number_series <- number_series * 1000

number_series
```

-   Use `rep()` to make a vector of 100 values containing 25 each of WT, KO1, KO2 and KO3

We can create a vector containing a single WT, KO1, KO2 and KO3 and check that that works first.

```{r}
c("WT", "KO1", "KO2", "KO3")
```

Then we can include that within the `rep()` function. We'll get different results depending on whether we use the `times` or `each` option.

```{r}
rep(c("WT", "KO1", "KO2", "KO3"), times = 25)
```

```{r}
rep(c("WT", "KO1", "KO2", "KO3"), each = 25)
```

## Statistical functions and vectors

Since R is a language built around data manipulation and statistics we can use some of the built in statistical functions.

We can use `rnorm` to generate a sampled set of values from a normal distribution

```{r}
normal_numbers <- rnorm(20)
```

Note that if you run this multiple times you'll get slightly different results.

We can now use the `t.test` function to test whether this vector of numbers has a mean which is significantly different from zero.

```{r}
t.test(normal_numbers)
```

Not surprisingly, it isn't significantly different.

If we do the same thing again but this time use a distribution with a mean of 1 we should see a difference in the statistical results we get.

```{r}
t.test(rnorm(20, mean=1))
```

This time the result is significant.

## (Optional) Median Centering

To median centre the values in `number_series` we simply calculate the median value and then subtract this from each of the values in the vector.

```{r}
number_series - median(number_series)
```

To calculate the mean and standard deviation of these median centred values, we can pass them to the `mean` and `sd` functions.

```{r}
median_centred <- number_series - median(number_series)

mean(median_centred)
sd(median_centred)
```

## (Optional) T-test power

If we are sampling from two distributions with only a 1% difference in their means, how many observations do we need to have before we can detect them as significantly changing?

Let's try a few different thresholds to see.

```{r echo=FALSE}
set.seed(1235)
```

### 100 Samples

```{r}
samples <- 100

data1 <- rnorm(samples, mean = 10, sd = 2)
data2 <- rnorm(samples, mean = 10.1, sd = 2) 

t.test(data1, data2)
```

### 500 Samples

```{r}
samples <- 500

data1 <- rnorm(samples, mean = 10, sd = 2)
data2 <- rnorm(samples, mean = 10.1, sd = 2)

t.test(data1, data2)
```

### 1000 Samples

```{r}
samples <- 1000

data1 <- rnorm(samples, mean = 10, sd = 2)
data2 <- rnorm(samples, mean = 10.1, sd = 2)

t.test(data1, data2)
```

### 5000 Samples

```{r}
samples <- 5000

data1 <- rnorm(samples, mean = 10, sd = 2) 
data2 <- rnorm(samples, mean = 10.1, sd = 2)

t.test(data1, data2)
```

It's only really when we get up close to 5000 samples that we can reliably detect such a small difference. The answers will be different every time since `rnorm` involves a random component.

# Exercise 3

We're going to read some data from a file straight into R. To do this we're going to use the tidyverse `read_` functions. We therefore need to load tidyverse into our script.

```{r}
library(tidyverse)
```

## Reading a small file

We'll start off by reading in a small file.

```{r}
small <- read_tsv("small_file.txt")

small
```

Note that the only relevant name from now on is `small` which is the name we saved the data under. The original file name is irrelevant after the data is loaded.

We can see that Sample and Category have the 'character' type because they are text. Length has the 'double' type because it is a number.

We want to find the median of the log2 transformed lengths.

To start with we need to extract the length column using the `$` notation.

```{r}
small$Length
```

Now we can log2 transform this.

```{r}
log2(small$Length)
```

Finally we can find the median of this.

```{r}
median(log2(small$Length))
```

## Reading a larger variants file

The second file we're going to read is a CSV file of variant data. We can still use `read_delim` to read it in.

```{r}
child <- read_delim("Child_Variants.csv") 

child
```

We can see that each row is a set of different parameters for a genetic variant.

We want to calculate the mean and standard deviation (sd) of the MutantReadPercent column.

We can again use the `$` notation to extract the relevant column into a vector. We can use `head` to just look at the first few values.

```{r}
head(child$MutantReadPercent, n = 20)
```

We can then pass this to the relevant functions to calculate the mean and sd.

```{r}
mean(child$MutantReadPercent)

sd(child$MutantReadPercent)
```

# Exercise 4

In this section we are going to use two tidyverse functions to cut down our data. We will use `select` to pick only certain columns of information to keep, and `filter` to apply functional selections to the rows.

## Filtering small data

Extract all of the data for Category A.

```{r}
filter(small, Category == "A")
```

Filter for all samples with a length of more than 80.

```{r}
filter(small, Length > 80)
```

Remove the Sample column. We do a negative column selection by using a minus before the name.

```{r}
select(small, -Sample)
```

We could have instead selected the columns we wanted to keep.

```{r}
select(small, Length, Category)
```

## Filtering child variants

Select data from the mitochondrion (called MT in this data).

```{r}
filter(child, CHR == "MT")
```

Select variants with a MutantReadPercent of greater than (or equal to) 70.

```{r}
filter(child, MutantReadPercent >= 70)
```

Select variants with a quality of 200.

```{r}
filter(child, QUAL == 200)
```

Select all variants of the IGFN1 gene.

```{r}
filter(child, GENE == "IGFN1")
```

Remove the ENST and dbSNP columns.

```{r}
select(child, -ENST, -dbSNP)
```

# Exercise 5

Here we are going to do some more complicated multi-step filtering. For this we will be using the `%>%` or `|>` pipe operator to link together different steps in the selection and filtering.

## Quality filtering

We want variants with a quality of 200, coverage of more than 50 and MutantReadPercent of more than 70.

```{r}
child %>%
  filter(QUAL == 200) %>%
  filter(COVERAGE > 50) %>%
  filter(MutantReadPercent > 70)
```

Same as above but using the base R pipe.

```{r}
child |>
  filter(QUAL == 200) |>
  filter(COVERAGE > 50) |>
  filter(MutantReadPercent > 70)
```

## Positional Filtering

We want to remove all variants on the X, Y and MT chromosomes. At the moment we will have to do this in three separate filtering steps, but there is actually a quicker way to do it in one step if we get a bit more advanced.

```{r}
child |>
  filter(CHR != "X") |>
  filter(CHR != "Y") |>
  filter(CHR != "MT")
```

## Annotation Filtering

We want to see the chromosome and position only for variants with a valid dbSNP id. We can get the valid dbSNP variants by excluding the variants with a dot as their ID.

```{r}
child |>
  filter(dbSNP != ".") |>
  select(CHR, POS)
```

## Transform - find deletions

Here we're going to do a more complicated filtering where we will transform the data within the filter statement. In this example we're going to find deletions. In a deletion the REF column will have more than one character in it, so that's what we're going to select for.

We can use the `nchar()` function as before, but this time we're going to use it in a filter condition.

```{r}
child |>
  filter(nchar(REF) > 1)
```

## Transform - Q genes

As a second example we're going to use another transformation along with a conventional selection to find genes whose name starts with a Q.

We will use the `substr()` function to find the genes which start with Q.

```{r}
child |>
  filter(substr(GENE, 1, 1) == "Q")
```

::: {.callout-tip collapse="true"}
### Extended explanation

To break this down.

```{r}
gene_names <- c("AK5", "NEXN", "QSER1", "GIPC2", "QPRT")
substr(gene_names, 1, 1)
```

```{r}
substr(gene_names, 1, 1) == "Q"
```

This output is a logical vector. Inside our filter functions we want logical tests, and these logical tests produce logical vectors that are used to filter the dataset.
:::

# Exercise 6

Here we're going to draw our first ggplot plot. We're going to use a new dataset so we need to load that first.

```{r, message = FALSE}
brain <- read_delim("brain_bodyweight.txt")

brain
```

## Simple scatterplot

We'll start out by doing a simple scatterplot (using `geom_point()`) of the brain and bodyweight values.

```{r}
brain |>
  ggplot(aes(x = log.brain, y = log.body)) +
  geom_point()
```

We can change the colour and size of the points by adding some options to the `geom_point()` call. These are fixed values rather than aesthetic mappings so we don't need to put them in a call to `aes()`.

```{r}
brain %>%
  ggplot(aes(x = log.brain, y = log.body)) +
  geom_point(colour = "blue", size = 3)
```

## Adding a colour aesthetic

Before we modified the plot by adding in a fixed colour parameter. Here we're going to add a third aesthetic mapping where we use the colour of the points to represent the Category that the animal falls into.

```{r}
brain %>%
  ggplot(aes(x = log.brain, y = log.body, colour = Category)) +
  geom_point(size = 3)
```

Now we can see the three groups coloured differently, and a legend explaining the colours appears on the right.

We can remove the extinct species from the plot by filtering the dataset before we pipe it into ggplot.

```{r}
brain |>
  filter(Category != "Extinct") |>
  ggplot(aes(x = log.brain, y = log.body, colour = Category)) +
  geom_point(size = 3)
```

# Exercise 7

We're going to get into more complex plots here, with more configuration of the static and mapped aesthetics, and the inclusion of more geometry layers.

## Gene Barplots

We want to see the lengths of all of the samples in category A in the small data set.

```{r}
small %>%
  filter(Category == "A") %>%
  ggplot(aes(x = Sample, y = Length)) +
  geom_col()
```

## Plotting Distributions

We're going to plot out the distribution of the MutantReadPercent values from child in a couple of different ways. In each case the plot itself is going to summarise the data for us.

We'll start with a histogram.

```{r, message=FALSE}
child %>%
  ggplot(aes(MutantReadPercent)) +
  geom_histogram()
```

We can pass options to change the colouring of the bars and the resolution of the histogram.

```{r}
child %>%
  ggplot(aes(MutantReadPercent)) +
  geom_histogram(binwidth = 5, fill = "seagreen", colour = "black")
```

We can also use `geom_density()` to give a continuous readout of the frequency of values rather than using the binning approach of `geom_histogram()`.

```{r}
child %>%
  ggplot(aes(MutantReadPercent)) +
  geom_density(fill = "seagreen", colour = "black")
```

## Frequency barplot

We can also construct a barplot which will summarise our data for us. In the previous barplot example we explicitly set the heights of the bars using an aesthetic mapping. Here we're going to just pass a set of values to barplot and get it to count them and plot out the frequency for each unique value.

We're going to do this for the REF bases, but only where we only have single letter REFs.

```{r}
child %>%
  filter(nchar(REF) == 1) %>%
  ggplot(aes(REF)) +
  geom_bar()
```

If we wanted to make this a bit more colourful we could colour the bars by the REF (so they'll also be different colours). This would normally create a legend to explain the colours, but because this is self-explanatory in this case we can use the `show.legend` parameter to suppress it.

```{r}
child %>%
  filter(nchar(REF) == 1) %>%
  ggplot(aes(REF, fill = REF)) +
  geom_bar(show.legend = FALSE)
```

## Additional Exercises

### Adding text

In a previous exercise we drew a plot of the log brainweight vs the log bodyweight. In the scatterplot we drew we couldn't actually see which animal was represented by which dot, so we're going to fix that now.

We can add the names to the animals by adding a `geom_text` layer to the plot. This will use the same `x` and `y` aesthetics as the points (so the labels go in the same place), but it will define a new aesthetic mapping of label with the names of the species.

We can add the new aesthetic mapping either in the original ggplot call, or in the geom_text call.

```{r}
brain %>%
  ggplot(aes(x = log.brain, y = log.body, color = Category, label = Species)) +
  geom_point() + 
  geom_text()
```

There are a few improvements we could make:

1.  Make the text black - we can add a static aesthetic in `geom_text()` to fix this.
2.  Make the text smaller. Again we can add a `size` parameter to `geom_text()` for this.
3.  Move the text off the point. We can do this with `hjust`.
4.  Create a suitable graph title.
5.  Create better axis labels.

```{r}
brain |>
  ggplot(aes(x = log.brain, y = log.body, color = Category, label = Species)) +
  geom_point(size = 3) + 
  geom_text(colour = "black", size = 2, hjust = 1.2) + 
  ggtitle("Brainweight vs Bodyweight") +
  xlab("Brain weight (log2 g)") +
  ylab("Body weight (log2 kg)")
```

It's still not perfect. There are other options but we'll stick with this for now.

## Calculating stats

We can see what appears to be a relationship between the weight of the brain and the body, but is it real? We can do a correlation using the `cor.test` function to check.

The function takes two vectors as arguments, so we need to extract the two relevant columns out of the brain data and pass them to the function.

```{r}
cor.test(brain$log.body, brain$log.brain)
```

We can see that the association **is** significant with p=2.44e-06. We can add this information to the graph.

## Building a linear model

We can also calculate a line of best fit for these two variables. We will do this with the `lm()` function which builds a linear model.

Building models is beyond the scope of this course, but it is different to the `cor.test()` in as much as you pass the whole of the data tibble and then say which parameters you want to predict which other parameter.

The model we will build for the `brain` dataset is: `log.body ~ log.brain` which will build a model to predict log bodyweight from log brainweight.

```{r}
lm(data = brain, formula = log.body ~ log.brain)
```

This tells us that the model line which it constructs has a slope of 1.215 and an intercept of -2.283. We can use the `geom_abline()` geometry to add this as a new layer to the existing plot.

```{r}
brain %>%
  ggplot(aes(x = log.brain, y = log.body, color = Category, label = Species)) +
  geom_point(size = 3) + 
  geom_text(color = "black", size = 3, hjust = 1.2) + 
  ggtitle("Brainweight vs Bodyweight (corr p=2.44e-06)") +
  xlab("Brain weight (log2 g)") +
  ylab("Body weight (log2 kg)") +
  geom_abline(slope = 1.215, intercept = -2.283)
```

We can see that the extinct species have their own little outgroup away from everything else. They have much smaller brains than their bodyweight would predict.

```{r}
brain |>
  ggplot(aes(x = log.brain, y = log.body, color = Category, label = Species)) +
  geom_point(size = 3) + 
  geom_text(colour = "black", size = 2, hjust = 1.2) + 
  ggtitle("Brainweight vs Bodyweight (corr p=2.44e-06)") +
  xlab("Brain weight (log2 g)") +
  ylab("Body weight (log2 kg)")
```