About R

If you haven’t learned the basics of R prior to attending this course, you should check out our R crash course for an overview of R’s syntax. It’s also a great refresher if you feel it has been a while since you last worked with R.

About the practicals for this workshop

• The traditional way to enter R commands is via opening a Terminal or, or using the console in RStudio (bottom-left panel when RStudio opens for first time).
• For this course we will instead be using a relatively new feature called R Notebooks.
• An R notebook mixes plain text written in markdown with “chunks” of R code.

Markdown is a very simple way of writing a template to produce a pdf, HTML or word document. For example, the compiled version of this document is available online and is more convenient to browse here.

• “Chunks” of R code can be added using the insert option from the tool bar, or the CTRL + ALT + I shortcut
• Each line of R can be executed by clicking on the line and pressing CTRL and ENTER
• Or you can execute the whole chunk by pressing CTRL + SHIFT + ENTER
• Or you can press the green triangle on the right-hand side to run everything in the chunk
• The code might have different options which dictate how the output is displayed in the compiled document (e.g. HTML)
  • e.g. you might see EVAL = FALSE or echo = FALSE
  • you don’t have to worry about this if stepping through the markdown interactively

print("Hello World")

## [1] "Hello World"

When viewing the R notebooks directly, not the compiled documents, sections may have additional characters so as to format them nicely when compiled. This is markdown, a simple language that formats the text whle, crucially, being still readable in the raw state.

For example:

This will be displayed in italic

This will be displayed in bold

• this
• is
• a
• list
  • this is a sub-list

You can also add hyperlinks, images and tables.

Lastly, you can even embed chunks of code written in other programming languages.

a='Wow python'
print(a.split())

## ['Wow', 'python']

More help is available through RStudio Help -> Markdown Quick Reference or you can view a cheat sheet here.

To create markdown files for your own analysis; File -> New File -> R Markdown…

About the Bioconductor project

Established in 2001, Bioconductor provided a convenient method to distribute tools for the analysis and comprehension of high-throughput genomic data in R. Initially focused on microarrays, Bioconductor now has packages (read: software) to process data obtained from most modern data sources.

• R is rarely used for the primary processing of modern data
  • R is far slower than many other programming languages due to it being an interpreted language (Interpreted vs Compiled)
  • R is extensively-used for visualisation, interpretation and inference once data has been parsed into a more manageable form, e.g., a csv.

On the Bioconductor website, you will find

• Installation instructions
• Course Materials
• Support forum
  • a means of communicating with developers and power-users
  • Stack Overflow tends to get faster/better responses
• Example datasets
• Annotation Resources
• Conferences

For this session, we will introduce the Bioconductor project as a means of analysing high-throughput data

Installing a package

All Bioconductor software packages are listed under

• bioconductor.org -> Install -> Packages -> Analysis software packages
  • Many thousands of packages have been added over the years, so I would suggest just googling “bioconductor [package_name]”
  • e.g. edgeR landing page
• installation instructions are given, which involves running the install command from the BiocManager package
  • this will install and update any additional dependencies
  • Running BiocManager::install forgoes needing to load the package manager into the work enviroment
• you only need to run this procedure once for each version of R

## You don't need to run this, edgeR should already be installed for the course
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("edgeR")

Once installed, a Bioconductor package can be loaded in the usual way with the library function. All packages are required to have a vignette which gives detailed instructions on how to use the package and the workflow of commands. Some packages such as edgeR have very comprehensive user guides with lots of use-cases.

library(edgeR)
vignette("edgeR")
edgeRUsersGuide()

Package documentation can also be accessed via the Help tab in RStudio, which can also be invoked in the console using “?”

?edgeR

Structures for data analysis

Complex data structures are used in Bioconductor to represent high-throughput data, but we often have simple functions that we can use to access the data. We will use some example data available via Bioconductor to demonstrate how high-throughput data can be represented, and also to review some basic concepts in data manipulation in R.

• This data is from a microarray experiment. We will be concentrating on more modern technologies in this class, but most of the R techniques required will be eaxctly the same or at least very similar.
• experimental data packages are available through Bioconductor, and can be installed in the way we just described
  • the package should already be installed on your computer, so you won’t need to run this.

## No need to run this - for reference only!
BiocManager::install("breastCancerVDX")

To make the dataset accessible in R, we first need to load the package. If we navigate to the documentation for breastCancerVDX in RStudio, we find that it provides an object called vdx which we load into R’s memory using the data function.

library(breastCancerVDX)
data(vdx)

The object vdx is a representation of breast cancer dataset that has been converted for use with standard Bioconductor tools. The package authors don’t envisage that we will want to view the entire dataset at once, so have provided a number of ways to interact with the data

• typing the name of the object provides a summary, e.g.,
  • how many genes in the dataset
  • how many samples

vdx

## ExpressionSet (storageMode: lockedEnvironment)
## assayData: 22283 features, 344 samples 
##   element names: exprs 
## protocolData: none
## phenoData
##   sampleNames: VDX_3 VDX_5 ... VDX_2038 (344 total)
##   varLabels: samplename dataset ... e.os (21 total)
##   varMetadata: labelDescription
## featureData
##   featureNames: 1007_s_at 1053_at ... AFFX-TrpnX-M_at (22283
##     total)
##   fvarLabels: probe Gene.title ... GO.Component.1 (22 total)
##   fvarMetadata: labelDescription
## experimentData: use 'experimentData(object)'
##   pubMedIds: 17420468 
## Annotation: hgu133a

Accessing expression values

The expression values can be obtained by the exprs function:

• remember, <- is used for assignment to create a new variable
• the data are stored in a matrix in R
  • it is a good idea to check the dimensions using dim, ncol, nrow etc.

eValues <- exprs(vdx) # also found at vdx@assayData$exprs
class(eValues)

## [1] "matrix"

dim(eValues)

## [1] 22283   344

ncol(eValues)

## [1] 344

nrow(eValues)

## [1] 22283

• the row names are the manufacturer-assigned ID for a particular probe
• the column names are the identifiers for each patient in the study
• each entry is a normalised log\(_2\) intensity value for a particular gene in a given sample
  • we won’t talk about normalisation here, but basically the data has been transformed so that samples and/or genes can be compared
• subsetting a matrix is done using the [row, column] notation
  • the function c is used to make a one-dimensional vector
  • the shortcut : can used to stand for a sequence of consecutive numbers

eValues[c(1,2,3),c(1,2,3,4)]

##               VDX_3     VDX_5     VDX_6     VDX_7
## 1007_s_at 11.965135 11.798593 11.777625 11.538577
## 1053_at    7.895424  7.885696  7.949535  7.481396
## 117_at     8.259272  7.052025  8.225930  8.382408

eValues[1:3,1:4]

##               VDX_3     VDX_5     VDX_6     VDX_7
## 1007_s_at 11.965135 11.798593 11.777625 11.538577
## 1053_at    7.895424  7.885696  7.949535  7.481396
## 117_at     8.259272  7.052025  8.225930  8.382408

• subsetting can be chained together
• we can also omit certain rows or columns from the output by prefixing the indices with a -

eValues[1:3,1:4][,2:3]

##               VDX_5     VDX_6
## 1007_s_at 11.798593 11.777625
## 1053_at    7.885696  7.949535
## 117_at     7.052025  8.225930

eValues[1:3,1:4][,-(2:3)]

##               VDX_3     VDX_7
## 1007_s_at 11.965135 11.538577
## 1053_at    7.895424  7.481396
## 117_at     8.259272  8.382408

Simple visualisations

The most basic plotting function in R is plot

• using the plot function with a vector will plot the values of that vector against the index
  • what do you think is displayed in the plot below?

plot(eValues[1,])

• one possible use is to compare the values in a vector with respect to a given factor
• but we don’t know the clinical variables in our dataset yet (to come later)
• a boxplot can also accept a matrix or data frame as an argument
  • what do you think the following plot shows?

boxplot(eValues,outline=FALSE)

Accessing patient data

The metadata, or phenotypic data, for the samples is retrieved using the pData function.

metadata <- pData(vdx) #vdx@phenoData@data
head(metadata)

• head prints only the first few rows/values.
  • you can specificy the exact number of rows to output, e.g., head(metadata,10)
• individual columns can be accessed using the $ notation.
• columns are returned as a vector, which can be fed into other standard plotting and analysis functions.
• auto-complete is available in RStudio; once you type the $, it should give you a list of possible options.

head(metadata$samplename,15)

##  [1] "VDX_3"  "VDX_5"  "VDX_6"  "VDX_7"  "VDX_8"  "VDX_9"  "VDX_11"
##  [8] "VDX_14" "VDX_15" "VDX_17" "VDX_18" "VDX_19" "VDX_20" "VDX_21"
## [15] "VDX_22"

#FILL IN SOLUTIONS HERE

table(metadata$grade == 3)

## 
## FALSE  TRUE 
##    49   148

table(metadata$grade == 3,useNA="ifany")

## 
## FALSE  TRUE  <NA> 
##    49   148   147

metadata[which(metadata$grade == 3),]

head(colnames(eValues))
head(rownames(metadata))
table(colnames(eValues) == rownames(metadata))

eValues[,which(metadata$grade==3)][1:10,1:10]

##               VDX_5     VDX_6     VDX_7     VDX_9    VDX_11    VDX_14
## 1007_s_at 11.798593 11.777625 11.538577 12.248698 12.259243 11.695620
## 1053_at    7.885696  7.949535  7.481396  7.708049  8.182891  7.970394
## 117_at     7.052025  8.225930  8.382408  7.731319  8.226412  7.710118
## 121_at    10.666845 10.888819 10.795472 10.886687 10.192909 11.275077
## 1255_g_at  5.452859  6.456149  6.147714  6.551516  6.376777  6.620586
## 1294_at    9.585901  9.422906  9.456970  9.831782  9.404290  9.292552
## 1316_at    7.644577  7.065012  8.121534  7.573647  7.068241  8.377644
## 1320_at    5.066089  4.661065  6.400879  4.153805  3.678072  4.209453
## 1405_i_at 10.043711  9.350276 11.520619 13.012869  9.896030  8.618386
## 1431_at    6.230741  6.263034  5.951868  6.183883  5.741467  6.797013
##              VDX_15    VDX_17    VDX_18    VDX_20
## 1007_s_at 11.577193 12.483665 11.899130 12.414157
## 1053_at    6.721099  8.005625  6.815063  7.067165
## 117_at     7.375908  7.488644  7.209453  9.775939
## 121_at    10.614526 10.664625 10.772892 10.623607
## 1255_g_at  6.044394  6.244126  5.311067  6.266787
## 1294_at    8.572511 10.301839  9.430871  8.914684
## 1316_at    7.302867  7.780048  8.203103  7.651769
## 1320_at    4.887525  6.456149  5.385431  6.230741
## 1405_i_at  8.355792  9.875289  7.210428  9.133656
## 1431_at    5.744161  6.407693  6.560715  5.617651

vdx[,which(metadata$grade==3)]

## ExpressionSet (storageMode: lockedEnvironment)
## assayData: 22283 features, 148 samples 
##   element names: exprs 
## protocolData: none
## phenoData
##   sampleNames: VDX_5 VDX_6 ... VDX_913 (148 total)
##   varLabels: samplename dataset ... e.os (21 total)
##   varMetadata: labelDescription
## featureData
##   featureNames: 1007_s_at 1053_at ... AFFX-TrpnX-M_at (22283
##     total)
##   fvarLabels: probe Gene.title ... GO.Component.1 (22 total)
##   fvarMetadata: labelDescription
## experimentData: use 'experimentData(object)'
##   pubMedIds: 17420468 
## Annotation: hgu133a

fac <- metadata$er
table(fac)

## fac
##   0   1 
## 135 209

boxplot(eValues[1,] ~ fac,
        xlab="ER Status",
        ylab="Expression Level",
        col=c("steelblue","goldenrod"))

t.test(eValues[1,]~fac)

## 
##  Welch Two Sample t-test
## 
## data:  eValues[1, ] by fac
## t = -1.7476, df = 244.15, p-value = 0.08179
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.2101311  0.0125555
## sample estimates:
## mean in group 0 mean in group 1 
##        11.72515        11.82394

Accessing feature (gene) information

We could be lucky, and the first row in the expression matrix could be our favourite gene of interest! However, this is unlikely to be the case and we will have to figure out which row we want to plot

• we can use another aspect of the vdx object; the feature data
• there is a handy fData function to access these data
• again, this gives us a data frame
• this is a pre-built table supplied with the dataset
  • later in the course we will look at using online services and databases for annotation and converting between identifiers

features <- fData(vdx) #try to figure out the location in memory of this data.frame via @ and $
class(features)

## [1] "data.frame"

head(features[,1:5])

colnames(features)

##  [1] "probe"                 "Gene.title"           
##  [3] "Gene.symbol"           "Gene.ID"              
##  [5] "EntrezGene.ID"         "UniGene.title"        
##  [7] "UniGene.symbol"        "UniGene.ID"           
##  [9] "Nucleotide.Title"      "GI"                   
## [11] "GenBank.Accession"     "Platform_CLONEID"     
## [13] "Platform_ORF"          "Platform_SPOTID"      
## [15] "Chromosome.location"   "Chromosome.annotation"
## [17] "GO.Function"           "GO.Process"           
## [19] "GO.Component"          "GO.Function.1"        
## [21] "GO.Process.1"          "GO.Component.1"

As we know, gene symbols (more-common gene names) can be accessed using the $ syntax; returning a vector

head(features$Gene.symbol)

## [1] "DDR1"   "RFC2"   "HSPA6"  "PAX8"   "GUCA1A" "UBA7"

We could inspect the data frame manually (e.g. using View(features)) and identify the row number corresponding to our gene of interest. However, as aspiring R programmers, there is a better way

• == to test for exact matching
• match will return the index of the first match
• grep can be used for partial matches
• each of the above will give an vector that can be used to subset the expression values
  • Why do we get different results from the various commands?

which(features$Gene.symbol == "BRCA1")

## [1]  4058 11245

match("BRCA1",features$Gene.symbol)

## [1] 4058

grep("BRCA1",features$Gene.symbol)

## [1]  4058 11245

grep("BRCA",features$Gene.symbol)

## [1]  4058  7869 11245 14103

#FILL IN SOLUTIONS HERE

Testing all genes for significance

Later in the course we will see how to execute a differential expression analysis for RNA-seq data, and discuss some of the issues surrounding this. For now we will perform a simple two-sample t-test for all genes in our study, and derive a results table

• firstly, we load the genefilter package which has a very convenient function for performing many t tests in parallel
• rowttests will run a t-test for each row in a given matrix and produce an output table
  • statistic; test statistic
  • dm; difference in means
  • p.value; the p-value
• rowttests expects a factor as the second argument, so we have to use as.factor
  • as usual, we can get help by doing ?rowttests

library(genefilter)
tstats <- rowttests(eValues, as.factor(metadata$er))
head(tstats)

hist(tstats$statistic)

The rows are ordered in the same way as the input matrix

• to change this to increasing significance we can use the order function
• when given a vector, order will return a vector of the same length giving the permutation that rearranges that vector into ascending or descending order
• The sort function shortcuts using the output of order to rearrange the original vector by returning the sorted vector

x <- c(9, 3, 4, 2, 1, 6,5, 10, 8, 7)
x

##  [1]  9  3  4  2  1  6  5 10  8  7

order(x)

##  [1]  5  4  2  3  7  6 10  9  1  8

x[order(x)]

##  [1]  1  2  3  4  5  6  7  8  9 10

sort(x)

##  [1]  1  2  3  4  5  6  7  8  9 10

• so if we want to order by p-value we first use order on the p-value vector
• this can then be used to re-arrange the rows of the table

head(tstats[order(tstats$p.value,decreasing = FALSE),])

However, the table we get isn’t immediately useful unless we can recognise the manufacturer probe IDs

• to provide extra annotation to the table, we can column bind (cbind) the columns of test statistic with those from the feature matrix
• be careful though, we can only do this in cases where the rows are in direct correspondence

table(rownames(tstats) == rownames(features))

## 
##  TRUE 
## 22283

tstats.annotated <- cbind(features[,c("Gene.symbol","EntrezGene.ID","Chromosome.location")], tstats)
head(tstats.annotated)

Now when we order by p-value, the extra columns that we just added allow us to interpret the results more easily

tstats.ordered <- tstats.annotated[order(tstats$p.value,decreasing = FALSE),]
head(tstats.ordered)

We can also query this table to look for our favourite genes of interest

• %in% is a simplified way to perform matches to multiple items in a vector

tstats.ordered[grep("ESR1",tstats.ordered$Gene.symbol),]

tstats.ordered[tstats.ordered$Gene.symbol %in% c("ESR1","GATA3","FOXA1"),]

#FILL IN SOLUTIONS HERE

Using your own data

The object vdx, which contains the data we’ve been using, is an ExpressionSet object (check using the class function). ExpressionSet’s are used as input in a variety of bioinformatic packages and pipelines. As we found above, an ExpressionSet usually contains the following data:

• expression data
• meta data (optional)
• feature data (optional)

Let’s generate some data and create our own ExpressionSet from scratch - rnorm generates normally distrubruted values - matrix takes a vector and arranges it into the specified number of rows and columns - rep is used to repeat a vector or its elements - sample samples the input vector, with or without replacement - AnnotatedDataFrame is a special container that conveniently stores and manipulates data - Remeber that you can use help (?) to understand unknown functions

eData <- matrix(rnorm(100000),nrow = 1000,ncol = 100)
mData <- AnnotatedDataFrame(data.frame(Sex=rep(c("M","F"),each=50),biomarker.X=rnorm(100)))
fData <- AnnotatedDataFrame(data.frame(Gene=1:1000,Chr=sample(1:23,1000,replace = T)))

exampleSet <-  ExpressionSet(
  assayData=eData,
  phenoData=mData,
  featureData=fData
)

#FILL IN SOLUTIONS HERE