TESTS: Add tests from examples

NEWS

@@ -2,6 +2,7 @@ Package: aroma.apd
 ==================
 
 Version: 0.6.0 [2015-02-24]
+o ROBUSTNESS: Added package tests (from examples).
 o ROBUSTNESS: Importing suggested functions using ::.
 o ROBUSTNESS: Declaring all S3 methods.
 o Bumped package dependencies.

tests/celToApd.R

@@ -0,0 +1,30 @@
+library("aroma.apd")
+library("R.utils")  ## Arguments
+
+# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+# 1. Scan for existing CEL files
+# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+# a) Scan for CEL files
+files <- list.files(pattern="[.](cel|CEL)$")
+files <- files[!file.info(files)$isdir]
+if (length(files) > 0 && require("affxparser")) {
+  cat("Create an optimal read map for CEL file:", files[1], "\n")
+  cdffile <- findCdf(readCelHeader(files[1])$chiptype)
+  res <- cdfToApdMap(cdffile)
+
+  cat("Converting CEL file to APD file:", files[1], "\n")
+  apdfile <- celToApd(files[1])
+  cat("Created APD file:", apdfile, "\n")
+  file.remove(apdfile)
+
+  cat("Converting CEL file to APD file with an optimized read map:", files[1], "\n")
+  apdfile <- celToApd(files[1], mapType=res$mapType)
+  cat("Created APD file:", apdfile, "\n")
+
+  writeMap <- invertMap(res$readMap)
+  for (file in files[-1]) {
+    cat("Converting CEL file to APD file with an optimized read map:", file, "\n")
+    apdfile <- celToApd(file, mapType=res$mapType, writeMap=writeMap)
+    cat("Created APD file:", apdfile, "\n")
+  }
+} # if (length(files) > 0)

tests/createApd.R

@@ -0,0 +1,78 @@
+library("aroma.apd")
+
+# Float precision
+.Machine$float.eps <- (2^((8-4)*8)*.Machine$double.eps)
+tol <- .Machine$float.eps ^ 0.5
+
+
+# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+# 1. Create an Affymetrix Probe Signal (APD) file for a
+#    'Mapping50K_Hind240' with 1600-by-1600 probes.
+# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+chipType <- "Mapping50K_Hind240"
+nbrOfCells <- 1600^2
+
+pathname <- paste(tempfile(), "apd", sep=".")
+createApd(pathname, nbrOfCells=nbrOfCells, chipType=chipType)
+
+# File size
+cat("File name:", pathname, "\n")
+cat("File size:", file.info(pathname)$size, "bytes\n")
+cat("APD header:\n")
+header <- readApdHeader(pathname)
+print(header)
+
+
+# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+# 2. Update the signals for a subset of probes
+# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+cells <- c(1, 1100:1120, 220:201, 998300:999302)
+signals <- log(cells+1, base=2)  # Fake signals
+updateApd(pathname, indices=cells, data=signals)
+
+
+# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+# 3. Re-read the signals for a subset of probes
+# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+apd <- readApd(pathname, indices=cells)
+
+# Signals in APD files are stored as floats (since this is
+# the precision in CEL files).
+stopifnot(all.equal(signals, apd$intensities, tolerance=tol))
+
+
+# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+# 4. Re-read the signals for a subset of probes
+# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+if (require("affxparser")) {
+  cdfFile <- findCdf(chipType)
+  if (length(cdfFile) > 0) {
+
+    apd <- readApdUnits(pathname, units=100:104)
+
+    # Sample new data (with one decimal precision)
+    apd2 <- lapply(apd, function(unit) {
+      lapply(unit, function(groups) {
+        n <- length(groups$intensities)
+        values <- as.integer(runif(n, max=655350))/10
+        list(intensities=values)
+      })
+    })
+
+    # Update APD file with new data
+    updateApdUnits(pathname, units=100:104, data=apd2)
+
+    # Re-read data to verify correctness
+    apd <- readApdUnits(pathname, units=100:104)
+
+    # Signals in APD files are stored as floats (since this is
+    # the precision in CEL files).
+    stopifnot(all.equal(apd2, apd, tolerance=tol))
+  } # if (length(cdfFile) > 0 ...)
+} # if (require("affxparser"))
+
+
+# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+# 4. Clean up
+# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+file.remove(pathname)

tests/gtypeCelToPQ.R

@@ -0,0 +1,12 @@
+library("aroma.apd")
+
+# Scan for CEL files
+files <- list.files(pattern="[.](cel|CEL)$")
+
+# Convert each to RAW file
+for (file in files) {
+  rawFile <- gsub("[.][^.]*$", ".raw", file)
+  file.remove(rawFile)
+  cel <- gtypeCelToPQ(file, verbose=TRUE)
+  write.table(file=rawFile, cel, sep="\t", quote=FALSE)
+}

tests/readApdRectangle.R

@@ -0,0 +1,24 @@
+library("aroma.apd")
+
+# Local functions
+rotate270 <- function(x, ...) {
+  x <- t(x)
+  nc <- ncol(x)
+  if (nc < 2) return(x)
+  x[,nc:1,drop=FALSE]
+}
+
+# Scan current directory for APD files
+files <- list.files(pattern="[.](apd|APD)$")
+files <- files[!file.info(files)$isdir]
+if (length(files) > 0) {
+  apdFile <- files[1]
+
+  # Read APD intensities in the upper left corner
+  apd <- readApdRectangle(apdFile, xrange=c(0,250), yrange=c(0,250))
+  z <- rotate270(apd$intensities)
+  sub <- paste("Chip type:", apd$header$chipType)
+  image(z, col=gray.colors(256), axes=FALSE, main=apdFile, sub=sub)
+  text(x=0, y=1, labels="(0,0)", adj=c(0,-0.7), cex=0.8, xpd=TRUE)
+  text(x=1, y=0, labels="(250,250)", adj=c(1,1.2), cex=0.8, xpd=TRUE)
+}

tests/readApdUnits.R

@@ -0,0 +1,54 @@
+library("aroma.apd")
+library("R.utils") # Arguments
+
+verbose <- Arguments$getVerbose(TRUE)
+
+# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+# 1. Scan for existing CEL files
+# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+# a) Scan current directory for CEL files
+files <- list.files(pattern="[.](cel|CEL)$")
+files <- files[!file.info(files)$isdir]
+
+if (length(files) > 0 && require("affxparser")) {
+  # b) Corresponding APD filenames
+  celNames <- files
+  apdNames <- gsub(pattern, ".apd", files)
+
+  # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+  # 1. Copy the probe intensities from a CEL to an APD file
+  # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+  for (kk in 1) {
+    verbose && enter(verbose, "Reading CEL file #", kk)
+    cel <- readCel(celNames[kk])
+    verbose && exit(verbose)
+
+    if (!file.exists(apdNames[kk])) {
+      verbose && enter(verbose, "Creating APD file #", kk)
+      chipType <- cel$header$chiptype
+      writeApd(apdNames[kk], data=cel$intensities, chipType=chipType)
+      verbose && exit(verbose)
+    }
+
+    verbose && enter(verbose, "Verifying APD file #", kk)
+    apd <- readApd(apdNames[kk])
+    verbose && exit(verbose)
+    stopifnot(identical(apd$intensities, cel$intensities))
+  }
+
+
+  # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+  # 2. Read a subset of the units
+  # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+  units <- c(1, 20:205)
+  cel <- readCelUnits(celNames[1], units=units)
+  apd <- readApdUnits(apdNames[1], units=units)
+  stopifnot(identical(apd, cel))
+
+
+  # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+  # 3. The same, but stratified on PMs and MMs
+  # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+  apd <- readApdUnits(apdNames[1], units=units, stratifyBy="pmmm",
+                                                addDimnames=TRUE)
+} # if (length(files) > 0)

tests/writeApdMap.R

@@ -0,0 +1,98 @@
+library("aroma.apd")
+library("R.utils") # Arguments
+
+verbose <- Arguments$getVerbose(TRUE)
+
+# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+# 1. Scan for existing CEL files
+# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+# a) Scan current directory for CEL files
+files <- list.files(pattern="[.](cel|CEL)$")
+files <- files[!file.info(files)$isdir]
+if (length(files) > 0 && require("affxparser")) {
+  # b) Corresponding APD filenames
+  celNames <- files
+  apdNames <- gsub(pattern, ".apd", files)
+
+
+  # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+  # 1. Copy the probe intensities from a CEL to an APD file
+  # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+  for (kk in 1) {
+    verbose && enter(verbose, "Reading CEL file #", kk)
+    cel <- readCel(celNames[kk])
+    verbose && exit(verbose)
+
+    chipType <- cel$header$chiptype
+    verbose && enter(verbose, "Getting read map for '", chipType, "'")
+    mapType <- chipType
+    mapFile <- paste(mapType, ".apm", sep="")
+    if (file.exists(mapFile)) {
+      verbose && enter(verbose, "Reading read map from APD map file")
+      readMap <- readApdMap(mapFile)$map
+      verbose && exit(verbose)
+    } else {
+      verbose && enter(verbose, "Generating read map from CDF file")
+      cdfFile <- findCdf(chipType)
+      readMap <- readCdfUnitsMap(cdfFile)
+      writeApdMap(mapFile, map=readMap)
+      verbose && exit(verbose)
+    }
+    verbose && exit(verbose)
+
+    if (!file.exists(apdNames[kk])) {
+      verbose && enter(verbose, "Calculating write map from read map")
+      writeMap <- invertMap(readMap)
+      verbose && exit(verbose)
+
+      verbose && enter(verbose, "Creating APD file #", kk)
+      writeApd(apdNames[kk], data=cel$intensities, chipType=chipType,
+                                     mapType=mapType, writeMap=writeMap)
+      verbose && exit(verbose)
+      rm(writeMap)
+    }
+
+    verbose && enter(verbose, "Verifying APD file #", kk)
+    apd <- readApd(apdNames[kk], readMap=readMap)
+    verbose && exit(verbose)
+    stopifnot(identical(apd$intensities, cel$intensities))
+
+    rm(cel, mapType, mapFile, apd)
+  }
+
+
+  # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+  # 2. Read a subset of the units
+  # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+  units <- c(1, 20:205)
+  cel <- readCelUnits(celNames[1], units=units)
+  apd <- readApdUnits(apdNames[1], units=units)
+  stopifnot(identical(apd, cel))
+
+
+  # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+  # 3. Benchmark reading with and without read map
+  # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+  # To minimize the overhead from searching for and reading CDF
+  # file when benchmarking, we read the CDF already here
+  cdfFile <- findCdf(chipType)
+  cdf <- readCdfCellIndices(cdfFile, units=units)
+
+  cat("Benchmarks for read APD file:\n")
+  t0 <- system.time({
+    # Approximately 10 times faster than without a read map
+    apd <- readApdUnits(apdNames[1], cdf=cdf, readMap=readMap)
+  })[3]
+  cat(sprintf("a) with read map by vector: %.3fs [ 1.00x]\n", t0))
+
+  t <- system.time({
+    apd <- readApdUnits(apdNames[1], cdf=cdf)
+  })[3]
+  cat(sprintf("b) with read map from file: %.3fs [%5.2fx]\n", t, t/t0))
+
+  t <- system.time({
+    # Force no read map to be used
+    apd <- readApdUnits(apdNames[1], cdf=cdf, readMap=NULL)
+  })[3]
+  cat(sprintf("c) without read map       : %.3fs [%5.2fx]\n", t, t/t0))
+} # if (length(files) > 0)