[BioC] error and warning message with blocked samples with Limma
Gordon K Smyth
smyth at wehi.EDU.AU
Sun Apr 7 01:41:23 CEST 2013
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.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 intend...{{dropped:4}}
More information about the Bioconductor
mailing list