[BioC] error and warning message with blocked samples with Limma
Richard Friedman
friedman at cancercenter.columbia.edu
Fri Apr 5 16:42:22 CEST 2013
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
More information about the Bioconductor
mailing list