80-prog.Rmd

# Additional programming concepts {#sec-prog}

**Learning Objectives**

Learn programming concepts, including

- how to handle conditions
- iteration of data structures
- good coding practice
- code re-use through functions

When the size and the complexity of the data increases, or the data
science question of interest becomes more complex, the data analysis
techniques as we have seen them so far need to be complemented with
programming techniques. From a data science point of view, there is no
clear delimitation between data analysis and (data) programming, both
morphing into each other[^dataprog].

[^dataprog]: A fundamental difference however is how data analysis and
    programming are taught. When it comes to researchers, and
    biomedical researchers in particular, teaching programming to
    analyse data isn't successful. Teaching data analysis to
    eventually programme with data, however, has proven a successful
    strategy.

This chapter will introduce some additional programming skills and
demonstrate how to use them in the context of high troughput omics
data.

## Writing clean code


Writing clean code means writing easily readable code, hence easily
understandable code and, eventually code with less bugs. One easy way to
achieve this is through consistency, i.e. stick to a **style
guide**. The issue is that there are several style guides available,
often with conflicting suggestions. Two widely used ones are the
[Bioconductor style
guide](http://bioconductor.org/developers/how-to/coding-style/) and
[Hadley Wickhams's R Style
Guide](http://r-pkgs.had.co.nz/style.html). The advice from for naming
variables seen in the first and second chapters are also relevant.

Here are a couple of suggestions:

- Use `<-` to assign variables. Use `=` is also valid, but make sure
  that you stick with one.
- Use either camel case (`camelCaseNaming`) of snake case
  (`snake_case_naming`), and avoid using dots (don't use
  `dot.variable.names`). These conventions apply to functions,
  variables and files (for the latter, a `-` instead of `_` is also
  acceptable).
- Always spell out `TRUE` and `FALSE`, and resist the temptation to
  use `T` and `F` instead[^TF].
- Use 4 spaces for indenting. No tabs.
- No lines longer than 80 characters.
- Use spaces around binary operators.
- No spaces between a function name and the opening parenthesis.

[^TF]: `TRUE` and `FALSE` are reserved words; one can't use them as
variable names. This however doesn't hold for `T` and `F`. One can use
`T` and `F` as regular variable names such, as for example,
`T <- FALSE` and `F <- TRUE`!

Another good advise is to avoid re-writing the same code many
times. We will see below two strategies to do that, namely iteration
and writing new functions. This firstly reduces the amount of code
typed and hence the number of bugs, but more importantly enables
consistency. If something in your code changes, there's only one
change and it applies everywhere, rather than doing that same change
repeatedly (at the first of adding bugs and to miss some updates).

When writing code, **keep it simple and short**[^KISS]. Ideally, the
code should be evident. But when the questions tackled are not
trivial, is becomes essential to add **comments** to clarify aspects
of the script/program. Make sure to use them to describe why
something is done, rather than explaining how things are done (which
is tyically best done by the code itself).

[^KISS]: The official motto is
[KISS](https://en.wikipedia.org/wiki/KISS_principle), *keep your
functions simple, stupid*, widely used in programming.

A general guideline to avoid bugs is to apply **defensive
programming**:

- making the code work in a predicable manner
- writing code that fails in a well-defined manner
- if something weird happens, either properly deal with it, or it fail
  quickly and loudly

Here are some examples of how to do this:

- use functions like `is.numeric(x)`, `is.character(x)`,
  `is.data.frame(x)`, ... but make sure that the variable you are
  going to use is of the expected type.

- Make sure the length or dimensions of what you are going to use are
  what you expect:

```{r, eval = FALSE}
length(x) > 0
nrow(x) > 0
ncol(x) > 0
```

- Failing fast and well! Wrap such test in `stopifnot()`, so that if
  they fail, you get immediately an error, rather than risking that
  you code fails later with obscure error messages or, worse, the code
  runs to completion but returns meaningless results.

```{r, eval = FALSE}
stopifnot(length(x) > 0)
stopifnot(dim(x) > 0) ## same as next one
stopifnot(ncol(x) > 0, nrow(x) > 0)
```

## Iteration

Iteration describes the situtation when a specific operation has to be
repeated many times on different inputs of the same type. For example,
if we have a vector of numeric `x` shown below,

```{r}
(x <- 1:10)
```

and we wanted to apply the logarithm operation on each element of `x`,
it wouldn't be convenient to type

```
log(1)
log(2)
log(3)
...
log(10)
```

The concept of iteration allows us to program the following command:

> Repeat `log` for each value of my input `x`.


or more formally

> Repeat `log(i)` where `i` takes in turn each value of my input `x`.


We will see different ways of implementing such an iteration.

### Using a `for` loop {-}

```{r}
for (i in x) {
    print(log(i))
}
```

The loop above only prints the results on screen. They aren't stored
and are lost for any further re-use, which would be very annoying if
it took much more time to perform all the calculations. In the code
chunk below, we are going to first inititalise a vector with the
appropriate number of `NA` values and, at each iteration, we then
store the result. We however now need to change our loop and iterate
of the indices of the input vector, so that we can re-use these
indices to save the results in the output vector.


```{r}
res_loop <- rep(NA, length(x))
for (i in seq_along(x)) {
    res_loop[i] <- log(x[i])
}
res_loop
```

### Using the `apply` function {-}

The apply family of functions implements our defintion of iteration
quite literally

> Repeat `log` for each value of my input `x`

is reformulated as

> For each value of my input `x`, apply `log`

and coded as


```{r}
res_apply <- sapply(x, log)
res_apply
```

There are three such function that *apply* a function iteratively:

- `sapply` iterates over a vector and returns a new vector of the same
  length as the input vector[^sapplyexception].

- `lapply` iterates over a linear input (a vector of a list) and
  returns a list of the same length as the input.

- `apply` iterates of the rows or the columns of a `data.frame` or
  `matrix` and returns a list or vector or approritate length (number
  of rows or columns). The dimension over which the iterations
  proceeds is defined by the second argument, where `1` defines rows
  and `2` defined columns[^arrayapply].


[^sapplyexception]: This is a simplification of how `sapply` works,
    that is partly defined by the the `simplify` argument and whether
    the result of applying the function on each element of the input
    can be returned as a vector. The alternative is to return a
    `list`, like `lapply`.

[^arrayapply]: This generalises to arrays with more than 2 dimensions.

When performing the same operations in different ways (using different
implementations), it is essential to verify that the results are
identical or, if not, at least compatible. The most direct way to do
the former is to use the `identical` function:

```{r}
identical(res_loop, res_apply)
```

The `purrr` package a set of `map` functions similar to the `apply`
set described above.

### Vectorisation {-}

We must not forget the obvious, which is vectorisation. Many R
functions work by default iteratively on every element of a vectors,
i.e. they work irrespectively whether the vector is of length
1[^scalar] or longer.

[^scalar]: A vector of length 1 would be called a scalar in other
    programming languages.


```{r}
res_vec <- log(x)
res_vec
```

And, as before, we check that we obtain identical results:

```{r}
identical(res_loop, res_vec)
```

### Which iteration to use? {-}

Even though they produce the same result, the iteration strategies
above aren't equal, and some should be preferred in different
situations.

- When a vectorised solution exists, this is the one that should be
  chosen. It is by far the fastest solution, but only applicable to
  existing *vectorised* functions. If you were to write your own
  function to iterate over, it is advisable to write a vectorised
  function.

- The apply functions are extremely convenient and concise, and hence
  widely used. They have a couple of advantages, including that there
  is no need to explicitly initialise a result variable and that they
  can easily be parallelised.

- For loops are the most generic solution for iteration, and they
  require to initialise the result variable that will be populated
  during the loop. As opposed to popular belief, they aren't slower
  than using apply functions. They are however the best solution if,
  during the iteration, one has the access another element than the
  one currently processed (typically, if `i` is the counter, accessing
  `i + 1` or `i - 1`).


`r msmbstyle::question_begin()`

```{r, message = FALSE, echo = FALSE}
library(rWSBIM1207)
fs <- expressions.csv()
```

Use the `expressions.csv()` function from the `rWSBIM1207` package to
get the path to `r length(fs)` csv files containing gene expression data
for a single gene and hundreds of patience. Read all data in and
combine them into a single dataframe.

Before starting, open at least two of these files to familiarise
yourself with their structure and identify how and with what function
they need to be combined.

`r msmbstyle::question_end()`

`r msmbstyle::solution_begin()`

Let's start by loading the required packages and check two file. We
will also use the *tidyverse* to work with the data.

```{r}
library("tidyverse")
library("rWSBIM1207")
fls <- expressions.csv()
read_csv(fls[1])
read_csv(fls[2])
```

We see that the data will have to be combined column-wise. We could do
it manually, adding news columns one by on to the first data, but this
would require to check that the observations are in the same
order. Instead, we will use a `full_join`, as seen in chapter
\@ref(sec-join).

```{r, message = FALSE}
## read the first one
res <- read_csv(fls[1])
## iterate over the other data
for (i in 2:length(fls)) {
    data_i <- read_csv(fls[i], progress = FALSE)
    res <- full_join(res, data_i)
}
res
```

`r msmbstyle::solution_end()`

<!-- Possibly also mention `while`. -->


## Conditionals

When an operation has to be executed when a condition is met, one
typically uses a `if` and `else` construct:

```
if (CONDITION) {
   ## DO SOMETHING
} else {
  ## DO SOMETHING ELSE
}
```

For examples

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

```{r}
(x <- rnorm(1))

if (x > 0) {
   print(log(x))
} else {
   print(log(-x))
}
```

But note that in the example above, it would be much better to
simplify the code and use the absolute value of `x` before taking the
log, which will generalise the calculation for positive and negative
values ...

```{r}
log(abs(x))
```

... and works for vectors of any length thanks to vectorisation.

```{r}
x <- rnorm(10)
log(abs(x))
```

There are also in-line, vectorised versions of the `if/else`
statements seen above: the `ifelse` function that ships with the
`base` package, and `if_else` from `dplyr`. Both work similarly (the
differences are beyond the scope of this course) and have the form:

```
if_else(condition, true, false)
```

- where `condition` is the condition to be tested, and will return
  either `TRUE` or `FALSE`;
- `true` is the expression that is executed if the condition is
  `TRUE`;
- `false` is the expression that is executed otherwise (i.e if the
  condition is `FALSE`).

Here is an example, that will return 1 if `x` is strictly positive,
and 0 otherwise:


```{r}
x <- 0.5
if_else(x > 0, 1, 0)
```

The function being vectorised, the condition can be a vector of length
greater than 1:


```{r}
x <- rnorm(10)
x
if_else(x > 0, 1, 0)
```

The vectorised conditional functions can directly be used in a
standard data analysis pipeline:

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

```{r}
x <- tibble(x = letters[1:5],
            y = rnorm(5))
x
mutate(x, z = if_else(y > 0, 1, 0))
```

## Writing new functions

A function is composed of a name, inputs (inside the parenthesis), a
body (between curly brackets) and an ouput (last statement or variably
inside the `return` statement).

```{r myfun}
my_fun <- function(x, y) {
    message("First input: ", x)
    message("Second input: ", y)
    z <- x * abs(y)
    return(z)
}
my_fun(2, -5)
```

`r msmbstyle::question_begin()`

```{r, echo = FALSE}
fun <- function(x, y) {
    res <- NA
    if (x > y) {
        res <- sample(x, y)
    } else { ## here x <= y
        res <- rnorm(1000, x, y)
    }
    return(sum(res))
}
set.seed(1L)
res1 <- fun(5, 15)
set.seed(1L)
res2 <- fun(15, 5)
```


Complete the following function. It is supposed to take two inputs,
`x` and `y` and, depending whether the `x > y` or `x <= y`, it generates
the permutation `sample(x, y)` in the first case or draws a sample
from `rnorm(1000, x, y)` in the second case. Finally, it returns the
sum of all values.


```{r, eval = FALSE}
fun <- function(x, y) {
    res <- NA
    if (   ) {
        res <- sample(, )
    } else {
        res <- rnorm(, , )
    }
    return()
}
```

To check your answer, run it with inputs 5, 15 and 15, 5, after
setting the random number generator to 1 each time with `set.seed(1)`
and you should get `r res1` and `r res2`.

`r msmbstyle::question_end()`

`r msmbstyle::solution_begin()`

```{r}
fun
set.seed(1L)
fun(5, 15)
set.seed(1L)
fun(15, 5)
```

`r msmbstyle::solution_end()`


## Pattern matching

One skill that comes handy more often than not is the ability to find
patterns in text and replace these. We are going to see two functions
out of many to perform such tasks. As an illustration, we are going to
mine a vector or peptides identified by mass spectrometry. The data
can be loaded from the `rWSBIM1207` data:

```{r}
library("rWSBIM1207")
data(peptides)
length(peptides)
head(peptides)
```

1. Among the `r length(peptides)` petides, which ones contain the pattern `"AAGE"`?

```{r}
pattern <- "AAGE"
## indices of matching peptides
grep(pattern, peptides)
## matching peptides
grep(pattern, peptides, value = TRUE)
## position of matching peptides
head(grepl(pattern, peptides))
```

2. Replace a pattern by another string.

```{r}
head(peptides)
head(sub("AE", "XX", peptides))
```

But careful that `sub` only replaced the first occurrence in a string!

```{r}
head(peptides)
head(sub("A", "X", peptides))
head(gsub("A", "X", peptides))
```

## Analysing data from multiple files

This exercise recapitulates the most important material that we have
seen in this course. We are going to analyse the student tests A and B
results using functionality from the `dplyr` package (seen in chapter
\@ref(sec-dplyr)). The student test files are available in the
`rWSBIM1207` package (the `interroA.csv` and `interroB.csv` functions
return their respective paths). We want to compare how the male and
female students in groups A and B have performed. To do this, we want
to calculate the mean scores and visualise the score distributions in
each groups.

`r msmbstyle::question_begin()`

Before writing any code to answer the questions above, take a couple
of minutes to think about what packages are needed to answer the
questions above and identify the individual steps needed to first
prepare the data needed to find the answers and then to produce the
answers.

`r msmbstyle::question_end()`


1. Start by loading the `rWSBIM1207` and `tidyverse` package, to get
   all the functionality we will need.


`r msmbstyle::solution_begin()`
```{r, message}
library("rWSBIM1207")
library("tidyverse")
```
`r msmbstyle::solution_end()`


2. Using the `interroA.csv()` and `interroB.csv()` functions from the
   `rWSBIM1207` package, create a vector containing the two csv file
   names.

`r msmbstyle::solution_begin()`

```{r}
fls <- c(interroA.csv(),
         interroB.csv())
fls
```

`r msmbstyle::solution_end()`


3. Iterate over the files and read each into an element of a
   list. Then, bind the two elements of the list into a single
   dataframe containing all results.

`r msmbstyle::solution_begin()`

```{r}
l <- as.list(fls)
for (i in seq_along(fls))
    l[[i]] <- read_csv(fls[i])
interros <- rbind(l[[1]], l[[2]])
interros
dim(interros)
```

or

```{r}
l <- lapply(fls, read_csv)
interros2 <- rbind(l[[1]], l[[2]])
identical(interros, interros2)
```

`r msmbstyle::solution_end()`

4. What are the mean scores for male and female students in each test
   (`interro1` to `interro4`).

`r msmbstyle::solution_begin()`

```{r}
interros %>%
    group_by(gender) %>%
    summarise(m1 = mean(interro1, na.rm = TRUE),
              m2 = mean(interro2, na.rm = TRUE),
              m3 = mean(interro3, na.rm = TRUE),
              m4 = mean(interro4, na.rm = TRUE))
```

`r msmbstyle::solution_end()`


5. What are the mean scores for male and female students in each test
   (`interro1` to `interro4`) and in each group (A and B).

`r msmbstyle::solution_begin()`

To answer this question, we need to be able to tell the two groups
apart. This information is however not available in the full
dataset. We could analyse the two test individually to answer this
particular question, but a preferred approach is the add this
information, as it will be needed later anyway.

One way would be to check the number of students that tool the two
tests, and add a new column based on these values. We still have the
two dataframes in the `l` list, which allows us to calculate the
respective number of rows and use this information to create a group
variables.

```{r}
sapply(l, nrow)
interros$group <- c(rep("A", nrow(l[[1]])),
                    rep("B", nrow(l[[2]])))
interros
```

The assumption here is that the order of the observations in the large
dataframe hasn't changed!

An alternative is to add the group information before combining the
dataframes.

```{r}
l[[1]]$group <- "A"
l[[2]]$group <- "B"
interros2 <- rbind(l[[1]], l[[2]])
identical(interros, interros2)
```

Now we can group by `gender` and `group` to perform our summary
calculations:

```{r}
interros %>%
    group_by(gender, group) %>%
    summarise(m1 = mean(interro1, na.rm = TRUE),
              m2 = mean(interro2, na.rm = TRUE),
              m3 = mean(interro3, na.rm = TRUE),
              m4 = mean(interro4, na.rm = TRUE))
```

`r msmbstyle::solution_end()`

6. Instead of comparing the means for the four tests in each group and
   gender, visualise their distributions using boxplots. In addition
   to the boxplots, overlay the actual score values (using the
   *jitter* geom) and colour the points based on the students height.


`r msmbstyle::solution_begin()`

The produce this visualisation, we first need to convert the data we
have in a long format using `pivot_longer()`. Below, we create two new
variables: `test`, that contains `interro1`, ..., `interro4`, and
`score`, that contains the students marks the actual test.

We obtain a new tibble that contains 4 times more rows that initially,
as the tests, that were in 4 different columns, have now been gathered
as a single new variable.

```{r}
interros2 <- pivot_longer(interros,
                          names_to = "test",
                          values_to = "score",
                          starts_with("interro"))
interros2
```

We can now map the new `test` and `score` variables to the *x* and *y*
axes and construct the visualisation.

```{r}
ggplot(interros2, aes(y = score, x = test)) +
    geom_boxplot() +
    geom_jitter(aes(colour = height)) +
    facet_grid(gender ~ group)

```

`r msmbstyle::solution_end()`


## Additional exercises

`r msmbstyle::question_begin()`

Complete the following function that converts meters to yards, knowing
that that 1 yard corresponds to 0.9144 meters

```{r, eval = FALSE}
meters_to_yards <- function(argument goes here) {
  yards <- _code goes here!_
  return(yards)
}

```
`r msmbstyle::question_end()`

`r msmbstyle::question_begin()`

Write functions that convert temperatures:

- from Fahrenheit to Kelvin, where $T_K = (T_F - 32) \times 5/9 + 273.15$
- fom Celsius to Kelvin, where $T_K = T_C + 273.15$
- from Kelvin to Celsis
- from Fahrenheit to Celsius

`r msmbstyle::question_end()`