Limma analysis - decomposition problem
1
0
Entering edit mode
6.7 years ago
iside ▴ 20

Dear all,

I am running a differential expression analysis with limma. My dataset includes unrelated individuals as well as twins, and I am aware that limma is not designed to deal with familiar data, but, as explained in the users' guide I may use the family as a blocking variable (in this way limma will estimate a within-family correlation from the expression data and will use it for the differential analysis).

Unfortunately, when running lmFit with the correlation factor I get this error:

Error in chol.default(V) : 
  the leading minor of order 2 is not positive definite
  
I understand that this depends on my data structure, and clearly I do not have a positive definite matrix, therefore it doesn't respect the assumptions of cholesky decomposition. But do you have any suggestion on how to overcome the problem?  
 
This is the code I am using:

formula <- ~ condition
design <- model.matrix(formula, data=info)
corfit <- duplicateCorrelation(eset, design, block=info$famID)
fit <- lmFit(eset, design, block=info$famID, correlation=corfit.corr$consensus)

Thanks in advance

Limma lmFit Cholesky duplicateCorrelation • 2.0k views
ADD COMMENT
0
Entering edit mode

I have the same problem with lmFit! Have you found any solution yet? I have log2 transformed data as well.

Thanks a lot!

ADD REPLY
1
Entering edit mode
6.7 years ago

Is your normalised eset in log2 space? It'd be useful to see the output of:

range(eset)

edit: It'd also be useful to see a pheno table so we can gauge the experimental design.

ADD COMMENT
0
Entering edit mode

Hi Andrew,

thanks for your reply. 

The data had already been normalised, and corrected for batch effect with the package Combat. The output you asked for is the following:

range(expression)

[1]  6.109361 16.553681.

My design matrix is composed by two columns: Intercept (136 rows set to 1) and  conditioncase, where control is set to 0 (68 rows) and case is set to 1 (68 rows).

 

ADD REPLY
0
Entering edit mode

Thanks, how many levels are there to your condition variable? i.e. How many conditions do you have?

ADD REPLY
0
Entering edit mode

I have only 2 conditions. Either case or control. 

ADD REPLY
0
Entering edit mode

ok, that explains your intercept model. What's the output of corfit?

ADD REPLY
0
Entering edit mode

I get only a matrix of NAs!

ADD REPLY
0
Entering edit mode

Ah, no errors when running duplicateCorrelation?

ADD REPLY

Login before adding your answer.

Traffic: 2063 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6