Entering edit mode
Milena Gongora
▴
80
@milena-gongora-2035
Last seen 10.2 years ago
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