[BioC] Problem in limma-voom
Gordon K Smyth
smyth at wehi.EDU.AU
Wed Feb 19 07:32:06 CET 2014
Hi Ming,
This has nothing to do with voom. It is the same message that you will
always receive from limma when you try to fit an over-parametrized model.
(By comparison, edgeR would have simply told you the design matrix was not
of full rank and stopped.)
You can only used a paired analysis (+Patient etc) when the treatments to
be compared are applied to the same patient. You can't compare different
tumour types (KRas and KRasP) which affect different patients.
Gordon
On Wed, 19 Feb 2014, Ming Yi wrote:
> Hi, Dr. Gordon and all:
>
>
>
> I used limma before, but not much on limma-voom package. I run into some
> issue of warnings and then errors, not sure what is. Below are my
> commands, only show essential parts:
>
>
>
>> show(head(targets))
>
> ID RasType Subject Type
> 1 T4626_01A KRas TCGA_38_4626 Tumor
> 2 N4626_11A KRas TCGA_38_4626 Normal
> 3 T2668_01A KRas TCGA_44_2668 Tumor
> 4 N2668_11A KRas TCGA_44_2668 Normal
> 5 T6145_01A KRas TCGA_44_6145 Tumor
>
> > RasTum<-paste(targets$RasType, targets$Type,sep=".")
>> RasTum<-factor(RasTum, levels=c("KRas.Tumor","KRas.Normal","KRasP.Tumor","KRasP.Normal"));
>> RasTum
> [1] KRas.Tumor KRas.Normal KRas.Tumor KRas.Normal KRas.Tumor
> [6] KRas.Normal KRas.Tumor KRas.Normal KRas.Tumor KRas.Normal
> [11] KRas.Tumor KRas.Normal KRas.Tumor KRas.Normal KRas.Tumor
> [16] KRas.Normal KRas.Tumor KRas.Normal KRas.Tumor KRas.Normal
> [21] KRas.Tumor KRas.Normal KRas.Tumor KRas.Normal KRasP.Tumor
> [26] KRasP.Normal KRasP.Tumor KRasP.Normal KRasP.Tumor KRasP.Normal
> [31] KRasP.Tumor KRasP.Normal KRasP.Tumor KRasP.Normal KRasP.Tumor
> [36] KRasP.Normal KRasP.Tumor KRasP.Normal KRasP.Tumor KRasP.Normal
> [41] KRasP.Tumor KRasP.Normal KRasP.Tumor KRasP.Normal KRasP.Tumor
> [46] KRasP.Normal KRasP.Tumor KRasP.Normal KRasP.Tumor KRasP.Normal
> [51] KRasP.Tumor KRasP.Normal KRasP.Tumor KRasP.Normal KRasP.Tumor
> [56] KRasP.Normal KRasP.Tumor KRasP.Normal
> Levels: KRas.Tumor KRas.Normal KRasP.Tumor KRasP.Normal
>
>> Patient<-factor(targets$Subject);
>
>> rawdata<-rawdata[,c("symbols","entrezID" ,"gene_id",targets$ID)];
>
>> y<-DGEList(count=rawdata[,4:ncol(rawdata)],genes=rawdata[,1:3]);
>
>> show(dim(y));
> [1] 20531 58
>
>> y3<-y[rowSums(y$counts)>=50,];
>
>> show(dim(y3));
> [1] 18215 58
>> y3$samples$lib.size <- colSums(y3$counts)
>> rownames(y3$counts) <- rownames(y3$genes) <- y3$genes$entrezID
>> y3 <- calcNormFactors(y3,method=c("TMM"));
>> show(y3$samples[1:6,])
> group lib.size norm.factors
> T4626_01A 1 97435132 0.9834587
> N4626_11A 1 62583672 0.9154837
> T2668_01A 1 77034746 0.9687860
> N2668_11A 1 58631311 0.8644352
> T6145_01A 1 31634419 1.0303445
> N6145_11A 1 29872761 0.9615021
>> design<-model.matrix(~Patient+RasTum);
>> show(design[1:6,c(1:10,30:32)])
> (Intercept) PatientTCGA_38_4626 PatientTCGA_38_4627 PatientTCGA_38_4632
> 1 1 1 0 0
> 2 1 1 0 0
> 3 1 0 0 0
> 4 1 0 0 0
> 5 1 0 0 0
> 6 1 0 0 0
> PatientTCGA_44_2655 PatientTCGA_44_2657 PatientTCGA_44_2661
> 1 0 0 0
> 2 0 0 0
> 3 0 0 0
> 4 0 0 0
> 5 0 0 0
> 6 0 0 0
> PatientTCGA_44_2662 PatientTCGA_44_2665 PatientTCGA_44_2668 RasTumKRas.Normal
> 1 0 0 0 0
> 2 0 0 0 1
> 3 0 0 1 0
> 4 0 0 1 1
> 5 0 0 0 0
> 6 0 0 0 1
> RasTumKRasP.Tumor RasTumKRasP.Normal
> 1 0 0
> 2 0 0
> 3 0 0
> 4 0 0
> 5 0 0
> 6 0 0
>> colnames(design)<-sub("Patient","",colnames(design));
>> colnames(design)<-sub("RasTum","",colnames(design));
>> colnames(design)[1] <- "Intercept"
>> show(design[1:6,c(1:10,30:32)])
> Intercept TCGA_38_4626 TCGA_38_4627 TCGA_38_4632 TCGA_44_2655 TCGA_44_2657
> 1 1 1 0 0 0 0
> 2 1 1 0 0 0 0
> 3 1 0 0 0 0 0
> 4 1 0 0 0 0 0
> 5 1 0 0 0 0 0
> 6 1 0 0 0 0 0
> TCGA_44_2661 TCGA_44_2662 TCGA_44_2665 TCGA_44_2668 KRas.Normal KRasP.Tumor
> 1 0 0 0 0 0 0
> 2 0 0 0 0 1 0
> 3 0 0 0 1 0 0
> 4 0 0 0 1 1 0
> 5 0 0 0 0 0 0
> 6 0 0 0 0 1 0
> KRasP.Normal
> 1 0
> 2 0
> 3 0
> 4 0
> 5 0
> 6 0
>> dim(design)
> [1] 58 32
>> con.matrix<-makeContrasts(
> + KRas.Tumor_KRas.Normal = -KRas.Normal,
> + KRasP.Tumor_KRasP.Normal = KRasP.Tumor-KRasP.Normal,
> + KRas.Tumor_KRasP.Tumor=-KRasP.Tumor,
> + levels=design)
>> v <- voom(y3,design);
> Coefficients not estimable: KRasP.Normal
> Warning message:
> Partial NA coefficients for 18215 probe(s)
>
>> fit<- lmFit(v,design)
> Coefficients not estimable: KRasP.Normal
> Warning message:
> Partial NA coefficients for 18215 probe(s)
>> fit1<-contrasts.fit(fit, con.matrix)
> Error in contrasts.fit(fit, con.matrix) :
> trying to take contrast of non-estimable coefficient
>
>
> I used a similar way before for array data, but not with limma-voom for RNAseq data. Any suggestions here?
>
>
>
> Thanks a lot in advance!
>
>
>
> Ming
>
>
>
> ATRF/NCI-Frederick
>
> Maryland, USA
>
>
>
> Here is the sessionInfo:
>
>> show(sessionInfo())
> R version 3.0.1 (2013-05-16)
> Platform: x86_64-unknown-linux-gnu (64-bit)
>
> locale:
> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
> [7] LC_PAPER=C LC_NAME=C
> [9] LC_ADDRESS=C LC_TELEPHONE=C
> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
>
> attached base packages:
> [1] stats graphics grDevices utils datasets methods base
>
> other attached packages:
> [1] edgeR_3.4.2 limma_3.16.8
>
>
>
>
>
>
>
>
>
>
______________________________________________________________________
The information in this email is confidential and intend...{{dropped:4}}
More information about the Bioconductor
mailing list