inst/testScripts/system/chipTypes/Mapping10K_Xba142/22.doASCRMAv2,extract.R

In aroma.affymetrix: Analysis of Large Affymetrix Microarray Data Sets

library("aroma.affymetrix")
verbose <- Arguments$getVerbose(-4, timestamp=TRUE)

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Setup
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
dataSet <- "GSE8605"
chipType <- "Mapping10K_Xba142"
csR <- AffymetrixCelSet$byName(dataSet, chipType=chipType, verbose=verbose)
print(csR)


# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# AS-CRMAv2
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
res <- doASCRMAv2(csR, drop=FALSE, verbose=verbose)
print(res)


# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Extraction test
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
ces <- res$cesN

# (a) Single array
ce <- ces[[1]]
Y <- extractMatrix(ce, verbose=verbose)
print(head(Y))

theta <- extractTheta(ce, verbose=verbose)
print(head(theta))

data <- extractDataFrame(ce, addNames=TRUE, verbose=verbose)
print(head(data))


# (b) Multiple arrays
Y <- extractMatrix(ces, verbose=verbose)
print(head(Y))

theta <- extractTheta(ces, verbose=verbose)
print(theta[1:6,,])

data <- extractDataFrame(ces, addNames=TRUE, verbose=verbose)
print(head(data))

# Extract unit and group name and indices together with chip effects
dataT <- extractDataFrame(ces, units=1:50, addNames=TRUE)
print(head(dataT))
# Sanity check
stopifnot(all.equal(dataT[1:50,], data[1:50,]))

# SNPs only
unf <- getUnitNamesFile(ces)
print(unf)

units <- indexOf(unf, pattern="^SNP")
str(units)

theta <- extractTheta(ces, units=units[1:10])
print(theta)


# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Plot theta:s
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
theta <- extractTheta(ces, units=units)
dimnames(theta)[[2]] <- c("A", "B")
ltheta <- log2(theta)
a <- (ltheta[,"A",]+ltheta[,"B",])/2
m <- ltheta[,"A",]-ltheta[,"B",]

toPNG(getFullName(ces), tags=c("LogDiffvsLogSum"), aspectRatio=0.7, {
  Alim <- c(5,14)
  Mlim <- diff(Alim)*c(-1,1)
  Alab <- expression(1/2%*%log[2](theta[A]*theta[B]))
  Mlab <- expression(log[2](theta[A]/theta[B]))

  I <- dim(theta)[3]
  layout(matrix(1:I, ncol=3, byrow=TRUE))
  par(mar=c(3.8,4,3,1)+0.1)

  for (ii in 1:I) {
    name <- dimnames(theta)[[3]][ii]
    plot(a[,ii], m[,ii], pch=".", xlim=Alim, ylim=Mlim, xlab=Alab, ylab=Mlab, main=name)
  }
})