Question: Limma analysis - decomposition problem
0
iside10 wrote:

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)

modified 4.4 years ago by andrew.j.skelton735.9k • written 4.4 years ago by iside10

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

Thanks a lot!

1
andrew.j.skelton735.9k wrote:

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,

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)

  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!