This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Limma analysis - decomposition problem

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

duplicatecorrelation lmfit cholesky limma

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

Thanks a lot!

1 answer

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.

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).

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

I have only 2 conditions. Either case or control.

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

I get only a matrix of NAs!

Ah, no errors when running duplicateCorrelation?

Log in to answer this question.