Hi,

 

I'm building an array analysis package using an MCMC process for each gene &
am writing the key values (posterior quantiles, mean, sd etc) for each gene
to disk using CEL files and the aroma.affymetrix architecture. However I've
noticed that on rare occasion, the incorrect values are written to the CEL
files. 

 

As an example of what goes wrong, if I have a vector of 8 values (e.g.
intensities=1:8) and am updating 8 cells (e.g. indices=221632:221639) with
these 8 values, occasionally the first value (i.e. intensities[1]) would be
incorrectly given to all 8 cells on the CEL file. This has happened in an
obvious manner only to 2 different genes in 2 of 3 datasets & I can't spot a
connection in the data structures that may cause the misfire. If I can see
these 2 obvious problems, I'm unsure as to what less obvious errors there
may be lurking in the dataset & hence my confidence in the results is
somewhat diminished.

 

The code from the relevant section of my function (where I write batches of
genes at a time) is:

monoUgcMap <- getUnitGroupCellMap(monoCdf, unit=batchUnits) # Get the
UnitGroupCellMap from a monocell cdf for the current batch units

indices <- monoUgcMap$cell

for (j in 1:9) { # The output after the MCMC process will always have 9
columns, which are written to separate CEL files

      updateCel(filenames[j], indices=indices, intensities=exp(output[,j]))
#update the CEL file for the current batch

}

 

I am wondering if the issue arises from the usage of 'seek' within the
function 'updateCel', as I am on a Windows 7 x64 system? The line I'm
referring to  from the 'updateCel' source code is:

 

seek(con, origin="start", where=offset, rw="read")

 

Reading the 'seek' help page didn't really shed much light (to my flawed
brain) on what this step does, but I did spot the warning about this
function under Windows.  I'm still a bit unclear as to what this step does,
so am looking for either a clarification as to the point of this line of
code, or any suggestions for a simple workaround. 

 

As it appears to be a random & rare misfire, I'll have to make changes & run
a few datasets (over multiple days) just to check it so may be a bit slow
with any responses...

 

Thanks in advance,

 

Steve


	[[alternative HTML version deleted]]

