# Example code for processing of Affymetrix data:

# install R on your computer, and put all the .CEL files into a directory on your computer.
# then run the following code in R:

source("http://www.bioconductor.org/biocLite.R")
biocLite()
biocLite("rat2302.db")
biocLite("siggenes")
# no longer downloading and installing rat2302cdf - this is automatic with : Data=ReadAffy() below

memory.limit(size=4095)

library(affy)
library(limma)
library(siggenes)
library(MASS)
library(rat2302.db)
rat2302()

# this must indicate the directory ON YOUR COMPUTER where the .CEL files are stored:
setwd("C:/Documents and Settings/Carlo/Desktop/2008.2")

Data=ReadAffy()
Data
colnames(exprs(Data))
dim(exprs(Data))
dim(pm(Data))

Groups=c(rep("OLD",4),rep("YNG",4))

cbind(colnames(exprs(Data)),Groups)

#########################
# normalization

RMA=rma(Data)

#########################
# annotation

AffyIDs=rownames(exprs(RMA))
GeneNames=unlist(mget(AffyIDs,rat2302GENENAME))
GeneSymbols=unlist(mget(AffyIDs,rat2302SYMBOL))

class(GeneNames)
length(GeneNames)
GeneNames[1:5]

class(GeneSymbols)
length(GeneSymbols)
GeneSymbols[1:5]

AffyIDanno=cbind(AffyIDs,GeneSymbols,GeneNames)
class(AffyIDanno)
AffyIDanno[1:5,]

#########################
# RNA degradation plot

deg <- AffyRNAdeg(Data)
jpeg(file="RNAdegradPlot.jpg",width=500,height=500)
plotAffyRNAdeg(deg)
dev.off()

#########################
# LogRatio vs. Intensity plots

jpeg(file="IntVSratio_Scatter_RAW.jpg",width=1500,height=1500)
MAplot(Data, pairs = TRUE)
dev.off()

###############
# RAW log2 intensitiy density plots

jpeg(file="IntDensRAWlog2.jpg",width=750,height=750)
plot(density(log2(pm(Data))),type="n",ylim=c(0,.45))
for(i in 1:dim(pm(Data))[2])
{
lines(density(log2(pm(Data)[,i])),col=i,lty=i)
}
ArrayIndex=c(1:dim(pm(Data))[2])
legend(x="topright", as.character(ArrayIndex), lt = 1:length(ArrayIndex), col = 1:length(ArrayIndex), cex = 0.75)
dev.off()

###############
# RMA intensitiy density plots

jpeg(file="IntDensRMA.jpg",width=750,height=750)
plot(density(exprs(RMA)),type="n")
for(i in 1:dim(exprs(RMA))[2])
{
lines(density(exprs(RMA)[,i]),col=i,lty=i)
}
ArrayIndex=c(1:dim(exprs(RMA))[2])
legend(x="topright", as.character(ArrayIndex), lt = 1:length(ArrayIndex), col = 1:length(ArrayIndex), cex = 0.75)
dev.off()

###############
# PCA

PCA<-prcomp(t(exprs(RMA)))

jpeg(file="PCAColorByGroup.jpg",width=750,height=750)
plot(PCA$x,main="PCA of Affy Experiments\nColored by Array Group: Red=YNG, BLK=OLD",xlab=paste("PC1: ",signif((((PCA$sdev)^2)[1]/(sum((PCA$sdev)^2))),2),sep=""),ylab=paste("PC2: ",signif((((PCA$sdev)^2)[2]/(sum((PCA$sdev)^2))),2),sep=""))
GroupsUni=unique(Groups)
for(i in 1:length(GroupsUni)){points(PCA$x[Groups==GroupsUni[i],],pch=20,col=i,cex=5)}
text(PCA$x,labels=1:dim(exprs(RMA))[2],col="white",font=2)
dev.off()

###############
# SAM

Groups
SAMclasses=as.numeric(Groups=="OLD")
cbind(Groups,SAMclasses)

SAM.OLD.v.YNG=sam(exprs(RMA),cl=SAMclasses,method="d.stat",q.version=1,gene.names=rownames(exprs(RMA)))
str(SAM.OLD.v.YNG)
SAM.OLD.v.YNG

indx.hi=SAM.OLD.v.YNG@d==max(SAM.OLD.v.YNG@d)
names(SAM.OLD.v.YNG@d)[indx.hi]

indx.lo=SAM.OLD.v.YNG@d==min(SAM.OLD.v.YNG@d)
names(SAM.OLD.v.YNG@d)[indx.lo]

findDelta(SAM.OLD.v.YNG,fdr=0.02)

jpeg(file="SAM.jpg",width=750,height=750)
plot(SAM.OLD.v.YNG,1.821499)
dev.off()

jpeg(file="SAMdensity.jpg",width=750,height=750)
plot(density(SAM.OLD.v.YNG@d),lwd=3,xlim=c(-6,6),ylim=c(0,.5),main="Distribution of SAM d-stats: Observed (solid), and Permuted (dashed)")
lines(density(SAM.OLD.v.YNG@d.bar),lty=2,lwd=3)
dev.off()

###############
# Check the differential expression of a single group

# grab a group of interest

indx=grep("Gst",GeneSymbols)
length(indx)
cbind(AffyIDanno[indx,],SAM.OLD.v.YNG@d[indx],SAM.OLD.v.YNG@q.value[indx])

# geneSetTest

GSTp=geneSetTest(indx,SAM.OLD.v.YNG@d)
GSTp

jpeg(file="GstGroupDensity.jpg",width=750,height=750)
plot(density(SAM.OLD.v.YNG@d),lwd=3,xlim=c(-8,8),ylim=c(0,.3),main="Distribution of SAM d-stats:\nAll Genes (solid), and Gst Genes (red dashed)")
lines(density(SAM.OLD.v.YNG@d[indx]),lty=2,col=2,lwd=3)
dev.off()

#########################################################
# save

save.image(file="ExampleAffyData.RData")