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
 

 

 

 

 
 		 	   		  
	[[alternative HTML version deleted]]

