Limma analysis - decomposition problem

Question

[iside]

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

[elenichri]

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

Thanks a lot!

Answer 1

[andrew.j.skelton73]

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.

[iside]

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

[andrew.j.skelton73]

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

[iside]

I have only 2 conditions. Either case or control.

[andrew.j.skelton73]

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

[iside]

I get only a matrix of NAs!

[andrew.j.skelton73]

Ah, no errors when running duplicateCorrelation?