[BioC] Problem in limma-voom

Yunshun Chen yuchen at wehi.EDU.AU
Wed Feb 26 01:19:56 CET 2014


Dear Ming,

Regarding to your first question, a much easier way is to create your design
matrix without the intercept.
> design <- model.matrix(~0+RasTum)

Then if you intend to pool both RasP class and RasP1Hit class together to
compare to RasPNot class, you can use the following contrast:
> con.matrix <- makeContrasts( RasP_RasPNot=0.5*(RasP1Hit.Tumor+
RasP.Tumor)-RasPNot.Tumor, levels=design )

Similarly, for your second question, you can avoid confusion by creating
your design matrix without the intercept.
Suppose the six columns of your design matrix correspond to: Ras(Tum),
Ras1H(Tum), RasNot(Tum),   Ras(Nor), Ras1H(Nor), RasNot(Nor).
Then you can make contrasts for any comparisons you are interested. 
For example:
> con.matrix <- cbind( Tum_Nor=1/3*c(1,1,1,-1,-1,-1) ,
Ras1H_RasNot=0.5*c(0,1,-1,0,1,-1), RasAllNor-RasNotNor=c(0,0,0,0.5,0.5,-1) )

Hope that helps.

Yunshun Chen


-------------------------------------------------------------------

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_Ra
sPNot.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   


______________________________________________________________________
The information in this email is confidential and intend...{{dropped:4}}



More information about the Bioconductor mailing list