Dear Somnath, You are asking questions about variances, rather than questions about means; and questions about differential variability rather than questions about differential expression. You cannot easily answer these questions using limma, or I think any other Bioconductor package. Have a look at the following paper for the sort of methods that are relevant for analysing variability: http://www.ncbi.nlm.nih.gov/pubmed/21655321 These are statistically advanced methods, and you will need to find a very experienced statistician that you can work with. Best wishes Gordon --------------------------------------------- Professor Gordon K Smyth, Bioinformatics Division, Walter and Eliza Hall Institute of Medical Research, 1G Royal Parade, Parkville, Vic 3052, Australia. Tel: (03) 9345 2326, Fax (03) 9347 0852, smyth at wehi.edu.au http://www.wehi.edu.au http://www.statsci.org/smyth On Wed, 8 Feb 2012, somnath bandyopadhyay wrote: > > Hi Gordon, > > I tried model.matrix(~gender+gender:donor+day). I am still running into the same error. > > In terms of what I want to get out of this analysis is as follows: > > 1. Does gene expression vary significantly across donors on the same > day? Is the variation higher in males than in females? What are the > least varying genes across donors? > 2. Does gene expression vary significantly across days for a given > donor? Is the variation higher in males than in females? What are the > least varying genes across days? > 3. What are the least varying genes across donors, days and gender all > taken together? (more like house keeping genes) > > Any help would be greatly appreciated. > > What I tried: > > R version 2.10.1 (2009-12-14) > Copyright (C) 2009 The R Foundation for Statistical Computing > ISBN 3-900051-07-0 > R is free software and comes with ABSOLUTELY NO WARRANTY. > You are welcome to redistribute it under certain conditions. > Type 'license()' or 'licence()' for distribution details. > Natural language support but running in an English locale > R is a collaborative project with many contributors. > Type 'contributors()' for more information and > 'citation()' on how to cite R or R packages in publications. > Type 'demo()' for some demos, 'help()' for on-line help, or > 'help.start()' for an HTML browser interface to help. > Type 'q()' to quit R. >> library("limma") >> rna <- read.table('qced_non_redundant_rma_data.txt', header=T, row.names=1, sep= "\t") >> annotation <- read.table('Test_Samples_Annotation.txt', header=T, row.names=1, sep= "\t") >> dim(rna) > [1] 54675 36 >> ## check that the annotations and data match >> if (sum(sort(rownames(annotation)) == sort(colnames(rna))) != length(rownames(annotation))) { > + stop ("ERROR: annotation rownames and RMA colnames do not match") > + } >> donor <- factor(annotation$donor) >> day <- factor(annotation$day) >> gender <- factor(annotation$gender) >> design <- model.matrix(~ donor + gender:donor + day, data=annotation) >> fit <- lmFit(rna,design) > Coefficients not estimable: donorF2:gendermale donorF3:gendermale donorF5:gendermale donorF7:gendermale donorF8:gendermale donorM1:gendermale donorM4:gendermale donorM6:gendermale > Warning message: > Partial NA coefficients for 54675 probe(s) >> fit2 <- eBayes(fit) > Warning message: > In ebayes(fit = fit, proportion = proportion, stdev.coef.lim = stdev.coef.lim) : > Estimation of var.prior failed - set to default value >> write.fit(fit,file="effects.txt",digits=30, adjust="BH", method="separate",sep="\t") > Warning message: > In ebayes(fit = fit, proportion = proportion, stdev.coef.lim = stdev.coef.lim) : > Estimation of var.prior failed - set to default value >> sessionInfo() > R version 2.10.1 (2009-12-14) > x86_64-unknown-linux-gnu > locale: > [1] LC_CTYPE=en_US.iso885915 LC_NUMERIC=C > [3] LC_TIME=en_US.iso885915 LC_COLLATE=en_US.iso885915 > [5] LC_MONETARY=C LC_MESSAGES=en_US.iso885915 > [7] LC_PAPER=en_US.iso885915 LC_NAME=C > [9] LC_ADDRESS=C LC_TELEPHONE=C > [11] LC_MEASUREMENT=en_US.iso885915 LC_IDENTIFICATION=C > attached base packages: > [1] stats graphics grDevices utils datasets methods base > other attached packages: > [1] limma_3.2.3 >> > > > >> Date: Wed, 8 Feb 2012 14:23:44 +1100 >> From: smyth at wehi.EDU.AU >> To: genome1976 at hotmail.com >> CC: naomi at stat.psu.edu; bioconductor at r-project.org >> Subject: Re: [BioC] Coefficients not estimable >> >> Or perhaps as a nested fixed model >> >> model.matrix(~gender+gender:donor+day) >> >> Whether to fit donor as a fixed term or as a block really depends on what >> questions you want to answer and, in particular, what you want to >> interpret from a donor "main effect". >> >> Best wishes >> Gordon >> >> --------------------------------------------- >> Professor Gordon K Smyth, >> Bioinformatics Division, >> Walter and Eliza Hall Institute of Medical Research, >> 1G Royal Parade, Parkville, Vic 3052, Australia. >> smyth at wehi.edu.au >> http://www.wehi.edu.au >> http://www.statsci.org/smyth >> >> On Tue, 7 Feb 2012, Naomi Altman wrote: >> >>> The problem is that donor is nested within gender. The only way I am aware >>> of to handle this in limma is to use donor as a block. >>> >>> Regards, >>> Naomi >>> >>> >>> At 11:24 AM 2/7/2012, somnath bandyopadhyay wrote: >>> >>> >>>> Dear Gordon and limma users, >>>> >>>> I am trying to fit a linear model to assess variability in healthy donors >>>> over days for both genders using Affymetrix microarray data. I am >>>> interested in assessing variability arising from geneder, day of blood >>>> donated and individual donors. When I tried to do a linear fit in using >>>> limma, I am not able to assess the coefficients for gender.System is >>>> computationally singular for the gender main effect. Is there any way to >>>> assess the variability arising from all three main effects? Any help would >>>> be greatly appreciated. >>>> >>>> Following is my target file: >>>> >>>> >>>> >>>> >>>> >>>> >>>> >>>> >>>> Array_NUID >>>> gender >>>> day >>>> donor >>>> >>>> NUID.0000.0133.2260 >>>> male >>>> d1 >>>> M1 >>>> >>>> NUID.0000.0133.2261 >>>> male >>>> d2 >>>> M1 >>>> >>>> NUID.0000.0133.2262 >>>> male >>>> d8 >>>> M1 >>>> >>>> NUID.0000.0133.2230 >>>> female >>>> d1 >>>> F2 >>>> >>>> NUID.0000.0133.6809 >>>> female >>>> d2 >>>> F2 >>>> >>>> NUID.0000.0133.2251 >>>> female >>>> d8 >>>> F2 >>>> >>>> NUID.0000.0133.2228 >>>> female >>>> d15 >>>> F2 >>>> >>>> NUID.0000.0133.2229 >>>> female >>>> d28 >>>> F2 >>>> >>>> NUID.0000.0133.2231 >>>> female >>>> d1 >>>> F3 >>>> >>>> NUID.0000.0133.2232 >>>> female >>>> d2 >>>> F3 >>>> >>>> NUID.0000.0133.2233 >>>> female >>>> d8 >>>> F3 >>>> >>>> NUID.0000.0133.2234 >>>> female >>>> d15 >>>> F3 >>>> >>>> NUID.0000.0133.2265 >>>> male >>>> d1 >>>> M4 >>>> >>>> NUID.0000.0133.2267 >>>> male >>>> d8 >>>> M4 >>>> >>>> NUID.0000.0133.2268 >>>> male >>>> d15 >>>> M4 >>>> >>>> NUID.0000.0133.2269 >>>> male >>>> d28 >>>> M4 >>>> >>>> NUID.0000.0133.2236 >>>> female >>>> d1 >>>> F5 >>>> >>>> NUID.0000.0133.2237 >>>> female >>>> d2 >>>> F5 >>>> >>>> NUID.0000.0133.2238 >>>> female >>>> d8 >>>> F5 >>>> >>>> NUID.0000.0133.6810 >>>> female >>>> d15 >>>> F5 >>>> >>>> NUID.0000.0133.2240 >>>> female >>>> d28 >>>> F5 >>>> >>>> NUID.0000.0133.2270 >>>> male >>>> d1 >>>> M6 >>>> >>>> NUID.0000.0133.2271 >>>> male >>>> d2 >>>> M6 >>>> >>>> NUID.0000.0133.2272 >>>> male >>>> d8 >>>> M6 >>>> >>>> NUID.0000.0133.2273 >>>> male >>>> d15 >>>> M6 >>>> >>>> NUID.0000.0133.2274 >>>> male >>>> d28 >>>> M6 >>>> >>>> NUID.0000.0133.2241 >>>> female >>>> d1 >>>> F7 >>>> >>>> NUID.0000.0133.2242 >>>> female >>>> d2 >>>> F7 >>>> >>>> NUID.0000.0133.2243 >>>> female >>>> d8 >>>> F7 >>>> >>>> NUID.0000.0133.2244 >>>> female >>>> d15 >>>> F7 >>>> >>>> NUID.0000.0133.2245 >>>> female >>>> d28 >>>> F7 >>>> >>>> NUID.0000.0133.2246 >>>> female >>>> d1 >>>> F8 >>>> >>>> NUID.0000.0133.2247 >>>> female >>>> d2 >>>> F8 >>>> >>>> NUID.0000.0133.2248 >>>> female >>>> d8 >>>> F8 >>>> >>>> NUID.0000.0133.2249 >>>> female >>>> d15 >>>> F8 >>>> >>>> NUID.0000.0133.2250 >>>> female >>>> d28 >>>> F8 >>>> >>>> >>>> Following is the code I am using: >>>> library("limma") >>>> rna <- read.table('qced_non_redundant_rma_data.txt', header=T, row.names=1, >>>> sep= "\t") >>>> annotation <- read.table('Test_Samples_Annotation.txt', header=T, >>>> row.names=1, sep= "\t") >>>> ## rearrange the annotations to conform to the RMA data >>>> annotation <- annotation[colnames(rna),] >>>> ## check that the annotations and data match >>>> if (sum(sort(rownames(annotation)) == sort(colnames(rna))) != >>>> length(rownames(annotation))) { >>>> stop ("ERROR: annotation rownames and RMA colnames do not match") >>>> } >>>> donor <- factor(annotation$donor) >>>> day <- factor(annotation$day) >>>> gender <- factor(annotation$gender) >>>> design <- model.matrix(~ donor + gender + day, data=annotation) >>>> fit <- lmFit(rna,design) >>>> fit2 <- eBayes(fit) >>>> write.fit(fit,file="effects.txt",digits=30, adjust="BH", >>>> method="separate",sep="\t") >>>> >>>> >>>> >>>> Following is my sessionInfo: >>>>> sessionInfo() >>>> R version 2.10.1 (2009-12-14) >>>> x86_64-unknown-linux-gnu >>>> locale: >>>> [1] LC_CTYPE=en_US.iso885915 LC_NUMERIC=C >>>> [3] LC_TIME=en_US.iso885915 LC_COLLATE=en_US.iso885915 >>>> [5] LC_MONETARY=C LC_MESSAGES=en_US.iso885915 >>>> [7] LC_PAPER=en_US.iso885915 LC_NAME=C >>>> [9] LC_ADDRESS=C LC_TELEPHONE=C >>>> [11] LC_MEASUREMENT=en_US.iso885915 LC_IDENTIFICATION=C >>>> attached base packages: >>>> [1] stats graphics grDevices utils datasets methods base >>>> other attached packages: >>>> [1] limma_3.2.3 >>>> loaded via a namespace (and not attached): >>>> [1] tools_2.10.1 >>>> [[alternative HTML version deleted]] >>>> >>>> _______________________________________________ >>>> Bioconductor mailing list >>>> Bioconductor at r-project.org >>>> https://stat.ethz.ch/mailman/listinfo/bioconductor >>>> Search the archives: >>>> http://news.gmane.org/gmane.science.biology.informatics.conductor >>> >>> >>> >>> >> >> ______________________________________________________________________ >> The information in this email is confidential and intended solely for the addressee. >> You must not disclose, forward, print or use it without the permission of the sender. >> ______________________________________________________________________ > ______________________________________________________________________ The information in this email is confidential and intend...{{dropped:4}}