[BioC] Problem with function limmaCtData in HTqPCR package: "leading minor of order 2 is not positive definite"
Heidi Dvinge
heidi at ebi.ac.uk
Fri Jul 23 11:11:47 CEST 2010
Hello Kevin,
On 21 Jul 2010, at 19:50, Bass, Kevin wrote:
> Hi,
>
> I am having a problem with using the function limmaCtData on a qPCRset
> object created with the package HTqPCR. When I try to execute
> limmaCtData, I get the following error:
>
> "Error in chol.default(V) :
> the leading minor of order 2 is not positive definite"
>
as your traceback() shows, in the first step the error comes from
lmFit from the limma package. As I recall, it means that one of
internal design matrices has become singular. I'm afraid I don't
know exactly why this is happening, however it can be caused by
trying to do to much with too few observations/replicates. Is it
possible to use a smaller design matrix? Looking at your design
matrix it would appear that you have no replicates of either of the 5
treatments you list there. Based on your description of the
experiment, I'm not really sure whether this is the case or not?
By the way, it looks like you have quite a complex plate/sample
combination design compared to a standard qPCR analysis - I can see
we you end up with an object called "raw_monster2" after all the
different rbind and cbind ;)
Cheers
\Heidi
> Below, I will describe the experimental design and the steps taken to
> create my qPCRset object. Then I will paste the commands used, and
> their results, in the steps leading up to running the limmaCtData
> function on my qPCRset object.
>
> We have 21 96-well plates. Each plate contains 5 experimental groups
> and 4 genes--2 target genes, and 2 endogenous controls. Each
> experimental group sampled all 4 genes, and there were 3 biological
> replicates per sample, for a total of 12 wells per experimental group.
>
> Every 7 plates among the total 21 plates constitutes a "set" of
> plates: they each contain the same 14 target genes. This means that
> each gene, in each experimental condition, has 3 samples among the 21
> plates--one sample per experimental condition for each 7-plate set.
>
> The goal is to compare the Ct values for each gene in each
> experimental group, to the Ct values for the same gene in every other
> experimental group.
>
> Using rbind (HTqPCR), I collated 7 of the data files into one file,
> so that all 14 genes could be analyzed simultaneously, at least among
> a single set of plates--once I had figured that part out, I had
> planned on combining the 3 sets.
>
> To give a clear idea what my data looks like--and how it was
> implemented in my qPCRset object--this is the Slot "history" and Slot
> "exprs" of my combined qPCRset object (with the data removed):
>
> Slot "exprs":
> 01_veh+FA 02_low+FA 03_mid+FA 04_high+FA
> 05_no_treatment
> PGES
> PGES
> PGES
> c-Fos
> c-Fos
> c-Fos
> SPP1
> SPP1
> SPP1
> CD200
> CD200
> CD200
> COX-1
> COX-1
> COX-1
> COX-2
> COX-2
> COX-2
> OX-42
> OX-42
> OX-42
> iBA-1
> iBA-1
> iBA-1
> IL-2
> IL-2
> IL-2
> IL-4
> IL-4
> IL-4
> IL-6
> IL-6
> IL-6
> IL-8
> IL-8
> IL-8
> IL-10
> IL-10
> IL-10
> CD4
> CD4
> CD4
>
> Slot "history":
> history
> 1 raw8: readCtData(files = "NS398_08b.txt", path = barrPath,
> n.features = 12,
> 2 flag = NULL, feature = 5, type = 7, position = 2, Ct = 6,
> 3 header = TRUE, n.data = 5)
> 4 raw9: readCtData(files = "NS398_09b.txt", path = barrPath,
> n.features = 12,
> 5 flag = NULL, feature = 5, type = 7, position = 2, Ct = 6,
> 6 header = TRUE, n.data = 5)
> 7 raw10: readCtData(files = "NS398_10b.txt", path = barrPath,
> n.features = 12,
> 8 flag = NULL, feature = 5, type = 7, position = 2, Ct = 6,
> 9 header = TRUE, n.data = 5)
> 10 raw11: readCtData(files = "NS398_11b.txt", path = barrPath,
> n.features = 12,
> 11 flag = NULL, feature = 5, type = 7, position = 2, Ct = 6,
> 12 header = TRUE, n.data = 5)
> 13 raw12: readCtData(files = "NS398_12b.txt", path = barrPath,
> n.features = 12,
> 14 flag = NULL, feature = 5, type = 7, position = 2, Ct = 6,
> 15 header = TRUE, n.data = 5)
> 16 raw13: readCtData(files = "NS398_13b.txt", path = barrPath,
> n.features = 12,
> 17 flag = NULL, feature = 5, type = 7, position = 2, Ct = 6,
> 18 header = TRUE, n.data = 5)
> 19 raw14: readCtData(files = "NS398_14b.txt", path = barrPath,
> n.features = 12,
> 20 flag = NULL, feature = 5, type = 7, position = 2, Ct = 6,
> 21 header = TRUE, n.data = 5)
> 22 rbind(deparse.level, ..1, ..2, ..3, ..4, ..5, ..6, ..7)
> 23 normalizeCtData(q = raw_monster2, norm = "deltaCt",
> deltaCt.genes = "GAPDH")
> 24 filterCtDataNew(q = d.raw2, remove.type = "Endogenous
> Control")
> 25 setCategory(q = fd.raw2, Ct.max = 100, Ct.min = 0,
> quantile = 0.9,
>
> So, then I prepared the matrix for analysis with limma:
>
>> design<-model.matrix(~0+sampleNames(test.d.raw2))
> Warning message:
> In model.matrix.default(~0 + sampleNames(test.d.raw2)) :
> variable 'sampleNames(test.d.raw2)' converted to a factor
>> colnames(design)<-c("VehFA","LowFA","MidFA","HighFA","NoTreat")
>> print(design)
> VehFA LowFA MidFA HighFA NoTreat
> 1 1 0 0 0 0
> 2 0 1 0 0 0
> 3 0 0 1 0 0
> 4 0 0 0 1 0
> 5 0 0 0 0 1
> attr(,"assign")
> [1] 1 1 1 1 1
> attr(,"contrasts")
> attr(,"contrasts")$`sampleNames(test.d.raw2)`
> [1] "contr.treatment"
>> contrasts<-makeContrasts(VehFA-LowFA, VehFA-MidFA, VehFA-HighFA,
> + VehFA-NoTreat, LowFA-MidFA, LowFA-HighFA, LowFA-NoTreat,
> + MidFA-HighFA, MidFA-NoTreat,HighFA-NoTreat, levels=design)
>> colnames(contrasts)<-c("V-L", "V-M", "V-H", "V-NT", "L-M", "L-H",
> + "L-NT", "M-H", "M-NT", "H-NT")
>> print(contrasts)
> Contrasts
> Levels V-L V-M V-H V-NT L-M L-H L-NT M-H M-NT H-NT
> VehFA 1 1 1 1 0 0 0 0 0 0
> LowFA -1 0 0 0 1 1 1 0 0 0
> MidFA 0 -1 0 0 -1 0 0 1 1 0
> HighFA 0 0 -1 0 0 -1 0 -1 0 1
> NoTreat 0 0 0 -1 0 0 -1 0 -1 -1
>> test.d.raw2b<-test.d.raw2[order(featureNames(test.d.raw2)), ]
> ======================================================================
> =
>> qDE.limma <- limmaCtData(test.d.raw2b,design=design,
> + contrasts=contrasts,ndups=3,spacing=1)
> Error in chol.default(V) :
> the leading minor of order 2 is not positive definite
> In addition: Warning message:
> In sqrt(dfitted.values) : NaNs produced
>> traceback()
> 6: .Call("La_chol", as.matrix(x), PACKAGE = "base")
> 5: chol.default(V)
> 4: chol(V)
> 3: gls.series(y$exprs, design = design, ndups = ndups,
> spacing = spacing, block = block, correlation = correlation,
> weights = weights, ...)
> 2: lmFit(data, design = design, ndups = ndups, spacing = spacing,
> correlation = dup.cor$consensus, ...)
> 1: limmaCtData(test.d.raw2b, design = design, contrasts = contrasts,
> ndups = 3, spacing = 1)
>
> Any ideas on why I am getting this error and what I might do to avoid
> it? If there is any other information needed, please let me know.
>
> Thanks,
> Kevin
> bassk1 at email.chop.edu
>
>
>
> =====
>
> Kevin Bass, Research Technician
> Barr Lab
> Children's Hospital of Philadelphia
> Abramson Research Center
> 3615 Civic Center Blvd, Suite 714
> Philadelphia PA 19104-4399
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives: http://news.gmane.org/
> gmane.science.biology.informatics.conductor
More information about the Bioconductor
mailing list