[BioC] limma/multtest comparison question

[Sean Davis]

Dear all,

I have a limma/multtest question.  I have performed what I think is the equivalent of a t-test (of course, moderated) in limma below.  Compare that with the output from mt.maxT in the multtest package.  Why are there such large discrepancies (ie., what am I doing wrong), particularly in p-values?  The data (abinorms3) are single-channel log2 intensities.  The commands and output are below.  We have the same samples hybridized on two-color arrays, also analyzed with limma, with p-values in mt.maxT and limma that agree with the mt.maxT result from below, so I think there are differentially-expressed genes.

Thanks,
Sean

> colnames(cma1) <- c("U","I")

> dma1 <- makeContrasts(U-I,levels=cma1)

> cma1

U I

1 1 0

2 0 1

3 0 1

4 1 0

5 0 1

6 1 0

8 0 1

9 1 0

10 0 1

> cl <- cma1[,2]

> tmp3 <- mt.maxT(abinorms3,cl)

We'll do complete enumerations

We're doing 126 complete permutations

b=1 b=2 b=3 b=4 b=5 b=6 b=7 b=8 

b=9 b=10 b=11 b=12 b=13 b=14 b=15 b=16 b=17 b=18 

b=19 b=20 b=21 b=22 b=23 b=24 b=25 b=26 b=27 b=28 

b=29 b=30 b=31 b=32 b=33 b=34 b=35 b=36 b=37 b=38 

b=39 b=40 b=41 b=42 b=43 b=44 b=45 b=46 b=47 b=48 

b=49 b=50 b=51 b=52 b=53 b=54 b=55 b=56 b=57 b=58 

b=59 b=60 b=61 b=62 b=63 b=64 b=65 b=66 b=67 b=68 

b=69 b=70 b=71 b=72 b=73 b=74 b=75 b=76 b=77 b=78 

b=79 b=80 b=81 b=82 b=83 b=84 b=85 b=86 b=87 b=88 

b=89 b=90 b=91 b=92 b=93 b=94 b=95 b=96 b=97 b=98 

b=99 b=100 b=101 b=102 b=103 b=104 b=105 b=106 b=107 b=108 

b=109 b=110 b=111 b=112 b=113 b=114 b=115 b=116 b=117 b=118 

b=119 b=120 b=121 b=122 b=123 b=124 b=125 b=126 

> tmp3[order(tmp3$rawp)[1:10],]

index teststat rawp adjp

5928 5928 -9.184354 0.007936508 0.3809524

8064 8064 7.712025 0.007936508 0.6190476

2878 2878 -7.574579 0.007936508 0.6825397

2503 2503 -7.491922 0.007936508 0.6984127

10283 10283 -7.415540 0.007936508 0.6984127

7927 7927 -7.213376 0.007936508 0.7222222

4638 4638 -7.091922 0.007936508 0.7301587

105 105 -7.076814 0.007936508 0.7460317

1074 1074 -7.061793 0.007936508 0.7539683

17870 17870 6.967126 0.007936508 0.7777778

> fit1s <- lmFit(abinorms3,cma1)

> fit2s <- contrasts.fit(fit1s,dma1)

> fit3s <- eBayes(fit2s)

> topTable(fit3s)

M t P.Value B

8064 -3.155503 -7.695828 0.2987386 1.52861663

105 1.523110 6.053780 1.0000000 0.46775473

2865 1.512303 5.764098 1.0000000 0.23825116

16035 3.125686 5.754778 1.0000000 0.23063006

1074 1.394530 5.602738 1.0000000 0.10416651

34 1.428812 5.595450 1.0000000 0.09800220

10283 1.316832 5.591068 1.0000000 0.09429127

4638 1.365565 5.530195 1.0000000 0.04239040

8571 -4.982818 -5.516468 1.0000000 0.03059528

8588 1.430188 5.472810 1.0000000 -0.00714332


	[[alternative HTML version deleted]]