### Load the appropriate library for reading affy data, and inspect the data.
library( 'affy' )
soy.ab <- ReadAffy( 'geo_data/GSM209576.CEL.gz',
                   'geo_data/GSM209585.CEL.gz',
                   'geo_data/GSM209594.CEL.gz',
                   'geo_data/GSM209577.CEL.gz',
                   'geo_data/GSM209586.CEL.gz',
                   'geo_data/GSM209595.CEL.gz',

                   ## we have gz files which R can read in.
                   compress=TRUE)



## Inspect the loaded data. This will make a network connection first time
## around and can be slow.
soy.ab

## Check out the names of the samples.
sampleNames( soy.ab )



## as the current sample names refer to the original files, we change this for
## something more, er, easy to remember.
new.sampleNames <- c('hr.a3.12','hr.b3.12','hr.c3.12',
                     'ts.a4.12','ts.b4.12','ts.c4.12')
sampleNames(soy.ab) <- new.sampleNames

## and check that it has worked
sampleNames( soy.ab )



##
## We are trying to do some subsetting because not all of the probes on the
## chip are from soy


## read in another data frame called Species.Affy.ID.
## this links species names to affy ids.
Species.Affy.ID <- read.table('SpeciesAffyID.txt', header = T, sep = "")
dim(Species.Affy.ID)



load( 'SoybeanCutObjects.RData' )

tv.for.glycine.max <- Species.Affy.ID$species == 'Glycine max'
table( tv.for.glycine.max )
listOutProbeSets <- Species.Affy.ID$affyID[ tv.for.glycine.max==FALSE ]

length( listOutProbeSets )
is.factor( listOutProbeSets )

## Create a character vector for listOutProbeSets
## One way: rename listOutProbeSets as a character vector
listOutProbeSets <- as.character(listOutProbeSets)

## Confirm that listOutProbeSets is a character vector
is.character(listOutProbeSets)

## check object
soy.ab


## this is the bit which actually removes the stuff we are not intereste
RemoveProbes(listOutProbes=NULL, listOutProbeSets, cdfpackagename, probepackagename)

## Check that the object has less IDs now. There should be 37444.
soy.ab

# Start preparation for phenoData slot in AffyBatch object
pd <- data.frame(population = c(1,1,1,2,2,2), replicate = c(1,2,3,1,2,3))

# Display contents of pd
pd

# Assign the sampleNames(soy.ab) to the rownames of pd
rownames(pd) <- sampleNames(soy.ab)

# Display contents of pd again, notice change in rownames
pd

## Continue preparation for phenoData slot
metaData <- data.frame(labelDescription = c( 'population', 'replicate' ))

## Establish new phenoData slot
phenoData(soy.ab) <- new( 'AnnotatedDataFrame', data = pd, varMetadata = metaData)

## Display pData(soy.ab)
pData(soy.ab)

## Display phenoData(soy.ab)
phenoData(soy.ab)

palette.gray <- c(rep(gray(0:10/10), times = seq(1,41, by = 4)))

library('RColorBrewer')
brewer.cols <- brewer.pal(6, 'Set1')


## Now we want to look at the images of the chips. Basically, we are looking
## for something unusual here; one chip too bright for example.

## Construct image plots for quality control assessment
## Note: that these plots appear dark; they are not log transformed
par(oma = c(1,1,3,1))
par(mfrow = c(2,3))


palette.gray <- c(rep(gray(0:10/10), times = seq(1,41, by = 4)))
image(soy.ab[,1], transfo = function(x) x, col = palette.gray)
image(soy.ab[,2], transfo = function(x) x, col = palette.gray)
image(soy.ab[,3], transfo = function(x) x, col = palette.gray)
image(soy.ab[,4], transfo = function(x) x, col = palette.gray)
image(soy.ab[,5], transfo = function(x) x, col = palette.gray)
image(soy.ab[,6], transfo = function(x) x, col = palette.gray)
mtext ("> image(soy.ab[,i=1:6], transfo = function(x) x, col = palette.gray)", side = 3, outer = T, cex = .8)


## Construct individual image plots using log intensities
## Note: that these plots appear lighter; they are log transformed
##
## This allows us to see a single microarray chip.
par(mfrow = c(1,1))

palette.gray <- c(rep(gray(0:10/10), times = seq(1,41, by = 4)))
image(soy.ab[,1],  col = palette.gray)
title(sub = "> image(soy.ab[,1],  col = palette.gray)")

## Repeat above commands with .. soy.ab[,i], i = 2-6

## Display all six log-transformed image plots in one plot
par(oma = c(3,1,3,1))
par(mfrow = c(2,3))
image(soy.ab[,1],  col = palette.gray)
image(soy.ab[,2],  col = palette.gray)
image(soy.ab[,3],  col = palette.gray)
image(soy.ab[,4],  col = palette.gray)
image(soy.ab[,5],  col = palette.gray)
image(soy.ab[,6],  col = palette.gray)
mtext("Image Plots - Hawaii/Resistant vs. Taiwan/Susceptible", side = 3, outer = T)
mtext("Top Row: Hawaii/Resistant - Bottom Row: Taiwan/Susceptible", side = 1, outer = T)

## Reset to single screen.
par(mfrow = c(1,1))