[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