diff --git a/NEWS b/NEWS index ad66d71..40386fa 100644 --- a/NEWS +++ b/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. diff --git a/tests/celToApd.R b/tests/celToApd.R new file mode 100644 index 0000000..1f29131 --- /dev/null +++ b/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) diff --git a/tests/createApd.R b/tests/createApd.R new file mode 100644 index 0000000..bf1904f --- /dev/null +++ b/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) diff --git a/tests/gtypeCelToPQ.R b/tests/gtypeCelToPQ.R new file mode 100644 index 0000000..8f41bcf --- /dev/null +++ b/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) +} diff --git a/tests/readApdRectangle.R b/tests/readApdRectangle.R new file mode 100644 index 0000000..27e8481 --- /dev/null +++ b/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) +} diff --git a/tests/readApdUnits.R b/tests/readApdUnits.R new file mode 100644 index 0000000..2659b29 --- /dev/null +++ b/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) diff --git a/tests/writeApdMap.R b/tests/writeApdMap.R new file mode 100644 index 0000000..b4456ff --- /dev/null +++ b/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)