[BioC] LIMMA - problem with block analysis

[Milena Gongora]

Dear LIMMA users,

Hello. I am analysing single colour microarrays in LIMMA and am facing 
some obstacles, which I am not sure have a solution at all. So wondering 
if someone can confirm this or illuminate me on why these errors.

I have a really bad experimental design. It is a simple pair-wise 
comparison between mutant and a control conditions. I have 7 samples, 3 
control samples (comming from independent individuals) and 4 mutant 
samples where 3 come from the same individual and 1 from its mother.

So in essence, I have a high level of correlation between the mutant 
samples, but independence in the control samples.

What I thought I could do to get the most realistic information out of 
this data, was to do a "blocked" analysis where I would look at the 
differences between mutant and control, but blocking the samples by 
patient first.
I tried doing this by following the example in the LIMMA guide for 
section 8.3 "paired samples" as below:

 > samples
  Sample Mutation Patient Related
1  Mu_1A       mu       1      R1
2  Mu_1B       mu       1      R1
3  Mu_1C       mu       1      R1
4    C_3     ctrl       3      R2
5    C_4     ctrl       4      R3
6    C_5     ctrl       5      R4
7   Mu_2       mu       2      R1

# extract the variables of interest as factors
 > mutation <- factor(samples$Mutation)
 > patient <- factor(samples$Patient)
 
 > design_c1 <- model.matrix(~patient+mutation)
 > fit <- lmFit(data_Qnorm_log2, design_c1)
Coefficients not estimable: mutationmu

My question here is: Why can the coefficient not be estimated?

So I tried a second approach. Instead of putting my blocking vairable in 
the design matrix. I tried fitting the model specifying the correlation 
using duplicateCorrelation() like in section 8.5 of LIMMA guide 
"Technical Replication". In this approach, the problem is that my 
correlation is negative and can't fit the model. See below:

 > block_p <- as.vector(patient)
 > corr_patient <- duplicateCorrelation(exprs(data_Qnorm_log2),design_1, 
block=block_p)
 > corr_patient$consensus.correlation
[1] -1
 > fit_corr_pat <- lmFit(data_Qnorm_log2, design_1, block=block_p, 
correlation=corr_patient$consensus)
Error in chol.default(V) :
  the leading minor of order 2 is not positive definite

Why is it producing a negative correlation? if I do cor() on the 
expression values, I get positive correlations as seen here

 >  cor(exprs(data_Qnorm_log2))
          Mu_1A     Mu_1B     Mu_1C       C_3       C_4       C_5      Mu_2
Mu_1A 1.0000000 0.9814710 0.9768791 0.9753479 0.9783894 0.9744542 0.9710319
Mu_1B 0.9814710 1.0000000 0.9785653 0.9765288 0.9783644 0.9755863 0.9706648
Mu_1C 0.9768791 0.9785653 1.0000000 0.9727115 0.9750562 0.9756139 0.9673867
C_3   0.9753479 0.9765288 0.9727115 1.0000000 0.9762357 0.9765884 0.9673692
C_4   0.9783894 0.9783644 0.9750562 0.9762357 1.0000000 0.9761812 0.9707479
C_5   0.9744542 0.9755863 0.9756139 0.9765884 0.9761812 1.0000000 0.9710843
Mu_2  0.9710319 0.9706648 0.9673867 0.9673692 0.9707479 0.9710843 1.0000000
 >

In summary is it possible to run this analysis or is there just too many 
constraints in this data? If possible, how could I do so?

I appreciate every insight and thanks very much for reading this long post!

Milena