Introduction to Solving Biological Problems Using R - Day 2

4. Programming in R

Motivation

From the previous exercise, you should see how we can easily adapt our markdown scripts:

• e.g. ESR1 versus GATA3
• But what if we want to analyse many genes?
• It would be tedious to create a new markdown document for every gene
• …and prone to error too

Introducing loops

• Many programming languages have ways of doing the same thing many times, perhaps changing some variable each time. This is called looping
• Loops are not used in R so often, because we can usually achieve the same thing using vector calculations
• For example, to add two vectors together, we do not need to add each pair of elements one by one, we can just add the vectors

x <- 1:10
y <- 11:20
x+y

• But there are some situations where R functions can not take vectors as input. For example, t.test() will only test one gene at a time
• What if we wanted to test multiple genes?

For completeness, we can re-run the R code to import the data

geneAnnotation    <- read.delim("gene.description.txt",stringsAsFactors = FALSE)
patientMetadata <- read.delim("cancer.patients.txt",stringsAsFactors = FALSE)
normalizedValues    <- read.delim("gene.expression.txt")

• We could run the following code to perform t-tests on the first two genes

t.test(as.numeric(normalizedValues[1,]) ~ factor(patientMetadata$er))
t.test(as.numeric(normalizedValues[2,]) ~ factor(patientMetadata$er))

• But for many genes this will be boring to type, difficult to change, and prone to error
• As we are doing the same thing multiple times, but with a different index each time, we can use a loop instead

Loops: Commands and flow control

• R has two basic types of loop
  • a for loop: run some code on every value in a vector
  • a while loop: run some code while some condition is true (hardly ever used!)

for

for(i in 1:10) {
  print(i)
  }

while

i <- 1
while(i <= 10 ) {
  print(i)
  i <- i + 1
  }

• Here’s how we might use a for loop to test the first 10 genes

for(i in 1:10) {
  
  t.test(as.numeric(normalizedValues[i,]) ~ factor(patientMetadata$er))
  
  }

• This is exactly the same as:

i <- 1
t.test(as.numeric(normalizedValues[i,]) ~ factor(patientMetadata$er))
i <- 2
t.test(as.numeric(normalizedValues[i,]) ~ factor(patientMetadata$er))
i <- 3
###....etc....####

Storing results

However, this for loop is doing the calculations but not storing the results

• The output of t.test() is an object with data placed in different slots
  • the names() of the object tells us what data we can retrieve, and what variable name to use

testResult <- t.test(as.numeric(normalizedValues[1,]) ~ factor(patientMetadata$er))
names(testResult)
testResult$statistic

• When using a loop, we often create an empty “dummy” variable
• This is used store the results at each stage of the loop

stats <- NULL
for(i in 1:10) {
  testResult <- t.test(as.numeric(normalizedValues[i,]) ~ factor(patientMetadata$er))
  stats[i] <- testResult$statistic
  }
stats

Practical application

Previously we have identified probes on chromosome 8

• Lets say that we want to do a t-test for each gene on chromosome 8

chr8Genes <- geneAnnotation[geneAnnotation$Chromosome=="chr8",]
head(chr8Genes)
chr8GenesOrd <- chr8Genes[order(chr8Genes$Start),]
head(chr8GenesOrd)

• The first step is to extract the expression values for chromosome 8 genes from our expression matrix, which has expression values for all genes
• We can use the match function to tell us which rows in the matrix correspond to chromosome 8 genes

match(chr8GenesOrd$probe, rownames(normalizedValues))
chr8Expression <- normalizedValues[match(chr8GenesOrd$probe, rownames(normalizedValues)),]
dim(chr8Expression)

We are now ready to write the for loop

Exercise:

• Create a for loop to perform to test if the expression level of each gene on chromosome 8 is significantly different between ER positive and negative samples
• Store the p-value from each individual test
• How many genes have a p-value < 0.05?
• N.B. Our code will be more robust if we store the number of chromosome 8 genes as a variable
  • if the data change, the code should still run

chr8Genes <- geneAnnotation[geneAnnotation$Chromosome=="chr8",]
chr8GenesOrd <- chr8Genes[order(chr8Genes$Start),]
chr8Expression <- normalizedValues[match(chr8GenesOrd$probe, rownames(normalizedValues)),]
### Your Answer Here ###

Conditional branching: Commands and flow control

• 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)
is.numeric("TEN")
is.character(10)

• is.na() is useful for seeing if an NA value is found
  • cannot use == NA!
• could be used to check if a gene symbol is found in the data before proceeding with statistical test

match("foo", geneAnnotation$HUGO.gene.symbol)
is.na(match("foo", geneAnnotation$HUGO.gene.symbol))

• Using the for loop we wrote before, we could add some code to plot the expression of each gene
  • a boxplot would be ideal
• However, we might only want plots for genes with a “significant” pvalue
• Here’s how we can use an if statement to test for this
  • for each iteration of the the loop:
    1. test if the p-value from the test is below 0.05 or not
    2. if the p-value is less than 0.05 make a boxplot
    3. if not, do nothing

pdf("Chromosome8Genes.pdf")
pvals <- NULL
for (i in 1:18) {
  testResult <- t.test(as.numeric(chr8Expression[i,]) ~ factor(patientMetadata$er))
  pvals[i] <- testResult$p.value
  if(testResult$p.value < 0.05){
    boxplot(as.numeric(chr8Expression[i,]) ~ factor(patientMetadata$er),
            main=chr8Genes$HUGO.gene.symbol[i])
  }
} 
pvals
dev.off()

Code formatting avoids bugs!

Compare:

f<-26
while(f!=0){
print(letters[f])
f<-f-1}

to:

f <- 26
while(f != 0 ){
  print(letters[f])
  f <- f-1
  }

• 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

Making a heatmap

• A heatmap is often used to visualise how the expression level of a set of genes vary between conditions
• Making the plot is actually quite straightforward
  • providing you have processed the data appropriately!
• Let’s take a list of “most-variable genes”
  • see below for how we identified such genes
• heatmap requires a matrix object rather than a data frame

genelist <- c("CLIC6","TFF3","PDZK1","SCUBE2","CYP2B6","HOXB13","NAT1","LY6D","SLC7A2")
probes   <- geneAnnotation$probe[match(genelist, geneAnnotation$HUGO.gene.symbol)]
probes

[1] "Contig46937_RC" "NM_003226"      "NM_002614"      "NM_020974"      "M29873"         "Contig28549"   
[7] "NM_000662"      "X82693"         "NM_003046"     

exprows  <- match(probes, rownames(normalizedValues))
heatmap(as.matrix(normalizedValues[exprows,]),labCol="")

Heatmap adjustments

• We can provide a colour legend for the samples
• Adjust colour of cells
• Label the rows according to gene name
• Don’t print the sample names as they are too cluttered

library(RColorBrewer)
sampcol <- rep("blue", ncol(normalizedValues))
sampcol[patientMetadata$er == 1 ] <- "yellow"
rbPal <- brewer.pal(10, "RdBu")
heatmap(as.matrix(normalizedValues[exprows,]), 
        ColSideColors = sampcol, 
        col=rbPal,
        labRow = genelist,labCol="")

• see also
  • heatmap.2 from library(gplots); example(heatmap.2)
  • heatmap.plus from library(heatmap.plus); example(heatmap.plus)

(Supplementary) Choosing the genes for the heatmap

Often when using R you will come across a convenient shortcut function that can save you many hours of coding and frustration.

• the genefilter package in Bioconductor contains many useful methods for filtering genomic datasets
• you can install this package with the following commands

source("http://www.bioconductor.org/biocLite.R")
biocLite("genefilter")

• the rowSds function will calculate the standard deviation for each row in a numeric matrix
• the output will be vector, with each element being the standard deviation for a corresponding gene

library(genefilter)
geneVar <- rowSds(normalizedValues)
geneVar[1]
sd(normalizedValues[1,])

• we can now order this matrix to get the subset with [] to get the indices of the most-variable genes (10 in this case).
• the same indices can be used to subset the gene annotation data frame
  • we can do this because the annotation data frame and expression matrix are in the same order

topVar <- order(geneVar,decreasing = TRUE)[1:10]
topVar
geneAnnotation[topVar,]