###################################################
### chunk number 1: 
###################################################
library(GEOquery)


###################################################
### chunk number 2: 
###################################################
gds <- getGEO('GDS10')


###################################################
### chunk number 3: 
###################################################
gsm <- getGEO('GSM3')


###################################################
### chunk number 4: 
###################################################
# Look at gsm metadata:
Meta(gsm)
# Look at data associated with the GSM:
# but restrict to only first 5 rows, for brevity
Table(gsm)[1:5,]
# Look at Column descriptions:
Columns(gsm)


###################################################
### chunk number 5: 
###################################################
Columns(gds)


###################################################
### chunk number 6: 
###################################################
gse <- getGEO("GSE462",GSEMatrix=FALSE)
Meta(gse)
# names of all the GSM objects contained in the GSE
names(GSMList(gse))
# and get the first GSM object on the list
GSMList(gse)[[1]]
# and the names of the GPLs represented
names(GPLList(gse))


###################################################
### chunk number 7: 
###################################################
gse2553 <- getGEO('GSE2553',GSEMatrix=TRUE)
show(gse2553)
show(pData(phenoData(gse2553[[1]]))[1:5,c(1,6,8)])


###################################################
### chunk number 8: 
###################################################
eset <- GDS2eSet(gds,do.log2=TRUE)


###################################################
### chunk number 9: 
###################################################
eset
pData(eset)


###################################################
### chunk number 10: 
###################################################
#get the platform from the GDS metadata
Meta(gds)$platform
#So use this information in a call to getGEO
gpl <- getGEO("GPL5")


###################################################
### chunk number 11: 
###################################################
MA <- GDS2MA(gds,GPL=gpl)
MA


###################################################
### chunk number 12: 
###################################################
gsmplatforms <- lapply(GSMList(gse),function(x) {Meta(x)$platform})
gsmplatforms


###################################################
### chunk number 13: 
###################################################
Table(GSMList(gse)[[1]])[1:5,]
# and get the column descriptions
Columns(GSMList(gse)[[1]])[1:5,]


###################################################
### chunk number 14: 
###################################################
# get the probeset ordering
probesets <- Table(GPLList(gse)[[1]])$ID
# make the data matrix from the VALUE columns from each GSM
# being careful to match the order of the probesets in the platform
# with those in the GSMs
data.matrix <- do.call('cbind',lapply(GSMList(gse),function(x) 
                                      {tab <- Table(x)
                                       mymatch <- match(probesets,tab$ID_REF)
                                       return(tab$VALUE[mymatch])
                                     }))
data.matrix <- apply(data.matrix,2,function(x) {as.numeric(as.character(x))})
data.matrix <- log2(data.matrix)
data.matrix[1:5,]


###################################################
### chunk number 15: 
###################################################
require(Biobase)
# go through the necessary steps to make a compliant ExpressionSet
rownames(data.matrix) <- probesets
colnames(data.matrix) <- names(GSMList(gse))
pdata <- data.frame(samples=names(GSMList(gse)))
rownames(pdata) <- names(GSMList(gse))
pheno <- as(pdata,"AnnotatedDataFrame")
eset2 <- new('ExpressionSet',exprs=data.matrix,phenoData=pheno)
eset2


###################################################
### chunk number 16: 
###################################################
toLatex(sessionInfo())