library(Biobase)		# the Bioconductor base package
library(affy)			# for dealing with Affymetrix arrays
library(limma)			# for analysis of differential expression
library(hgu133plus2)	# reference data for the affy array
library(annaffy)		# affymetrix annotations
library(GOstats)		# GO statistics


#------------------------------------------------------------------
# Loading data
# let's load some data into an affybatch
# -----------------------------------------------------------------

data.pd	<- read.phenoData("phenodata.txt")
data 	<- ReadAffy("LN_1.CEL", "LN_2.CEL", "C4_1.CEL", "C4_2.CEL",
					phenoData=data.pd)


# -----------------------------------------------------------------
# QC section
# we're not actually going to do this due to memory limits
# -----------------------------------------------------------------

# Pset <- fitPLM(data)	# do the QC fitting
# image(Pset, which=2)	# draw the chip pseudoimage


# -----------------------------------------------------------------
# Background correction and normalization
# ...this will take a while...
# -----------------------------------------------------------------

data.ex <- expresso(data, 
					bgcorrect.method="rma",
					normalize.method="quantiles",
					pmcorrect.method="pmonly",
					summary.method="avgdiff")


# -----------------------------------------------------------------
# Analysis
# what kinds of tests can we run???
# 	* differential expression tests
#	* clustering tests
#	* GO categorization
# -----------------------------------------------------------------

# differential expression
design 	<- cbind(WT=c(1,1,0,0), MU=c(0,0,1,1))	# experimental design
fit		<- eBayes(lmFit(data.ex, design))		# do the fitting
t50		<- topTable(fit, number=50)				# find top 50 genes

# get english descriptions of top 50 genes
t50.descs <- aafDescription(t50$ID, "hgu133plus2")

# clustering
cutoff 		<- t50$P.Value[50]		# highest cutoff in the top 50
map.genes 	<- data.ex[which(fit$p.value[,2] <= cutoff),] 
heatmap(exprs(map.genes))

# GO analysis of Biological Process (BP)
ll.ids <- unique(as.numeric(aafLocusLink(t50$ID, "hgu133plus2")))
go.analysis <- GOHyperG(ll.ids, lib="hgu133plus2", what="BP")
go.names <- names(go.analysis$p.values)