reordering in extractMatrix.AffymetrixCelSet?
Mark Robinson
When I call readUnits on an AffymetrixCelSet object, I get data
returned in a different order than when I call extractMatrix on the
same object. Is there a reason for that? (I'm creating some plots
where I rely on the order to be the same as in the CDF file)
An example:
> x<-readUnits(csN,units=1)
> log2(x[[1]][[1]]$intensities)[,1:5]
[,1] [,2] [,3] [,4] [,5]
[1,] 1.5622211 2.6349930 1.6494993 0.6173638 0.8673267
[2,] 4.8310593 0.6491851 0.4831402 2.3608759 1.5856072
[3,] 1.6537622 4.4431581 3.8837151 2.5679079 0.9464137
[4,] 0.8082923 1.3786186 1.6937647 0.8430860 0.7713563
> log2(x[[1]][[2]]$intensities)[,1:5]
[,1] [,2] [,3] [,4] [,5]
[1,] 1.8642625 1.1082062 1.1517775 1.0505186 2.1497801
[2,] 0.5102162 0.5758794 0.8986854 0.7551015 0.5601526
[3,] 0.7713563 1.0885455 0.6978568 0.6915522 0.8386120
[4,] 0.9673580 1.3029830 1.2212527 0.6657550 0.5758794
> ind<-unlist(getCellIndices(cdf,units=1),use.names=F)
> xx<-extractMatrix(csN,cells=ind)
> log2(xx[1:8,1:5])
10_5N 11_6T 12_6N 13_7T 14_7N
[1,] 1.5622211 2.6349930 1.6494993 0.6173638 0.8673267
[2,] 0.8082923 1.3786186 1.6937647 0.8430860 0.7713563
[3,] 1.6537622 4.4431581 3.8837151 2.5679079 0.9464137
[4,] 1.8642625 1.1082062 1.1517775 1.0505186 2.1497801
[5,] 0.9673580 1.3029830 1.2212527 0.6657550 0.5758794
[6,] 0.5102162 0.5758794 0.8986854 0.7551015 0.5601526
[7,] 4.8310593 0.6491851 0.4831402 2.3608759 1.5856072
[8,] 0.7713563 1.0885455 0.6978568 0.6915522 0.8386120
> sessionInfo()
R version 2.7.1 (2008-06-23)
i386-apple-darwin8.10.1
locale:
en_AU.UTF-8/en_AU.UTF-8/C/C/en_AU.UTF-8/en_AU.UTF-8
attached base packages:
[1] grid tools stats graphics grDevices utils datasets
[8] methods base
other attached packages:
[1] GenomeGraphs_1.0.0 biomaRt_1.14.0 RCurl_0.9-3
[4] preprocessCore_1.2.0 affyio_1.8.0 Biobase_2.1.7
[7] aroma.affymetrix_0.9.3 aroma.apd_0.1.3 R.huge_0.1.5
[10] affy_1.18.2 affxparser_1.12.2 aroma.core_0.9.3
[13] sfit_0.1.5 aroma.light_1.8.1 digest_0.3.1
[16] matrixStats_0.1.2 R.rsp_0.3.4 R.cache_0.1.7
[19] R.utils_1.0.2 R.oo_1.4.3 R.methodsS3_1.0.1
loaded via a namespace (and not attached):
[1] XML_1.95-2
I've made a modified version of extractMatrix.AffymetrixCelSet below,
with a line commented out, that seems to fix the problem. Just wanted
to know if there is a reason to have it that way ... or, am I doing
something wrong?
THanks,
Mark
extractMatrix.AffymetrixCelSet<-
function (this, cells = NULL, ..., field = c("intensities", "stdvs",
"pixels"), verbose = FALSE)
{
cdf <- getCdf(this)
if (is.null(cells)) {
ncells <- nbrOfCells(cdf)
}
else {
cells <- Arguments$getIndices(cells, range = c(1,
nbrOfCells(cdf)))
ncells <- length(cells)
}
field <- match.arg(field)
verbose <- Arguments$getVerbose(verbose)
if (verbose) {
pushState(verbose)
on.exit(popState(verbose))
}
settings <- getOption("aroma.affymetrix.settings")
gcArrayFrequency <- settings$memory$gcArrayFrequency
if (is.null(gcArrayFrequency))
gcArrayFrequency <- 10
verbose && enter(verbose, "Getting data for the array set")
verbose && enter(verbose, "Allocating matrix")
arrayNames <- getNames(this)
nbrOfArrays <- length(arrayNames)
df <- matrix(NA, nrow = ncells, ncol = nbrOfArrays)
colnames(df) <- arrayNames
verbose && str(verbose, df)
verbose && printf(verbose, "RAM: %.2fMB\n", object.size(df)/1024^2)
verbose && exit(verbose)
if (!is.null(cells)) {
verbose && enter(verbose, "Optimize reading order")
o <- order(cells)
cells <- cells[o]
#o <- order(o) #### Mark R added this HERE ####
verbose && exit(verbose)
}
else {
o <- seq(length = ncells)
}
gc <- gc()
verbose && print(verbose, gc)
verbose && enter(verbose, "Retrieving data")
for (aa in seq_len(nbrOfArrays)) {
verbose && printf(verbose, "Array %d,\n", aa)
cf <- getFile(this, aa)
df[o, aa] <- getData(cf, indices = cells, fields = field,
verbose = less(verbose))[[field]]
if (aa%%gcArrayFrequency == 0) {
gc <- gc()
verbose && print(verbose, gc)
}
}
verbose && exit(verbose)
verbose && exit(verbose)
df
}
------------------------------
Mark Robinson
Epigenetics Laboratory, Garvan
Bioinformatics Division, WEHI
e: m.rob...@garvan.org.au
e: mrob...@wehi.edu.au
p: +61 (0)3 9345 2628
f: +61 (0)3 9347 0852
------------------------------
Henrik Bengtsson
yes, the extractMatrix(..., cells) of AffymetrixCelSet should return
the data in the order of the cells you've requested. It looks like
you're nailed the cause of the problem.
FYI, the extractMatrix() of AffymetrixCelSet has been used much less
than the corresponding one for ChipEffectSet - the quality of the
latter(s) should be higher - I hope ;)
Thanks for the troubleshooting and the bug fix. It will be available
in the next release. In meantime, use the follow patch:
library(aroma.affymetrix);
downloadPackagePatch("aroma.affymetrix");
/Henrik
> e: m.robin...@garvan.org.au
> e: mrobin...@wehi.edu.au