[BioC] Questions about model matrix, logFC and adjusted P.value of t-test
Mingkwan.Nipitwattanaphon at unil.ch
Mingkwan.Nipitwattanaphon at unil.ch
Wed Oct 27 00:52:22 CEST 2010
Dear Bio-C users,
My samples are queen ants. I use two-color spotted
microarrays, hybridisation
against a
reference, no dye swaps. I would like to compare samples
with different
genotypes, developmental stages and social form.
There are 3 genotypes: 1. BB (D=dominant), 2. Bb
(H=heterozygous) and 3. bb
(R=recessive).
Three developmental stages: 1. young (2d) virgin queen, 2.
mature (11d) virgin
queen, 3. mated/mother queen (mom)
Two social forms: 1. Monogyne (M), 2, Polygyne (P).
>From these 3 genotypes, 3 developmental stages and 2 types
of social form, I can
group my 99 slides into 9 different categories. These slides
contain two batches
(I and J series). I treated batch effect as a fixed effect.
My model matrix is:
design <-
model.matrix(~0+factor(targets$Cy3)+factor(targets$batch))
colnames(design) <- c("P2dD", "P2dH", "P11dD", "P11dH",
"P11dR", "M2dD", "M11dD",
"MomD", "MomH", "batch")
design
P2dD P2dH P11dD P11dH P11dR M2dD M11dD MomD MomH batch
1 0 0 0 0 0 1 0 0 0 0
2 0 0 0 0 1 0 0 0 0 0
3 0 0 0 0 0 1 0 0 0 0
4 0 0 0 0 1 0 0 0 0 0
5 0 0 0 0 0 1 0 0 0 0
6 0 0 0 0 1 0 0 0 0 0
7 0 0 0 0 0 1 0 0 0 0
8 0 0 0 0 1 0 0 0 0 0
9 0 0 0 0 0 1 0 0 0 0
10 0 0 0 0 1 0 0 0 0 0
11 0 0 0 0 0 1 0 0 0 0
12 0 0 0 0 1 0 0 0 0 0
13 0 0 0 0 0 1 0 0 0 0
14 0 0 0 0 1 0 0 0 0 0
15 0 0 0 0 0 0 1 0 0 0
16 0 0 0 0 0 0 1 0 0 0
17 0 0 0 0 0 0 1 0 0 0
18 0 0 0 0 0 0 1 0 0 0
19 0 1 0 0 0 0 0 0 0 0
20 0 1 0 0 0 0 0 0 0 0
21 0 1 0 0 0 0 0 0 0 0
22 0 1 0 0 0 0 0 0 0 0
23 1 0 0 0 0 0 0 0 0 0
24 1 0 0 0 0 0 0 0 0 0
25 1 0 0 0 0 0 0 0 0 0
26 1 0 0 0 0 0 0 0 0 0
27 0 0 0 0 0 0 0 1 0 0
28 0 0 0 0 0 0 0 0 1 0
29 0 0 0 0 0 0 0 1 0 0
30 0 0 0 0 0 0 0 0 1 0
31 0 0 0 0 0 0 0 1 0 0
32 0 0 0 0 0 0 0 0 1 0
33 0 0 0 0 0 0 0 1 0 0
34 0 0 0 0 0 0 0 0 1 0
35 0 0 0 0 0 0 0 1 0 0
36 0 0 0 0 0 0 0 0 1 0
37 0 0 0 0 0 0 0 1 0 0
38 0 0 0 0 0 0 0 0 1 0
39 0 0 0 0 0 0 0 1 0 0
40 0 0 0 0 0 0 0 0 1 0
41 0 0 0 0 0 0 0 1 0 0
42 0 0 0 0 0 0 0 0 1 0
43 0 0 0 0 0 0 0 1 0 1
44 0 0 0 0 0 0 0 0 1 1
45 0 0 0 0 0 1 0 0 0 1
46 0 0 0 0 1 0 0 0 0 1
47 0 0 0 0 0 0 0 0 1 1
48 0 0 0 0 0 0 0 1 0 1
49 0 0 0 0 1 0 0 0 0 1
50 0 0 0 0 0 1 0 0 0 1
51 0 1 0 0 0 0 0 0 0 1
52 1 0 0 0 0 0 0 0 0 1
53 0 1 0 0 0 0 0 0 0 1
54 1 0 0 0 0 0 0 0 0 1
55 0 1 0 0 0 0 0 0 0 1
56 1 0 0 0 0 0 0 0 0 1
57 0 0 0 0 0 0 0 0 1 1
58 0 0 0 0 0 0 0 1 0 1
59 0 0 0 0 0 1 0 0 0 1
60 0 0 0 0 1 0 0 0 0 1
61 0 0 0 0 0 0 0 1 0 1
62 0 0 0 0 0 0 0 0 1 1
63 0 0 0 0 0 1 0 0 0 1
64 0 0 0 0 1 0 0 0 0 1
65 0 0 0 0 0 0 0 0 1 1
66 0 0 0 0 0 0 0 1 0 1
67 0 0 0 0 0 1 0 0 0 1
68 0 0 0 0 1 0 0 0 0 1
69 0 0 0 0 0 0 0 0 1 1
70 0 0 0 0 0 0 0 1 0 1
71 0 0 0 0 0 1 0 0 0 1
72 0 0 0 0 1 0 0 0 0 1
73 0 0 0 0 0 0 0 1 0 1
74 0 0 0 0 0 0 0 0 1 1
75 0 0 0 0 0 1 0 0 0 1
76 0 0 0 0 1 0 0 0 0 1
77 0 0 0 0 0 0 0 1 0 1
78 0 0 0 0 0 0 0 0 1 1
79 0 0 0 0 0 0 1 0 0 1
80 0 0 0 0 0 1 0 0 0 1
81 0 0 0 0 1 0 0 0 0 1
82 0 1 0 0 0 0 0 0 0 1
83 1 0 0 0 0 0 0 0 0 1
84 0 0 1 0 0 0 0 0 0 1
85 0 0 1 0 0 0 0 0 0 1
86 0 0 1 0 0 0 0 0 0 1
87 0 0 1 0 0 0 0 0 0 1
88 0 0 1 0 0 0 0 0 0 1
89 0 0 1 0 0 0 0 0 0 1
90 0 0 1 0 0 0 0 0 0 1
91 0 0 1 0 0 0 0 0 0 1
92 0 0 0 1 0 0 0 0 0 1
93 0 0 0 1 0 0 0 0 0 1
94 0 0 0 1 0 0 0 0 0 1
95 0 0 0 1 0 0 0 0 0 1
96 0 0 0 1 0 0 0 0 0 1
97 0 0 0 1 0 0 0 0 0 1
98 0 0 0 1 0 0 0 0 0 1
99 0 0 0 1 0 0 0 0 0 1
fitMA <- lmFit(normalizedMA, design)
There are 36 possible tests that I can make but I am only
interested in 16 tests
below.
contrastALL16 <- makeContrasts(P2dD-P2dH, P11dD-P11dH,
P11dH-P11dR, P11dD-P11dR,
M2dD-P2dD, M11dD-P11dD, MomD-MomH, P2dD-P11dD, P11dD-MomD,
P2dD-MomD, P2dH-P11dH,
P11dH-MomH, P2dH-MomH, M2dD-M11dD, M11dD-MomD, M2dD-MomD,
levels=design)
fitContrast.MA <- contrasts.fit(fitMA, contrastALL16)
fit_eBayes_MA <- eBayes(fitContrast.MA)
write.table(fit_eBayes_MA, file="Q16contrasts.txt",
sep="\t")
My result file contains coefficients of all the 16 contrasts
that I asked for and
the p-value of each contrast (each t-test) but NOT the
adjusted p-value. It also
gives me the F value and the p-value of the F-test and again
NOT the adjusted p-
value of the F-test.
I can get the adjusted p-value from the F-test by using the
command "p.adjust"
but not with the t-test. When I used the command "topTable"
with coeff=1 (until
16 each time for all of my 16 contrasts), I can get the
adjusted p-value of each
contrast.
My questions are:
1. Why does not the command "eBayes" give adjusted p-value?
Is there an easier or
more direct way to get adjusted p-value of the t-test?
2. How does the logFC calculate from?
If I took the M-values for a single spot, after
normalisation between arrays from
slides that are belonged to one of my contrasts (8 slides of
momD vs 8 slides of
momH, because this case all slides are from the same batch)
M-value of momD slide 1 to 8 = 3.00, 3.26, 2.73, 3,32, 2.93,
2.81, 2.55, 2.85
M-value of momH slide 1 to 8 = -0.44, -0.54, 0.03, -0.38,
0.49, -0.56, 0.07, 0.37
The mean M-value is -0.12 for momH and 2.93 for momD. I'd
expect that the
relative expression level in momD compared to momH would be:
(2^(2.93))/(2^(-
0.12)) = 8.29
The logFC should simply be: log2(8.29) = 3.05 However, the
logFC given by
limma is 3.37.
3. If I want to do more complicated model like:
~1 + Age + Geno type *nested within* Social form + Age :
Genotype *nested within*
Social form, (fixed factor = Batch)
Is it possible to do? How can I do it in limma?
I am really sorry for writing such a long email because I
want to make everything
clear. I really appreciate your help.
best regards,
Mingkwan
More information about the Bioconductor
mailing list