##################
## LOAD THE DATA #
##################

## Set the "working directory" to the directory which contains our ".CEL" files.
setwd("/home/nextgen/CEL_files")

## Load the "affy" library. This is an R library which provides
convenient pre-written functions that allow
## us to easily analyse Affymetrix microarray data.
library("affy")

## Read the raw microarray data. The ReadAffy() function is provided
by the "affy" library that we have loaded above.
## By default ReadAffy() will attempt to load all ".cel" files in the
current working directory.
## We will store the loaded data in an object called "rawData".
rawData <- ReadAffy()

## "rawData" is an object of class "affyBatch". This object contains
the information from the ".cel" files
## in a structured format that allows us to conveniently query and
process the data quickly.
rawData

## the following lines will print a matrix containing intensity levels
for the individual probes on the microarray
## rows of this matrix correspond to the individual probes, and
columns to the different samples
exprs(rawData)

## As you can see, the matrix is too big to print, the following will print the
## dimensions of the matrix (the number of probes in each array x the
number of arrays).
## The 2nd line will print the contents of the top corner (first 5
rows and columns) of the matrix
dim(exprs(rawData))
exprs(rawData)[1:5, 1:5]



###################################
## ASSESS THE QUALITY OF THE DATA #
###################################

## A boxplot provides a convenient means by which to visualize the
distribution of the data in these matrices
## and allows us to easily compare expression levels between the 12
different samples.
## If the data is okay, expression intensities should be roughly
comparable between samples
boxplot(exprs(rawData))

## Because of the distribution of the data (most values small, but
some very large), it is extremely difficult to read the previous
boxplot.
## We can "log transform" the data to make it easier to compare the boxplots
## We also add the "las=2" parameter to the "boxplot()" function,
which rotates the lables on the X-axis so we can see them properly
## typing "?boxplot()" will allow you to investigate other possible
parameters to a function.
boxplot(log2(exprs(rawData)), las=2)

## The previous boxplot shows us that the samples "N7S.CEL" tends to
have higher expression intensities than the other samples.
## This may be problematic.
## We can further investigate if there may be a problem with "N7S.CEL"
using "Principal Components Analysis" (PCA)
## which is a mathematical technique that allows us to visualize which
samples tend to have similar patterns of expression.
## We would expect that samples of the same phenotype will have more
similar expression profiles and hence
## cluster together on the PCA plot.
pca_proc <- prcomp(t(exprs(rawData))) #preform "Principal Components
Anaysis" calculations on this matrix and return data as object of
class "prcomp"

## Plot the results of PCA
## Here we introduce some more complex plotting functionality, whereby
we have use different colors for arrays of each of our two phenotypes
## i.e. P1S.CEL ... P6S.CEL are in red while N9S.CEL ... N12S.CEL are in black

## We have also used the "legend()" function to indicate which samples
are which.
mycol <- c(rep(1, 6), rep(2,6))
mypch <- seq(1,12)
plot(pca_proc$x, xlab="Component 1", ylab="Component 2", col=mycol,
pch=mypch, main="PCA of Preprocessed Expression set" )
legend("right", sampleNames(rawData), pch = mypch, col = mycol)


#########################################################
## BACKGROUND CORRECT, NORMALIZE AND SUMMARIZE THE DATA #
#########################################################

## But we may be able to correct this using "normalization", which
will attempt to adjust for scaling differences between samples
## We will use a method called "Robust Multi-array Average" (RMA).
## The rma() function in fact performs background correction,
summarization (of individual probe level expression estimates into
gene level expression estimates),
## log2 transforms the data and quantiles normalizes the data in 1 easy step.
normalizedData <- rma(rawData)

## The object normalizedData (of class "expressionSet") is an object
of normalized data for the 12 microarrays
normalizedData
exprs(normalizedData)[1:5, 1:5] #expression matrix row names are now
gene names (as Affymetrix gene identifiers)
dim(exprs(normalizedData)) #there are now much less rows in the
expression matrix, as probe intensities have been summarized to gene
level estimates

## View a boxplot of the normlized data (as you can see this seems to
corrects the scaling issues for sample "N7S.CEL")
boxplot(exprs(normalizedData), las=2)

## We can check if the array now clusteres on the PCA plot
pca_proc_norm <- prcomp(t(exprs(normalizedData))) #preform "Principal
Components Anaysis" calculations on this matrix and return data as
object of class "prcomp"
plot(pca_proc_norm$x, xlab="Component 1", ylab="Component 2",
col=mycol, pch=mypch, main="PCA of Preprocessed Expression set" )
legend("right", sampleNames(normalizedData), pch = mypch, col = mycol)