Hi Jim and Martin,

Thanks for looking into this.  That explanation makes a lot of sense and
would certainly cause the new version to behave differently.  In fact, it
sounds like the former version could give different results depending on
the order of the arrays.

Reproducible results is a big part of what GenePattern is trying to do, so
it's good to know about that.  I'll be sure to document the change for our
users.

Thanks again for the help. I really appreciate it!

Regards,
David

-- 
David Eby
Consultant
Cancer Informatics Development
Broad Institute of MIT and Harvard
7 Cambridge Center, Cambridge, MA 02142, USAhttp://www.broadinstitute.org/cancer



On Mon, Jan 28, 2013 at 11:08 AM, James W. MacDonald <jmacdon@uw.edu> wrote:

> Hi David,
>
>
> On 1/28/2013 10:14 AM, Martin Morgan wrote:
>
>> On 01/24/2013 01:41 PM, David Eby [guest] wrote:
>>
>>>
>>> Hello,
>>>
>>> I am in the midst of updating some inherited legacy code from R 2.8 and
>>> affy_1.20.2 to R 2.15.2 and affy_1.36.0, trying to fix some longstanding
>>> bugs.
>>>
>>>
>> That's a lot of water under the bridge (4 1/2 years? do you still have a
>> working installation?)...
>>
>>  Everything is going well except for some very different results coming
>>> out of an expresso call in one of my regression tests, specifically from:
>>>     eset <- expresso(afbatch, normalize.method="**invariantset",
>>> bg.correct=FALSE, pmcorrect.method="pmonly",**summary.method="liwong",
>>> verbose=TRUE)
>>>
>>
>> ...there is this at least in svn
>>
>> affy$ svn log -r57229
>> ------------------------------**------------------------------**
>> ------------
>> r57229 | jmacdon | 2011-08-04 12:18:56 -0700 (Thu, 04 Aug 2011) | 1 line
>>
>> fixed bug in normalize.AffyBatch.**invariantset
>> ------------------------------**------------------------------**
>> ------------
>>
>> and I've added Jim to the email for any insights.
>>
>
> Thanks for the detective work, Martin. This was on my to-do list for
> today, and now I can just press the easy button.
>
> I don't recall how this came up, but there was definitely a bug. The basic
> idea is that we want to say which of the arrays is the reference, based on
> one of four criteria. The default is to base this on the array that has the
> median overall intensity (where overall intensity is defined as the mean
> intensity).
>
> Prior to the fix we weren't correctly figuring out which array had the
> median overall intensity, and were instead simply picking the trunc(N/2)
> array, where N = # arrays. So if you had 8 arrays, this function would
> choose the 4th array as reference regardless.
>
> Now it actually chooses the array with the median overall intensity, so
> the difference between your old results and the new is caused by switching
> the array that is used as reference.
>
> Best,
>
> Jim
>
>>
>> Martin
>>
>>
>>
>>  It's reporting expression values approx. 0.4-0.7x the values in the
>>> earlier setup (examples below).  I take for granted that many things have
>>> changed both at the R level and at the package level given the large leap
>>> in versions, but the difference seemed a bit odd since the results for
>>> other methods (RMA, GCRMA, MAS5) remained reasonably consistent in the
>>> update.
>>>
>>> It's quite possible that there is a package conflict or that I've
>>> accidentally broken something in the code base, though the legacy code
>>> itself is largely unchanged and the other regression tests seem to hold up
>>> pretty well.  I have looked at a number of things already, but before
>>> heading further down that path I wanted to check whether the implementation
>>> itself might have changed in a way where these differences would be
>>> expected.  Has there been a major change to liwong/invariantset somewhere
>>> along the way between 1.20.2 and 1.36.0?
>>>
>>> Checking with other members of the team, we're actually completely OK
>>> with the differences if they are in line with known changes to the affy
>>> package and can be explained to our users (this is a GenePattern module).
>>>
>>> Here is some example output for reference.  This is for a cut-down data
>>> set; I've seen similar results with 20 samples.  I can provide more if it
>>> would be helpful.
>>>
>>>> From the original setup (using write.table(exprs(eset)):
>>>>
>>> "CL20030502207AA.CEL" "CL20030502208AA.CEL" "CL20030502307AA.CEL"
>>> "CL20030502308AA.CEL"
>>> "1007_s_at" 228.013212425507 214.877883873475 287.677963272274
>>> 306.997621485651
>>> "1053_at" 193.206766296132 168.017787218035 169.430151901596
>>> 157.819950438341
>>> <...snip...>
>>>
>>>  From the new setup:
>>>>
>>> "CL20030502207AA.CEL" "CL20030502208AA.CEL" "CL20030502307AA.CEL"
>>> "CL20030502308AA.CEL"
>>> "1007_s_at" 500.169414439674 461.734001304198 704.700664579873
>>> 735.43106514776
>>> "1053_at" 321.921661422307 251.464504650157 261.491260842992
>>> 228.793486504634
>>> <...snip...>
>>>
>>> Thanks in advance!
>>> Regards,
>>> David
>>>
>>>
>>
>>
> --
> James W. MacDonald, M.S.
> Biostatistician
> University of Washington
> Environmental and Occupational Health Sciences
> 4225 Roosevelt Way NE, # 100
> Seattle WA 98105-6099
>
>

	[[alternative HTML version deleted]]

