Search
0
7.7 years ago by
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 <- eBayesfitContrast.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