[BioC] pm summarization method

James W. MacDonald jmacdon at uw.edu
Wed Apr 18 15:10:15 CEST 2012


Hi Assa,

On 4/18/2012 2:33 AM, Assa Yeroslaviz wrote:
> Hi Jim,
>
> On Tue, Apr 17, 2012 at 20:04, James W. MacDonald <jmacdon at uw.edu 
> <mailto:jmacdon at uw.edu>> wrote:
>
>     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> <mailto: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?
>
>
> I basically wanted to show the raw values to the biologists in the 
> lab. The problem was, that they suggested that I somehow might have 
> switched the way of analyzing the data. I mean they thought, that 
> instead of doing Treatment- Control I did Control-Treatment, which 
> would have easily explained the results of the higher expression of 
> this gene in the control.
>
> By showing them the raw data, I can give them a proof, that the switch 
> wasn't made by me, but the arrays were switched, or that the gene is 
> somehow truly down-regulated in the treatment.
>
>
>
>
>         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
>         <tel:36%20%20%20%2015%20%20%20%2022%20%20%20%202%20%20%20%208%20%20%20%202.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.
>
>
> Sorry,but this here I don't understand.
> I didn't present any probes from my data, but just the results of _one 
> gene_ after running the expresso command with _different options_.
> How can you say that: "I only see a few probes in your raw data that 
> have very different expression levels between the sample types"?
> I mean after analyzing the data , I think you are also right, it just 
> wonder me how could you have known that.

What are the data you show just below? You said they were the PM and MM 
probe data for the gene in question.

>
> I was just wondering here, why the medianpolish option shows basically 
> no differences at all, while the others three with the mas or 
> subtractmm options shows such a huge difference.
> Can it be explained?

First, these aren't huge differences. Sure, if you do fold changes they 
are, but remember that these data range from 0 to 2^14 or so 
(hypothetically 2^16, but it is rare to get much over 2^14). A 
difference of 40 at the very low end of that range is indistinguishable 
from noise. Second, both mas5 and avgdiff suffer from the same problems; 
they assume MM probes accurately measure just background (not true), and 
they ignore probe-specific binding, which adds nuisance data to the 
result. Once you account for probe-specific binding, the differences 
largely disappear.

Best,

Jim


>
> Thanks,
> Assa
>
>
>
>     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>
>         <mailto: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
>
>

-- 
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