[BioC] pm summarization method

James W. MacDonald jmacdon at uw.edu
Tue Apr 17 20:04:25 CEST 2012


Hi Assa,

On 4/17/2012 3:09 AM, Assa Yeroslaviz wrote:
> Hi Jim and bioC-Users,
>
>
> On Mon, Apr 16, 2012 at 17:54, James W. MacDonald <jmacdon at uw.edu 
> <mailto:jmacdon at uw.edu>> wrote:
>
>     Hi Assa,
>
>
>     On 4/16/2012 9:05 AM, Assa Yeroslaviz wrote:
>
>         Hi everybody,
>
>         I have a question about the behavior of the expresso command when
>         extracting the raw data from an affyBatch.
>
>         I wanted to evaluate the raw intensities values of a specific
>         gene from my
>         data set and tried to extract it like that:
>         rawdata<- expresso(totalExpressionData, bg.correct=FALSE,
>                                normalize=FALSE,
>                                pmcorrect.method="pmonly",
>         summary.method="avgdiff",
>                                verbose=TRUE)
>
>         I've got the result I wanted:
>                      wt1    wt2    wt3    treat1    treat2    treat3
>         gene_at    125.5    101    123.5    52.5    63.5    58
>
>         The problem was that i expected them to be the other way around.
>
>
>     Why did you expect it to be the other way around? What exactly did
>     you expect this call to expresso() to do?
>
>
> Sorry I didn't explain a bit more. The gene I was investigating is the 
> 'white' gene in Drosophila, which sets the eye color of the flies. It 
> was used in our experiment to test the flies if they have a certain 
> knockout/mutant/treatment. By 'expect it to be the other way around' I 
> meant that only the treated flies should'veh ad a high level of the 
> white gene. In this data it is the other way around.
> We are not sure, if the error is in the data or maybe it was a error 
> made in the lab who analyzed the arrays and made an erroneous 
> identification symbols for the arrays.
>
> That was a little background information about the experiment. I hope 
> it clarify the problem better.
>
>     Since affy is primarily based on S4 methods, it can be a bit
>     difficult to figure out what a given function is going to do, so I
>     can understand you not knowing what this call to expresso() is
>     going to end up doing. However, what you are doing is pretty
>     weird, no? The avgdiff method implies pm-mm, so what do you expect
>     to happen if you then specify pmonly?
>
>     Given that the pmcorrect.method controls how we correct the PM
>     probes, and there is a subtractmm option, one would normally
>     assume that the 'difference' part of avgdiff might happen in that
>     step. But you said not to compute that, so all you are left with
>     is the 'avg' part of avgdiff.
>
>
> Now for the point of my command. To see if I did wrongly classify the 
> experiments or it was given to me as such, I wanted to extract the raw 
> data from the arrays. As I said before, I know one can extract the PM 
> and MM values for each of the probes of a probe-set. What I really 
> would like to have is _one value for each probe-set and array each_.
>
> For that reason I used the summarization of the values with the 
> expresso command. I will try to answer your questions as I understand 
> them.
> I used the pmonly correction method, because if I understood it the 
> right way, this is what the RMA algorithm is using when running the 
> normalization procedure on the data. Am I wrong about that?
>
> Which way is the best way to extract the 'real' raw values from the 
> arrays?
> I know that asking such question get you never a straightforward 
> answer, but a near one will be good :-)

Why do you need a raw value from the array to see if you mixed samples? 
To put the question a different way, if normalized, background corrected 
expression values are good enough for doing statistics, why are they not 
good enough to determine if you mixed samples up?

>
> If my logical thinking is wrong, will this cmmand will present me with 
> a more accurate results?
> rawdatasubtractmm <- expresso(totalExpressionData, bg.correct=FALSE,
>                                 normalize=FALSE,
>                                 pmcorrect.method="subtractmm", 
> summary.method="avgdiff",
>                                 verbose=TRUE)
>
> I will run it now and look for the differences.
>
>
>     But let's set the logical assumptions aside and look at the actual
>     code. Going through the code to expresso() is a bit like following
>     Alice down the rabbit hole, so I will cut to the chase. In the
>     end, you will be calling two pieces of code that will handle the
>     pm adjustment and the summary statistic calculation. In your call
>     to expresso() these will be (respectively):
>
>     > pmcorrect.pmonly
>     function (object)
>     {
>        return(pm(object))
>     }
>
>     and
>
>     > generateExprVal.method.avgdiff
>     function (probes, ...)
>     {
>        list(exprs = apply(probes, 2, median), se.exprs = apply(probes,
>            2, sd)/sqrt(nrow(probes)))
>     }
>
>     So the pmonly will just give you the pm probes, and then avgdiff
>     will give you the column medians of the pm data. Therefore you
>     have basically told expresso() that you want the median value for
>     the non-background corrected, unnormalized pm probes.
>
>
> In the vignette from affy it said to give back the average of the 
> values, so I expected it to give me back the mean, not the median. May 
> be this is why I didn't really understood the results. Now it make 
> more sense.
>
>     Was that your intention?
>
>
> As I said earlier, I wanted to extract the 'real' values for the 
> probe-sets and I thought that extracting the PM values with no 
> correction or normalization will be the best method to achieve it.
> I ran now several other possibilities:
> pm.correction    summary.method    gene    wt1    wt2    wt3    
> treat1    treat2    treat3
> mas    mas    gene_at    52.01982773    34.16503645    46.87094917    
> 4.081420059    9.505471139    3.538333392
> pmonly    avgdiff    gene_at    125.5    101    123.5    52.5    
> 63.5    58
> pmonly    medianpolish    gene_at    6.862247387    6.564551864    
> 6.926547933    6.093298307    6.355810999    6.105348101
> subtractmm    avgdiff    gene_at    36    15    22    2    8    2.5
>
> This are the results.
> The mas and the subtractmm gives me the differences between PM and MM. 
> Though the results are quite different, both of them gives me the same 
> behavior. the expression in higher in the wild-type.
> The medianpolish gives me the log2 values of the PM values.  This 
> results I don't understand, a I can clearly see a difference between 
> the wt and the treatemnt, but the log-scale values show almost no 
> difference at all.

Here is where our interpretations of the raw data diverge. I only see a 
few probes in your raw data that have very different expression levels 
between the sample types. The only probes that look to be binding at any 
appreciable level are the first two, and there really isn't any 
difference between them. In addition, if you believe that MM probes 
measure background, it doesn't really look like you have much if any of 
this transcript being expressed in either sample.

Best,

Jim


>
> I hope it make some more sense and I would be happy to hear some more 
> suggestions/comments on that topic
>
> thanks a lot,
> Assa
>
>
>     Best,
>
>     Jim
>
>
>
>         So I decided to look into the specific probe values of the
>         probes for this
>         probe-set.
>         This are the values I've got from the PM and MM respectively:
>             wt1    wt2    wt3    treat1    treat2    treat3
>         probe1    403    379    220    420    530    316
>         probe2    117    84    104    52    57    54
>         probe3    49    49    73    38    58    52
>         probe4    87    67    110    55    43    49
>         probe5    66    61    51    46    72    62
>         probe6    118    100    104    69    87    74
>         probe7    180    142    170    45    46    45
>         probe8    133    102    137    95    132    81
>         probe9    80    71    65    52    54    46
>         probe10    63    45    56    53    53    54
>         probe11    293    321    260    444    618    408
>         probe12    171    167    169    49    75    72
>         probe13    198    197    307    40    67    50
>         probe14    247    265    348    53    60    62
>
>         probe1    533    519    294    507    739    404
>         probe2    1789    1271    1468    1430    1666    1552
>         probe3    56    66    59    51    45    48
>         probe4    49    52    64    47    45    47
>         probe5    54    47    33    49    65    55
>         probe6    84    72    90    53    92    71
>         probe7    73    72    65    40    53    54
>         probe8    83    108    115    81    111    94
>         probe9    49    56    43    52    41    53
>         probe10    56    46    62    68    77    57
>         probe11    54    83    55    47    64    46
>         probe12    106    98    76    52    66    53
>         probe13    43    48    37    36    52    39
>         probe14    94    92    99    43    43    49
>
>         When I calculate the average of these two tables for each
>         array I don't get
>         the same values as presented in the top table.
>         I would like to understand how from the values on the last two
>         tables I
>         come to a summarized value I get. Even if I ignore the MM values
>         completely, which I think it does, I still don't see how it
>         comes to these
>         values. The two probes (Nr. 1 and 11 of the PM values) are
>         strongly differ
>         from the rest of the probes for this probe-set. Are they being
>         ignored in
>         the summarization?
>
>         Thanks in advance
>
>         Assa
>
>                [[alternative HTML version deleted]]
>
>         _______________________________________________
>         Bioconductor mailing list
>         Bioconductor at r-project.org <mailto:Bioconductor at r-project.org>
>         https://stat.ethz.ch/mailman/listinfo/bioconductor
>         Search the
>         archives:http://news.gmane.org/gmane.science.biology.informatics.conductor
>
>
>     -- 
>     James W. MacDonald, M.S.
>     Biostatistician
>     University of Washington
>     Environmental and Occupational Health Sciences
>     4225 Roosevelt Way NE, # 100
>     Seattle WA 98105-6099
>
>

-- 
James W. MacDonald, M.S.
Biostatistician
University of Washington
Environmental and Occupational Health Sciences
4225 Roosevelt Way NE, # 100
Seattle WA 98105-6099



More information about the Bioconductor mailing list