[BioC] error and warning message with blocked samples with Limma I AM STILL HAVING DIFFICULTY
Richard Friedman
friedman at cancercenter.columbia.edu
Mon Apr 8 18:57:35 CEST 2013
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=probeids)
>> 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=probeids)
>> 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.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}}
More information about the Bioconductor
mailing list