[BioC] limma, voom, weights: factorial versus nested interaction results
Reiner Schulz
reiner.schulz at kcl.ac.uk
Wed Jan 9 09:38:17 CET 2013
Dear Gordon,
Thanks again. That's very helpful. I am working w/ iCLIP data. In terms
of size, probably similar to SAGE: only 10^5 to 10^6 tags/reads for a
typical human sample. Have applied both edgeR and voom, and the results
do not seem to be widely different; though not what I expected, but
that's another matter.
Best, Reiner
On 06/01/13 05:08, Gordon K Smyth wrote:
> Dear Reiner,
>
> Well, I don't honestly know yet which of edgeR or voom are better in
> what circumstances. I feel that I can't know this until I have pushed
> both methods as far as they can go, something that I am still doing, and
> then have tried them both out on a wide variety of datasets.
>
> I started the edgeR project way back in 2004, before RNA-Seq technology
> even existed. My thought then was that non-normal methods would be
> essential. And I think that it remains true that edgeR has a clear edge
> over other methods for the SAGE data that was available then.
>
> edgeR performs very well, but has not solved all our data analysis
> problems for RNA-Seq data, not yet anyway.
>
> voom has surprised me by performing much better than I expected. It
> solves the biggest problem with normal-based methods by estimating the
> mean-variance relationship adaptively. It holds its size correctly in
> almost any circumstance, scales to very large datasets, can adaptively
> estimate the amount of smoothing necessary, and gives immediate access
> to all the downstream limma pipeline including gene set testing.
>
> My feeling the moment is that edgeR is superior for small counts and but
> that voom is safer and more reliable for noisy heterogeneous data. Only
> edgeR can estimate the biological coefficient of variation (as we
> defined this in our 2012 paper). But we are actively working on both
> methods, and are open to what we find.
>
> Best wishes
> Gordon
>
> ---------------------------------------------
> Professor Gordon K Smyth,
> Bioinformatics Division,
> Walter and Eliza Hall Institute of Medical Research,
> 1G Royal Parade, Parkville, Vic 3052, Australia.
> http://www.statsci.org/smyth
>
> On Fri, 4 Jan 2013, Reiner Schulz wrote:
>
>> Thank you very much, Gordon.
>>
>> Out of curiosity, why did you develop the voom method when there is
>> edgeR? Are there good, perhaps data-dependent, reasons for using one
>> over the other?
>>
>> Best, Reiner
>>
>> On 03/01/13 23:59, Gordon K Smyth wrote:
>>> Dear Reiner,
>>>
>>> You don't mention contrasts.fit(), but I assume that you are using that
>>> function to get equivalent quantities from the two design matrices.
>>>
>>> If I am guessing your problem correctly, then the issue is that
>>> contrasts.fit() gives results that are not exact in the presence of voom
>>> weights. See the warning note on the help page for contrasts.fit.
>>>
>>> You will get the exact correct t-statistic if the contrast you want is a
>>> coefficient of your design matrix. You should also get exact correct
>>> results when using the one way layout design formula ~0+cond, but you
>>> may get an approximate test statistic if you use a contrast from the
>>> nested formula. The fold changes are exact in all cases.
>>>
>>> There was a thread about this a few years ago, see:
>>>
>>> https://www.stat.math.ethz.ch/pipermail/bioconductor/2009-December/030947.html
>>>
>>>
>>>
>>> Basically I have promised to change contrasts.fit() to give exact
>>> results in all cases, even though this will make it much slower in some
>>> cases. I have not so far found the time to do it, but will definitely
>>> do so by the time we publish the voom method.
>>>
>>> Best wishes
>>> Gordon
>>>
>>> ---------------------------------------------
>>> Professor Gordon K Smyth,
>>> Bioinformatics Division,
>>> Walter and Eliza Hall Institute of Medical Research,
>>> 1G Royal Parade, Parkville, Vic 3052, Australia.
>>> http://www.statsci.org/smyth
>>>
>>>
>>>
>>>> Date: Wed, 2 Jan 2013 13:31:58 +0100
>>>> From: Reiner Schulz <reiner.schulz at kcl.ac.uk>
>>>> To: <bioconductor at r-project.org>
>>>> Subject: [BioC] limma, voom, weights: factorial versus nested
>>>> interaction results
>>>>
>>>> Dear all,
>>>>
>>>> Am analysing an RNA-seq experiment with a two-way factorial design, and
>>>> am observing different results for two contrasts of interest depending
>>>> on whether I use a factorial model matrix or a nested interaction model
>>>> matrix. I only observe this if the data are pre-processed w/ voom:
>>>>> counts.voom= voom( ...
>>>> i.e., when there are observational level weights passed to lmFit. If I
>>>> only pass counts.voom$E (no observational level weights) to lmFit, the
>>>> results for the two contrasts of interest are identical for the
>>>> factorial and nested interaction model matrices.
>>>> Is this to be expected, and if so, what is the explanation?
>>>>
>>>> Much appreciated, Reiner
>>>>
>>>> PS: some more details ...
>>>>
>>>>> strain
>>>> [1] MUT MUT MUT MUT MUT MUT MUT MUT MUT MUT MUT MUT MUT WT WT WT WT
>>>> WT WT
>>>> [20] WT WT WT WT WT WT WT
>>>> Levels: WT MUT
>>>>> treat
>>>> [1] U T2 T1 U T2 T1 U U T2 T1 U T2 T1 U T2 T1 U T2 T1 U U T2
>>>> T1 U T2
>>>> [26] T1
>>>> Levels: U T1 T2
>>>>
>>>>> cond= factor( paste( sep='.', strain, treat), lev=c('WT.U', 'MUT.U',
>>>> 'WT.T1', 'MUT.T1', 'WT.T2', 'MUT.T2'))
>>>>> design.matrix= model.matrix( ~0 + cond)
>>>>> design.matrix
>>>> WT.U MUT.U WT.T1 MUT.T1 WT.T2 MUT.T2
>>>> 1 0 1 0 0 0 0
>>>> 2 0 0 0 0 0 1
>>>> 3 0 0 0 1 0 0
>>>> 4 0 1 0 0 0 0
>>>> 5 0 0 0 0 0 1
>>>> 6 0 0 0 1 0 0
>>>> 7 0 1 0 0 0 0
>>>> 8 0 1 0 0 0 0
>>>> 9 0 0 0 0 0 1
>>>> 10 0 0 0 1 0 0
>>>> 11 0 1 0 0 0 0
>>>> 12 0 0 0 0 0 1
>>>> 13 0 0 0 1 0 0
>>>> 14 1 0 0 0 0 0
>>>> 15 0 0 0 0 1 0
>>>> 16 0 0 1 0 0 0
>>>> 17 1 0 0 0 0 0
>>>> 18 0 0 0 0 1 0
>>>> 19 0 0 1 0 0 0
>>>> 20 1 0 0 0 0 0
>>>> 21 1 0 0 0 0 0
>>>> 22 0 0 0 0 1 0
>>>> 23 0 0 1 0 0 0
>>>> 24 1 0 0 0 0 0
>>>> 25 0 0 0 0 1 0
>>>> 26 0 0 1 0 0 0
>>>>
>>>>> design.matrix.alt= model.matrix( ~strain + strain:treat)
>>>>> design.matrix.alt
>>>> Intercept MUT WT:T1 MUT:T1 WT:T2 MUT:T2
>>>> 1 1 1 0 0 0 0
>>>> 2 1 1 0 0 0 1
>>>> 3 1 1 0 1 0 0
>>>> 4 1 1 0 0 0 0
>>>> 5 1 1 0 0 0 1
>>>> 6 1 1 0 1 0 0
>>>> 7 1 1 0 0 0 0
>>>> 8 1 1 0 0 0 0
>>>> 9 1 1 0 0 0 1
>>>> 10 1 1 0 1 0 0
>>>> 11 1 1 0 0 0 0
>>>> 12 1 1 0 0 0 1
>>>> 13 1 1 0 1 0 0
>>>> 14 1 0 0 0 0 0
>>>> 15 1 0 0 0 1 0
>>>> 16 1 0 1 0 0 0
>>>> 17 1 0 0 0 0 0
>>>> 18 1 0 0 0 1 0
>>>> 19 1 0 1 0 0 0
>>>> 20 1 0 0 0 0 0
>>>> 21 1 0 0 0 0 0
>>>> 22 1 0 0 0 1 0
>>>> 23 1 0 1 0 0 0
>>>> 24 1 0 0 0 0 0
>>>> 25 1 0 0 0 1 0
>>>> 26 1 0 1 0 0 0
>>>>
>>>> Observed equalities (as expected):
>>>> MUT == MUT.U - WT.U
>>>> WT:T1 == WT.T1 - WT.U
>>>> MUT:T1 == MUT.T1 - MUT.U
>>>> WT:T2 == WT.T2 - WT.U
>>>> MUT:T2 == MUT.T2 - MUT.U
>>>> MUT:T1 - WT:T1 == (MUT.T1 - MUT.U) - (WT.T1 - WT.U)
>>>> MUT:T2 - WT:T2 == (MUT.T2 - MUT.U) - (WT.T2 - WT.U)
>>>>
>>>> Observed inequalities (only w/ voom pre-processed data; unexpected):
>>>> WT:T1 - WT:T2 != (WT.T1 - WT.U) - (WT.T2 - WT.U)
>>>> MUT:T1 - MUT:T2 != (MUT.T1 - MUT.U) - (MUT.T2 - MUT.U)
>>>>
>>>> --
>>>>
>>>
>>
>> --
>> Reiner Schulz, Ph.D.
>> Senior Lecturer in Bioinformatics and Epigenomics
>> Dept. of Medical and Molecular Genetics
>> King's College London
>> Guy's Hospital, 8th floor Tower Wing
>> London SE1 9RT
>> https://atlas.genetics.kcl.ac.uk/~rschulz
>>
>>
>
> ______________________________________________________________________
> The information in this email is confidential and intended solely for
> the addressee.
> You must not disclose, forward, print or use it without the permission
> of the sender.
> ______________________________________________________________________
>
--
Reiner Schulz, Ph.D.
Senior Lecturer in Bioinformatics and Epigenomics
Dept. of Medical and Molecular Genetics
King's College London
Guy's Hospital, 8th floor Tower Wing
London SE1 9RT
https://atlas.genetics.kcl.ac.uk/~rschulz
--
Reiner Schulz, Ph.D.
Senior Lecturer in Bioinformatics and Epigenomics
Dept. of Medical and Molecular Genetics
King's College London
Guy's Hospital, 8th floor Tower Wing
London SE1 9RT
https://atlas.genetics.kcl.ac.uk/~rschulz
More information about the Bioconductor
mailing list