[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