Lmfit Method Of Limma Package Throwing Error While Calculating P Values

Question

[Curious Mind]

Hi,

I am trying to find the p-values for golub data. I have loaded the data as follows:

library(multtest)
library(annotate)
library(limma)

dat<-golub[,]
dat<-as.data.frame(dat)

The first 27 columns are controls and the remaining 11 are tests. I want to calculate the p-values using the eBayes(0 method in the limma package. The following is my code:

design<-as.data.frame(cbind(dat[,1:27], dat[,28:38]))
fit<-lmFit(dat,design)

However, I am getting the following error:

Error in lm.fit(design, t(M)) : incompatible dimensions

When I checked, the dim of design and dat, they are of same dimension as shown below:

dim(design)
[1] 3051   38
dim(dat)
[1] 3051   38

Thanks

Answer 1

[Sean Davis]

Your design matrix should not have the same dimensions as your data. I find it useful to create both design and contrast matrices and to create the design matrix using a factor (or a character vector, as below, that will be converted to factor).

design = model.matrix(~0+c(rep('control',27),rep('test',11)))
colnames(design) = c('control','test')
cm = makeContrasts(test-control,levels=design)
fit = lmFit(dat,design)
fit2 = contrasts.fit(fit,cm)
fit3 = eBayes(fit2)

[Curious Mind]

Hello Sean, Thanks for the reply. Now the error is gone. When I tried to get the p-Values using .p.Values<-fit3$p.value[,2] , I am getting the following error: Error in fit3$p.value[, 2] : subscript out of bounds Thanks

[Curious Mind]

I have found the p values as follows. p.Values<-fit3$p.value Thanks

[Sean Davis]

You may want to use topTable() for this type of thing.

Answer 2

[ahmedc3.ri]

You will need to define the whole data set first then you define your dimensions:

For example I have a data of 17 samples 2 groups.

7 sample in the first group and 10 in the second.

I can say

des=cbind(c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1),c(1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2))

so cbind is the total samples and the c will specify how many for each group.

Then you can fit your model

fit=lmFit(exprs.eset,design=des)
gene.ps=rownames(exprs.eset)