# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Code chunk #01
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
library(aroma.affymetrix)
cdf <- AffymetrixCdfFile$fromChipType("Mapping50K_Hind240")
print(cdf)

getPathname(cdf)
getPath(cdf)
getFilename(cdf)
getFileSize(cdf)



# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Code chunk #02
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
library(aroma.affymetrix)
cdf <- AffymetrixCdfFile$fromChipType("Mapping50K_Hind240")
nbrOfUnits(cdf)
nbrOfCells(cdf)
nbrOfRows(cdf)
nbrOfColumns(cdf)
AffymetrixCdfFile


unitNames <- getUnitNames(cdf)
unitNames <- getUnitNames(cdf, units=101:120)
unit <- indexOf(cdf, "SNP_A-1683825")

data <- readUnits(cdf, units=121)
names(data)
names(data[[1]]$groups)
idxs <- getCellIndices(cdf, units=121)


# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Code chunk #03
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
cs <- AffymetrixCelSet$byName("HapMap100K-bioc2007", 
                                 chipType="Mapping50K_Hind240")
print(cs)

getName(cs)
getPath(cs)
getNames(cs)
getFullNames(cs)
getFileSize(cs)

cdf <- getCdf(cs)
AffymetrixCelSet



# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Code chunk #04
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
as.list(cs)
cf <- getFile(cs, 1)
print(cf)


# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Code chunk #05
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
library(aroma.affymetrix)
cs <- AffymetrixCelSet$byName("HapMap100K-bioc2007", 
                              chipType="Mapping50K_Hind240")
plotDensity(cs, stratifyBy="pm")  # ~2s/array



# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Code chunk #06                          [Requires: EBImage]
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
library(aroma.affymetrix)
cs <- AffymetrixCelSet$byName("HapMap100K-bioc2007", 
                              chipType="Mapping50K_Hind240")
cf <- getFile(cs, 1)
print(cf)
plotImage(cf)


# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Code chunk #07                          [Requires: EBImage]
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
library(aroma.affymetrix)
cs <- AffymetrixCelSet$byName("HapMap100K-bioc2007", 
                              chipType="Mapping50K_Hind240")
ae <- ArrayExplorer(cs)
setColorMaps(ae, "sqrt,yellow")
process(ae, verbose=-3)
display(ae)



# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Code chunk #08
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
library(aroma.affymetrix)
cs <- AffymetrixCelSet$byName("HapMap100K-bioc2007", 
                              chipType="Mapping50K_Hind240")
qn <- QuantileNormalization(cs)
print(qn)

csQN <- process(qn, verbose=TRUE)

plotDensity(csQN, stratifyBy="pm")  # ~2s/array

gc()



# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Code chunk #09
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
quit("no")

# Restart R
library(aroma.affymetrix)
cs <- AffymetrixCelSet$byName("HapMap100K-bioc2007", 
                              chipType="Mapping50K_Hind240")
qn <- QuantileNormalization(cs)
csQN <- process(qn, verbose=TRUE)
print(csQN)




# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Code chunk #10
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
library(aroma.affymetrix)
cs <- AffymetrixCelSet$byName("HapMap100K-bioc2007", 
                              chipType="Mapping50K_Hind240")
qn <- QuantileNormalization(cs)
csQN <- process(qn, verbose=TRUE)

plm <- RmaCnPlm(csQN, combineAlleles=TRUE, mergeStrands=TRUE)

ces <- getChipEffectSet(plm)
theta <- readUnits(ces, units=10)
print(theta)

fit(plm, units=10)
theta <- readUnits(ces, units=10)
print(theta)



# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Code chunk #11
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
fit(plm, units=1:1000, verbose=TRUE)
fit(plm, units=1:1000, verbose=TRUE)



# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Code chunk #12
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
library(aroma.affymetrix)
cs <- AffymetrixCelSet$byName("HapMap100K-bioc2007", 
                              chipType="Mapping50K_Hind240")
qn <- QuantileNormalization(cs)
csQN <- process(qn, verbose=TRUE)
plm <- RmaCnPlm(csQN, combineAlleles=TRUE, mergeStrands=TRUE)
fit(plm, verbose=TRUE)




# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Code chunk #13
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
cdf <- AffymetrixCdfFile$fromChipType("Mapping50K_Hind240")

# dChip genome information files
gi <- getGenomeInformation(cdf)
print(gi)
units <- getUnitsOnChromosome(gi, 16)
x <- getPositions(gi, units)
plot(density(x/1e6, from=0, adjust=0.2))


# dChip SNP information files
si <- getSnpInformation(cdf)
print(si)



# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Code chunk #14
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
library(aroma.affymetrix)
cs <- AffymetrixCelSet$byName("HapMap100K-bioc2007", 
                              chipType="Mapping50K_Hind240")
qn <- QuantileNormalization(cs)
csQN <- process(qn, verbose=TRUE)
ces <- getChipEffectSet(plm)

theta <- extractMatrix(ces)
thetaAvg <- rowMedians(theta)
M <- log2(theta/thetaAvg)

cdf <- getCdf(ces)
si <- getSnpInformation(cdf)
fl <- getFragmentLengths(si)

layout(matrix(1:6, ncol=3, byrow=TRUE))
for (kk in seq(ces)) {
  smoothScatter(fl, M[,kk], ylim=c(-1.5,1.5), main=kk)
  lines(smooth.spline(fl, M[,kk]), col="red", lwd=3)
}


# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Code chunk #15
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
fln <- FragmentLengthNormalization(ces)
cesN <- process(fln, verbose=TRUE)

theta <- extractMatrix(cesN)
thetaAvg <- rowMedians(theta)
M <- log2(theta/thetaAvg)

layout(matrix(1:6, ncol=3, byrow=TRUE))
for (kk in seq(ces)) {
  smoothScatter(fl, M[,kk], ylim=c(-1.5,1.5), main=kk)
  lines(smooth.spline(fl, M[,kk]), col="orange", lwd=3)
}




# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Code chunk #16
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
glad <- GladModel(cesN)
fit(glad, arrays=1, chromosomes=19, verbose=TRUE)
plot(glad, arrays=1, chromosomes=19)

plot(glad, arrays=1, chromosomes=18)



# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Code chunk #17
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
glad <- GladModel(cesN)
ce <- ChromosomeExplorer(glad)
print(ce)
process(ce, arrays=1, chromosomes=1:4, verbose=TRUE)
display(ce)


# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Code chunk #18
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
library(aroma.affymetrix)
cs <- AffymetrixCelSet$byName("HapMap100K-bioc2007", 
                              chipType="Mapping50K_Hind240")
qn <- QuantileNormalization(cs)
csQN <- process(qn, verbose=TRUE)
plm <- RmaCnPlm(csQN, combineAlleles=TRUE, mergeStrands=TRUE)
fit(plm, verbose=TRUE)
ces <- getChipEffectSet(plm)
fln <- FragmentLengthNormalization(ces)
cesN <- process(fln, verbose=TRUE)
glad <- GladModel(cesN)
ce <- ChromosomeExplorer(glad)
process(ce, verbose=-3)




# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Code chunk #19
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
library(aroma.affymetrix)
chipTypes <- c("Mapping50K_Hind240", "Mapping50K_Xba240")
cesList <- list()
for (chipType in chipTypes) {
  cs <- AffymetrixCelSet$byName("HapMap100K-bioc2007", chipType=chipType)
  print(cs)
  qn <- QuantileNormalization(cs)
  csQN <- process(qn, verbose=TRUE)
  plm <- RmaCnPlm(csQN, combineAlleles=TRUE, mergeStrands=TRUE)
  fit(plm, verbose=TRUE)
  ces <- getChipEffectSet(plm)
  fln <- FragmentLengthNormalization(ces)
  cesN <- process(fln, verbose=TRUE)
  cesList[[chipType]] <- cesN
}
glad <- GladModel(cesList)
ce <- ChromosomeExplorer(glad)
process(ce, verbose=-3)


# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Code chunk #20                          [Requires: EBImage]
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
library(aroma.affymetrix)
cs <- AffymetrixCelSet$byName("HapMap100K-bioc2007", 
                              chipType="Mapping50K_Hind240")
qn <- QuantileNormalization(cs)
csQN <- process(qn, verbose=TRUE)
plm <- RmaPlm(csQN)
fit(plm, verbose=TRUE)

rs <- calculateResidualSet(plm, verbose=TRUE)  # ~5 sec/array

rf <- getFile(rs, 1)
plotImage(rf)

ae <- ArrayExplorer(rs)
setColorMaps(ae, "log2,log2abs,rainbow")
process(ae, verbose=-3)