###################################################
### chunk number 1: startup
###################################################
library(Biostrings)
library(hgu95av2probe)
library(hgu95av2cdf)


###################################################
### chunk number 2: matchprobes
###################################################
pm <- DNAStringSet(hgu95av2probe$sequence)
dict <- pm[3801:4000]
pdict <- PDict(dict)
m <- vcountPDict(pdict, pm)
dim(m)
table(rowSums(m))
which(rowSums(m) == 3)
ii <- which(m[77, ] != 0)
pm[ii]


###################################################
### chunk number 3: basecontent
###################################################
bcpm <- alphabetFrequency(pm, baseOnly=TRUE)
head(bcpm)
alphabetFrequency(pm, baseOnly=TRUE, collapse=TRUE)


###################################################
### chunk number 4: hgu95av2dim
###################################################
nr = hgu95av2dim$NROW
nc = hgu95av2dim$NCOL


###################################################
### chunk number 5: 
###################################################
library(affy)
abseq = rep(as.character(NA), nr*nc) 
ipm = with(hgu95av2probe, xy2indices(x, y, nr=nr))
any(duplicated(ipm))  # just a sanity check
abseq[ipm] = hgu95av2probe$sequence
table(is.na(abseq))


###################################################
### chunk number 6: pm2mm
###################################################
pm2mm <- function(probes)
{
    parts <- threebands(DNAStringSet(probes), start=13, end=13)
    xscat(parts$left, complement(parts$middle), parts$right)
}
mm <- pm2mm(pm)
cat(as.character(pm[[1]]), as.character(mm[[1]]), sep="\n")


###################################################
### chunk number 7: mismatchSeq
###################################################
imm = with(hgu95av2probe, xy2indices(x, y+1, nr=nr))
intersect(ipm, imm)  # just a sanity check
abseq[imm] = as.character(mm)
table(is.na(abseq))


###################################################
### chunk number 8: bc
###################################################
freqs <- alphabetFrequency(DNAStringSet(abseq[!is.na(abseq)]), baseOnly=TRUE)
bc <- matrix(nrow=length(abseq), ncol=5)
colnames(bc) <- colnames(freqs)
bc[!is.na(abseq), ] <- freqs
head(na.omit(bc))


###################################################
### chunk number 9: GC
###################################################
GC = ordered(bc[,"G"] + bc[,"C"])
colores = rainbow(nlevels(GC))


###################################################
### chunk number 10: abatch
###################################################
library(affydata)
f <- system.file("extracelfiles", "CL2001032020AA.cel", package="affydata")
pd <- new("AnnotatedDataFrame", data=data.frame(fromFile=I(f), row.names="f"))
abatch <- read.affybatch(filenames=f, compress=TRUE, phenoData=pd)


###################################################
### chunk number 11: bap
###################################################
barplot(table(GC), col=colores, xlab="GC", ylab="number")


###################################################
### chunk number 12: bxp
###################################################
boxplot(log2(exprs(abatch)[,1]) ~ GC, outline=FALSE,
  col=colores, , xlab="GC", ylab=expression(log[2]~intensity))


###################################################
### chunk number 13: p2p
###################################################
png("matchprobes-p2p.png", width=900, height=480)
plot(exprs(abatch)[ipm,1], exprs(abatch)[imm,1], asp=1, pch=".", log="xy",
     xlab="PM", ylab="MM", col=colores[GC[ipm]])
abline(a=0, b=1, col="#404040", lty=3)
dev.off()