inst/testScripts/system/chipTypes/HuEx-1_0-st-v2/21.doFIRMA,QC.R

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

library("aroma.affymetrix")
library("matrixStats"); # rowMedians()
verbose <- Arguments$getVerbose(-3, timestamp=TRUE)


# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Setup data set
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
dataSet <- "Affymetrix-HeartBrain"
chipType <- "HuEx-1_0-st-v2"
cdf <- AffymetrixCdfFile$byChipType(chipType, tags="coreR3,A20071112,EP")
print(cdf)

# Setup CEL set using the core CDF.
csR <- AffymetrixCelSet$byName(dataSet, cdf=cdf)


# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Process only cerebellum and heart
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
types <- c("cerebellum", "heart")
csR <- csR[indexOf(csR, patterns=types)]

setFullName(csR, sprintf("%s,%s", dataSet, paste(types, collapse="+")))
print(csR)


# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# FIRMA
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
res <- doFIRMA(csR, drop=FALSE, verbose=verbose)
print(res)


# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# PLMs
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
plmList <- list(
  merge   = res$plm,                               # all exons together
  noMerge = ExonRmaPlm(res$csN, mergeGroups=FALSE) # each exon separately
)
print(plmList)

# Fit the per-exon PLM
fit(plmList$noMerge, verbose=verbose)

# Chip effects
cesList <- lapply(plmList, FUN=getChipEffectSet)
print(cesList)


# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# RLE/NUSE for first 100 units
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# The reason for observing a small difference is scores below,
# is that when saving to file the estimates are rounded of to
# floats, whereas the estimates calculated in memory are kept
# in full precision.
tol <- 1e-5

units <- 1:100

for (key in names(cesList)) {
  verbose && enter(verbose, key)

  ces <- cesList[[key]]

  # Assert correctness of unit-specific RLE scores
  theta <- extractMatrix(ces, field="theta", units=units)
  thetaR <- 2^rowMedians(log2(theta), na.rm=TRUE)
  rle0 <- log2(theta/thetaR)
  rle1 <- extractMatrix(ces, field="RLE", units=units, verbose=verbose)
  stopifnot(all.equal(rle1, rle0, tolerance=tol))

  # Assert correctness of boxplot statistics of RLE scores
  stats0 <- boxplot.stats(rle0[,1])
  stats1 <- boxplotStats(ces, type="RLE", arrays=1:2, subset=units)
  stopifnot(all.equal(stats1[[1]], stats0, tolerance=tol))

  # Assert correctness of unit-specific NUSE scores
  se <- extractMatrix(ces, field="sdTheta", units=units)
  seR <- 2^rowMedians(log2(se), na.rm=TRUE)
  nuse0 <- log2(se)/log2(seR)
  nuse1 <- extractMatrix(ces, field="NUSE", units=units, verbose=verbose)
  stopifnot(all.equal(nuse1, nuse0, tolerance=tol))

  # Assert correctness of boxplot statistics of NUSE scores
  stats0 <- boxplot.stats(nuse0[,1])
  stats1 <- boxplotStats(ces, type="NUSE", arrays=1:2, subset=units)
  stopifnot(all.equal(stats1[[1]], stats0, tolerance=tol))

  verbose && exit(verbose)
} # for (key in ...)