Error Message in nbinomLRT in DESeq2

[ToddB]

Hello all. I am working on some code for my lab and I have run into an error message I can seem to get around. My goal is to test each variable in the model given all other variables are present. Here is the code and the error that are tripping me up.

DESrun=nbinomLRT(DESrun, useOptim=FALSE, maxit = 500 reduced=as.formula(paste("~",paste(Vars[-i],collapse="+"))))

Error in optim(betaRow, objectiveFn, method = "L-BFGS-B", lower = -large, :
L-BFGS-B needs finite values of 'fn'


Vars contain my variables of interest, for example Vars = c("Age","Gender","Disorder"). And DESrun has my full model, dispersion estimates, size factor estimates, and data. I am using DESeq2 version 1.0.18.

If I am still getting this error with the useOptim set to false and maxit set to 500 anyone have a suggestion?

Thanks,
Todd

[Yohann]

Hi,

I have a pretty similar problem :

Error in optim(betaRow, objectiveFn, method = "L-BFGS-B", lower = -large, :
non-finite finite-difference value [3]

What is weird is that i'm pretty sure it worked with an older version of DESeq2 as i'm re-running an old script.

Did you find a solution ?

[dpryan]

Can you post the relevant portion of your script (including the design data.frame)? I suspect that something is off in your fit function.

[ToddB]

No solution yet. But, I have been in contact with Mike Love and he is tracking it down. He sent me the following to see which rows are causing an issue.


if you do:

debug(DESeq2:::fitNbinomGLMs)

this will show which rows are crashing.

nbinomLRT first runs fitNbinomGLMs on the full formula which is not
throwing the error. So you can hit enter through the function once.
The second time it is fitting the reduced formula. On the second time,
after this code chunk, you can identify the rows causing the problem:

# switch based on whether we should also use optim
# on rows which did not converge
if (useOptim) {
rowsForOptim <- which(!betaConv | !rowStable | !rowVarPositive)
} else {
rowsForOptim <- which(!rowStable | !rowVarPositive)
}

The rowsForOptim is the subset of rows which are causing problems.

My current plan is to use this to find the rows causing a problem and remove them, then rerun DESeq2.

[Michael Love]

I now am committing DESeq2 version 1.0.19 to Bioc which uses Nelder-Mead for this back-up optimization rather than L-BFGS-B. This appears to avoid this halting error.

Furthermore, for Todd's data, it is necessary, for now, to also standardize the numeric predictors, e.g.:

# scale the numeric predictors
for (var in c("age","weight")) {
colData(dds)[[var]] <- as.numeric(scale(colData(dds)[[var]]))
}

Then in the end, the log2 fold changes and lfcSE can be multiplied by sd(age) to obtain what would have been the original log2 fold change expected per year for the age variable. (Or you can leave them as they are and they are interpretable as log2 fold changes expected per standard deviation)

I am working on improving the fitting algorithm in the devel branch, so that this manual scaling will not be necessary.

[Yohann]

Hi,

I was using a continuous variable that was the number of bacteria as a covariate.
This value scales from ~ 1e4 to 1e7
If i use this scale, i have the error message : "Error in optim(betaRow, objectiveFn, method = "L-BFGS-B", lower = -large, : non-finite finite-difference value [3]"
If i instead divide my values by 1e6 (it's now millions of bacteria), it works.

[Simon Anders]

Are you sure that you don't want to rather regress on the logarithm of your predictor, given its huge dynamic range?

[Yohann]

Originally Posted by Simon Anders

Are you sure that you don't want to rather regress on the logarithm of your predictor, given its huge dynamic range?

That's a good idea, thanks!
I'll check that !

[Michael Love]

I have implemented a change to the fitting algorithm in DESeq2 version >= 1.1.32 which should eliminate the need to standardize predictors as described in my post above.

[Yohann]

Hi Michael,

Here is a part of my code :

colData <- data.frame(
	AFcomp = sub_individuals$AFcomp,
	bacteria_count = sub_individuals$Contact,
	batch = as.factor(paste0("F", sub_individuals$Flowcell)))

dds <- DESeqDataSetFromMatrix(
	countData = sub_reads,
	colData = colData,
	design = ~ AFcomp + bacteria_count + batch)

dds <- estimateSizeFactors(dds)
dds <- estimateDispersions(dds, maxit = 500, quiet = TRUE)
dds <- nbinomLRT(dds, full = design(dds), maxit = 500, reduced = ~ bacteria_count + batch)

Where :

• AFcomp is a continuous variable scaling from 0 to 1.
• bacteria_count is a continuous variable from 1e6 to 8e6.
• batch is a discrete variable (batch name).

if i use the bacteria_count like that, i get this error :

Error in solve.default(xtwx + ridge) :
  system is computationally singular: reciprocal condition number = 9.67189e-20

If i divide the bacteria_count by 1e6 (now looking at millions of bacteria), it works.

I'm using the 1.1.33 version from the dev branch of bioconductor.

Do you have an idea?

Thanks!

[Michael Love]

hi Yohann,

I see the problem and will look into it.

[Michael Love]

hi Yohann,

The problem was coming from some code which goes over rows which did not converge in the standard GLM fitting steps using R's optimization function. So I had not updated this code to account for very large columns in the design matrix. I have now updated this bit of code in version 1.1.35 and added a unit test to confirm that the coefficients and standard errors are identical with or without this extra optimization loop.

If you could check if this solves your error that would be great.

[Yohann]

hi Michael, thanks for your answer.

It looks like i can only retrieve the 1.1.34 version on http://bioconductor.org/packages/dev...ml/DESeq2.html.

Is there an other way to get the 1.1.35 ?

Thanks,

[Michael Love]

It should show up on the Bioc site around this time tomorrow.