Dear BioC users, I am using 2 color-spotted microarrays. My samples are queen ants with different genotypes (BB=D, Bb=H, bb=R), social forms (Monogyne, Polygyne) and ages (young virgin queen=2d, mature virgin queen=11d, and mated/reproductive/ mother queen=mom). They are hybridized to the reference RNA. I would like to analyze my data by making model matrix like y~ Age + Genotype *nested within* Social form + Age : Genotype *nested within* Social form + fixed factor (Batch) I have tried to analyze my data as show below and I got many questions. 1. Differences between different commands of model.matrix. From the examples below, I know that ~0 or ~1 is just whether I want the intercept or not but when there is an intercept, one of the groups is disappeared from the model (in this case, batchI or M11dD is disappeared depending on which factor comes first). Would this affect the result? Is the order of the factors in the model important? > design1 <- model.matrix(~0+factor(targets$Batch)+factor(targets$Age) %in%factor(targets$Genotype)%in%factor(targets$SocialForm)) > colnames(design1) [1] "factor(targets$Batch)I" [2] "factor(targets$Batch)J" [3] "factor(targets$Age)11d:factor(targets$Genotype)D:factor(targets $SocialForm)M" [4] "factor(targets$Age)2d:factor(targets$Genotype)D:factor(targets $SocialForm)M" [5] "factor(targets$Age)momQ:factor(targets$Genotype)D:factor(targets $SocialForm)M" [6] "factor(targets$Age)11d:factor(targets$Genotype)H:factor(targets $SocialForm)M" [7] "factor(targets$Age)2d:factor(targets$Genotype)H:factor(targets $SocialForm)M" [8] "factor(targets$Age)momQ:factor(targets$Genotype)H:factor(targets $SocialForm)M" [9] "factor(targets$Age)11d:factor(targets$Genotype)R:factor(targets $SocialForm)M" [10] "factor(targets$Age)2d:factor(targets$Genotype)R:factor(targets $SocialForm)M" [11] "factor(targets$Age)momQ:factor(targets$Genotype)R:factor(targets $SocialForm)M" [12] "factor(targets$Age)11d:factor(targets$Genotype)D:factor(targets $SocialForm)P" [13] "factor(targets$Age)2d:factor(targets$Genotype)D:factor(targets $SocialForm)P" [14] "factor(targets$Age)momQ:factor(targets$Genotype)D:factor(targets $SocialForm)P" [15] "factor(targets$Age)11d:factor(targets$Genotype)H:factor(targets $SocialForm)P" [16] "factor(targets$Age)2d:factor(targets$Genotype)H:factor(targets $SocialForm)P" [17] "factor(targets$Age)momQ:factor(targets$Genotype)H:factor(targets $SocialForm)P" [18] "factor(targets$Age)11d:factor(targets$Genotype)R:factor(targets $SocialForm)P" [19] "factor(targets$Age)2d:factor(targets$Genotype)R:factor(targets $SocialForm)P" [20] "factor(targets$Age)momQ:factor(targets$Genotype)R:factor(targets $SocialForm)P" > design1.1 <- model.matrix(~1+factor(targets$Batch)+factor(targets $Age)%in%factor(targets$Genotype)%in%factor(targets$SocialForm)) > colnames(design1.1) [1] "(Intercept)" [2] "factor(targets$Batch)J" [3] "factor(targets$Age)11d:factor(targets$Genotype)D:factor(targets $SocialForm)M" [4] "factor(targets$Age)2d:factor(targets$Genotype)D:factor(targets $SocialForm)M" [5] "factor(targets$Age)momQ:factor(targets$Genotype)D:factor(targets $SocialForm)M" [6] "factor(targets$Age)11d:factor(targets$Genotype)H:factor(targets $SocialForm)M" [7] "factor(targets$Age)2d:factor(targets$Genotype)H:factor(targets $SocialForm)M" [8] "factor(targets$Age)momQ:factor(targets$Genotype)H:factor(targets $SocialForm)M" [9] "factor(targets$Age)11d:factor(targets$Genotype)R:factor(targets $SocialForm)M" [10] "factor(targets$Age)2d:factor(targets$Genotype)R:factor(targets $SocialForm)M" [11] "factor(targets$Age)momQ:factor(targets$Genotype)R:factor(targets $SocialForm)M" [12] "factor(targets$Age)11d:factor(targets$Genotype)D:factor(targets $SocialForm)P" [13] "factor(targets$Age)2d:factor(targets$Genotype)D:factor(targets $SocialForm)P" [14] "factor(targets$Age)momQ:factor(targets$Genotype)D:factor(targets $SocialForm)P" [15] "factor(targets$Age)11d:factor(targets$Genotype)H:factor(targets $SocialForm)P" [16] "factor(targets$Age)2d:factor(targets$Genotype)H:factor(targets $SocialForm)P" [17] "factor(targets$Age)momQ:factor(targets$Genotype)H:factor(targets $SocialForm)P" [18] "factor(targets$Age)11d:factor(targets$Genotype)R:factor(targets $SocialForm)P" [19] "factor(targets$Age)2d:factor(targets$Genotype)R:factor(targets $SocialForm)P" [20] "factor(targets$Age)momQ:factor(targets$Genotype)R:factor(targets $SocialForm)P" > designOld <- model.matrix(~0+factor(targets$Cy3)+factor(targets $Batch)) > colnames(designOld) [1] "factor(targets$Cy3)M11dD" "factor(targets$Cy3)M2dD" "factor(targets$Cy3)MomD" "factor(targets$Cy3)MomH" [5] "factor(targets$Cy3)P11dD" "factor(targets$Cy3)P11dH" "factor(targets$Cy3)P11dR" "factor(targets$Cy3)P2dD" [9] "factor(targets$Cy3)P2dH" "factor(targets$Batch)J" > > designOld1 <- model.matrix(~1+factor(targets$Cy3)+factor(targets $Batch)) > colnames(designOld1) [1] "(Intercept)" "factor(targets$Cy3)M2dD" "factor(targets$Cy3)MomD" "factor(targets$Cy3)MomH" [5] "factor(targets$Cy3)P11dD" "factor(targets$Cy3)P11dH" "factor(targets$Cy3)P11dR" "factor(targets$Cy3)P2dD" [9] "factor(targets$Cy3)P2dH" "factor(targets$Batch)J" > designOldx <- model.matrix(~0+factor(targets$Batch)+factor(targets $Cy3)) > colnames(designOldx) [1] "factor(targets$Batch)I" "factor(targets$Batch)J" "factor(targets$Cy3)M2dD" "factor(targets$Cy3)MomD" [5] "factor(targets$Cy3)MomH" "factor(targets$Cy3)P11dD" "factor(targets$Cy3)P11dH" "factor(targets$Cy3)P11dR" [9] "factor(targets$Cy3)P2dD" "factor(targets$Cy3)P2dH" > designOld1 <- model.matrix(~1+factor(targets$Batch)+factor(targets $Cy3)) > colnames(designOld1) [1] "(Intercept)" "factor(targets$Batch)J" "factor(targets$Cy3)M2dD" "factor(targets$Cy3)MomD" [5] "factor(targets$Cy3)MomH" "factor(targets$Cy3)P11dD" "factor(targets$Cy3)P11dH" "factor(targets$Cy3)P11dR" [9] "factor(targets$Cy3)P2dD" "factor(targets$Cy3)P2dH" 2. Although I have 3 different genotypes, 2 social forms and 3 time points, some combinations are biologically impossible, e.g. M11dH, M2dH, MmomH, M11dR, M2dR, MmomR, PmomD, PmomR. When I use the command below, limma gives me all of the combinations but only some do exist. So, the ones with do not exist, have only 0 value in all rows. These groups will result in NA values for all the spots after the model fit. Is it OK to have NA values or it is better to remove these groups from the model matrix before fitting the model? design1 <- model.matrix(~0+factor(targets$Batch)+factor(targets$Age)%in %factor(targets$Genotype)%in%factor(targets$SocialForm)) > design1 BatchI BatchJ M11dD M2dD MmomD M11dH M2dH MmomH M11dR M2dR MmomR P11dD P2dD PmomD P11dH P2dH PmomH P11dR P2dR PmomR 1 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 . . . 99 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 attr(,"assign") [1] 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 attr(,"contrasts") attr(,"contrasts")$`factor(targets$Batch)` [1] "contr.treatment" attr(,"contrasts")$`factor(targets$Age)` [1] "contr.treatment" attr(,"contrasts")$`factor(targets$Genotype)` [1] "contr.treatment" attr(,"contrasts")$`factor(targets$SocialForm)` [1] "contr.treatment" 3. Problem with model fit After fitting the model. I got NA values for all spots of all the groups that do not exist. This could be normal, but the group "P11dR" do exist. Why did I also get NAs from this group? isGene <- MA$genes$Status == "cDNA" > fitDesign1 <- lmFit(MA[isGene,], design1) Coefficients not estimable: M11dH M2dH MmomH M11dR M2dR MmomR PmomD P11dR P2dR PmomR Warning message: Partial NA coefficients for 23523 probe(s) > fitDesign1 An object of class "MArrayLM" $coefficients BatchI BatchJ M11dD M2dD MmomD M11dH M2dH MmomH M11dR M2dR MmomR P11dD P2dD [1,] 0.1905989 0.0814195 -0.01294777 -0.26631768 1.85216328 NA NA NA NA NA NA 0.25491128 0.201785033 [2,] 0.6172982 0.8268361 -0.04929359 0.01953275 -0.31916109 NA NA NA NA NA NA -0.01204587 -0.008850553 [3,] 2.1798376 2.0604288 0.51830689 0.32399634 -0.03122528 NA NA NA NA NA NA -0.32814253 0.042235511 [4,] 1.0773968 0.9888733 0.12482602 0.20639996 -0.91445161 NA NA NA NA NA NA -0.41316719 -0.050187922 [5,] 1.0516807 0.4855826 -0.62717591 -0.21835907 -2.16438471 NA NA NA NA NA NA -0.75978080 -0.231660818 PmomD P11dH P2dH PmomH P11dR P2dR PmomR [1,] NA -0.12924261 -0.10818149 2.0982355 NA NA NA [2,] NA -0.18572498 -0.04461951 -0.5751055 NA NA NA [3,] NA 0.01181522 0.37474136 -0.2172751 NA NA NA [4,] NA -0.21312225 -0.07978014 -0.3549516 NA NA NA [5,] NA -0.38056336 -0.16362561 -1.7943584 NA NA NA 23526 more rows ... $stdev.unscaled BatchI BatchJ M11dD M2dD MmomD M11dH M2dH MmomH M11dR M2dR MmomR P11dD P2dD PmomD [1,] 0.4494704 0.4820557 0.5740672 0.5740672 0.6125642 NA NA NA NA NA NA 0.5218120 0.5198694 NA [2,] 0.5831241 0.6001124 0.6782392 0.6782392 0.7092194 NA NA NA NA NA NA 0.6534186 0.6520187 NA [3,] 0.4494973 0.4824567 0.5741145 0.5741145 0.6128798 NA NA NA NA NA NA 0.5253914 0.5199107 NA [4,] 0.5829695 0.5995112 0.6782060 0.6782060 0.7087107 NA NA NA NA NA NA 0.6533627 0.6458128 NA [5,] 0.7071068 0.7580744 0.8022900 0.8102921 0.8610130 NA NA NA NA NA NA 0.7911144 0.7870251 NA P11dH P2dH PmomH P11dR P2dR PmomR [1,] 0.5253396 0.5167718 0.5978108 NA NA NA [2,] 0.6412859 0.6424903 0.7484216 NA NA NA [3,] 0.5253914 0.5168244 0.6320057 NA NA NA [4,] 0.6412651 0.6374888 0.7479396 NA NA NA [5,] 0.7786772 0.7786772 0.9528956 NA NA NA 23526 more rows ... $sigma [1] 0.4750674 0.2759988 0.7108862 0.3199147 0.2901134 23526 more elements ... $df.residual [1] 86 69 83 71 56 23526 more elements ... $cov.coefficients BatchI BatchJ M11dD M2dD MmomD P11dD P2dD P11dH P2dH PmomH BatchI 0.2019737 0.1921053 -0.1970395 -0.1970395 -0.1921053 -0.1967105 -0.1970395 -0.1967105 -0.1970395 -0.1921053 BatchJ 0.1921053 0.2315789 -0.2118421 -0.2118421 -0.2315789 -0.2131579 -0.2118421 -0.2131579 -0.2118421 -0.2315789 M11dD -0.1970395 -0.2118421 0.3294408 0.2044408 0.2118421 0.2049342 0.2044408 0.2049342 0.2044408 0.2118421 M2dD -0.1970395 -0.2118421 0.2044408 0.3294408 0.2118421 0.2049342 0.2044408 0.2049342 0.2044408 0.2118421 MmomD -0.1921053 -0.2315789 0.2118421 0.2118421 0.3565789 0.2131579 0.2118421 0.2131579 0.2118421 0.2315789 P11dD -0.1967105 -0.2131579 0.2049342 0.2049342 0.2131579 0.2721491 0.2049342 0.2054825 0.2049342 0.2131579 P2dD -0.1970395 -0.2118421 0.2044408 0.2044408 0.2118421 0.2049342 0.2669408 0.2049342 0.2044408 0.2118421 P11dH -0.1967105 -0.2131579 0.2049342 0.2049342 0.2131579 0.2054825 0.2049342 0.2721491 0.2049342 0.2131579 P2dH -0.1970395 -0.2118421 0.2044408 0.2044408 0.2118421 0.2049342 0.2044408 0.2049342 0.2669408 0.2118421 PmomH -0.1921053 -0.2315789 0.2118421 0.2118421 0.2315789 0.2131579 0.2118421 0.2131579 0.2118421 0.3565789 $pivot [1] 1 2 3 4 5 12 13 15 16 17 6 7 8 9 10 11 14 18 19 20 $genes Block Row Column ID Name Status 9 1 1 9 SiJWA06AAB_pcr_2 NA cDNA 10 1 1 10 SiJWA12AAB_pcr_2 NA cDNA 11 1 1 11 SiJWC06AAB_pcr_2 NA cDNA 12 1 1 12 SiJWC12AAB_pcr_2 NA cDNA 13 1 1 13 SiJWE06AAB_pcr_2 NA cDNA 23526 more rows ... $Amean [1] 9.465811 8.056152 8.586644 8.151375 7.580591 23526 more elements ... $method [1] "ls" $design BatchI BatchJ M11dD M2dD MmomD M11dH M2dH MmomH M11dR M2dR MmomR P11dD P2dD PmomD P11dH P2dH PmomH P11dR P2dR PmomR 1 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 5 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 94 more rows ... 4. After fitting the model, do I need to make contrasts between interesting comparisons? Although all slides are compared to reference, with this complicated linear model, I am not sure that if I do not make contrast, the limma will give me the coefficient of each group compared to reference or something else. I truly appreciate all the answers in advance! best regards, Mingkwan Nipitwattanaphon [[alternative HTML version deleted]]