[BioC] [limma] Strange results in contrasts with dye-swap

Gordon K Smyth smyth at wehi.EDU.AU
Fri Aug 5 00:42:59 CEST 2011


Dear Maciej,

Thanks for including detailed R code.  I'm going to delete most of the 
code though in my reply, so that we can get quickly to the question and 
answer at the bottom.

> Date: Thu, 4 Aug 2011 09:27:51 GMT
> From: mjonczyk at biol.uw.edu.pl
> To: bioconductor at r-project.org
> Subject: [BioC] [limma] Strange results in contrasts with dye-swap
> Message-ID: <201108040927.p749Rnjm001752 at woreczko.biol.uw.edu.pl>
>
> Dear List Members,
>
> I've spent last day searching archive to find out how to account for dye 
> effect in dye swap common-reference design. I've found a lot of useful 
> information but there are few things I still don't understand.
>
> I have time-course experiment with 7 time points, 2 conditions (c, k),
> 4 biological replications (with dye-swap).
>
> Here is my targets frame (powt stands for biol. replication):
>
>> trg
>               FileName Cy3 Cy5 powt
> 1     _024_1_vs_7ex.gpr  k0  k2    1
> 2   _228_k0_vs_16ex.gpr  k0 c12    1
> 3    _229_k0_vs_8ex.gpr  k0  c4    1
> 4   _230_k0_vs_15ex.gpr  k0 k10    1
> 5   _235_k0_vs_19ex.gpr  k0 k14    1
> 6   _237_k0_vs_14ex.gpr  k0 c10    1
> 7   _238_k0_vs_10ex.gpr  k0  c6    1
> 8   _239_k0_vs_12ex.gpr  k0  c8    1
> 9   _240_k0_vs_13ex.gpr  k0  k8    1
> 10  _248_k0_vs_17ex.gpr  k0 k12    1
> 11   _249_k0_vs_9ex.gpr  k0  k4    1
> 12  _253_k0_vs_11ex.gpr  k0  k6    1
> 13   _254_k0_vs_6ex.gpr  k0  c2    1
> 14  _256_k0_vs_18ex.gpr  k0 c14    1
> 15  _189_K8_vs_K0ex.gpr  k8  k0    2
> 16  _190_K4_vs_K0ex.gpr  k4  k0    2
> 17  _191_K6_vs_K0ex.gpr  k6  k0    2
> 18  _192_C8_vs_K0ex.gpr  c8  k0    2
> 19 _193_C10_vs_K0ex.gpr c10  k0    2
> 20  _194_C2_vs_K0ex.gpr  c2  k0    2
> 21 _198_C14_vs_K0ex.gpr c14  k0    2
> 22 _200_K12_vs_K0ex.gpr k12  k0    2
> 23 _202_K10_vs_K0ex.gpr k10  k0    2
> 24  _206_C6_vs_K0ex.gpr  c6  k0    2
> 25 _207_C12_vs_K0ex.gpr c12  k0    2
> 26  _208_K2_vs_K0ex.gpr  k2  k0    2
> 27 _209_K14_vs_K0ex.gpr k14  k0    2
> 28  _210_C4_vs_K0ex.gpr  c4  k0    2
> 29  _201_k0_vs_C2ex.gpr  k0  c2    3
> 30  _203_k0_vs_K6ex.gpr  k0  k6    3
> 31 _205_k0_vs_C12ex.gpr  k0 c12    3
> 32 _217_k0_vs_K14ex.gpr  k0 k14    3
> 33  _219_k0_vs_C4ex.gpr  k0  c4    3
> 34 _231_k0_vs_C10ex.gpr  k0 c10    3
> 35  _232_k0_vs_K8ex.gpr  k0  k8    3
> 36  _233_k0_vs_C8ex.gpr  k0  c8    3
> 37 _234_k0_vs_C14ex.gpr  k0 c14    3
> 38  _235_k0_vs_C6ex.gpr  k0  c6    3
> 39 _236_k0_vs_K10ex.gpr  k0 k10    3
> 40  _237_k0_vs_K2ex.gpr  k0  k2    3
> 41  _238_k0_vs_K4ex.gpr  k0  k4    3
> 42 _239_k0_vs_K12ex.gpr  k0 k12    3
> 43    _004_8_vs_1ex.gpr  c8  k0    4
> 44   _005_10_vs_1ex.gpr c10  k0    4
> 45    _006_2_vs_1ex.gpr  c2  k0    4
> 46   _007_15_vs_1ex.gpr k14  k0    4
> 47    _008_4_vs_1ex.gpr  c4  k0    4
> 48    _010_3_vs_1ex.gpr  k2  k0    4
> 49    _011_9_vs_1ex.gpr  k8  k0    4
> 50   _012_14_vs_1ex.gpr c14  k0    4
> 51   _014_12_vs_1ex.gpr c12  k0    4
> 52   _026_13_vs_1ex.gpr k12  k0    4
> 53   _036_11_vs_1ex.gpr k10  k0    4
> 54    _037_6_vs_1ex.gpr  c6  k0    4
> 55    _043_7_vs_1ex.gpr  k6  k0    4
> 56    _044_5_vs_1ex.gpr  k4  k0    4
>
> I'm interested in contrasts between "c" and "k" in each separate time point.
> I also want to account for dye-effect.

[code deleted]

> contrast.matrix.d2=makeContrasts(DyeEffect,p2="c2-k2",p4="c4-k4",p6="c6-k6",p8="c8-k8",p10="c10-k10",p12="c12-k12",p14="c14-k14",levels=design_dye)

[more code deleted]

>> test.ug.d2=decideTests(fit2.ug.d2,method="global",adjust.method="BH",p.value=0.05)
>
> Results:
>
> TEST WITHOUT DYE EFFECT
>> summary(test.ug)
>   c2 - k2 c4 - k4 c6 - k6 c8 - k8 c10 - k10 c12 - k12 c14 - k14
> -1     115      96     377     141       263       175       265
> 0    43082   43048   42326   42973     42694     42842     42752
> 1      196     249     690     279       436       376       376
>
>
> TEST WITH DYE EFFECT
>> summary(test.ug.d)
>      p2    p4    p6    p8   p10   p12   p14
> -1   246   217   686   317   530   316   526
> 0  42797 42755 41594 42530 42103 42446 42239
> 1    350   421  1113   546   760   631   628
>
> TEST WITH DYE EFFECT AND CONTRAST FOR DYE EFFECT
>> summary(test.ug.d2)
>   DyeEffect    p2    p4    p6    p8   p10   p12   p14
> -1      6660   452   402  1103   593   858   539   883
> 0      30803 42353 42261 40659 41853 41359 41888 41564
> 1       5930   588   730  1631   947  1176   966   946
>
> *QUESTION*
> While I understand that second test (dye-effect included in model) can give more
> significant genes (adjustment for dye-effect) than model without dye-effect, I
> can't understand why the third test (with contrast for dye-effect) gives even
> more significant genes. In the third analysis more contrasts have been performed
> so there should be lower p-value cut-off and consequently one could suppose that
> there will be less significant genes.

Your reasoning would be correct for ordinary hypothesis testing, however 
false discovery rates do not follow the same rules.  FDRs are extensible 
in the sense that you can add more tests without reducing the proportion 
that are called significant.

You are using method="global".  When you have one contrast (in this case 
DyeEffect) that contains many hugely DE genes, it will pull up the other 
contrasts as well, causing them to find more DE genes.  This is because 
the large number of DE genes will allow the BH algorithm to go down to 
lower t-statistics, and this affects all the contrasts.

The lesson here is that you should not include DyeEffect in decideTests() 
with your other contrasts.  When using method="global", you should only 
include contrasts that are closely comparable to one another, and about 
which you will be making conclusions as a group.

> *OTHER QUESTIONS*
> 1. Is second model (WITH DYE EFFECT) correct?

Fine.

> 2. So my anova model is y=Treatment x + DyeEffect, right?

You code looks fine, but I don't have time to check it through line by line.

> 3. How Array effect is taken into account in analysis?

The analysis is in terms of log-ratio (M-values), so any array effect is 
subtracted out right at the start.

> 4. Should I include biol-replication effect in analysis (as block)?

If the biol replicates seem to vary randomly, and are only slightly 
different, then I would suggest that you enter them as a random effect 
using block instead.  If there are large differences between the biol 
reps, in particular if one rep is different to the others, then including 
them in the design matrix as you have done is safer and probably better.

Best wishes
Gordon

> Thanks for any help,
> Best Regards,
> Maciej
>
> ________________________________________________________
> Maciej Jończyk, MSc
> Department of Plant Molecular Ecophysiology
> Institute of Plant Experimental Biology
> Faculty of Biology, University of Warsaw
> 02-096 Warszawa, Miecznikowa
> 1

______________________________________________________________________
The information in this email is confidential and intend...{{dropped:4}}



More information about the Bioconductor mailing list