infinite solver loop when growing on computed draft medium

[mbiomeai]

Hello @jotech and @Waschina, first many thanks for the continued development of gapseq and the performance improvements that have been implemented. I have been following the batch_alignmet_tests branch and running the April 9 gapseq's "doall" on our collection of gut genomes. This has generated a usable SBML for the vast majority (99.8%), but several handfull got into extreme running time and had actually entered into what looks like an infinite loop of the GLPK solver or the R code calling it:

running Rscript /bifx/conda/envs/cuda-gapseq2/gapseq-2a/src/gf.suite.R -m /bifx/db/gems/gutmost/gm2308-gapseq2a/8SCYL354/./8SCYL354-draft.RDS -n /bifx/db/gems/gutmost/gm2308-gapseq2a/8SCYL354/./8SCYL354-medium.csv -b 100 -f /bifx/db/gems/gutmost/gm2308-gapseq2a/8SCYL354/.
LP solver: glpk
Loading model file /bifx/db/gems/gutmost/gm2308-gapseq2a/8SCYL354/./8SCYL354-draft.RDS
using media file /bifx/db/gems/gutmost/gm2308-gapseq2a/8SCYL354/./8SCYL354-medium.csv


1. Initial gapfilling: Make model grow on given media using all reactions
Warning: numerical instability (primal simplex, phase II)
Warning: numerical instability (primal simplex, phase II)
Warning: numerical instability (primal simplex, phase II)
Warning: numerical instability (primal simplex, phase II)
Warning: numerical instability (primal simplex, phase II)
...

So I was wondering if this is solver specific and what can be done about it to solve the issue. I have not tested the very latest development version. I have joined an archive with an example genome failing the initial gapfilling. Below is my gapseq test output.

8SCYL354.tgz

gapseq version: 1.4.0 95c8e31d
linux-gnu
#1 SMP PREEMPT_DYNAMIC Mon Feb 10 05:23:56 EST 2025


#######################
#Checking dependencies#
#######################
ldconfig (GNU libc) 2.34
libsbml.so.5 -> libsbml.so.5.20.4
libglpk.so.40 -> libglpk.so.40.3.1
GNU Awk 5.3.1, API 4.0, PMA Avon 8-g1, (GNU MPFR 4.2.1, GNU MP 6.3.0)
sed (GNU sed) 4.8
grep (GNU grep) 3.11
This is perl 5, version 32, subversion 1 (v5.32.1) built for x86_64-linux-thread-multi
tblastn: 2.16.0+
exonerate from exonerate version 2.4.0
bedtools v2.31.1
barrnap 0.9 - rapid ribosomal RNA prediction
R version 4.3.3 (2024-02-29) -- "Angel Food Cake"
git version 2.47.0
GNU parallel 20240922
HMMER 3.4 (Aug 2023); http://hmmer.org/
bc 1.07.1

Missing dependencies: 0


#####################
#Checking R packages#
#####################
data.table 1.15.4
stringr 1.5.1
cobrar 0.2.0
getopt 1.20.4
R.utils 2.12.3
stringi 1.8.4
BiocManager 1.30.25
Biostrings 2.70.1
jsonlite 1.8.9
httr 1.4.7

Missing R packages:  0


##############################
#Checking basic functionality#
##############################
Optimization test: OK
Building full model: OK
Blast test: OK

Passed tests: 3/3

[Waschina]

Hi @MRMHmdeleeuw

Thank you for reporting the issue and providing the model files. I can reproduce the error when using glpk as solver.

The issue is most likely due to a bug in GLPK, as this infinite loop of the warning is already reported in 2018 (!) here and here. The bug seems to persist to this day.

I will check if we can adjust something related to the constraint matrix scaling before optimization in GLPK.

[mbiomeai]

Ok thank you @Waschina for the feedback. If it cannot be solved I'd appreciate to know with which solvers available to you the model does solve.

[Waschina]

I usually use IBM's CPLEX as solver. This requires the installation of the R-package cobrarCPLEX, and of course IBM ILOG CPLEX itself.

But I think I solved the issue with GLPK now, too; by a slight modification in cobrar. I'll run some tests and if they work, I push a new version of cobrar.

[Waschina]

The latest cobrar development version now has the bug fix. You can update your cobrar version by running e.g. remotes::install_github("Waschina/cobrar") .

[mbiomeai]

Hello @Waschina, many thanks for the fix. I updated gapseq to the latest development branch commit 265f6ce and cobrar as indicated above. This made it possible to complete 41/42 previously failing GEMs, the last one, joined for reference, is still failing with the glpk infinite loop.

LP solver: glpk
Loading model file /mnt/bifx1/db/gems/gutmost/gm2308-gapseq2504/1KFYEC4J/./1KFYEC4J-draft.RDS
using media file /mnt/bifx1/db/gems/gutmost/gm2308-gapseq2504/1KFYEC4J/./1KFYEC4J-medium.csv


1. Initial gapfilling: Make model grow on given media using all reactions
Warning: numerical instability (primal simplex, phase II)
Warning: numerical instability (primal simplex, phase II)
...

1KFYEC4J.tgz

[Waschina]

Hi @MRMHmdeleeuw ,

I had a closer look at this again. The special case with the 1KFYEC4J.tgz model is that the medium prediction added N₂ and no Ammonium to the gapfill medium, because gapseq predicted a nitrogenase. That alone is not an issue, but in combination with the result that no energy source (like glucose or H₂ was detected for the gapfill medium, the solver has problems finding a solution for the LP... For some reason, GLPK enters an infinite loop in the specific example...

For gapseq, we could try to detect cases like this, where the predicted medium likely does not contain sufficient resources to enable a gapfilling+flux distribution solution; and maybe add glucose as a fallback... I'll discuss this with @jotech .

[mbiomeai]

Hello @Waschina, many thanks for the detailed follow-up. The GEM was by no means crucial to us and we have been able to make very good use of our 27k+ Gapseq GEMS. Good to see the explanation and if the case can be used to make Gapseq even more robust, all the better.