error and warning message with blocked samples with Limma
2
0
Entering edit mode
@richard-friedman-513
Last seen 9.6 years ago
Dear Bioconductor List, I want to do a linear model analysis of a blocked dataset in Limma. I get an error message if I closely follow the manual, or a warning message if I do it a little differentially. I understand from previous posts that the warning message may not be serious, but I would like to make sure that I am doing this correctly. My dataset is: ileName Name mouse condition QIU01.CEL AA.het.con.01.01 1 AA.het.con QIU02.CEL AB.het.sta.01.02 1 AB.het.sta QIU03.CEL AA.het.con.02.03 2 AA.het.con QIU04.CEL AB.het.sta.02.04 2 AB.het.sta QIU05.CEL BA.mut.con.03.05 3 BA.mut.con QIU06.CEL BB.mut.sta.03.06 3 BB.mut.sta QIU07.CEL BA.mut.con.04.07 4 BA.mut.con QIU08.CEL BB.mutsta.04.08 4 BB.mut.sta QIU09.CEL AA.het.con.05.09 5 AA.het.con QIU10.CEL AB.het.sta.05.10 5 AB.het.sta QIU11.CEL AA.het.con.06.11 6 AA.het.con QIU12.CEL AB.het.sta.06.12 6 AB.het.sta QIU13.CEL BA.mut.con.07.13 7 BA.mut.con QIU14.CEL BB.mut.sta.07.14 7 BB.mut.sta QIU15.CEL BA.mut.con.08.15 8 BA.mut.con QIU16.CEL BB.mut.sta.08.16 8 BB.mut.sta QIU17.CEL BA.mut.con.09.17 9 BA.mut.con QIU18.CEL BB.mut.sta.09.18 9 BB.mut.sta QIU19.CEL BA.mut.con.10.19 10 BA.mut.con QIU20.CEL BB.mut.sta.10.20 10 BB.mut.sta QIU21.CEL AA.het.con.11.21 11 AA.het.con QIU22.CEL AB.het.sta.11.22 11 AB.het.sta QIU23.CEL AA.het.con.12.23 12 AA.het.con QIU24.CEL AB.het.sta.12.24 12 AB.het.sta Arrays from the same mouse are paired but they are treated differently. I need the following 3 comparisons: AB.het.sta-AA.het.con BB.mut.sta-BA.mut.con ((BB.mut.sta-BA.mut.con)-(AB.het.sta-AA.het.con)) Here is my environment: > sessionInfo() R version 2.15.2 (2012-10-26) 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] mouse4302.db_2.8.1 mouse4302probe_2.11.0 mouse4302cdf_2.11.0 [4] org.Mm.eg.db_2.8.0 limma_3.14.3 gcrma_2.30.0 [7] annaffy_1.30.0 KEGG.db_2.8.0 GO.db_2.8.0 [10] RSQLite_0.11.2 DBI_0.2-5 AnnotationDbi_1.20.3 [13] affy_1.36.0 Biobase_2.18.0 BiocGenerics_0.4.0 loaded via a namespace (and not attached): [1] affyio_1.26.0 BiocInstaller_1.8.3 Biostrings_2.26.2 [4] IRanges_1.16.4 parallel_2.15.2 preprocessCore_1.20.0 [7] splines_2.15.2 stats4_2.15.2 tools_2.15.2 [10] zlibbioc_1.4.0 (I have not had the chance to switch to R3.0 yet and am waiting a while until any bugs get worked out). Here is my script: library(affy) library(annaffy) library(gcrma) library(limma) library(mouse4302.db) raw<-ReadAffy() gcrmanm<-gcrma(raw) write.exprs(gcrmanm,"gloria2gcexp.txt") targets<-readTargets("targets.txt") targets condition<-factor(targets$condition,levels=c("AA.het.con","AB.het.sta" ,"BA.mut.con","BB.mut.sta")) condition mouse<-factor(targets$mouse) mouse design<-model.matrix(~mouse+condition) design colnames(design)<-c("mouse2","mouse3","mouse4", "mouse5","mouse6","mouse7","mouse8", "mouse9","mouse10","mouse11","mouse12", "AB.het.sta","BA.mut.con","BB.mut.sta") fit<-lmFit(gcrmanm,design) contrast.matrix<-makeContrasts(AB.het.sta-AA.het.con,BB.mut.sta- BA.mut.con, ((BB.mut.sta-BA.mut.con)-(AB.het.sta-AA.het.con)),levels=design) contrast.matrix fit<-contrasts.fit(fit,contrast.matrix) fit<-eBayes(fit) comp1<-topTable(fit,coef=1,number=45101,adjust.method="BH",sort.by="B" ) names(comp1) names(comp1)[1]="Probe" probeids=comp1$Probe aaf.handler() anncols<-aaf.handler()[c(1:5,11:12)] anncols comp1table<-aafTableFrame(comp1,colnames=names(comp1),probeids=probeid s) anntable<-aafTableAnn(probeids, "mouse4302.db", colnames=anncols) comp1ann<-merge(comp1table,anntable) saveText(comp1ann,"AB.het.staVsAA.het.con.ann.txt") comp2<-topTable(fit,coef=2,number=45101,adjust.method="BH",sort.by="B" ) names(comp2) names(comp2)[1]="Probe" probeids=comp2$Probe aaf.handler() anncols<-aaf.handler()[c(1:5,11:12)] anncols comp2table<-aafTableFrame(comp2,colnames=names(comp2),probeids=probeid s) anntable<-aafTableAnn(probeids, "mouse4302.db", colnames=anncols) comp2ann<-merge(comp2table,anntable) saveText(comp1ann,"BB.mut.staVsBA.mut.con.ann.txt") comp3<-topTable(fit,coef=3,number=45101,adjust.method="BH",sort.by="B" ) names(comp3) names(comp3)[1]="Probe" probeids=comp3$Probe aaf.handler() anncols<-aaf.handler()[c(1:5,11:12)] anncols comp3table<-aafTableFrame(comp3,names(comp3),probeids=probeids) anntable<-aafTableAnn(probeids, "mouse4302.db", colnames=anncols) comp3ann<-merge(comp3table,anntable) saveText(comp3ann,"starvationresponse.mutVshet.ann.txt") I get an error message at the point at which I am naming the columns: > colnames(design)<-c("mouse2","mouse3","mouse4", + "mouse5","mouse6","mouse7","mouse8", + "mouse9","mouse10","mouse11","mouse12", + "AB.het.sta","BA.mut.con","BB.mut.sta") Error in `colnames<-`(`*tmp*`, value = c("mouse2", "mouse3", "mouse4", : length of 'dimnames' [2] not equal to array extent Calls: colnames<- -> colnames<- Execution halted I then tried the following revised script: library(affy) library(annaffy) library(gcrma) library(limma) library(mouse4302.db) raw<-ReadAffy() gcrmanm<-gcrma(raw) write.exprs(gcrmanm,"gloria2gcexp.txt") targets<-readTargets("targets.txt") targets condition<-factor(targets$condition,levels=c("AA.het.con","AB.het.sta" ,"BA.mut.con","BB.mut.sta")) condition mouse<-factor(targets$mouse) mouse design<-model.matrix(~0+condition+mouse) design colnames(design)<-c("AA.het.con","AB.het.sta","BA.mut.con","BB.mut.sta ", "mouse2","mouse3","mouse4", "mouse5","mouse6","mouse7","mouse8", "mouse9","mouse10","mouse11","mouse12") fit<-lmFit(gcrmanm,design) (everything else the same). I get the following warning message: > colnames(design)<-c("AA.het.con","AB.het.sta","BA.mut.con","BB.mut.s ta", + "mouse2","mouse3","mouse4", + "mouse5","mouse6","mouse7","mouse8", + "mouse9","mouse10","mouse11","mouse12") > fit<-lmFit(gcrmanm,design) Coefficients not estimable: mouse10 Warning message: Partial NA coefficients for 45101 probe(s) THE SCRIPT RAN TO COMPLETION MY QUESTIONS: 1. Should I worry about the warning message? 2. Is there any way in which I can do this analysis and not get the warning message. Thanks and best wishes, Rich Richard A. Friedman, PhD Associate Research Scientist, Biomedical Informatics Shared Resource Herbert Irving Comprehensive Cancer Center (HICCC) Lecturer, Department of Biomedical Informatics (DBMI) Educational Coordinator, Center for Computational Biology and Bioinformatics (C2B2)/ National Center for Multiscale Analysis of Genomic Networks (MAGNet)/ Columbia Initiative in Systems Biology Room 824 Irving Cancer Research Center Columbia University 1130 St. Nicholas Ave New York, NY 10032 (212)851-4765 (voice) friedman at cancercenter.columbia.edu http://cancercenter.columbia.edu/~friedman/ Fritz Lang. Didn't he do "Star Trek". -Rose Friedman, age 16
GO Cancer mouse4302 limma GO Cancer mouse4302 limma • 1.8k views
ADD COMMENT
0
Entering edit mode
@gordon-smyth
Last seen 10 hours ago
WEHI, Melbourne, Australia
Dear Richard, The error is coming when you set colnames for a matrix. When you try to set colnames for matrix, you must give the correct number of names as there are columns in the matrix, and the error is saying that you are not doing so. Looking at your code, it would seem that you have ommited the Intercept from your list of names. Even if you fixed this, the code would still fail later on because you use "AA.het.con" in makeContrasts() but there is no column by this name in your design matrix. Your second analysis runs because you give 15 column names instead of 14. You can get the three contrasts you want by design <- model.matrix(~mouse+condition) contrast.matrix <- makeContrasts( "AB.het.sta-AA.het.con" = conditionAB.het.sta, "BB.mut.sta-BA.mut.con" = conditionBB.mut.sta-conditionBA.mut.con, "Diff" = (conditionBB.mut.sta- conditionBA.mut.con)-conditionAB.het.sta, levels=design) Best wishes Gordon > Date: Fri, 5 Apr 2013 10:42:22 -0400 > From: Richard Friedman <friedman at="" cancercenter.columbia.edu=""> > To: "bioconductor at r-project.org list" <bioconductor at="" r-project.org=""> > Subject: [BioC] error and warning message with blocked samples with > Limma > > Dear Bioconductor List, > > I want to do a linear model analysis of a blocked dataset > in Limma. I get an error message if I closely follow the manual, > or a warning message if I do it a little differentially. I understand > from previous posts that the warning message may not be serious, > but I would like to make sure that I am doing this correctly. > > My dataset is: > > ileName Name mouse condition > QIU01.CEL AA.het.con.01.01 1 AA.het.con > QIU02.CEL AB.het.sta.01.02 1 AB.het.sta > QIU03.CEL AA.het.con.02.03 2 AA.het.con > QIU04.CEL AB.het.sta.02.04 2 AB.het.sta > QIU05.CEL BA.mut.con.03.05 3 BA.mut.con > QIU06.CEL BB.mut.sta.03.06 3 BB.mut.sta > QIU07.CEL BA.mut.con.04.07 4 BA.mut.con > QIU08.CEL BB.mutsta.04.08 4 BB.mut.sta > QIU09.CEL AA.het.con.05.09 5 AA.het.con > QIU10.CEL AB.het.sta.05.10 5 AB.het.sta > QIU11.CEL AA.het.con.06.11 6 AA.het.con > QIU12.CEL AB.het.sta.06.12 6 AB.het.sta > QIU13.CEL BA.mut.con.07.13 7 BA.mut.con > QIU14.CEL BB.mut.sta.07.14 7 BB.mut.sta > QIU15.CEL BA.mut.con.08.15 8 BA.mut.con > QIU16.CEL BB.mut.sta.08.16 8 BB.mut.sta > QIU17.CEL BA.mut.con.09.17 9 BA.mut.con > QIU18.CEL BB.mut.sta.09.18 9 BB.mut.sta > QIU19.CEL BA.mut.con.10.19 10 BA.mut.con > QIU20.CEL BB.mut.sta.10.20 10 BB.mut.sta > QIU21.CEL AA.het.con.11.21 11 AA.het.con > QIU22.CEL AB.het.sta.11.22 11 AB.het.sta > QIU23.CEL AA.het.con.12.23 12 AA.het.con > QIU24.CEL AB.het.sta.12.24 12 AB.het.sta > > Arrays from the same mouse are paired but they are > treated differently. I need the following 3 comparisons: > > AB.het.sta-AA.het.con > BB.mut.sta-BA.mut.con > ((BB.mut.sta-BA.mut.con)-(AB.het.sta-AA.het.con)) > > Here is my environment: > >> sessionInfo() > R version 2.15.2 (2012-10-26) > 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] mouse4302.db_2.8.1 mouse4302probe_2.11.0 mouse4302cdf_2.11.0 > [4] org.Mm.eg.db_2.8.0 limma_3.14.3 gcrma_2.30.0 > [7] annaffy_1.30.0 KEGG.db_2.8.0 GO.db_2.8.0 > [10] RSQLite_0.11.2 DBI_0.2-5 AnnotationDbi_1.20.3 > [13] affy_1.36.0 Biobase_2.18.0 BiocGenerics_0.4.0 > > loaded via a namespace (and not attached): > [1] affyio_1.26.0 BiocInstaller_1.8.3 Biostrings_2.26.2 > [4] IRanges_1.16.4 parallel_2.15.2 preprocessCore_1.20.0 > [7] splines_2.15.2 stats4_2.15.2 tools_2.15.2 > [10] zlibbioc_1.4.0 > > > (I have not had the chance to switch to R3.0 yet and am waiting a while until > any bugs get worked out). > > Here is my script: > > library(affy) > library(annaffy) > library(gcrma) > library(limma) > library(mouse4302.db) > raw<-ReadAffy() > gcrmanm<-gcrma(raw) > write.exprs(gcrmanm,"gloria2gcexp.txt") > > > targets<-readTargets("targets.txt") > targets > condition<-factor(targets$condition,levels=c("AA.het.con","AB.het.st a","BA.mut.con","BB.mut.sta")) > condition > mouse<-factor(targets$mouse) > mouse > design<-model.matrix(~mouse+condition) > design > colnames(design)<-c("mouse2","mouse3","mouse4", > "mouse5","mouse6","mouse7","mouse8", > "mouse9","mouse10","mouse11","mouse12", > "AB.het.sta","BA.mut.con","BB.mut.sta") > > fit<-lmFit(gcrmanm,design) > > contrast.matrix<-makeContrasts(AB.het.sta-AA.het.con,BB.mut.sta- BA.mut.con, > ((BB.mut.sta-BA.mut.con)-(AB.het.sta-AA.het.con)),levels=design) > contrast.matrix > fit<-contrasts.fit(fit,contrast.matrix) > fit<-eBayes(fit) > > comp1<-topTable(fit,coef=1,number=45101,adjust.method="BH",sort.by=" B") > > > names(comp1) > names(comp1)[1]="Probe" > probeids=comp1$Probe > > aaf.handler() > anncols<-aaf.handler()[c(1:5,11:12)] > anncols > > comp1table<-aafTableFrame(comp1,colnames=names(comp1),probeids=probe ids) > anntable<-aafTableAnn(probeids, "mouse4302.db", colnames=anncols) > comp1ann<-merge(comp1table,anntable) > saveText(comp1ann,"AB.het.staVsAA.het.con.ann.txt") > > > comp2<-topTable(fit,coef=2,number=45101,adjust.method="BH",sort.by=" B") > > names(comp2) > names(comp2)[1]="Probe" > probeids=comp2$Probe > > aaf.handler() > anncols<-aaf.handler()[c(1:5,11:12)] > anncols > > comp2table<-aafTableFrame(comp2,colnames=names(comp2),probeids=probe ids) > anntable<-aafTableAnn(probeids, "mouse4302.db", colnames=anncols) > comp2ann<-merge(comp2table,anntable) > saveText(comp1ann,"BB.mut.staVsBA.mut.con.ann.txt") > > comp3<-topTable(fit,coef=3,number=45101,adjust.method="BH",sort.by=" B") > > names(comp3) > names(comp3)[1]="Probe" > probeids=comp3$Probe > > aaf.handler() > anncols<-aaf.handler()[c(1:5,11:12)] > anncols > > comp3table<-aafTableFrame(comp3,names(comp3),probeids=probeids) > anntable<-aafTableAnn(probeids, "mouse4302.db", colnames=anncols) > comp3ann<-merge(comp3table,anntable) > saveText(comp3ann,"starvationresponse.mutVshet.ann.txt") > > > I get an error message at the point at which I am naming the columns: > >> colnames(design)<-c("mouse2","mouse3","mouse4", > + "mouse5","mouse6","mouse7","mouse8", > + "mouse9","mouse10","mouse11","mouse12", > + "AB.het.sta","BA.mut.con","BB.mut.sta") > Error in `colnames<-`(`*tmp*`, value = c("mouse2", "mouse3", "mouse4", : > length of 'dimnames' [2] not equal to array extent > Calls: colnames<- -> colnames<- > Execution halted > > I then tried the following revised script: > > library(affy) > library(annaffy) > library(gcrma) > library(limma) > library(mouse4302.db) > raw<-ReadAffy() > gcrmanm<-gcrma(raw) > write.exprs(gcrmanm,"gloria2gcexp.txt") > > > targets<-readTargets("targets.txt") > targets > condition<-factor(targets$condition,levels=c("AA.het.con","AB.het.st a","BA.mut.con","BB.mut.sta")) > condition > mouse<-factor(targets$mouse) > mouse > design<-model.matrix(~0+condition+mouse) > design > colnames(design)<-c("AA.het.con","AB.het.sta","BA.mut.con","BB.mut.s ta", > "mouse2","mouse3","mouse4", > "mouse5","mouse6","mouse7","mouse8", > "mouse9","mouse10","mouse11","mouse12") > fit<-lmFit(gcrmanm,design) > > (everything else the same). > > I get the following warning message: > >> colnames(design)<-c("AA.het.con","AB.het.sta","BA.mut.con","BB.mut. sta", > + "mouse2","mouse3","mouse4", > + "mouse5","mouse6","mouse7","mouse8", > + "mouse9","mouse10","mouse11","mouse12") >> fit<-lmFit(gcrmanm,design) > Coefficients not estimable: mouse10 > Warning message: > Partial NA coefficients for 45101 probe(s) > > THE SCRIPT RAN TO COMPLETION > > MY QUESTIONS: > 1. Should I worry about the warning message? > 2. Is there any way in which I can do this analysis and not get the warning message. > > Thanks and best wishes, > Rich > Richard A. Friedman, PhD > Associate Research Scientist, > Biomedical Informatics Shared Resource > Herbert Irving Comprehensive Cancer Center (HICCC) > Lecturer, > Department of Biomedical Informatics (DBMI) > Educational Coordinator, > Center for Computational Biology and Bioinformatics (C2B2)/ > National Center for Multiscale Analysis of Genomic Networks (MAGNet)/ > Columbia Initiative in Systems Biology > Room 824 > Irving Cancer Research Center > Columbia University > 1130 St. Nicholas Ave > New York, NY 10032 > (212)851-4765 (voice) > friedman at cancercenter.columbia.edu > http://cancercenter.columbia.edu/~friedman/ > > Fritz Lang. Didn't he do "Star Trek". > -Rose Friedman, age 16 > ______________________________________________________________________ The information in this email is confidential and intend...{{dropped:4}}
ADD COMMENT
0
Entering edit mode
On Apr 6, 2013, at 7:41 PM, Gordon K Smyth wrote: > Dear Richard, > > The error is coming when you set colnames for a matrix. When you try to set colnames for matrix, you must give the correct number of names as there are columns in the matrix, and the error is saying that you are not doing so. Looking at your code, it would seem that you have ommited the Intercept from your list of names. > > Even if you fixed this, the code would still fail later on because you use "AA.het.con" in makeContrasts() but there is no column by this name in your design matrix. > > Your second analysis runs because you give 15 column names instead of 14. > > You can get the three contrasts you want by > > design <- model.matrix(~mouse+condition) > contrast.matrix <- makeContrasts( > "AB.het.sta-AA.het.con" = conditionAB.het.sta, > "BB.mut.sta-BA.mut.con" = conditionBB.mut.sta-conditionBA.mut.con, > "Diff" = (conditionBB.mut.sta- conditionBA.mut.con)-conditionAB.het.sta, > levels=design) > > > Best wishes > Gordon Dear Gordon, Thank you for your reply. I am still having difficulty. I have tried both putting the code you recommended in and I also tried with the intercept, but I still get error messages. Is the script I wrote with design<-model.matrix(~0+condition+mouse) that ran all the way through reliable? Here is an excerpt from the output from just putting the code above in: ################################################################## > design<-model.matrix(~mouse+condition) > design (Intercept) mouse2 mouse3 mouse4 mouse5 mouse6 mouse7 mouse8 mouse9 mouse10 1 1 0 0 0 0 0 0 0 0 0 2 1 0 0 0 0 0 0 0 0 0 3 1 1 0 0 0 0 0 0 0 0 4 1 1 0 0 0 0 0 0 0 0 5 1 0 1 0 0 0 0 0 0 0 6 1 0 1 0 0 0 0 0 0 0 7 1 0 0 1 0 0 0 0 0 0 8 1 0 0 1 0 0 0 0 0 0 9 1 0 0 0 1 0 0 0 0 0 10 1 0 0 0 1 0 0 0 0 0 11 1 0 0 0 0 1 0 0 0 0 12 1 0 0 0 0 1 0 0 0 0 13 1 0 0 0 0 0 1 0 0 0 14 1 0 0 0 0 0 1 0 0 0 15 1 0 0 0 0 0 0 1 0 0 16 1 0 0 0 0 0 0 1 0 0 17 1 0 0 0 0 0 0 0 1 0 18 1 0 0 0 0 0 0 0 1 0 19 1 0 0 0 0 0 0 0 0 1 20 1 0 0 0 0 0 0 0 0 1 21 1 0 0 0 0 0 0 0 0 0 22 1 0 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 mouse11 mouse12 conditionAB.het.sta conditionBA.mut.con conditionBB.mut.sta 1 0 0 0 0 0 2 0 0 1 0 0 3 0 0 0 0 0 4 0 0 1 0 0 5 0 0 0 1 0 6 0 0 0 0 1 7 0 0 0 1 0 8 0 0 0 0 1 9 0 0 0 0 0 10 0 0 1 0 0 11 0 0 0 0 0 12 0 0 1 0 0 13 0 0 0 1 0 14 0 0 0 0 1 15 0 0 0 1 0 16 0 0 0 0 1 17 0 0 0 1 0 18 0 0 0 0 1 19 0 0 0 1 0 20 0 0 0 0 1 21 1 0 0 0 0 22 1 0 1 0 0 23 0 1 0 0 0 24 0 1 1 0 0 attr(,"assign") [1] 0 1 1 1 1 1 1 1 1 1 1 1 2 2 2 attr(,"contrasts") attr(,"contrasts")$mouse [1] "contr.treatment" attr(,"contrasts")$condition [1] "contr.treatment" > > fit<-lmFit(gcrmanm,design) Coefficients not estimable: conditionBB.mut.sta Warning message: Partial NA coefficients for 45101 probe(s) > > contrast.matrix <- makeContrasts( + "AB.het.sta-AA.het.con" = conditionAB.het.sta, + "BB.mut.sta-BA.mut.con" = conditionBB.mut.sta-conditionBA.mut.con, + "Diff" = (conditionBB.mut.sta- conditionBA.mut.con)-conditionAB.het.sta, + levels=design) Warning message: In makeContrasts(`AB.het.sta-AA.het.con` = conditionAB.het.sta, : Renaming (Intercept) to Intercept > > > contrast.matrix Contrasts Levels AB.het.sta-AA.het.con BB.mut.sta-BA.mut.con Diff Intercept 0 0 0 mouse2 0 0 0 mouse3 0 0 0 mouse4 0 0 0 mouse5 0 0 0 mouse6 0 0 0 mouse7 0 0 0 mouse8 0 0 0 mouse9 0 0 0 mouse10 0 0 0 mouse11 0 0 0 mouse12 0 0 0 conditionAB.het.sta 1 0 -1 conditionBA.mut.con 0 -1 -1 conditionBB.mut.sta 0 1 1 > fit<-contrasts.fit(fit,contrast.matrix) Error in contrasts.fit(fit, contrast.matrix) : trying to take contrast of non-estimable coefficient In addition: Warning message: In contrasts.fit(fit, contrast.matrix) : row names of contrasts don't match col names of coefficients Execution halted ########################################################### Here is an excerpt from the output from my explicitly naming the intercept and naming the columns: Levels: 1 2 3 4 5 6 7 8 9 10 11 12 > design<-model.matrix(~mouse+condition) > design (Intercept) mouse2 mouse3 mouse4 mouse5 mouse6 mouse7 mouse8 mouse9 mouse10 1 1 0 0 0 0 0 0 0 0 0 2 1 0 0 0 0 0 0 0 0 0 3 1 1 0 0 0 0 0 0 0 0 4 1 1 0 0 0 0 0 0 0 0 5 1 0 1 0 0 0 0 0 0 0 6 1 0 1 0 0 0 0 0 0 0 7 1 0 0 1 0 0 0 0 0 0 8 1 0 0 1 0 0 0 0 0 0 9 1 0 0 0 1 0 0 0 0 0 10 1 0 0 0 1 0 0 0 0 0 11 1 0 0 0 0 1 0 0 0 0 12 1 0 0 0 0 1 0 0 0 0 13 1 0 0 0 0 0 1 0 0 0 14 1 0 0 0 0 0 1 0 0 0 15 1 0 0 0 0 0 0 1 0 0 16 1 0 0 0 0 0 0 1 0 0 17 1 0 0 0 0 0 0 0 1 0 18 1 0 0 0 0 0 0 0 1 0 19 1 0 0 0 0 0 0 0 0 1 20 1 0 0 0 0 0 0 0 0 1 21 1 0 0 0 0 0 0 0 0 0 22 1 0 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 mouse11 mouse12 conditionAB.het.sta conditionBA.mut.con conditionBB.mut.sta 1 0 0 0 0 0 2 0 0 1 0 0 3 0 0 0 0 0 4 0 0 1 0 0 5 0 0 0 1 0 6 0 0 0 0 1 7 0 0 0 1 0 8 0 0 0 0 1 9 0 0 0 0 0 10 0 0 1 0 0 11 0 0 0 0 0 12 0 0 1 0 0 13 0 0 0 1 0 14 0 0 0 0 1 15 0 0 0 1 0 16 0 0 0 0 1 17 0 0 0 1 0 18 0 0 0 0 1 19 0 0 0 1 0 20 0 0 0 0 1 21 1 0 0 0 0 22 1 0 1 0 0 23 0 1 0 0 0 24 0 1 1 0 0 attr(,"assign") [1] 0 1 1 1 1 1 1 1 1 1 1 1 2 2 2 attr(,"contrasts") attr(,"contrasts")$mouse [1] "contr.treatment" attr(,"contrasts")$condition [1] "contr.treatment" > > colnames(design)<-c("Intercept","mouse2","mouse3","mouse4", + "mouse5","mouse6","mouse7","mouse8", + "mouse9","mouse10","mouse11","mouse12", + "AB.het.sta","BA.mut.con","BB.mut.sta") > > > fit<-lmFit(gcrmanm,design) Coefficients not estimable: BB.mut.sta Warning message: Partial NA coefficients for 45101 probe(s) > > #contrast.matrix <- makeContrasts( > #"AB.het.sta-AA.het.con" = conditionAB.het.sta, > #"BB.mut.sta-BA.mut.con" = conditionBB.mut.sta-conditionBA.mut.con, > #"Diff" = (conditionBB.mut.sta- conditionBA.mut.con)-conditionAB.het.sta, > #levels=design) > > contrast.matrix <- makeContrasts(AB.het.sta,BB.mut.sta-BA.mut.con, + (BB.mut.sta-BA.mut.con)-AB.het.sta,levels=design) > > > contrast.matrix Contrasts Levels AB.het.sta BB.mut.sta - BA.mut.con Intercept 0 0 mouse2 0 0 mouse3 0 0 mouse4 0 0 mouse5 0 0 mouse6 0 0 mouse7 0 0 mouse8 0 0 mouse9 0 0 mouse10 0 0 mouse11 0 0 mouse12 0 0 AB.het.sta 1 0 BA.mut.con 0 -1 BB.mut.sta 0 1 Contrasts Levels (BB.mut.sta - BA.mut.con) - AB.het.sta Intercept 0 mouse2 0 mouse3 0 mouse4 0 mouse5 0 mouse6 0 mouse7 0 mouse8 0 mouse9 0 mouse10 0 mouse11 0 mouse12 0 AB.het.sta -1 BA.mut.con -1 BB.mut.sta 1 > fit<-contrasts.fit(fit,contrast.matrix) Error in contrasts.fit(fit, contrast.matrix) : trying to take contrast of non-estimable coefficient Execution halted ################################################# I am clearly missing something. Again, I would appreciate any suggestions. Thanks and best wishes, Rich > >> Date: Fri, 5 Apr 2013 10:42:22 -0400 >> From: Richard Friedman <friedman at="" cancercenter.columbia.edu=""> >> To: "bioconductor at r-project.org list" <bioconductor at="" r-project.org=""> >> Subject: [BioC] error and warning message with blocked samples with >> Limma >> >> Dear Bioconductor List, >> >> I want to do a linear model analysis of a blocked dataset >> in Limma. I get an error message if I closely follow the manual, >> or a warning message if I do it a little differentially. I understand >> from previous posts that the warning message may not be serious, >> but I would like to make sure that I am doing this correctly. >> >> My dataset is: >> >> ileName Name mouse condition >> QIU01.CEL AA.het.con.01.01 1 AA.het.con >> QIU02.CEL AB.het.sta.01.02 1 AB.het.sta >> QIU03.CEL AA.het.con.02.03 2 AA.het.con >> QIU04.CEL AB.het.sta.02.04 2 AB.het.sta >> QIU05.CEL BA.mut.con.03.05 3 BA.mut.con >> QIU06.CEL BB.mut.sta.03.06 3 BB.mut.sta >> QIU07.CEL BA.mut.con.04.07 4 BA.mut.con >> QIU08.CEL BB.mutsta.04.08 4 BB.mut.sta >> QIU09.CEL AA.het.con.05.09 5 AA.het.con >> QIU10.CEL AB.het.sta.05.10 5 AB.het.sta >> QIU11.CEL AA.het.con.06.11 6 AA.het.con >> QIU12.CEL AB.het.sta.06.12 6 AB.het.sta >> QIU13.CEL BA.mut.con.07.13 7 BA.mut.con >> QIU14.CEL BB.mut.sta.07.14 7 BB.mut.sta >> QIU15.CEL BA.mut.con.08.15 8 BA.mut.con >> QIU16.CEL BB.mut.sta.08.16 8 BB.mut.sta >> QIU17.CEL BA.mut.con.09.17 9 BA.mut.con >> QIU18.CEL BB.mut.sta.09.18 9 BB.mut.sta >> QIU19.CEL BA.mut.con.10.19 10 BA.mut.con >> QIU20.CEL BB.mut.sta.10.20 10 BB.mut.sta >> QIU21.CEL AA.het.con.11.21 11 AA.het.con >> QIU22.CEL AB.het.sta.11.22 11 AB.het.sta >> QIU23.CEL AA.het.con.12.23 12 AA.het.con >> QIU24.CEL AB.het.sta.12.24 12 AB.het.sta >> >> Arrays from the same mouse are paired but they are >> treated differently. I need the following 3 comparisons: >> >> AB.het.sta-AA.het.con >> BB.mut.sta-BA.mut.con >> ((BB.mut.sta-BA.mut.con)-(AB.het.sta-AA.het.con)) >> >> Here is my environment: >> >>> sessionInfo() >> R version 2.15.2 (2012-10-26) >> 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] mouse4302.db_2.8.1 mouse4302probe_2.11.0 mouse4302cdf_2.11.0 >> [4] org.Mm.eg.db_2.8.0 limma_3.14.3 gcrma_2.30.0 >> [7] annaffy_1.30.0 KEGG.db_2.8.0 GO.db_2.8.0 >> [10] RSQLite_0.11.2 DBI_0.2-5 AnnotationDbi_1.20.3 >> [13] affy_1.36.0 Biobase_2.18.0 BiocGenerics_0.4.0 >> >> loaded via a namespace (and not attached): >> [1] affyio_1.26.0 BiocInstaller_1.8.3 Biostrings_2.26.2 >> [4] IRanges_1.16.4 parallel_2.15.2 preprocessCore_1.20.0 >> [7] splines_2.15.2 stats4_2.15.2 tools_2.15.2 >> [10] zlibbioc_1.4.0 >> >> >> (I have not had the chance to switch to R3.0 yet and am waiting a while until >> any bugs get worked out). >> >> Here is my script: >> >> library(affy) >> library(annaffy) >> library(gcrma) >> library(limma) >> library(mouse4302.db) >> raw<-ReadAffy() >> gcrmanm<-gcrma(raw) >> write.exprs(gcrmanm,"gloria2gcexp.txt") >> >> >> targets<-readTargets("targets.txt") >> targets >> condition<-factor(targets$condition,levels=c("AA.het.con","AB.het.s ta","BA.mut.con","BB.mut.sta")) >> condition >> mouse<-factor(targets$mouse) >> mouse >> design<-model.matrix(~mouse+condition) >> design >> colnames(design)<-c("mouse2","mouse3","mouse4", >> "mouse5","mouse6","mouse7","mouse8", >> "mouse9","mouse10","mouse11","mouse12", >> "AB.het.sta","BA.mut.con","BB.mut.sta") >> >> fit<-lmFit(gcrmanm,design) >> >> contrast.matrix<-makeContrasts(AB.het.sta-AA.het.con,BB.mut.sta- BA.mut.con, >> ((BB.mut.sta-BA.mut.con)-(AB.het.sta-AA.het.con)),levels=design) >> contrast.matrix >> fit<-contrasts.fit(fit,contrast.matrix) >> fit<-eBayes(fit) >> >> comp1<-topTable(fit,coef=1,number=45101,adjust.method="BH",sort.by= "B") >> >> >> names(comp1) >> names(comp1)[1]="Probe" >> probeids=comp1$Probe >> >> aaf.handler() >> anncols<-aaf.handler()[c(1:5,11:12)] >> anncols >> >> comp1table<-aafTableFrame(comp1,colnames=names(comp1),probeids=prob eids) >> anntable<-aafTableAnn(probeids, "mouse4302.db", colnames=anncols) >> comp1ann<-merge(comp1table,anntable) >> saveText(comp1ann,"AB.het.staVsAA.het.con.ann.txt") >> >> >> comp2<-topTable(fit,coef=2,number=45101,adjust.method="BH",sort.by= "B") >> >> names(comp2) >> names(comp2)[1]="Probe" >> probeids=comp2$Probe >> >> aaf.handler() >> anncols<-aaf.handler()[c(1:5,11:12)] >> anncols >> >> comp2table<-aafTableFrame(comp2,colnames=names(comp2),probeids=prob eids) >> anntable<-aafTableAnn(probeids, "mouse4302.db", colnames=anncols) >> comp2ann<-merge(comp2table,anntable) >> saveText(comp1ann,"BB.mut.staVsBA.mut.con.ann.txt") >> >> comp3<-topTable(fit,coef=3,number=45101,adjust.method="BH",sort.by= "B") >> >> names(comp3) >> names(comp3)[1]="Probe" >> probeids=comp3$Probe >> >> aaf.handler() >> anncols<-aaf.handler()[c(1:5,11:12)] >> anncols >> >> comp3table<-aafTableFrame(comp3,names(comp3),probeids=probeids) >> anntable<-aafTableAnn(probeids, "mouse4302.db", colnames=anncols) >> comp3ann<-merge(comp3table,anntable) >> saveText(comp3ann,"starvationresponse.mutVshet.ann.txt") >> >> >> I get an error message at the point at which I am naming the columns: >> >>> colnames(design)<-c("mouse2","mouse3","mouse4", >> + "mouse5","mouse6","mouse7","mouse8", >> + "mouse9","mouse10","mouse11","mouse12", >> + "AB.het.sta","BA.mut.con","BB.mut.sta") >> Error in `colnames<-`(`*tmp*`, value = c("mouse2", "mouse3", "mouse4", : >> length of 'dimnames' [2] not equal to array extent >> Calls: colnames<- -> colnames<- >> Execution halted >> >> I then tried the following revised script: >> >> library(affy) >> library(annaffy) >> library(gcrma) >> library(limma) >> library(mouse4302.db) >> raw<-ReadAffy() >> gcrmanm<-gcrma(raw) >> write.exprs(gcrmanm,"gloria2gcexp.txt") >> >> >> targets<-readTargets("targets.txt") >> targets >> condition<-factor(targets$condition,levels=c("AA.het.con","AB.het.s ta","BA.mut.con","BB.mut.sta")) >> condition >> mouse<-factor(targets$mouse) >> mouse >> design<-model.matrix(~0+condition+mouse) >> design >> colnames(design)<-c("AA.het.con","AB.het.sta","BA.mut.con","BB.mut. sta", >> "mouse2","mouse3","mouse4", >> "mouse5","mouse6","mouse7","mouse8", >> "mouse9","mouse10","mouse11","mouse12") >> fit<-lmFit(gcrmanm,design) >> >> (everything else the same). >> >> I get the following warning message: >> >>> colnames(design)<-c("AA.het.con","AB.het.sta","BA.mut.con","BB.mut .sta", >> + "mouse2","mouse3","mouse4", >> + "mouse5","mouse6","mouse7","mouse8", >> + "mouse9","mouse10","mouse11","mouse12") >>> fit<-lmFit(gcrmanm,design) >> Coefficients not estimable: mouse10 >> Warning message: >> Partial NA coefficients for 45101 probe(s) >> >> THE SCRIPT RAN TO COMPLETION >> >> MY QUESTIONS: >> 1. Should I worry about the warning message? >> 2. Is there any way in which I can do this analysis and not get the warning message. >> >> Thanks and best wishes, >> Rich >> Richard A. Friedman, PhD >> Associate Research Scientist, >> Biomedical Informatics Shared Resource >> Herbert Irving Comprehensive Cancer Center (HICCC) >> Lecturer, >> Department of Biomedical Informatics (DBMI) >> Educational Coordinator, >> Center for Computational Biology and Bioinformatics (C2B2)/ >> National Center for Multiscale Analysis of Genomic Networks (MAGNet)/ >> Columbia Initiative in Systems Biology >> Room 824 >> Irving Cancer Research Center >> Columbia University >> 1130 St. Nicholas Ave >> New York, NY 10032 >> (212)851-4765 (voice) >> friedman at cancercenter.columbia.edu >> http://cancercenter.columbia.edu/~friedman/ >> >> Fritz Lang. Didn't he do "Star Trek". >> -Rose Friedman, age 16 >> > > ______________________________________________________________________ > The information in this email is confidential and inte...{{dropped:6}}
ADD REPLY
0
Entering edit mode
@gordon-smyth
Last seen 10 hours ago
WEHI, Melbourne, Australia
Dear Richard, I wrote a little quickly the first time. The warnings are serious in this context. Looking more closely at the experimental design, I see now that only two comparisons are being made at the within-mouse level. For six mice, the AA.het.con vs AB.het.sta comparison is made. For six other mice, the BA.mut.con vs BB.mut.sta comparison made. Hence it is only possible to estimate two coefficients between conditions in the paired design. Hence the model ~mouse+condition will lead to a warning because it sets up three coefficients for differences between the conditions whereas only two are estimable. There is no automatic way in R to create the design matrix for a special design like this. You have to do it by hand. First setup the mouse effects: design <- model.matrix(~mouse) Then add the two estimable comparisons: het.sta.vs.con <- (condition=="AB.het.sta") mut.sta.vs.con <- (condition=="BB.mut.sta") design <- cbind(design,het.sta.vs.con,mut.sta.vs.con) Then the contrasts will be: contrast.matrix <- makeContrasts( het.sta.vs.con, mut.sta.vs.con, mut.sta.vs.con-het.sta.vs.con, levels=design) Best wishes Gordon > Date: Mon, 8 Apr 2013 12:57:35 -0400 > From: Richard Friedman <friedman at="" cancercenter.columbia.edu=""> > To: Gordon K Smyth <smyth at="" wehi.edu.au=""> > Cc: Bioconductor mailing list <bioconductor at="" r-project.org=""> > Subject: Re: [BioC] error and warning message with blocked samples > with Limma I AM STILL HAVING DIFFICULTY > > On Apr 6, 2013, at 7:41 PM, Gordon K Smyth wrote: > >> Dear Richard, >> >> The error is coming when you set colnames for a matrix. When you try >> to set colnames for matrix, you must give the correct number of names >> as there are columns in the matrix, and the error is saying that you >> are not doing so. Looking at your code, it would seem that you have >> ommited the Intercept from your list of names. >> >> Even if you fixed this, the code would still fail later on because you >> use "AA.het.con" in makeContrasts() but there is no column by this name >> in your design matrix. >> >> Your second analysis runs because you give 15 column names instead of >> 14. >> >> You can get the three contrasts you want by >> >> design <- model.matrix(~mouse+condition) >> contrast.matrix <- makeContrasts( >> "AB.het.sta-AA.het.con" = conditionAB.het.sta, >> "BB.mut.sta-BA.mut.con" = conditionBB.mut.sta-conditionBA.mut.con, >> "Diff" = (conditionBB.mut.sta- conditionBA.mut.con)-conditionAB.het.sta, >> levels=design) >> >> >> Best wishes >> Gordon > > Dear Gordon, > > Thank you for your reply. I am still having difficulty. I have tried > both putting the code you recommended in and I also tried with the intercept, > but I still get error messages. Is the script I wrote with > design<-model.matrix(~0+condition+mouse) > that ran all the way through reliable? > > Here is an excerpt from the output from just putting the code above in: > > ################################################################## > > >> design<-model.matrix(~mouse+condition) >> design > (Intercept) mouse2 mouse3 mouse4 mouse5 mouse6 mouse7 mouse8 mouse9 mouse10 > 1 1 0 0 0 0 0 0 0 0 0 > 2 1 0 0 0 0 0 0 0 0 0 > 3 1 1 0 0 0 0 0 0 0 0 > 4 1 1 0 0 0 0 0 0 0 0 > 5 1 0 1 0 0 0 0 0 0 0 > 6 1 0 1 0 0 0 0 0 0 0 > 7 1 0 0 1 0 0 0 0 0 0 > 8 1 0 0 1 0 0 0 0 0 0 > 9 1 0 0 0 1 0 0 0 0 0 > 10 1 0 0 0 1 0 0 0 0 0 > 11 1 0 0 0 0 1 0 0 0 0 > 12 1 0 0 0 0 1 0 0 0 0 > 13 1 0 0 0 0 0 1 0 0 0 > 14 1 0 0 0 0 0 1 0 0 0 > 15 1 0 0 0 0 0 0 1 0 0 > 16 1 0 0 0 0 0 0 1 0 0 > 17 1 0 0 0 0 0 0 0 1 0 > 18 1 0 0 0 0 0 0 0 1 0 > 19 1 0 0 0 0 0 0 0 0 1 > 20 1 0 0 0 0 0 0 0 0 1 > 21 1 0 0 0 0 0 0 0 0 0 > 22 1 0 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 > mouse11 mouse12 conditionAB.het.sta conditionBA.mut.con conditionBB.mut.sta > 1 0 0 0 0 0 > 2 0 0 1 0 0 > 3 0 0 0 0 0 > 4 0 0 1 0 0 > 5 0 0 0 1 0 > 6 0 0 0 0 1 > 7 0 0 0 1 0 > 8 0 0 0 0 1 > 9 0 0 0 0 0 > 10 0 0 1 0 0 > 11 0 0 0 0 0 > 12 0 0 1 0 0 > 13 0 0 0 1 0 > 14 0 0 0 0 1 > 15 0 0 0 1 0 > 16 0 0 0 0 1 > 17 0 0 0 1 0 > 18 0 0 0 0 1 > 19 0 0 0 1 0 > 20 0 0 0 0 1 > 21 1 0 0 0 0 > 22 1 0 1 0 0 > 23 0 1 0 0 0 > 24 0 1 1 0 0 > attr(,"assign") > [1] 0 1 1 1 1 1 1 1 1 1 1 1 2 2 2 > attr(,"contrasts") > attr(,"contrasts")$mouse > [1] "contr.treatment" > > attr(,"contrasts")$condition > [1] "contr.treatment" > >> >> fit<-lmFit(gcrmanm,design) > Coefficients not estimable: conditionBB.mut.sta > Warning message: > Partial NA coefficients for 45101 probe(s) >> >> contrast.matrix <- makeContrasts( > + "AB.het.sta-AA.het.con" = conditionAB.het.sta, > + "BB.mut.sta-BA.mut.con" = conditionBB.mut.sta-conditionBA.mut.con, > + "Diff" = (conditionBB.mut.sta- conditionBA.mut.con)-conditionAB.het.sta, > + levels=design) > Warning message: > In makeContrasts(`AB.het.sta-AA.het.con` = conditionAB.het.sta, : > Renaming (Intercept) to Intercept >> >> >> contrast.matrix > Contrasts > Levels AB.het.sta-AA.het.con BB.mut.sta-BA.mut.con Diff > Intercept 0 0 0 > mouse2 0 0 0 > mouse3 0 0 0 > mouse4 0 0 0 > mouse5 0 0 0 > mouse6 0 0 0 > mouse7 0 0 0 > mouse8 0 0 0 > mouse9 0 0 0 > mouse10 0 0 0 > mouse11 0 0 0 > mouse12 0 0 0 > conditionAB.het.sta 1 0 -1 > conditionBA.mut.con 0 -1 -1 > conditionBB.mut.sta 0 1 1 >> fit<-contrasts.fit(fit,contrast.matrix) > Error in contrasts.fit(fit, contrast.matrix) : > trying to take contrast of non-estimable coefficient > In addition: Warning message: > In contrasts.fit(fit, contrast.matrix) : > row names of contrasts don't match col names of coefficients > Execution halted > > ########################################################### > > Here is an excerpt from the output from my explicitly naming the intercept > and naming the columns: > > Levels: 1 2 3 4 5 6 7 8 9 10 11 12 >> design<-model.matrix(~mouse+condition) >> design > (Intercept) mouse2 mouse3 mouse4 mouse5 mouse6 mouse7 mouse8 mouse9 mouse10 > 1 1 0 0 0 0 0 0 0 0 0 > 2 1 0 0 0 0 0 0 0 0 0 > 3 1 1 0 0 0 0 0 0 0 0 > 4 1 1 0 0 0 0 0 0 0 0 > 5 1 0 1 0 0 0 0 0 0 0 > 6 1 0 1 0 0 0 0 0 0 0 > 7 1 0 0 1 0 0 0 0 0 0 > 8 1 0 0 1 0 0 0 0 0 0 > 9 1 0 0 0 1 0 0 0 0 0 > 10 1 0 0 0 1 0 0 0 0 0 > 11 1 0 0 0 0 1 0 0 0 0 > 12 1 0 0 0 0 1 0 0 0 0 > 13 1 0 0 0 0 0 1 0 0 0 > 14 1 0 0 0 0 0 1 0 0 0 > 15 1 0 0 0 0 0 0 1 0 0 > 16 1 0 0 0 0 0 0 1 0 0 > 17 1 0 0 0 0 0 0 0 1 0 > 18 1 0 0 0 0 0 0 0 1 0 > 19 1 0 0 0 0 0 0 0 0 1 > 20 1 0 0 0 0 0 0 0 0 1 > 21 1 0 0 0 0 0 0 0 0 0 > 22 1 0 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 > mouse11 mouse12 conditionAB.het.sta conditionBA.mut.con conditionBB.mut.sta > 1 0 0 0 0 0 > 2 0 0 1 0 0 > 3 0 0 0 0 0 > 4 0 0 1 0 0 > 5 0 0 0 1 0 > 6 0 0 0 0 1 > 7 0 0 0 1 0 > 8 0 0 0 0 1 > 9 0 0 0 0 0 > 10 0 0 1 0 0 > 11 0 0 0 0 0 > 12 0 0 1 0 0 > 13 0 0 0 1 0 > 14 0 0 0 0 1 > 15 0 0 0 1 0 > 16 0 0 0 0 1 > 17 0 0 0 1 0 > 18 0 0 0 0 1 > 19 0 0 0 1 0 > 20 0 0 0 0 1 > 21 1 0 0 0 0 > 22 1 0 1 0 0 > 23 0 1 0 0 0 > 24 0 1 1 0 0 > attr(,"assign") > [1] 0 1 1 1 1 1 1 1 1 1 1 1 2 2 2 > attr(,"contrasts") > attr(,"contrasts")$mouse > [1] "contr.treatment" > > attr(,"contrasts")$condition > [1] "contr.treatment" > >> >> colnames(design)<-c("Intercept","mouse2","mouse3","mouse4", > + "mouse5","mouse6","mouse7","mouse8", > + "mouse9","mouse10","mouse11","mouse12", > + "AB.het.sta","BA.mut.con","BB.mut.sta") >> >> >> fit<-lmFit(gcrmanm,design) > Coefficients not estimable: BB.mut.sta > Warning message: > Partial NA coefficients for 45101 probe(s) >> >> #contrast.matrix <- makeContrasts( >> #"AB.het.sta-AA.het.con" = conditionAB.het.sta, >> #"BB.mut.sta-BA.mut.con" = conditionBB.mut.sta-conditionBA.mut.con, >> #"Diff" = (conditionBB.mut.sta- conditionBA.mut.con)-conditionAB.het.sta, >> #levels=design) >> >> contrast.matrix <- makeContrasts(AB.het.sta,BB.mut.sta-BA.mut.con, > + (BB.mut.sta-BA.mut.con)-AB.het.sta,levels=design) >> >> >> contrast.matrix > Contrasts > Levels AB.het.sta BB.mut.sta - BA.mut.con > Intercept 0 0 > mouse2 0 0 > mouse3 0 0 > mouse4 0 0 > mouse5 0 0 > mouse6 0 0 > mouse7 0 0 > mouse8 0 0 > mouse9 0 0 > mouse10 0 0 > mouse11 0 0 > mouse12 0 0 > AB.het.sta 1 0 > BA.mut.con 0 -1 > BB.mut.sta 0 1 > Contrasts > Levels (BB.mut.sta - BA.mut.con) - AB.het.sta > Intercept 0 > mouse2 0 > mouse3 0 > mouse4 0 > mouse5 0 > mouse6 0 > mouse7 0 > mouse8 0 > mouse9 0 > mouse10 0 > mouse11 0 > mouse12 0 > AB.het.sta -1 > BA.mut.con -1 > BB.mut.sta 1 >> fit<-contrasts.fit(fit,contrast.matrix) > Error in contrasts.fit(fit, contrast.matrix) : > trying to take contrast of non-estimable coefficient > Execution halted > > ################################################# > > I am clearly missing something. > Again, I would appreciate any suggestions. > > Thanks and best wishes, > Rich > > > >> >>> Date: Fri, 5 Apr 2013 10:42:22 -0400 >>> From: Richard Friedman <friedman at="" cancercenter.columbia.edu=""> >>> To: "bioconductor at r-project.org list" <bioconductor at="" r-project.org=""> >>> Subject: [BioC] error and warning message with blocked samples with >>> Limma >>> >>> Dear Bioconductor List, >>> >>> I want to do a linear model analysis of a blocked dataset >>> in Limma. I get an error message if I closely follow the manual, >>> or a warning message if I do it a little differentially. I understand >>> from previous posts that the warning message may not be serious, >>> but I would like to make sure that I am doing this correctly. >>> >>> My dataset is: >>> >>> ileName Name mouse condition >>> QIU01.CEL AA.het.con.01.01 1 AA.het.con >>> QIU02.CEL AB.het.sta.01.02 1 AB.het.sta >>> QIU03.CEL AA.het.con.02.03 2 AA.het.con >>> QIU04.CEL AB.het.sta.02.04 2 AB.het.sta >>> QIU05.CEL BA.mut.con.03.05 3 BA.mut.con >>> QIU06.CEL BB.mut.sta.03.06 3 BB.mut.sta >>> QIU07.CEL BA.mut.con.04.07 4 BA.mut.con >>> QIU08.CEL BB.mutsta.04.08 4 BB.mut.sta >>> QIU09.CEL AA.het.con.05.09 5 AA.het.con >>> QIU10.CEL AB.het.sta.05.10 5 AB.het.sta >>> QIU11.CEL AA.het.con.06.11 6 AA.het.con >>> QIU12.CEL AB.het.sta.06.12 6 AB.het.sta >>> QIU13.CEL BA.mut.con.07.13 7 BA.mut.con >>> QIU14.CEL BB.mut.sta.07.14 7 BB.mut.sta >>> QIU15.CEL BA.mut.con.08.15 8 BA.mut.con >>> QIU16.CEL BB.mut.sta.08.16 8 BB.mut.sta >>> QIU17.CEL BA.mut.con.09.17 9 BA.mut.con >>> QIU18.CEL BB.mut.sta.09.18 9 BB.mut.sta >>> QIU19.CEL BA.mut.con.10.19 10 BA.mut.con >>> QIU20.CEL BB.mut.sta.10.20 10 BB.mut.sta >>> QIU21.CEL AA.het.con.11.21 11 AA.het.con >>> QIU22.CEL AB.het.sta.11.22 11 AB.het.sta >>> QIU23.CEL AA.het.con.12.23 12 AA.het.con >>> QIU24.CEL AB.het.sta.12.24 12 AB.het.sta >>> >>> Arrays from the same mouse are paired but they are >>> treated differently. I need the following 3 comparisons: >>> >>> AB.het.sta-AA.het.con >>> BB.mut.sta-BA.mut.con >>> ((BB.mut.sta-BA.mut.con)-(AB.het.sta-AA.het.con)) >>> >>> Here is my environment: >>> >>>> sessionInfo() >>> R version 2.15.2 (2012-10-26) >>> 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] mouse4302.db_2.8.1 mouse4302probe_2.11.0 mouse4302cdf_2.11.0 >>> [4] org.Mm.eg.db_2.8.0 limma_3.14.3 gcrma_2.30.0 >>> [7] annaffy_1.30.0 KEGG.db_2.8.0 GO.db_2.8.0 >>> [10] RSQLite_0.11.2 DBI_0.2-5 AnnotationDbi_1.20.3 >>> [13] affy_1.36.0 Biobase_2.18.0 BiocGenerics_0.4.0 >>> >>> loaded via a namespace (and not attached): >>> [1] affyio_1.26.0 BiocInstaller_1.8.3 Biostrings_2.26.2 >>> [4] IRanges_1.16.4 parallel_2.15.2 preprocessCore_1.20.0 >>> [7] splines_2.15.2 stats4_2.15.2 tools_2.15.2 >>> [10] zlibbioc_1.4.0 >>> >>> >>> (I have not had the chance to switch to R3.0 yet and am waiting a while until >>> any bugs get worked out). >>> >>> Here is my script: >>> >>> library(affy) >>> library(annaffy) >>> library(gcrma) >>> library(limma) >>> library(mouse4302.db) >>> raw<-ReadAffy() >>> gcrmanm<-gcrma(raw) >>> write.exprs(gcrmanm,"gloria2gcexp.txt") >>> >>> >>> targets<-readTargets("targets.txt") >>> targets >>> condition<-factor(targets$condition,levels=c("AA.het.con","AB.het. sta","BA.mut.con","BB.mut.sta")) >>> condition >>> mouse<-factor(targets$mouse) >>> mouse >>> design<-model.matrix(~mouse+condition) >>> design >>> colnames(design)<-c("mouse2","mouse3","mouse4", >>> "mouse5","mouse6","mouse7","mouse8", >>> "mouse9","mouse10","mouse11","mouse12", >>> "AB.het.sta","BA.mut.con","BB.mut.sta") >>> >>> fit<-lmFit(gcrmanm,design) >>> >>> contrast.matrix<-makeContrasts(AB.het.sta-AA.het.con,BB.mut.sta- BA.mut.con, >>> ((BB.mut.sta-BA.mut.con)-(AB.het.sta-AA.het.con)),levels=design) >>> contrast.matrix >>> fit<-contrasts.fit(fit,contrast.matrix) >>> fit<-eBayes(fit) >>> >>> comp1<-topTable(fit,coef=1,number=45101,adjust.method="BH",sort.by ="B") >>> >>> >>> names(comp1) >>> names(comp1)[1]="Probe" >>> probeids=comp1$Probe >>> >>> aaf.handler() >>> anncols<-aaf.handler()[c(1:5,11:12)] >>> anncols >>> >>> comp1table<-aafTableFrame(comp1,colnames=names(comp1),probeids=pro beids) >>> anntable<-aafTableAnn(probeids, "mouse4302.db", colnames=anncols) >>> comp1ann<-merge(comp1table,anntable) >>> saveText(comp1ann,"AB.het.staVsAA.het.con.ann.txt") >>> >>> >>> comp2<-topTable(fit,coef=2,number=45101,adjust.method="BH",sort.by ="B") >>> >>> names(comp2) >>> names(comp2)[1]="Probe" >>> probeids=comp2$Probe >>> >>> aaf.handler() >>> anncols<-aaf.handler()[c(1:5,11:12)] >>> anncols >>> >>> comp2table<-aafTableFrame(comp2,colnames=names(comp2),probeids=pro beids) >>> anntable<-aafTableAnn(probeids, "mouse4302.db", colnames=anncols) >>> comp2ann<-merge(comp2table,anntable) >>> saveText(comp1ann,"BB.mut.staVsBA.mut.con.ann.txt") >>> >>> comp3<-topTable(fit,coef=3,number=45101,adjust.method="BH",sort.by ="B") >>> >>> names(comp3) >>> names(comp3)[1]="Probe" >>> probeids=comp3$Probe >>> >>> aaf.handler() >>> anncols<-aaf.handler()[c(1:5,11:12)] >>> anncols >>> >>> comp3table<-aafTableFrame(comp3,names(comp3),probeids=probeids) >>> anntable<-aafTableAnn(probeids, "mouse4302.db", colnames=anncols) >>> comp3ann<-merge(comp3table,anntable) >>> saveText(comp3ann,"starvationresponse.mutVshet.ann.txt") >>> >>> >>> I get an error message at the point at which I am naming the columns: >>> >>>> colnames(design)<-c("mouse2","mouse3","mouse4", >>> + "mouse5","mouse6","mouse7","mouse8", >>> + "mouse9","mouse10","mouse11","mouse12", >>> + "AB.het.sta","BA.mut.con","BB.mut.sta") >>> Error in `colnames<-`(`*tmp*`, value = c("mouse2", "mouse3", "mouse4", : >>> length of 'dimnames' [2] not equal to array extent >>> Calls: colnames<- -> colnames<- >>> Execution halted >>> >>> I then tried the following revised script: >>> >>> library(affy) >>> library(annaffy) >>> library(gcrma) >>> library(limma) >>> library(mouse4302.db) >>> raw<-ReadAffy() >>> gcrmanm<-gcrma(raw) >>> write.exprs(gcrmanm,"gloria2gcexp.txt") >>> >>> >>> targets<-readTargets("targets.txt") >>> targets >>> condition<-factor(targets$condition,levels=c("AA.het.con","AB.het. sta","BA.mut.con","BB.mut.sta")) >>> condition >>> mouse<-factor(targets$mouse) >>> mouse >>> design<-model.matrix(~0+condition+mouse) >>> design >>> colnames(design)<-c("AA.het.con","AB.het.sta","BA.mut.con","BB.mut .sta", >>> "mouse2","mouse3","mouse4", >>> "mouse5","mouse6","mouse7","mouse8", >>> "mouse9","mouse10","mouse11","mouse12") >>> fit<-lmFit(gcrmanm,design) >>> >>> (everything else the same). >>> >>> I get the following warning message: >>> >>>> colnames(design)<-c("AA.het.con","AB.het.sta","BA.mut.con","BB.mu t.sta", >>> + "mouse2","mouse3","mouse4", >>> + "mouse5","mouse6","mouse7","mouse8", >>> + "mouse9","mouse10","mouse11","mouse12") >>>> fit<-lmFit(gcrmanm,design) >>> Coefficients not estimable: mouse10 >>> Warning message: >>> Partial NA coefficients for 45101 probe(s) >>> >>> THE SCRIPT RAN TO COMPLETION >>> >>> MY QUESTIONS: >>> 1. Should I worry about the warning message? >>> 2. Is there any way in which I can do this analysis and not get the warning message. >>> >>> Thanks and best wishes, >>> Rich >>> Richard A. Friedman, PhD >>> Associate Research Scientist, >>> Biomedical Informatics Shared Resource >>> Herbert Irving Comprehensive Cancer Center (HICCC) >>> Lecturer, >>> Department of Biomedical Informatics (DBMI) >>> Educational Coordinator, >>> Center for Computational Biology and Bioinformatics (C2B2)/ >>> National Center for Multiscale Analysis of Genomic Networks (MAGNet)/ >>> Columbia Initiative in Systems Biology >>> Room 824 >>> Irving Cancer Research Center >>> Columbia University >>> 1130 St. Nicholas Ave >>> New York, NY 10032 >>> (212)851-4765 (voice) >>> friedman at cancercenter.columbia.edu >>> http://cancercenter.columbia.edu/~friedman/ >>> >>> Fritz Lang. Didn't he do "Star Trek". >>> -Rose Friedman, age 16 >>> ______________________________________________________________________ The information in this email is confidential and intend...{{dropped:4}}
ADD COMMENT

Login before adding your answer.

Traffic: 801 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6