[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