Hi Jim, Thank you so much for your time and input. I reckon we should try going with a single channel analysis... Hopefully we can still get out of the "corner". Best, Joaquin On Thu, Mar 8, 2012 at 11:25 AM, James W. MacDonald <jmacdon@uw.edu> wrote: > Hi Joaquin, > > I think you may have painted yourself into a corner with this experiment. > You have done something similar to the technical replication examples in > section 8.2, but the problem is that you haven't done enough replicates to > analyze that way. Note that in the examples in section 8.2, there are at > least two identical (or dye swap) slides for each comparison, whereas you > only have one. > > The experiment at the bottom of p.40 only has mu and wt samples, with > three technical replicates each, and they used 18 slides! To do something > similar, you would need 54 slides, not 21. > > If this were my experiment to analyze, I would probably go with a single > channel analysis with either a blocking variable or duplicateCorrelation() > to account for the technical replicates. This is a non-trivial analysis and > I would highly recommend getting some help from a local statistician, > preferably one with microarray experience. > > Best, > > Jim > > > > > On 3/7/2012 6:26 PM, Joaquin Martinez wrote: > >> Hi Jim, >> >> Once again, thank you for your helpful advice. Maybe I was not very >> clear in my original email. In the targets file that I sent, we have both >> biological replication (e.g. Puninf.1, Puninf.2 and Puninf.3) and technical >> replication (e.g. Puninf.1 present in slides29, 30 and 31) instead of >> including a common reference sample in all hybridizations. Ours is more >> like the "Direct Two-Color Designs" in section 7.4 including "Technical >> Replication" as in section 8.2. Based on this, does your previous advice >> still hold? Should we remove the numbers from the sample names? Or are >> there any further considerations we should take into account in our >> analysis? >> >> I truly appreciate your help. >> >> Best, >> Joaquin >> >> On Wed, Mar 7, 2012 at 3:07 PM, James W. MacDonald <jmacdon@uw.edu<mailto:>> jmacdon@uw.edu>> wrote: >> >> Hi Joaquin, >> >> >> On 3/7/2012 2:12 PM, Joaquin Martinez wrote: >> >> Hi Jim, >> >> Thank you for your quick reply and advice. Actually we had >> considered doing exactly what you suggest and remove the >> numbers from the sample names. My concern then is what is then >> the reference or baseline?, i.e. Puninf.1 or an average value >> of Puninf.1, Puninf.2 and Puninf.3. Also, how does removing >> numbers take into consideration slight value differences among >> biological replicate samples? >> >> >> The baseline is the average of the Puninf samples. Which is, I >> believe, what you want. To make the comparisons you request, limma >> is going to compute a t-statistic, which simply put is >> >> difference in means between two samples/ measure of the >> variability of the means >> >> So the numerator of the statistic tells you how different the two >> samples are. The denominator (which is a measure of the slight >> differences among biological replicates) is used to determine if >> the differences between samples is 'large' or not. If there is a >> lot of variability within sample types, the numerator has to be >> pretty big. Conversely, if there isn't much variability within >> sample types then a smaller numerator may be considered >> significantly large. >> >> Note that you aren't forced to use Puninf as the reference. That's >> just the easiest way to use modelMatrix(). You can also pass in a >> matrix for the 'param' argument, selecting direct comparisons to >> make. Or you can just leave Puninf as the reference, and make >> further comparisons with a contrast matrix, which is probably simpler. >> >> Note here that the design matrix with Puninf as reference is >> computing implicit comparisons (e.g., the M163 coefficient is >> really log2(M163) - log2(Puninf)). So if you want to know if there >> is a significant difference between M163 and M86, you can just add >> a contrast matrix that looks like >> >> 1 >> -1 >> 0 >> 0 >> >> Which simply means that we will be computing log2(M163) - >> log2(Puninf) - (log2(M86) - log2(Puninf)). Algebraically, the two >> log2(Puninf) will cancel out, and you end up with log2(M163) - >> log2(M86). You can create any contrast matrix you want using >> makeContrasts(). >> >> >> >> In the example we are following in the limma user's guide >> there are 3 wild-type mice (wt1,wt2,wt3) and 3 mutant mice >> (mu1, mu2, mu3), where wt1 is used as the reference. I think >> we have the same type of study but with 6 "mouse types". >> >> >> Unless I am mistaken, you are looking at the wrong example. The >> example you reference involves technical replication, where RNA >> from each mouse was hybridized to more than one slide. From what >> you said originally, each of the samples in your experiment is a >> biological replicate. Your experiment is more similar to the >> common reference design experiment on p34 of the limma User's Guide. >> >> Best, >> >> Jim >> >> >> >> Thanks, >> Joaquin >> >> >> >> >> >> On Wed, Mar 7, 2012 at 12:01 PM, James W. MacDonald >> <jmacdon@uw.edu <mailto:jmacdon@uw.edu=""> <mailto:jmacdon@uw.edu>> >> <mailto:jmacdon@uw.edu>>> wrote: >> >> Hi Joaquin, >> >> >> On 3/7/2012 11:08 AM, Joaquin Martinez wrote: >> >> Dear All, >> >> My colleague Ben and I have fitted a linear model to >> microarray data using >> the limma package but we get "coefficients not >> estimable" for >> the last two >> columns in the design matrix during the fit. As a >> consequence >> when trying >> to fit our contrast matrix we get the following error >> message, >> "Error in >> contrasts.fit(fit, contrasts = A) : trying to take >> contrast of >> non-estimable coefficient." We have discovered that if we >> simply reorder >> the design matrix columns and contrasts matrix rows, we >> encounter the same >> error. We would very much appreciate if someone could >> explain >> to us why the >> not-estimable coefficients arise, and how to correct >> the problem. >> >> We are following the guide (section 8.2) for limma 3.10.0 >> using this >> version of R 2.14.0 (sessionInfo at end of email). We >> have 6 >> different >> types of samples (Puninf, P86, P163, Muninf, M86, >> M163), and 3 >> biological >> replicates of each (1,2,3; 4,5,6; 7,8,9; 10,11,12; >> 13,14,15; >> 16,17,18). >> >> >> Yes, but your targets file has numbers appended to the sample >> types, which makes R think they are different. >> Hypothetically you >> would have noticed this when using the modelMatrix() function: >> >> >> > modelMatrix(targets, ref = "Puninf.1") >> Found unique target names: >> M163.16 M163.17 M163.18 M86.13 M86.14 M86.15 Muninf.10 >> Muninf.11 >> Muninf.12 P163.7 P163.8 P163.9 P86.4 P86.5 P86.6 Puninf.1 >> Puninf.2 >> Puninf.3 >> >> This indicates that limma thinks all of your replicates are >> different sample types. This won't work out (as you have >> already >> found) because the model will be over specified, which simply >> means that you cannot estimate all the parameters you have >> in the >> model. This is the same as saying you can't determine the >> values >> of x and y in the formula x + y = 3 (but if you also know >> that x >> - y = -1, then you can solve for both x and y). >> >> Now back to the problem at hand. If you remove the numbers from >> your sample types, >> >> targets$Cy3 <- sapply(strsplit(targets$Cy3, "\\."), >> function(x) x[1]) >> targets$Cy5 <- sapply(strsplit(targets$Cy5, "\\."), >> function(x) x[1]) >> >> then you will create the correct model matrix >> >> > modelMatrix(targets, ref="Puninf") >> Found unique target names: >> M163 M86 Muninf P163 P86 Puninf >> M163 M86 Muninf P163 P86 >> [1,] 0 0 1 0 0 >> [2,] 0 0 0 0 1 >> [3,] 0 0 0 -1 0 >> [4,] 0 -1 1 0 0 >> [5,] 1 0 -1 0 0 >> [6,] 0 1 0 0 -1 >> [7,] -1 0 0 1 0 >> [8,] 0 0 -1 0 0 >> [9,] 0 0 0 0 -1 >> [10,] 0 0 0 1 0 >> [11,] 0 1 -1 0 0 >> [12,] -1 0 1 0 0 >> [13,] 0 -1 0 0 1 >> [14,] 1 0 0 -1 0 >> [15,] 0 0 1 0 0 >> [16,] 0 0 0 0 1 >> [17,] 0 0 0 -1 0 >> [18,] 0 -1 1 0 0 >> [19,] 1 0 -1 0 0 >> [20,] 0 1 0 0 -1 >> [21,] -1 0 0 1 0 >> >> Which will compare all samples to Puninf. >> >> Best, >> >> Jim >> >> >> >> >> >> >> # The targets file looks like this >> >> targets<- structure(list(SlideNumber = c(29L, 30L, 31L, >> 32L, >> 33L, 34L, >> 35L, 36L, 37L, 38L, 39L, 40L, 41L, 42L, 43L, 44L, >> 51L, 46L, >> 47L, >> 48L, 50L), FileName = c("BE T22h slide 29", "BE T22h >> slide 30", >> "BE T22h slide 31", "BE T22h slide 32", "BE T22h >> slide 33", >> "BE T22h >> slide 34", >> "BE T22h slide 35", "BE T22h slide 36", "BE T22h >> slide 37", >> "BE T22h >> slide 38", >> "BE T22h slide 39", "BE T22h slide 40", "BE T22h >> slide 41", >> "BE T22h >> slide 42", >> "BE T22h slide 43", "BE T22h slide 44", "BE T22h >> slide 51", >> "BE T22h >> slide 46", >> "BE T22h slide 47", "BE T22h slide 48", "BE T22h >> slide 50"), >> Cy3 = c("Puninf.1", "Puninf.1", "P163.7", "M86.13", >> "Muninf.10", >> "P86.4", "M163.16", "Muninf.11", "P86.5", "Puninf.2", >> "Muninf.11", >> "M163.17", "M86.14", "P163.8", "Puninf.3", "Puninf.3", >> "P163.9", >> "M86.15", "Muninf.12", "P86.6", "M163.18"), Cy5 = >> c("Muninf.10", >> "P86.4", "Puninf.1", "Muninf.10", "M163.16", "M86.13", >> "P163.7", >> "Puninf.2", "Puninf.2", "P163.8", "M86.14", >> "Muninf.11", >> "P86.5", "M163.17", "Muninf.12", "P86.6", "Puninf.3", >> "Muninf.12", >> "M163.18", "M86.15", "P163.9")), .Names = >> c("SlideNumber", >> "FileName", "Cy3", "Cy5"), class = "data.frame", >> row.names >> = c(NA, >> -21L)) >> >> # for readability we also paste it in here... >> #> targets >> # SlideNumber FileName Cy3 Cy5 >> # 1 29 BE T22h slide 29 Puninf.1 Muninf.10 >> # 2 30 BE T22h slide 30 Puninf.1 P86.4 >> # 3 31 BE T22h slide 31 P163.7 Puninf.1 >> # 4 32 BE T22h slide 32 M86.13 Muninf.10 >> # 5 33 BE T22h slide 33 Muninf.10 M163.16 >> # 6 34 BE T22h slide 34 P86.4 M86.13 >> # 7 35 BE T22h slide 35 M163.16 P163.7 >> # 8 36 BE T22h slide 36 Muninf.11 Puninf.2 >> # 9 37 BE T22h slide 37 P86.5 Puninf.2 >> # 10 38 BE T22h slide 38 Puninf.2 P163.8 >> # 11 39 BE T22h slide 39 Muninf.11 M86.14 >> # 12 40 BE T22h slide 40 M163.17 Muninf.11 >> # 13 41 BE T22h slide 41 M86.14 P86.5 >> # 14 42 BE T22h slide 42 P163.8 M163.17 >> # 15 43 BE T22h slide 43 Puninf.3 Muninf.12 >> # 16 44 BE T22h slide 44 Puninf.3 P86.6 >> # 17 51 BE T22h slide 51 P163.9 Puninf.3 >> # 18 46 BE T22h slide 46 M86.15 Muninf.12 >> # 19 47 BE T22h slide 47 Muninf.12 M163.18 >> # 20 48 BE T22h slide 48 P86.6 M86.15 >> # 21 50 BE T22h slide 50 M163.18 P163.9 >> >> # The design matrix is computed like this design<- >> modelMatrix(targets, >> ref = "Puninf.1") >> # but here is the structure if it is helpful >> design<- structure(c(0, 0, 0, 0, 1, 0, -1, 0, 0, 0, 0, >> 0, 0, >> 0, 0, 0, >> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, >> 0, 1, 0, >> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, >> 0, 0, 0, 0, >> 0, 0, 0, 1, 0, -1, 0, 0, 0, -1, 0, 1, 0, 0, 0, 0, 0, >> 0, 0, 0, >> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, >> 1, 0, -1, >> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, >> 0, 0, 0, 0, >> 0, 0, 0, 0, -1, 0, 1, 0, 1, 0, 0, 1, -1, 0, 0, 0, 0, >> 0, 0, 0, >> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, >> 0, 0, -1, >> 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, >> 0, 0, 0, 0, >> 0, 0, 0, 1, 0, 0, 1, -1, 0, 0, 0, 0, -1, 0, 0, 0, 1, >> 0, 0, 0, >> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, >> 0, 0, 0, 1, >> 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, >> 0, 0, 0, >> 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 1, 0, 1, 0, 0, 0, >> -1, 0, 0, >> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, >> 0, 0, 0, 0, >> -1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, >> 0, 0, 0, >> 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, -1, 0, 0, 0, 0, >> 0, 0, 0, >> 0, 1, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, >> 0, 0, 0, >> 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, -1, 1, 0, 0, 0, 0), >> .Dim = >> c(21L, >> 17L), .Dimnames = list(c("slide 29", "slide 30", >> "slide 31", >> "slide 32", "slide 33", "slide 34", "slide 35", >> "slide 36", >> "slide 37", >> "slide 38", "slide 39", "slide 40", "slide 41", >> "slide 42", >> "slide 43", >> "slide 44", "slide 51", "slide 46", "slide 47", >> "slide 48", >> "slide 50" >> ), c("M163.16", "M163.17", "M163.18", "M86.13", >> "M86.14", >> "M86.15", >> "Muninf.10", "Muninf.11", "Muninf.12", "P163.7", >> "P163.8", >> "P163.9", >> "P86.4", "P86.5", "P86.6", "Puninf.2", "Puninf.3"))) >> >> >> # Fit linear model - commented out becuase we can't really >> send all the >> data in MA >> # fit<- lmFit(MA, design = design) >> # Coefficients not estimable: Puninf.2 Puninf.3 >> >> # Make a contrast matrix: >> A<- >> makeContrasts(PuninfvsMuninf=(** >> Muninf.10+Muninf.11+Muninf.12-**Puninf.2-Puninf.3)/3, >> levels = design) >> >> # and now fit the contrast >> con.fit<- contrasts.fit(fit, contrasts = A) >> Error in contrasts.fit(fit, contrasts = A) : trying to take >> contrast of >> non-estimable coefficient >> >> # now change the column order in the design matrix and row >> order of >> contrast matrix >> # we create a new sort index - arbitrarily placing the last >> two columns in >> the original >> # design matrix in the middle of the new one. Ditto >> for the >> rows of the >> contrast matrix. >> n<- nrow(A) >> newIndex<- c(1:5, n-1, n, 6:(n-2)) >> designB<- design[,newIndex] >> B<- A[newIndex,] >> attributes(B)<- attributes(A) >> >> # create a new fit, but we get the same knd of error, >> the last >> two columns >> are >> # identified as non-estimable >> fitB<- lmFit(X$MA, design = designB) >> #Coefficients not estimable: P86.5 P86.6 >> #Warning message: >> #Partial NA coefficients for 16278 probe(s) >> >> >> As far as we can tell, it's the trailing columns that get >> singled out, but >> we don't know why. Hopefully the information about is >> detailed >> enough, but >> we can provide more information if necessary to help >> troubleshooting. >> >> Thanks in advance, >> >> Joaquin >> >> sessionInfo() >> >> R version 2.14.0 (2011-10-31) >> Platform: i386-apple-darwin9.8.0/i386 (32-bit) >> >> locale: >> [1] >> en_US.UTF-8/en_US.UTF-8/en_US.**UTF-8/C/en_US.UTF-8/en_US.UTF-**8 >> >> attached base packages: >> [1] stats graphics grDevices utils datasets >> methods >> base >> >> other attached packages: >> [1] limma_3.10.0 >> >> loaded via a namespace (and not attached): >> [1] tools_2.14.0 >> >> [[alternative HTML version deleted]] >> >> ______________________________**_________________ >> Bioconductor mailing list >> Bioconductor@r-project.org <mailto:bioconductor@r-**project.org<bioconductor@r-project.org> >> > >> <mailto:bioconductor@r-**project.org <bioconductor@r-project.org=""> >> >> <mailto:bioconductor@r-**project.org <bioconductor@r-project.org=""> >> >> >> >> https://stat.ethz.ch/mailman/**listinfo/bioconductor<https: stat.ethz.ch="" mailman="" listinfo="" bioconductor=""> >> Search the archives: >> http://news.gmane.org/gmane.**science.biology.informatics.** >> conductor<http: news.gmane.org="" gmane.science.biology.informatics.c="" onductor=""> >> >> >> -- James W. MacDonald, M.S. >> Biostatistician >> University of Washington >> Environmental and Occupational Health Sciences >> 4225 Roosevelt Way NE, # 100 >> Seattle WA 98105-6099 >> >> >> >> -- James W. MacDonald, M.S. >> Biostatistician >> University of Washington >> Environmental and Occupational Health Sciences >> 4225 Roosevelt Way NE, # 100 >> Seattle WA 98105-6099 >> >> >> > -- > James W. MacDonald, M.S. > Biostatistician > University of Washington > Environmental and Occupational Health Sciences > 4225 Roosevelt Way NE, # 100 > Seattle WA 98105-6099 > > -- Joaquin Martinez Martinez, Ph.D Post Doctoral Research Scientist Bigelow Laboratory for Ocean Sciences 180 McKown Point P.O. Box 475 West Boothbay Harbor, Maine 04575-0475 Office phone: 202-747-3255 EXT-403 e-mail:jmartinez@bigelow.org [[alternative HTML version deleted]]