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}}
