Hi, Dr. Gordon and all:

 

I run into some limma design issue not sure what they are:

 

I had three tumor types want to compare with each other: RasP, RasP1Hit, RasPNot

 

here are the critical commands:

>targets<-targets[targets$RasType!="RasOnly" & targets$Type=="Tumor",]

> show(head(targets));
        Subject SampleName  Type RasPType RasType
1  TCGA_38_4625  T4625_01A Tumor     RasP    RasP
5  TCGA_38_4627  T4627_01A Tumor     RasP    RasP
7  TCGA_38_4632  T4632_01A Tumor     RasP    RasP
9  TCGA_44_2655  T2655_01A Tumor     RasP    RasP
11 TCGA_44_2657  T2657_01A Tumor     RasP    RasP
13 TCGA_44_2661  T2661_01A Tumor     RasP    RasP

 

> RasTum<-paste(targets$RasType, targets$Type,sep=".")
> RasTum<-factor(RasTum, levels=c("RasP1Hit.Tumor","RasP.Tumor","RasPNot.Tumor"));
> show(RasTum);
 [1] RasP.Tumor     RasP.Tumor     RasP.Tumor     RasP.Tumor     RasP.Tumor
 [6] RasP.Tumor     RasP.Tumor     RasP.Tumor     RasP1Hit.Tumor RasP1Hit.Tumor
[11] RasP.Tumor     RasPNot.Tumor  RasP1Hit.Tumor RasP.Tumor     RasP1Hit.Tumor
[16] RasP1Hit.Tumor RasP.Tumor     RasP.Tumor     RasPNot.Tumor  RasP1Hit.Tumor
[21] RasP1Hit.Tumor RasP1Hit.Tumor RasPNot.Tumor  RasP.Tumor     RasPNot.Tumor
[26] RasPNot.Tumor  RasPNot.Tumor  RasPNot.Tumor  RasPNot.Tumor  RasPNot.Tumor
[31] RasPNot.Tumor  RasPNot.Tumor  RasPNot.Tumor  RasPNot.Tumor  RasPNot.Tumor
[36] RasPNot.Tumor  RasPNot.Tumor  RasPNot.Tumor  RasPNot.Tumor  RasP.Tumor
[41] RasP.Tumor     RasP.Tumor     RasP.Tumor     RasPNot.Tumor  RasPNot.Tumor
Levels: RasP1Hit.Tumor RasP.Tumor RasPNot.Tumor


> design<-model.matrix(~RasTum);
> colnames(design)<-sub("RasTum","",colnames(design));
> colnames(design)[1] <- "Intercept"
> rownames(design)<-targets$Sample
> show(design[1:5,])
          Intercept RasP.Tumor RasPNot.Tumor
T4625_01A         1          1             0
T4627_01A         1          1             0
T4632_01A         1          1             0
T2655_01A         1          1             0
T2657_01A         1          1             0
...

> con.matrix<-makeContrasts(
+    RasP.Tumor_RasPNot.Tumor=RasP.Tumor-RasPNot.Tumor,RasPRasP1Hit.Tumor_RasPNot.Tumor=RasP.Tumor-RasPNot.Tumor,
+    RasP1Hit.Tumor_RasPNot.Tumor=-RasPNot.Tumor,
+ levels=design)


> show(con.matrix);
               Contrasts
Levels          RasP.Tumor_RasPNot.Tumor RasPRasP1Hit.Tumor_RasPNot.Tumor
  Intercept                            0                                0
  RasP.Tumor                           1                                1
  RasPNot.Tumor                       -1                               -1
               Contrasts
Levels          RasP1Hit.Tumor_RasPNot.Tumor
  Intercept                                0
  RasP.Tumor                               0
  RasPNot.Tumor                           -1


 

My first question is: as shown in show(RasTum), there are three subtype of tumors I want to compare between:

Levels: RasP1Hit.Tumor, RasP.Tumor and RasPNot.Tumor

shown in show(design[1:5,]), everything is relative to RasP1Hit.Tumor. So in my con.matrix<-makeContrasts command:

RasP1Hit.Tumor_RasPNot.Tumor=-RasPNot.Tumor

RasP.Tumor_RasPNot.Tumor=RasP.Tumor-RasPNot.Tumor

 

in particular, RasP1Hit.Tumor_RasPNot.Tumor=-RasPNot.Tumor since everything is relative to RasP1Hit.Tumor. 

However, for contrast: RasPRasP1Hit.Tumor_RasPNot.Tumor, I intended to pool both RasP class and RasP1Hit class together to compare to RasPNot class,  I end up with: 

RasPRasP1Hit.Tumor_RasPNot.Tumor=RasP.Tumor-RasPNot.Tumor

the right side is the same as 

RasP.Tumor_RasPNot.Tumor=RasP.Tumor-RasPNot.Tumor

because everything is relative to RasP1Hit.Tumor. So what is the right way to set up contrast RasPRasP1Hit.Tumor_RasPNot.Tumor in makeContrasts command? kind of confused here?

 

Because of the setting, the results are as below:

> show(summary(de <- decideTests(fit)));
   RasP.Tumor_RasPNot.Tumor RasPRasP1Hit.Tumor_RasPNot.Tumor
-1                        0                                0
0                     17764                            17764
1                         0                                0
   RasP1Hit.Tumor_RasPNot.Tumor
-1                            0
0                         17764
1                             0

 

 

My 2nd question is:  I did a similar design with including all normal samples, since at the same model I want to get tumor vs normal contrasts as well.

> con.matrix<-makeContrasts(
+    RasP.Tumor_RasP.Normal = RasP.Tumor-RasP.Normal,RasPRasP1Hit.Tumor_RasPRasP1Hit.Normal = RasP.Tumor-RasP.Normal-RasP1Hit.Normal,
+    RasPNot.Tumor_RasPNot.Normal=RasPNot.Tumor-RasPNot.Normal,
+    RasP.Tumor_RasPNot.Tumor=RasP.Tumor-RasPNot.Tumor,RasPRasP1Hit.Tumor_RasPNot.Tumor=RasP.Tumor-RasPNot.Tumor,
+    RasP.Normal_RasPNot.Normal=RasP.Normal-RasPNot.Normal,RasPRasP1Hit.Normal_RasPNot.Normal=RasP1Hit.Normal+RasP.Normal-RasPNot.Normal,
+    Tumor_Normal=RasP.Tumor+ RasPNot.Tumor+-RasP.Normal-RasPNot.Normal,
+    Tumor_Normal_2=RasP.Tumor+ RasPNot.Tumor-RasP.Normal-RasPNot.Normal-RasP1Hit.Normal,
+ levels=design)

 

here everything is relative to RasP1Hit.Tumor. and the result as below: 

  
> show(summary(de <- decideTests(fit)));
   RasP.Tumor_RasP.Normal RasPRasP1Hit.Tumor_RasPRasP1Hit.Normal
-1                   3947                                   4522
0                    9912                                   8830
1                    4274                                   4781
   RasPNot.Tumor_RasPNot.Normal RasP.Tumor_RasPNot.Tumor
-1                         4809                       24
0                          8461                    18102
1                          4863                        7
   RasPRasP1Hit.Tumor_RasPNot.Tumor RasP.Normal_RasPNot.Normal
-1                               24                        608
0                             18102                      17158
1                                 7                        367
   RasPRasP1Hit.Normal_RasPNot.Normal Tumor_Normal Tumor_Normal_2
-1                               2001         5611           5580
0                               13934         6728           6651
1                                2198         5794           5902

 

RasP.Tumor_RasPNot.Tumor contrast has small number of DEGs (24+7) at FDR 5% compared to using tumor data alone for the model shown earlier above, this might be cuased by filtering as well and I am fine with that small difference. What is bothering me is RasP.Normal_RasPNot.Normal contrast has about 608+367=~1k DEGs, which is way more than RasP.Tumor_RasPNot.Tumor contrast. I did plot plotMDS and see those normals are mixed together (I did see two subgroups seems having batch effects, but for each batch (if there is), there are RasP.Normal and RasPNot.Normal mixed there), all seem not supposed to having 1k DEGs either. Plus, the tumor contrast RasP.Tumor_RasPNot.Tumor has such small set of DEGs. Biologically, it might be true (population, ethics group, batch effect etc), but I just want to make sure what I did is reasonable. 

 

Any advice would be highly appreciated!

Thanks a lot for your help!

 

Best

Ming   

 

 

 

  

 

 

 

 

  



 

 

 

 

> Date: Wed, 19 Feb 2014 17:32:06 +1100
> From: smyth@wehi.EDU.AU
> To: yi02@hotmail.com
> CC: bioconductor@r-project.org
> Subject: Re: Problem in limma-voom
> 
> 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 inte...{{dropped:9}}

