# Microarray pipeline - MA570 

# This practical uses public microarray data from an experiment studying the effect of spaceflight 
# on Mesenchymal Stem Cells (MSCs)

# There are 8 samples in total:
# - 2 MSC samples from the experiment and 2 controls
# - 2 Osteo-induced MSC samples and controls

############################
#                          #
# Step 1: loading the data #
#                          #
############################

# Set the working directory to the folder containing the raw microarray .CEL files
setwd('/home/nextgen/data/microarray/MSC')

# Next, load BioConductor package "affy" 
# This is a library with pre-written functions for microarray analysis

library(affy)

# Use the "ReadAffy" function to read in the raw data

rawData<- ReadAffy()

# This creates an "Affy batch" object which contains the raw data, as well as information about the data
# (Type of microarray, number of genes on the array etc. )
# To see this type "rawData"

rawData

# To see the raw probe values for the data, use

head(exprs(rawData))

# This will show the first few probes.


###########################
#                         #
# Step 2: Quality control #
#                         #
###########################

# Generate a pseudo-image of the microarray chip 

image(rawData)

# Create a heatmap for the first 100 probes on the microarray

# eset<- exprs(rawData)
# heatmap(eset[1:100,])

# A boxplot is a convienient way to visualise the distribution of probe intensities
# This will allow us to compare distributions accross all samples. 
# If the data is ok, the distributions should be relaively comparable between samples

boxplot(exprs(rawData), las=2)

# Because of the distribution of the data, where values are small but some are relatively large,
# It is difficult to properly visualise the data on the boxplot.
# To overcome this, we are going to generate another boxplot,
# But this time we will transform the data to the log2 scale

boxplot(log2(exprs(rawData)), las=2)

# The next step is to conduct a principle components analysis.
# This is a mathematical technique that will allow us to see which samples tend to have similar patterns of expression.
# We expect that samples of the same phenotype should have similar expression patterns,
# and should therefore cluster together on the PCA plot.

pca_proc<- prcomp(t(exprs(rawData)))

# Next, we want to make the PCA plot, but also we want to be able to easily tell samples and phenotypes appart.
# We will do this by using different colors and point characters (pch) on the plot.

# In this experiment we have four phenotypes, Two MSC groups, one for the flight expreiment and the other control,
# and the other Osteo induced MSC, control and flight.

mycol<- c(1,1,2,2,3,3,4,4)
mypch<- seq(1,8)

# Finally, plot the data

plot(pca_proc$x, col=mycol, pch=mypch, xlab='Component 1', ylab='Component 2', main='PCA of raw data')

# And add a legend to the plot,

legend('right', sampleNames(rawData), col=mycol, pch=mypch)


#######################################################################
#                                                                     #
# Step 3: Background correction, normalization and summarize the data #
#                                                                     #
#######################################################################

# The next step is to process the array data using the "Robust Multiarray Average" (RMA) algorithm
# This algorithm carries out three steps:
# 1.) Background correction to filter out the 'noise' on the microarray.
# 2.) Quantile normalization.
# 3.) Summarization using Median Polish

# All this is done using a single function 'rma()'

normData<- rma(rawData)

# We can again check the heatmap, boxplot and PCA of the data, to see what effect this has had on the data

# eset<- exprs(normData)
# heatmap(eset[1:100,])

boxplot(exprs(normData), las=2)

# Note, RMA log transforms the data, so the log2() argument is not used here. 

pca_proc<- prcomp(t(exprs(normData)))
mycol<- c(1,1,2,2,3,3,4,4)
mypch<- seq(1,8)
plot(pca_proc$x, col=mycol, pch=mypch, xlab='Component 1', ylab='Component 2', main='PCA of normalized data')
legend('right', sampleNames(normData), col=mycol, pch=mypch)

##########################################
#                                        #
# Finding differentially expressed genes #
#                                        #
##########################################

# All of the above commands use the 'affy' package, 
# The following commands use the 'limma' package which stands for 'linear models for micrarray data'
# To load the limma package:

library(limma)

# Make a design matrix, which will place each sample into its corresponding groups, we can see which columns belong together using the head() command

head(exprs(normData))

# The column names should indicate which phenotype a sample falls under. 
# These samples represent four phenotypes, MSC flight and control, OI flight and control

design<- model.matrix(~-1+factor(c(1,1,2,2,3,3,4,4)))
colnames(design)<- c('MSC_control', 'MSC_flight', 'OI_control', 'OI_flight')

# Next, create a contrasts matrix, which will tell limma which comparisons we wish to make.
# Here we are going to do two comparrions, MSC_flight vs control, OI_flight vs control

contrasts<- makeContrasts(MSC_flight - MSC_control, OI_flight - OI_control, levels=design)

# Next, fit the data to a linear model

fit<- lmFit(normData, design)
fit2<- contrasts.fit(fit, contrasts)
fit3<- eBayes(fit2)

# Finally make a data frame, summarizing the data, for each of out comparisons

# This function will also carry out correction for multiple testing. 
# In this case it is using the 'Benjamini - Hochberg' method, shown here as 'BH',
# This option can be changed to a different method if desired (such as Bonferroni correction) 

top.MSC<- topTable(fit3, coef='MSC_flight - MSC_control', adjust='BH', number=10000, sort.by='P')
top.OI<- topTable(fit3, coef='OI_flight - OI_control', adjust='BH', number=10000, sort.by='P')

# Have a look at the data frames to see what sort of data it contains.
head(top.MSC)

# We can see that one of the columns contain the adjusted P values,
# We can use this to find statistically significant genes using conditional indexing

# To fing genes that have an adjusted p value < 0.05:

sig.MSC<- top.MSC[top.MSC$adj.P.Val<0.05,]
sig.OI<- top.OI[top.OI$adj.P.Val<0.05,]

# See how many genes we have using the dim() function

dim(sig.MSC)
dim(sig.OI)

# Write out the results to a file, for later analysis usige DAVID
# For DAVID, we only need the gene id (shown in the ID column)

MSC.ids<- as.vector(sig.MSC$ID)
OI.ids<- as.vector(sig.OI$ID)

write(MSC.ids, file='MSC.significant.txt')
write(OI.ids, file='OI.significant.txt')