[BioC] Problem with A mean in Limma

Gordon Smyth smyth at wehi.edu.au
Fri Feb 24 01:27:48 CET 2006


The component Amean in the linear model fit output contains the mean 
A-values for all the arrays in the original linear model fit. It does 
not change when a contrast is applied to the fit object, in the same 
way that other quantities such as sigma do not change. This is the 
intended behavior and is not a bug.

This issue has been discussed a few times on the list before. Try 
searching the archives for "Amean" or similar.

I personally do not see any good reason why Amean should change with 
the contrast, but I remain open to a good argument.

Best wishes
Gordon

>Date: Wed, 22 Feb 2006 14:45:07 +0100
>From: david hot <David.Hot at pasteur-lille.fr>
>Subject: [BioC] Problem with A mean in Limma
>To: Bioconductor Mailing list <bioconductor at stat.math.ethz.ch>
>Message-ID: <43FC6AE2.D7BF9597 at pasteur-lille.fr>
>Content-Type: text/plain; charset=iso-8859-2
>
>Dear All,
>
>I have posted a mail to the list 2 weeks ago with no answer so far. 
>My first message was probably not very clear so I try again.
>
>I am using Limma 2.4.7 on R 2.2.1 to analyse two-color arrays and I 
>have noticed something strange in the way the mean 'A' value is 
>calculated when one use a contrast matrix.
>In my experiment I want to compare the effect of a stimulation on 
>Tcells. This stimulation was applied on Tcells during 6h or during 
>24h. The 6h stimulated Tcells and 24h stimulated Tcells were 
>compared to the same non stimulated cells (= reference). I have got 
>the following Targets Frame :
>
>FileName Cy3 Cy5
>CtrlCy5vsTh2_24hCy3.dav 24h Ref
>CtrlCy5vsTh2_6hCy3.dav 6h Ref
>Th2_24hCy5vsCtrlCy3.dav Ref 24h
>Th2_6hCy5vsCtrlCy3.dav Ref 6h
>
>I would like to know which genes are modulated after 6h of 
>stimulation, which are modulated after 24h and genes modulated 
>between 6h and 24h. After normalisation of the data I use the 
>following script to create a matrix and a contrast matrix:
>
> > design <- cbind("6h"=c(0,-1,0,1),"24h"=c(-1,0,1,0))
> > targets <- readTargets()
> > rownames(design) <- removeExt(targets$FileName)
> > design
>                     6h 24h
>CtrlCy5vsTh2_24hCy3  0  -1
>CtrlCy5vsTh2_6hCy3  -1   0
>Th2_24hCy5vsCtrlCy3  0   1
>Th2_6hCy5vsCtrlCy3   1   0
> > fit <- lmFit(MA, design)
> > cont.matrix <- cbind("24h-6h"=c(-1, 1), "6h"=c(1,0), "24h"=c(0,1))
> > rownames(cont.matrix) <- cbind("6h","24h")
> > cont.matrix
>     24h-6h 6h 24h
>6h      -1  1   0
>24h      1  0   1
> > fitcont <- contrasts.fit(fit, cont.matrix)
> > fitcontbayes<-eBayes(fitcont)
>
>If one look at the results of the 3 coef of the contrast matrix (: 
>24h-6h, 6h and 24h) everything seem OK :
>
> > topTable(fitcontbayes, coef=1, n=10, adjust="fdr")
>       Block Column Row     ID      Name 
> Status         M        A         t    P.Value        B
>21897    39      9   1  46861      PRAC   Gene  2.876222 
>8.631290  8.433385 0.03855747 3.506137
>20473    36      1  14  28828     PAPD1   Gene -2.383738 6.927029 
>-7.473090 0.03855747 2.783536
>19040    34      8   2  78497    CLDN23   Gene -3.139618 5.859508 
>-8.299092 0.03855747 2.386795
>6221     11      5  20  84193 LOC115749   Gene  3.010418 
>5.792223  7.957571 0.03855747 2.193344
>19497    34      9  21 125848    CRISP3   Gene -2.376080 6.925017 
>-6.742931 0.07172051 2.143568
>21099    37      3  16 147050    HCG2P7   Gene -2.724596 6.706554 
>-7.202047 0.07172051 1.713956
>6196     11      4  19  64892    ETR101   Gene  1.954074 
>8.169356  6.048570 0.13538457 1.453035
>6199     11      7  19 149357       #NA   Gene  1.898843 
>8.060294  5.954056 0.14110973 1.352464
>21132    37     12  17 155265 LOC440590   Gene -2.527218 6.638719 
>-6.680307 0.10272218 1.336678
>20556    36     12  17 160410       #NA   Gene  2.082918 
>7.656239  5.881841 0.14362662 1.274528
> > topTable(fitcontbayes, coef=2, n=10, adjust="fdr")
>       Block Column Row     ID         Name 
> Status         M         A         t      P.Value        B
>2603      5     11  13 
>128320          CA2   Gene  2.312298  8.323775 12.690864 0.0004900596 7.446470
>23250    41     18   9  83486       CYP1B1   Gene  1.365898 
>9.351358  7.529896 0.0168728265 4.118037
>14667    26      3  12  89981          TPO   Gene 
>-1.370724  9.825983 -7.441705 0.0168728265 4.018976
>17543    31     23  11  95124          TNC   Gene  1.362982 
>11.478777  7.450669 0.0176675806 3.633118
>27313    48      1  11  59209         CTSC   Gene  1.337229 
>10.259389  6.727203 0.0276792498 3.169911
>2542      5     22  10  91105        FABP5   Gene  1.180159 
>10.683417  6.708199 0.0276792498 3.146168
>19497    34      9  21 
>125848       CRISP3   Gene  2.254023  6.925017  7.834144 0.0168728265 3.042473
>4862      9     14  11 
>158129       CAPN14   Gene  2.252979  5.170231  8.422209 0.0168728265 3.002382
>18922    33     10  21  49700 
>DKFZp547F072   Gene  2.241397  6.368605  8.378910 0.0168728265 2.975528
>20210    36      2   3 
>129217     MGC10854   Gene  2.155348  6.019866  8.057238 0.0176289711 2.768940
> > topTable(fitcontbayes, coef=3, n=10, adjust="fdr")
>       Block Column Row     ID        Name 
> Status         M         A          t      P.Value        B
>23257    41      1  10 
>119072       CDH26   Gene  2.045561  8.618405  10.977306 0.0007520237 8.012605
>2550      5      6  11 
>112799    SERPINB3   Gene  1.905378  9.951108  10.146438 0.0007520237 7.297240
>21897    39      9   1  46861        PRAC   Gene  3.036560  8.631290 
>  10.904529 0.0007520237 6.084763
>6221     11      5  20  84193   LOC115749   Gene  3.157852  5.792223 
>  11.804853 0.0007520237 5.857994
>20473    36      1  14  28828       PAPD1   Gene -2.693871  6.927029 
>-10.343418 0.0009737256 5.746556
>27313    48      1  11  59209        CTSC   Gene  1.677554 
>10.259389   8.439275 0.0036026318 5.575635
>17543    31     23  11  95124         TNC   Gene  2.432964 
>11.478777   9.404284 0.0022965417 5.110162
>18097    32      1  11 
>104580        CTSC   Gene  1.414608  9.663292   8.014200 0.0056176380 5.087315
>22151    39     23  11 
>104788      COL6A3   Gene  1.507607  9.470752   7.775592 0.0069709347 4.801786
>8275     15     19   9  53347 TXN delta 3   Gene  1.340699 
>11.341639   7.522172 0.0090096504 4.489250
>
>
>But when I sort the table according to the 'A' value one can see 
>that the genes have exactly the same 'A' values for the 1st the 2d 
>and the 3d coef :
>
> > topTable(fitcontbayes, coef=1, n=10, adjust="fdr", sort.by="A")
>       Block Column Row     ID      Name Status  M        A  t P.Value  B
>4716      9     12   5  17228     CPSF4   Gene NA 15.16203 NA      NA NA
>2944      6     16   3  96063      RBM4   Gene NA 15.12503 NA      NA NA
>772       2      4   9   5198    KCNAB3   Gene NA 15.06048 NA      NA NA
>18041    32     17   8 101204    GABBR1   Gene NA 15.03191 NA      NA NA
>16817    30     17   5  15058    PRKAA1   Gene NA 14.98546 NA      NA NA
>7642     14     10   7  90324     GSTA3   Gene NA 14.85254 NA      NA NA
>5105      9     17  21 149752 LOC199800   Gene NA 14.83062 NA      NA NA
>20702    36     14  23 104223       #NA   Gene NA 14.81930 NA      NA NA
>4270      8     22  10  94545  SERPINA6   Gene NA 14.79895 NA      NA NA
>6429     12     21   4  96944    MAD1L1   Gene NA 14.79082 NA      NA NA
> > topTable(fitcontbayes, coef=2, n=10, adjust="fdr", sort.by="A")
>       Block Column Row     ID      Name Status  M        A  t P.Value  B
>4716      9     12   5  17228     CPSF4   Gene NA 15.16203 NA      NA NA
>2944      6     16   3  96063      RBM4   Gene NA 15.12503 NA      NA NA
>772       2      4   9   5198    KCNAB3   Gene NA 15.06048 NA      NA NA
>18041    32     17   8 101204    GABBR1   Gene NA 15.03191 NA      NA NA
>16817    30     17   5  15058    PRKAA1   Gene NA 14.98546 NA      NA NA
>7642     14     10   7  90324     GSTA3   Gene NA 14.85254 NA      NA NA
>5105      9     17  21 149752 LOC199800   Gene NA 14.83062 NA      NA NA
>20702    36     14  23 104223       #NA   Gene NA 14.81930 NA      NA NA
>4270      8     22  10  94545  SERPINA6   Gene NA 14.79895 NA      NA NA
>6429     12     21   4  96944    MAD1L1   Gene NA 14.79082 NA      NA NA
> > topTable(fitcontbayes, coef=3, n=10, adjust="fdr", sort.by="A")
>       Block Column Row     ID      Name Status  M        A  t P.Value  B
>4716      9     12   5  17228     CPSF4   Gene NA 15.16203 NA      NA NA
>2944      6     16   3  96063      RBM4   Gene NA 15.12503 NA      NA NA
>772       2      4   9   5198    KCNAB3   Gene NA 15.06048 NA      NA NA
>18041    32     17   8 101204    GABBR1   Gene NA 15.03191 NA      NA NA
>16817    30     17   5  15058    PRKAA1   Gene NA 14.98546 NA      NA NA
>7642     14     10   7  90324     GSTA3   Gene NA 14.85254 NA      NA NA
>5105      9     17  21 149752 LOC199800   Gene NA 14.83062 NA      NA NA
>20702    36     14  23 104223       #NA   Gene NA 14.81930 NA      NA NA
>4270      8     22  10  94545  SERPINA6   Gene NA 14.79895 NA      NA NA
>6429     12     21   4  96944    MAD1L1   Gene NA 14.79082 NA      NA NA
>
>In all the 3 tables the mean 'A' value is the mean of the 'A' values 
>from the 4 slides where I expected for the second and third table 
>(coef=2 and coef=3) a mean between only 2 slides (slides 2 and 4 for 
>6h (coef=2) and slides 1 and 3 for 24h (coef=3)). I first though it 
>was just a problem in the report of 'A' values in topTable but it 
>seems to have an effect also on the P.value and B value.
>
>Is that a mistake in 'lmFit' or 'contrasts.fit' functions or am I 
>doing something wrong in the design of my matrix and contrast matrix?
>
>Many thanks for any help or comments.
>
>David
>
>
>######################################
>David HOT, PhD.
>Laboratoire d'Etudes Transcriptomiques et G?nomiques Appliqu?es - 
>Plate-forme Biopuces
>UMR8161 - IFR142
>
>Institut Pasteur de Lille
>1, rue du professeur Calmette
>59 019 LILLE, Cedex.
>FRANCE
>T?l : +33 (0)3 20 87 72 09
>Fax : +33 (0)3 20 87 73 11



More information about the Bioconductor mailing list