Introduction to Solving Biological Problems Using R - Week 4

Miriam Yeung

Modified version of material from the University of Cambridge Bioinformatics Training Unit.

4. Automation in R

Motivation

Typically when we perform an analysis we need to perform the same steps multiple times, for example, over multiple genes, samples, datasets.

With the basic analyses and graphs that we have created in the previous sessions it is possible to copy and paste the relevant section of code and adjust the code to produce multiple graphs.

However, take the scenarios where you want to:

1. Count the number of genes expressed in the patients for multiple chromosomes
2. Create boxplots of expression in ER -ve versus ER +ve patients for multiple genes
3. Creating boxplots for multiple genes if expression is different in ER +ve versus ER-ve

First lets load in the libraries we’ll need.

library(dplyr)
library(ggplot2)
library(RColorBrewer)

Lets also load in the data.

patients <-  read.delim("updated.patient.txt", stringsAsFactors = FALSE)
exprsAnno <- read.delim("anno.gene.expression.txt", stringsAsFactors = FALSE)

You can check the dimensions of the objects, to check the numbers of rows and columns they contain.

dim(patients)

[1] 168   6

dim(exprsAnno)

[1] 342 171

Take a look at the files to check what they contain and as a sanity check that all looks ok.

View(patients)
View(exprsAnno)

You can check the column names.

colnames(patients)

[1] "samplename" "age"        "er"         "grade"      "her2"      
[6] "pr"        

colnames(exprsAnno)

  [1] "NKI_154"          "NKI_199"          "NKI_35"          
  [4] "NKI_99"           "NKI_201"          "NKI_138"         
  [7] "NKI_256"          "NKI_71"           "NKI_378"         
 [10] "NKI_229"          "NKI_172"          "NKI_273"         
 [13] "NKI_383"          "NKI_276"          "NKI_321"         
 [16] "NKI_260"          "NKI_125"          "NKI_226"         
 [19] "NKI_294"          "NKI_324"          "NKI_120"         
 [22] "NKI_393"          "NKI_78"           "NKI_44"          
 [25] "NKI_297"          "NKI_214"          "NKI_140"         
 [28] "NKI_347"          "NKI_236"          "NKI_34"          
 [31] "NKI_72"           "NKI_380"          "NKI_23"          
 [34] "NKI_193"          "NKI_197"          "NKI_288"         
 [37] "NKI_27"           "NKI_29"           "NKI_221"         
 [40] "NKI_85"           "NKI_327"          "NKI_286"         
 [43] "NKI_390"          "NKI_355"          "NKI_354"         
 [46] "NKI_4"            "NKI_341"          "NKI_402"         
 [49] "NKI_51"           "NKI_280"          "NKI_157"         
 [52] "NKI_94"           "NKI_70"           "NKI_391"         
 [55] "NKI_239"          "NKI_333"          "NKI_150"         
 [58] "NKI_217"          "NKI_109"          "NKI_180"         
 [61] "NKI_404"          "NKI_135"          "NKI_302"         
 [64] "NKI_48"           "NKI_147"          "NKI_7"           
 [67] "NKI_243"          "NKI_222"          "NKI_365"         
 [70] "NKI_188"          "NKI_122"          "NKI_346"         
 [73] "NKI_304"          "NKI_31"           "NKI_60"          
 [76] "NKI_169"          "NKI_205"          "NKI_107"         
 [79] "NKI_148"          "NKI_39"           "NKI_131"         
 [82] "NKI_241"          "NKI_318"          "NKI_371"         
 [85] "NKI_79"           "NKI_287"          "NKI_271"         
 [88] "NKI_314"          "NKI_149"          "NKI_113"         
 [91] "NKI_265"          "NKI_61"           "NKI_91"          
 [94] "NKI_335"          "NKI_358"          "NKI_369"         
 [97] "NKI_98"           "NKI_307"          "NKI_284"         
[100] "NKI_374"          "NKI_401"          "NKI_123"         
[103] "NKI_106"          "NKI_336"          "NKI_215"         
[106] "NKI_110"          "NKI_38"           "NKI_86"          
[109] "NKI_292"          "NKI_181"          "NKI_363"         
[112] "NKI_366"          "NKI_62"           "NKI_183"         
[115] "NKI_348"          "NKI_126"          "NKI_184"         
[118] "NKI_190"          "NKI_325"          "NKI_129"         
[121] "NKI_359"          "NKI_381"          "NKI_403"         
[124] "NKI_357"          "NKI_161"          "NKI_128"         
[127] "NKI_104"          "NKI_338"          "NKI_136"         
[130] "NKI_296"          "NKI_59"           "NKI_218"         
[133] "NKI_245"          "NKI_203"          "NKI_319"         
[136] "NKI_254"          "NKI_251"          "NKI_295"         
[139] "NKI_283"          "NKI_146"          "NKI_166"         
[142] "NKI_385"          "NKI_281"          "NKI_95"          
[145] "NKI_356"          "NKI_308"          "NKI_97"          
[148] "NKI_322"          "NKI_186"          "NKI_328"         
[151] "NKI_275"          "NKI_270"          "NKI_11"          
[154] "NKI_114"          "NKI_185"          "NKI_195"         
[157] "NKI_272"          "NKI_361"          "NKI_264"         
[160] "NKI_345"          "NKI_360"          "NKI_237"         
[163] "NKI_213"          "NKI_191"          "NKI_323"         
[166] "NKI_337"          "NKI_187"          "NKI_298"         
[169] "HUGO.gene.symbol" "Chromosome"       "Start"           

patients has 6 columns for samplename, age, er, grade, her2, pr, for 168 patients. We can use this object to get ER/PR/HER2 status of the patients.

exprsAnno has the expression information for the 168 patients (who all have an ID that starts with “NKI”), followed by columns for gene symbol, chromosome and start position (HUGO.gene.symbol, Chromosome, Start).

Automation Example 1 - Counting genes for multiple chromosomes

We want to identify the number of the genes in the dataset for Chromosomes 1, 5, 6, 8.

First lets identify the number on chromosome 1. We can use dplyr’s filter and summarise functions. We can use summarise(n()) to get a count of rows.

exprsAnno %>%
    filter(Chromosome == "chr1") %>%
    summarise(n())

Note that tally() is a convenient wrapper function you can use instead of summarise(n())

exprsAnno %>%
    filter(Chromosome == "chr1") %>%
    tally()

That gives us the count for one chromosome (chr1). Now how do we identify the number of genes for multiple chromosomes, chromosomes 1, 5, 6 and 8?

We could copy and paste the code like below.

exprsAnno %>%
    filter(Chromosome == "chr1") %>%
    summarise(n())

exprsAnno %>%
    filter(Chromosome == "chr5") %>%
    summarise(n())

exprsAnno %>%
    filter(Chromosome == "chr6") %>%
    summarise(n())

exprsAnno %>%
    filter(Chromosome == "chr8") %>%
    summarise(n())

But as you can imagine, the more lines of code that need to be copied, pasted and edited, the more likely it is that errors will arise.

This method of copying, pasting and editing can be:

1. Tedious
2. Erroneous 

Automating Commands: Loops and flow control

• Many programming languages have ways of doing the same thing many times, perhaps changing some variable each time. This is called looping. It is a way to automate tasks.

• As we are doing the same thing multiple times, but with a different chromosome each time, we can use a loop instead
• R has two basic types of loop
  • a while loop: run some code while some condition is true (hardly ever used! Therefore will not be covered)
  • a for loop: run some code on every value in a vector

for The basic structure of a for loop:

for (element in vector){
    ... do this ...
}

• Therefore the code is only run for as many elements are in the vector. Therefore we can predict the number of times the code is run.

• Note: element and vector are just variable names and can therefore be named whatever you want them to be, as long as they satisfy the constraints on varible naming.

• Here’s how we might use a for loop to find out the number of genes on each chromosome. We store the count in a variable called numRows so we can then print it out.

chrom <- c("chr1", "chr5", "chr6", "chr8")
for (chr in chrom){
    numRows <- exprsAnno %>%
        filter(Chromosome == chr) %>%
        summarise(n())
    s <- paste("The number of genes on", chr, "is", numRows)
    print(s)
}

[1] "The number of genes on chr1 is 35"
[1] "The number of genes on chr5 is 16"
[1] "The number of genes on chr6 is 20"
[1] "The number of genes on chr8 is 12"

• The above for loop finds out the same information that we identified earlier by copying/pasting/editing
• To more accurately depict the operations of a for loop, the commands being run are similar to the following:

chr <- "chr1"
numRows <- exprsAnno %>% filter(Chromosome == chr) %>% summarise(n())
chr <- "chr5"
numRows <- exprsAnno %>% filter(Chromosome == chr) %>% summarise(n())
chr <- "chr6"
numRows <- exprsAnno %>% filter(Chromosome == chr) %>% summarise(n())
chr <- "chr8"
numRows <- exprsAnno %>% filter(Chromosome == chr) %>% summarise(n())

Exercise

Using the patients object, find the number of patients who are aged above 20, 30, 40, and print the answer to the console.

ages <- c(20, 30, 40)
### Insert your code here ###

Automation Example 2 - Creating boxplots for multiple genes

Whilst the first example only contains 3 lines of code, and therefore is less likely to lead to errors, imagine that you want to create multiple boxplots for the expression values of the following genes in ER -ve and ER +ve patients.

- FOXE1  
- TECTA  
- BAX  
- MAP3K8  
- GEMIN8

First lets create a boxplot for one gene, FOXE1, in ER -ve and ER +ve patients.

We need to select the rows where HUGO.gene.symbol == “FOXE1” and the columns with the expression values for the patients. We select the columns starting with the patient ids “NKI” as we don’t want the columns HUGO.gene.symbol, Chromosome, Start.

exprsAnno %>% 
    filter(HUGO.gene.symbol == "FOXE1") %>%
    select(starts_with("NKI")) %>%
    View()

This gives us the expression values all in one row. We need to tranpose this to plot. With dplyr we can pipe (%>%) the result into R’s t() function to easily achieve this.

exprsAnno %>% 
    filter(HUGO.gene.symbol == "FOXE1") %>%
    select(starts_with("NKI")) %>%
    t() %>%
    View()

This is a matrix (which we can see from the “V1” header in View or if we use class() around the code above). We need to convert this into a data frame, which we can do by adding another pipe.

exprsAnno %>% 
    filter(HUGO.gene.symbol == "FOXE1") %>%
    select(starts_with("NKI")) %>%
    t() %>%
    data.frame() %>%
    View()

Next we’ll save the data frame as an object called filtered.

filtered <- exprsAnno %>% 
    filter(HUGO.gene.symbol == "FOXE1") %>% 
    select(starts_with("NKI")) %>%
    t() %>%
    data.frame()

And add the gene name as a header.

colnames(filtered) <- "FOXE1"

Now we’re ready to make a boxplot. We have a data frame with all the patients values for FOXE1. We want to group these by ER status, to visualise the distribution of expression values for the ER+ve and ER-ve patients.

We can identify the ER status of the patients from the patients dataframe.

patients[, "er"]

  [1] 1 0 1 1 1 1 1 1 1 1 1 0 0 1 1 1 1 1 1 0 1 1 1 1 0 0 1 1 0 0 1 0 0 1 0
 [36] 1 0 1 1 1 1 1 1 1 1 0 0 0 1 1 1 0 1 0 1 1 1 1 1 1 1 1 0 1 0 1 0 1 1 1
 [71] 1 1 0 0 1 1 1 1 1 1 0 0 1 1 0 0 1 0 0 1 0 1 1 1 1 1 1 1 0 0 1 1 1 1 0
[106] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 0 1 1 1 0 0 1 1 0 0 1 1
[141] 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 0 1 1

A 1 means a patient is ER +ve, a 0 means they’re ER -ve. We can then tell ggplot to group the patients by the ER status column.

ggplot(filtered, aes(y = FOXE1, x = factor(patients[, "er"]))) +
    geom_boxplot()

We can also colour by ER status.

ggplot(filtered, aes(y = FOXE1, x = factor(patients[, "er"]), fill = factor(patients[, "er"]))) +
    geom_boxplot()

These are the default ggplot colours. We can also specify our own colours with scale_fill_manual. For example, we could colour the ER -ve goldenrod (yellow) and ER +ve dodgerblue and add labels.

ggplot(filtered, aes(y = FOXE1, x = factor(patients[, "er"]), fill = factor(patients[, "er"]))) +
    geom_boxplot() + 
    scale_fill_manual(values = c("goldenrod", "dodgerblue"),
                      name = "ER status",
                    labels = c("ER -ve", "ER +ve"))

We can label the axes with labs, and scale_x_discrete(labels = NULL) can be used to remove the “1” and “0” from the x axis ticks.

g <- ggplot(filtered, aes(y = FOXE1, x = factor(patients[, "er"]), fill = factor(patients[, "er"]))) +
    geom_boxplot() + 
    scale_fill_manual(values = c("goldenrod", "dodgerblue"),
                      name = "ER status",
                      labels = c("ER -ve", "ER +ve")) +
    labs(title = "Expression of FOXE1 ~ Estrogen Receptor Status",
         x = "Estrogen Receptor Status",
         y = "Expression values of FOXE1") +
    scale_x_discrete(labels = NULL)
print(g)

To centre the title theme(plot.title = element_text(hjust = 0.5)) can be used.

ggplot(filtered, aes(y = FOXE1, x = factor(patients[, "er"]), fill = factor(patients[, "er"]))) +
    geom_boxplot() + 
    scale_fill_manual(values = c("goldenrod", "dodgerblue"),
                      name = "ER status",
                      labels = c("ER -ve", "ER +ve")) +
    labs(title = "Expression of FOXE1 ~ Estrogen Receptor Status",
         x = "Estrogen Receptor Status",
         y = "Expression values of FOXE1") +
    scale_x_discrete(labels = NULL) + 
    theme(plot.title = element_text(hjust = 0.5))

We can save the plot in an object, lets call it g.

g <- ggplot(filtered, aes(y = FOXE1, x = factor(patients[, "er"]), fill = factor(patients[, "er"]))) +
    geom_boxplot() + 
    scale_fill_manual(values = c("goldenrod", "dodgerblue"),
                      name = "ER status",
                      labels = c("ER -ve", "ER +ve")) +
    labs(title = "Expression of FOXE1 ~ Estrogen Receptor Status",
         x = "Estrogen Receptor Status",
         y = "Expression values of FOXE1") +
    scale_x_discrete(labels = NULL) + 
    theme(plot.title = element_text(hjust = 0.5))

If we want to produce the plot from g we can use print().

print(g)

Now we’ve got a nice plot showing the expression of FOXE1 in ER -ve and ER +ve patients.

Next, to make the boxplots for the other genes, TECTA, BAX, MAP3K8, GEMIN8, we could copy and paste the code for the FOXE1 plot over and over like below.

## TECTA gene ##

filtered <- exprsAnno %>% 
    filter(HUGO.gene.symbol == "TECTA") %>% 
    select(starts_with("NKI")) %>%
    t() %>%
    data.frame()
colnames(filtered) <- "TECTA"
g <- ggplot(filtered, aes(y = TECTA, x = factor(patients[, "er"]), fill = factor(patients[, "er"]))) +
    geom_boxplot() + 
    scale_fill_manual(values = c("goldenrod", "dodgerblue"),
                      name = "ER status",
                      labels = c("ER -ve", "ER +ve")) +
    labs(title = "Expression of TECTA ~ Estrogen Receptor Status",
         x = "Estrogen Receptor Status",
         y = "Expression values of TECTA") +
    scale_x_discrete(labels = NULL) + 
    theme(plot.title = element_text(hjust = 0.5))
print(g)

## BAX gene ##

filtered <- exprsAnno %>% 
    filter(HUGO.gene.symbol == "BAX") %>% 
    select(starts_with("NKI")) %>%
    t() %>%
    data.frame()
colnames(filtered) <- "BAX"
g <- ggplot(filtered, aes(y = BAX, x = factor(patients[, "er"]), fill = factor(patients[, "er"]))) +
    geom_boxplot() + 
    scale_fill_manual(values = c("goldenrod", "dodgerblue"),
                      name = "ER status",
                      labels = c("ER -ve", "ER +ve")) +
    labs(title = "Expression of BAX ~ Estrogen Receptor Status",
         x = "Estrogen Receptor Status",
         y = "Expression values of BAX") +
    scale_x_discrete(labels = NULL) + 
    theme(plot.title = element_text(hjust = 0.5))
print(g)

## MAP3K8 gene ##

# ETC....

However, as already stated, this is tedious and error-prone.

So instead, lets automate it!

Exercise

Convert the code used to generate the boxplots in the second example of copy/paste/edit into a for loop.

Which variables are you looping over?

Hint: Consider using the function paste() to help change the y-axis label and title to reflect the gene that is being plotted.

Saving multiple plots

You can save all the plots produced by a for loop in a PDF by using the pdf() and dev.off() functions outside the for loop. This will create one file with all the plots.

pdf("myplots.pdf")
for (element in vector) {   
    # plot code
}
dev.off()

Otherwise, another option is to create a PDF (or PNG/JPEG) for each graph inside the for loop.

for (element in vector) {
    pdf(element)
    # plot code   
    dev.off()
}

Storing results

Note that this for loop is helping us identify the number of genes in each chromosome but the result is not stored. Thus, we can not access the results at a later time.

• When storing the results from a loop, we often create an empty variable before starting the for loop
• This is used store the result at each iteration of the loop

numGenes <- NULL
chrom <- c("chr1", "chr5", "chr6", "chr8")
for(chr in chrom) {
  numRows <- exprsAnno %>%
      filter(Chromosome == chr) %>%
      summarise(n())
  numGenes[chr] <- numRows
  }
numGenes

$chr1
[1] 35

$chr5
[1] 16

$chr6
[1] 20

$chr8
[1] 12

Exercise

Identify the number of patients who are positive for ER/PR/HER2 when considering each status individually. Store the results in a vector.

status <- c("er", "pr", "her2")
### Your answer here ###

Conditional branching: Commands and flow control

What if we only wanted to create boxplots for genes if they met certain criteria, for example, if the expression of the gene differed between ER -ve and ER +ve. To do that we could use an if statement.

• Use an if statement for any kind of condition testing
• Different outcomes can be selected based on a condition within brackets

if (condition) {
  ... do this ...
  } else {
    ... do something else ...
    }

• condition is any logical value, and can contain multiple conditions.
  • e.g. (a == 2 & b < 5), this is a compound conditional argument
• The condition should return a single value of TRUE or FALSE

Other conditional tests

• There are various tests that can check the type of data stored in a variable
  • these tend to be called is...().
    • try tab-complete on is.

is.numeric(10)

[1] TRUE

is.numeric("TEN")

[1] FALSE

is.character(10)

[1] FALSE

• is.na() is useful for seeing if an NA value is found
  • cannot use == NA!

x <- c(1, 2, NA)
mean(x)

[1] NA

is.na(x)

[1] FALSE FALSE  TRUE

Example

Checking if a gene of interest is in the dataset - There are 2 methods that could be used to achieve this 1. Make uses of the functions all() and any() 2. Make use of the function sum()

## Method 1
if (any(exprsAnno$HUGO.gene.symbol == "PIK3CA")){
    print("PIK3CA is in the dataset")
}else{
    print("PIK3CA is not in the dataset")
}

[1] "PIK3CA is not in the dataset"

## Method 2
if (sum(exprsAnno$HUGO.gene.symbol == "PIK3CA") == 1){
    print("PIK3CA is in the dataset")
}else{
    print("PIK3CA is not in the dataset")
}

[1] "PIK3CA is not in the dataset"

Exercise

Write an if else statement to check if all of the following genes are in the dataset: BCL2, HOXA9, MAPK1, ARID1A, GATA3, ESR1

Hint: Make use of %in%

Other useful conditionals

file.exists(), dir.exists() and dir.create() are also functions that are useful in if else statements. Particularly, when writing scripts that take inputs from the commandline.

Combining Loops and Conditional branching

Using the for loop we wrote before, we could add an if else branch to identify if the difference between the mean expression for a gene differs between the groups ER positive/negative.

Automation Example 3 - Creating boxplots if expression differs between ER -ve and ER +ve

We are interested in finding out if the following genes:
AMPD3, TECTA, TRPV4, CD244, ABHD10, GEMIN8, MAP1A, SMAD7 have difference between the ER status of greater than 0.04.

• Here’s how we can combine a for loop and an if statement to test for this for each iteration of the loop: 1. Identify the expression values associated with the gene of interest 2. Group the samples based on there ER status, and find the mean for each group 3. if the absolute difference is greater than 0.04, print a statement that informs us of this 4. else, do nothing

mygenes <- c("AMPD3", "TECTA", "TRPV4", "CD244", "ABHD10", "GEMIN8", "MAP1A", "SMAD7")
for (gen in mygenes){
    filtered <- exprsAnno %>%
        filter(HUGO.gene.symbol == gen) %>%
        select(starts_with("NKI")) %>%
        t() %>%
        data.frame()
    colnames(filtered) <- "gene"
    x <- filtered %>% 
        group_by(factor(patients[, "er"])) %>%
        summarise("exprsAv" = mean(gene, na.rm = TRUE))
    diff <- x$exprsAv[1] - x$exprsAv[2]
    if (abs(diff) > 0.04){
        print(paste("The difference between ER status for", gen, "is", abs(diff), "which is greater than 0.04"))
    }else{
        print(paste("The difference between ER status for", gen, "is not greater than 0.04"))
    }
}

[1] "The difference between ER status for AMPD3 is 0.046240625 which is greater than 0.04"
[1] "The difference between ER status for TECTA is 0.046746875 which is greater than 0.04"
[1] "The difference between ER status for TRPV4 is 0.0914234375 which is greater than 0.04"
[1] "The difference between ER status for CD244 is not greater than 0.04"
[1] "The difference between ER status for ABHD10 is 0.0404546875 which is greater than 0.04"
[1] "The difference between ER status for GEMIN8 is 0.0495082677165354 which is greater than 0.04"
[1] "The difference between ER status for MAP1A is not greater than 0.04"
[1] "The difference between ER status for SMAD7 is not greater than 0.04"

Whilst it is slightly informative to know that the difference between means of ER status for each gene. It would be more interesting to plot the expression of the genes that have an absolute difference of greater than 0.04.

Therefore we want to take the same steps as above but change from printing a statement to plotting a graph - Here’s how we can use an if statement to test for this + for each iteration of the the loop: 1. Identify the expression values associated with the gene of interest 2. Group the samples based on there ER status, and find the mean for each group 3. if the absolute difference is greater than 0.04, produce a boxplot to depict the difference 4. else, do nothing

mygenes <- c("AMPD3", "TECTA", "TRPV4", "CD244", "ABHD10", "GEMIN8", "MAP1A", "SMAD7")
for (gen in mygenes){
    filtered <- exprsAnno %>%
        filter(HUGO.gene.symbol == gen) %>%
        select(starts_with("NKI")) %>%
        t() %>%
        data.frame()
    colnames(filtered) <- "gene"
    x <- filtered %>% 
        group_by(factor(patients[, "er"])) %>%
        summarise("exprsAv" = mean(gene, na.rm = TRUE))
    diff <- x$exprsAv[1] - x$exprsAv[2]
    if (abs(diff) > 0.04){
        g <- ggplot(filtered, aes(y = gene, x = factor(patients[, "er"]), fill = factor(patients[, "er"]))) +
            geom_boxplot()
        print(g)
    }else{
        next
    }
}

Exercise

Add additional layers to the above plot so that the labels are more informative. See plot below, replace “x” in the labs for the gene that we are investigating:
Also, change the colours associated with each ER status.

Code formatting avoids bugs!

Compare:

for (a in ages){
numRows <- patients %>%
filter(age > a) %>%
summarise(n())
if(numRows < 10){
s <- paste("There are fewer than 10 patients younger than", a)
print(s)
}else{
s <- paste("There are more than 10 patients younger than", a)
print(s)
}
}

to:

for (a in ages) {
    numRows <- patients %>%
        filter(age > a) %>%
        summarise(n())
    if (numRows < 10) {
        s <- paste("There are fewer than 10 patients younger than", a)
        print(s)
    } else{
        s <- paste("There are more than 10 patients younger than", a)
        print(s)
    }
}

• The code between brackets {} always is indented, this clearly separates what is executed once, and what is run multiple times
• Trailing bracket } always alone on the line at the same indentation level as the initial bracket {
• Use white spaces to divide the horizontal space between units of your code, e.g. around assignments, comparisons