A running limma analysis of a dataset, and I got an error in the final step, eBayes(fit2)
This is my code
i <- pData(eset)$Sample_title
Part <- ifelse(pData(eset)$Part == "Tn", "nT", "T")
Grade <- ifelse(pData(eset)$Grade == "Low", "nH", "H")
sample_byGrade <- factor(paste(Grade, Part, sep="."))
design_Grade <- model.matrix(~ 0 + sample_byGrade)
colnames(design_Grade) <- levels(sample_byGrade)
corfit <- duplicateCorrelation(eset,design_Grade,block = i)
contrast_matrix_Grade<- makeContrasts(HighVsLowinT = H.T-nH.T, HighVsLowinnT = H.nT-H.nT,
TVsnTinHigh = H.T-H.nT, TVsnTinLow = nH.T-nH.nT,
levels = design_Grade)
fit <- lmFit(eset,design = design_Grade,block = i, correlation = corfit$consensus)
fit2<- contrasts.fit(fit,contrast_matrix_Grade)
fit2<- eBayes(fit2)
The error message:
Error in eigen(cor.matrix, symmetric = TRUE) :
infinite or missing values in 'x'
Besides: Warning messages:
1: In ebayes(fit = fit, proportion = proportion, stdev.coef.lim = stdev.coef.lim, :
Estimation of var.prior failed - set to default value
2: In cov2cor(object$cov.coefficients) :
diag(.) had 0 or NA entries; non-finite result is doubtful
I see previous reply and some said that it is possible due to NA in the data. I do anyNA(exprs(eset)) and it is FALSE, and I did do filtering before running LIMMMA so I believe my data got no NA values.
But some of my column in pData(eset) did have NA value, but these columns are not used during this analysis. Where is my problem?
Thanks everyone.
Please show us pData(eset).
Is it possible that Sample_title takes a different value for each sample? In that case, setting block=i wouldn't make any sense and would lead to an error such as you have.
Have you checked that corfit$consensus is between -1 and 1 (not equal to -1 or 1)?
The pData is here. Not sure if this work out.
https://www.dropbox.com/s/35hqpm11umynx19/eset.txt?dl=0
Sample_title, i, is the subject name
i
[1] "E005" "E005" "E012" "E012" "E021" "E021" "E028" "E028" "E033" "E033" "E041" "E041" "E042" "E042"
[15] "E045" "E045" "E048" "E048" "E052" "E052" "E058" "E058" "E064" "E064" "E073" "E073" "E082" "E082"
[29] "E086" "E086" "E088" "E088" "E094" "E094" "E109" "E109" "E119" "E119" "E120" "E120" "E123" "E123"
[43] "E130" "E130" "E133" "E133" "E136" "E136" "E138" "E138" "E139" "E139" "E152" "E152" "E162" "E162"
[57] "E165" "E165" "E178" "E178" "E179" "E179" "E181" "E181" "E182" "E182" "E186" "E186" "E203" "E203"
[71] "E208" "E208" "E209" "E209" "E218" "E218" "E220" "E220" "E229" "E229" "E234" "E234" "E245" "E245"
[85] "E246" "E246" "E249" "E249" "E250" "E250"
> corfit$consensus
[1] 0.05972559
The pData file on drop box doesn't agree with the output you give for Sample_title, i.
In the pData file, the Sample_title values are "E005-Tc8", "E005-Tn", "E012-Tn" etc. In the pData file, every sample has a different title.
Ok.
I modified i <- pData(eset)$Sample_title into this:
i <- sub("-.*", "", pData(eset)$Sample_title), then the data become as what i show up there.
Sorry for wasting your time, but I make a stupid mistake.
HighVsLowinnT = H.nT-H.nT should be HighVsLowinnT = H.nT-nH.nT.
Thank you for helping me.