From dd503a039a83f2fcdc30ee68261ecf52b87ff86d Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 30 Sep 2026 21:05:06 +0200 Subject: [PATCH 1/3] fix: gurobi without Parallel Computing Toolbox, report missing WSL (#753) * Fix gurobi solving without Parallel Computing Toolbox * Report missing WSL when HMMER check fails on Windows * Update documentation * Diagnose why the HMMER check fails on Windows * Update documentation * Restore and update documentation * Skip getCurrentTask when Parallel Computing Toolbox is not licensed --- doc/installation/checkInstallation.html | 193 +++++++----- doc/solver/optimizeProb.html | 381 ++++++++++++------------ installation/checkInstallation.m | 33 +- solver/optimizeProb.m | 11 +- 4 files changed, 347 insertions(+), 271 deletions(-) diff --git a/doc/installation/checkInstallation.html b/doc/installation/checkInstallation.html index 3ee042608..a36271d48 100644 --- a/doc/installation/checkInstallation.html +++ b/doc/installation/checkInstallation.html @@ -381,92 +381,121 @@

SOURCE CODE ^'hmmerTests.m')"); 0310 res=interpretResults(res); 0311 if res==false -0312 fprintf([' This is essential to run getKEGGModelFromHomology()\n'... -0313 ' when using a FASTA file as input\n']) -0314 end -0315 -0316 if developMode -0317 fprintf('\n=== Development binary executables ===\n'); -0318 fprintf('NOTE: These binaries are only required when using KEGG FTP dump files in getKEGGModelForOrganism\n'); -0319 -0320 fprintf(myStr(' > Checking CD-HIT',40)) -0321 [~,res]=evalc("runtests('cdhitTests.m');"); -0322 interpretResults(res); -0323 -0324 fprintf(myStr(' > Checking MAFFT',40)) -0325 [~,res]=evalc("runtests('mafftTests.m');"); -0326 interpretResults(res); -0327 end -0328 end -0329 -0330 fprintf('\n=== Compatibility ===\n'); -0331 fprintf(myStr(' > Checking function uniqueness',40)) -0332 checkFunctionUniqueness(); -0333 -0334 fprintf('\n*** checkInstallation complete ***\n'); -0335 end -0336 -0337 function res = interpretResults(results) -0338 if results.Failed==0 && results.Incomplete==0 -0339 fprintf('Pass\n'); -0340 res=true; -0341 else -0342 printOrange('Fail\n') -0343 fprintf(' Download/compile the binary and rerun checkInstallation\n'); -0344 res=false; -0345 end -0346 end -0347 -0348 function str = myStr(InputStr,len) -0349 str=InputStr; -0350 lenDiff = len - length(str); -0351 if lenDiff < 0 -0352 warning('String too long'); -0353 else -0354 str = [str blanks(lenDiff)]; -0355 end -0356 end -0357 -0358 function status = makeBinaryExecutable(ravenDir) -0359 % This function is required to run when RAVEN is downloaded as MATLAB -0360 % Add-On, in which case the file permissions are not correctly set -0361 if ispc -0362 status = 0; % No need to run on Windows -0363 return; +0312 fprintf([' HMMER is only required by getKEGGModelForOrganism() when\n'... +0313 ' using a FASTA file as input. If you do not use this function,\n'... +0314 ' this failure can be ignored.\n']) +0315 % On Windows, HMMER runs through Windows Subsystem for Linux +0316 if ispc +0317 [status,~]=system('where wsl'); +0318 if status~=0 +0319 printOrange([' MATLAB cannot find wsl.exe. On Windows, HMMER runs through WSL.\n'... +0320 ' Install it with "wsl --install" in a Command Prompt, then restart\n'... +0321 ' MATLAB (it inherits the system PATH at startup) and rerun checkInstallation\n']) +0322 else +0323 [status,~]=system('wsl -e true'); +0324 if status~=0 +0325 printOrange([' MATLAB finds wsl.exe, but cannot start a Linux distribution\n'... +0326 ' through it. Run "wsl -l -v" in a Command Prompt to check that a\n'... +0327 ' distribution is installed and set up, then restart MATLAB and\n'... +0328 ' rerun checkInstallation\n']) +0329 else +0330 try +0331 [status,~]=system(['wsl "' getWSLpath(fullfile(ravenDir,'software','hmmer','hmmsearch')) '" -h']); +0332 catch +0333 status=1; +0334 end +0335 if status~=0 +0336 printOrange([' WSL works, but the RAVEN HMMER binary cannot be run in it.\n'... +0337 ' Check that the files in software/hmmer are present and that\n'... +0338 ' the drive holding RAVEN is accessible from WSL\n']) +0339 end +0340 end +0341 end +0342 end +0343 end +0344 +0345 if developMode +0346 fprintf('\n=== Development binary executables ===\n'); +0347 fprintf('NOTE: These binaries are only required when using KEGG FTP dump files in getKEGGModelForOrganism\n'); +0348 +0349 fprintf(myStr(' > Checking CD-HIT',40)) +0350 [~,res]=evalc("runtests('cdhitTests.m');"); +0351 interpretResults(res); +0352 +0353 fprintf(myStr(' > Checking MAFFT',40)) +0354 [~,res]=evalc("runtests('mafftTests.m');"); +0355 interpretResults(res); +0356 end +0357 end +0358 +0359 fprintf('\n=== Compatibility ===\n'); +0360 fprintf(myStr(' > Checking function uniqueness',40)) +0361 checkFunctionUniqueness(); +0362 +0363 fprintf('\n*** checkInstallation complete ***\n'); 0364 end -0365 binDir = fullfile(ravenDir,'software'); -0366 -0367 binList = {fullfile(binDir,'blast+','blastp'); fullfile(binDir,'blast+','blastp.mac'); -0368 fullfile(binDir,'blast+','makeblastdb'); fullfile(binDir,'blast+','makeblastdb.mac'); -0369 fullfile(binDir,'cd-hit','cd-hit'); fullfile(binDir,'cd-hit','cd-hit.mac'); -0370 fullfile(binDir,'diamond','diamond'); fullfile(binDir,'diamond','diamond.mac'); -0371 fullfile(binDir,'hmmer','hmmbuild'); fullfile(binDir,'hmmer','hmmbuild.mac'); -0372 fullfile(binDir,'hmmer','hmmsearch'); fullfile(binDir,'hmmer','hmmsearch.mac'); -0373 fullfile(binDir,'GLPKmex','glpkcc.mexa64'); fullfile(binDir,'GLPKmex','glpkcc.mexglx'); fullfile(binDir,'GLPKmex','glpkcc.mexmaci64'); fullfile(binDir,'GLPKmex','glpkcc.mexmaca64'); -0374 fullfile(binDir,'libSBML','TranslateSBML_RAVEN.mexa64'); fullfile(binDir,'libSBML','TranslateSBML_RAVEN.mexglx'); fullfile(binDir,'libSBML','TranslateSBML_RAVEN.mexmaci64'); fullfile(binDir,'libSBML','TranslateSBML_RAVEN.mexmaca64'); -0375 fullfile(binDir,'libSBML','OutputSBML_RAVEN.mexa64'); fullfile(binDir,'libSBML','OutputSBML_RAVEN.mexglx'); fullfile(binDir,'libSBML','OutputSBML_RAVEN.mexmaci64'); fullfile(binDir,'libSBML','OutputSBML_RAVEN.mexmaca64'); -0376 fullfile(binDir,'mafft','mafft-linux64','mafft.bat'); -0377 fullfile(binDir,'mafft','mafft-mac','mafft.bat');}; -0378 -0379 for i=1:numel(binList) -0380 [status,cmdout] = system(['chmod +x "' binList{i} '"']); -0381 if status ~= 0 -0382 warning('Failed to make %s executable: %s ',binList{i},strip(cmdout)) -0383 end +0365 +0366 function res = interpretResults(results) +0367 if results.Failed==0 && results.Incomplete==0 +0368 fprintf('Pass\n'); +0369 res=true; +0370 else +0371 printOrange('Fail\n') +0372 fprintf(' Download/compile the binary and rerun checkInstallation\n'); +0373 res=false; +0374 end +0375 end +0376 +0377 function str = myStr(InputStr,len) +0378 str=InputStr; +0379 lenDiff = len - length(str); +0380 if lenDiff < 0 +0381 warning('String too long'); +0382 else +0383 str = [str blanks(lenDiff)]; 0384 end 0385 end 0386 -0387 function printOrange(stringToPrint) -0388 % printOrange -0389 % Duplicate of RAVEN/core/printOrange is also kept here, as this function -0390 % should be able to run before adding RAVEN to the MATLAB path. -0391 try useDesktop = usejava('desktop'); catch, useDesktop = false; end -0392 if useDesktop -0393 fprintf(['[\b' stringToPrint,']\b']) -0394 else -0395 fprintf(stringToPrint) -0396 end -0397 end +0387 function status = makeBinaryExecutable(ravenDir) +0388 % This function is required to run when RAVEN is downloaded as MATLAB +0389 % Add-On, in which case the file permissions are not correctly set +0390 if ispc +0391 status = 0; % No need to run on Windows +0392 return; +0393 end +0394 binDir = fullfile(ravenDir,'software'); +0395 +0396 binList = {fullfile(binDir,'blast+','blastp'); fullfile(binDir,'blast+','blastp.mac'); +0397 fullfile(binDir,'blast+','makeblastdb'); fullfile(binDir,'blast+','makeblastdb.mac'); +0398 fullfile(binDir,'cd-hit','cd-hit'); fullfile(binDir,'cd-hit','cd-hit.mac'); +0399 fullfile(binDir,'diamond','diamond'); fullfile(binDir,'diamond','diamond.mac'); +0400 fullfile(binDir,'hmmer','hmmbuild'); fullfile(binDir,'hmmer','hmmbuild.mac'); +0401 fullfile(binDir,'hmmer','hmmsearch'); fullfile(binDir,'hmmer','hmmsearch.mac'); +0402 fullfile(binDir,'GLPKmex','glpkcc.mexa64'); fullfile(binDir,'GLPKmex','glpkcc.mexglx'); fullfile(binDir,'GLPKmex','glpkcc.mexmaci64'); fullfile(binDir,'GLPKmex','glpkcc.mexmaca64'); +0403 fullfile(binDir,'libSBML','TranslateSBML_RAVEN.mexa64'); fullfile(binDir,'libSBML','TranslateSBML_RAVEN.mexglx'); fullfile(binDir,'libSBML','TranslateSBML_RAVEN.mexmaci64'); fullfile(binDir,'libSBML','TranslateSBML_RAVEN.mexmaca64'); +0404 fullfile(binDir,'libSBML','OutputSBML_RAVEN.mexa64'); fullfile(binDir,'libSBML','OutputSBML_RAVEN.mexglx'); fullfile(binDir,'libSBML','OutputSBML_RAVEN.mexmaci64'); fullfile(binDir,'libSBML','OutputSBML_RAVEN.mexmaca64'); +0405 fullfile(binDir,'mafft','mafft-linux64','mafft.bat'); +0406 fullfile(binDir,'mafft','mafft-mac','mafft.bat');}; +0407 +0408 for i=1:numel(binList) +0409 [status,cmdout] = system(['chmod +x "' binList{i} '"']); +0410 if status ~= 0 +0411 warning('Failed to make %s executable: %s ',binList{i},strip(cmdout)) +0412 end +0413 end +0414 end +0415 +0416 function printOrange(stringToPrint) +0417 % printOrange +0418 % Duplicate of RAVEN/core/printOrange is also kept here, as this function +0419 % should be able to run before adding RAVEN to the MATLAB path. +0420 try useDesktop = usejava('desktop'); catch, useDesktop = false; end +0421 if useDesktop +0422 fprintf(['[\b' stringToPrint,']\b']) +0423 else +0424 fprintf(stringToPrint) +0425 end +0426 end
Generated by m2html © 2005
\ No newline at end of file diff --git a/doc/solver/optimizeProb.html b/doc/solver/optimizeProb.html index 5fd16c557..328b92faa 100644 --- a/doc/solver/optimizeProb.html +++ b/doc/solver/optimizeProb.html @@ -144,194 +144,203 @@

SOURCE CODE ^if ~isempty(getCurrentTask) % If run in parallel, then one thread per gurobi -0094 solverparams.Threads=1; -0095 end -0096 solverparams = structUpdate(solverparams,params); -0097 -0098 % Restructering problem according to gurobi format -0099 if isfield(prob, 'csense') -0100 prob.sense = renameparams(prob.csense, {'L','G','E'}, {'<','>','='}); -0101 prob = rmfield(prob, {'csense'}); -0102 end -0103 if isfield(prob, 'osense') -0104 osense = prob.osense; -0105 prob.modelsense = renameparams(num2str(prob.osense), {'1','-1'}, {'min','max'}); -0106 prob = rmfield(prob, {'osense'}); -0107 end -0108 [prob.obj, prob.rhs, prob.vtype] = deal(prob.c, prob.b, prob.vartype); -0109 prob = rmfield(prob, {'c','b','vartype'}); -0110 -0111 resG = gurobi(prob,solverparams); -0112 -0113 try -0114 % Name output fields the same as COBRA does -0115 res.full = resG.x; -0116 res.obj = resG.objval; -0117 res.origStat = resG.status; -0118 if isfield(resG,{'pi','rc'}) -0119 res.dual = -resG.pi*osense; -0120 res.rcost = -resG.rc*osense; -0121 end -0122 if milp && strcmp(resG.status, 'TIME_LIMIT') -0123 % If res has the objval field, it succeeded, regardless of -0124 % time_limit status -0125 resG.status = 'OPTIMAL'; -0126 end -0127 switch resG.status -0128 case 'OPTIMAL' -0129 res.stat = 1; -0130 case 'UNBOUNDED' -0131 res.stat = 2; -0132 otherwise -0133 res.stat = 0; -0134 end -0135 if ~milp -0136 res.vbasis = resG.vbasis; -0137 res.cbasis = resG.cbasis; -0138 else -0139 res.mipgap = resG.mipgap; -0140 end -0141 catch -0142 res.stat = 0; -0143 res.origStat = resG.status; % useful information to have -0144 end -0145 %% Use GLPK using RAVEN-provided binary -0146 case 'glpk' -0147 solverparams.scale = 1; % Auto scaling -0148 %solverparams.tmlim = defaultparams.timeLimit; -0149 solverparams.tolbnd = defaultparams.feasTol; -0150 solverparams.toldj = defaultparams.optTol; -0151 solverparams.tolint = defaultparams.intTol; -0152 solverparams.tolobj = defaultparams.objTol; -0153 solverparams.msglev = 0; % Level of verbosity -0154 solverparams = structUpdate(solverparams,params); -0155 -0156 prob.csense = renameparams(prob.csense, {'L','G','E'}, {'U','L','S'}); -0157 -0158 if milp -0159 solverparams.tmlim = solverparams.tmlim*10; -0160 solverparams.msglev = 1; % Level of verbosity -0161 disp('Issues have been observed when using GLPK for MILP solving. Be advised to carefully observe the results, or us another solver.') -0162 end -0163 -0164 % Ensure that RAVEN glpk binary is used, return to original -0165 % directory afterwards -0166 [ravenDir,currDir]=findRAVENroot(); -0167 cd(fullfile(ravenDir,'software','GLPKmex')) -0168 [xopt, fmin, errnum, extra] = glpk(prob.c, prob.A, prob.b, prob.lb, prob.ub, prob.csense, prob.vartype, prob.osense, solverparams); -0169 cd(currDir) -0170 -0171 switch errnum % 1 = undefined; 2 = feasible; 3 = infeasible; 4 = no feasible solution; 5 = optimal; 6 = no unbounded solution -0172 case 5 -0173 res.stat = 1; % Optimal -0174 case 2 -0175 res.stat = 2; % Feasible, but not optimal -0176 otherwise -0177 res.stat = 0; -0178 end -0179 res.origStat = errnum; -0180 res.full = xopt; -0181 res.obj = fmin; -0182 res.dual = -extra.lambda*prob.osense; -0183 res.rcost = -extra.redcosts*prob.osense; -0184 %% Use scip -0185 case {'soplex','scip'} % Old 'soplex' option also allowed -0186 [xopt,fval,exitflag] = scip([], prob.c, prob.A,-prob.b, prob.b, prob.lb, prob.ub, prob.vartype); -0187 -0188 % [x,fval,exitflag,stats] = scip(H, f, A, rl, ru, lb, ub, xtype, sos, qc, nl, x0, opts) -0189 % -0190 % Input arguments*: -0191 % H - quadratic objective matrix (sparse, optional [NOT TRIL / TRIU]) -0192 % f - linear objective vector -0193 % A - linear constraint matrix (sparse) -0194 % rl - linear constraint lhs -0195 % ru - linear constraint rhs -0196 % lb - decision variable lower bounds -0197 % ub - decision variable upper bounds -0198 % xtype - string of variable integrality ('c' continuous, 'i' integer, 'b' binary) -0199 % sos - SOS structure with fields type, index and weight (see below) -0200 % qc - Quadratic Constraints structure with fields Q, l, qrl and qru (see below) -0201 % nl - Nonlinear Objective and Constraints structure (see below) -0202 % x0 - primal solution -0203 % opts - solver options (see below) -0204 % -0205 % Return arguments: -0206 % x - solution vector -0207 % fval - objective value at the solution -0208 % exitflag - exit status (see below) -0209 % stats - statistics structure -0210 % -0211 % Option Fields (all optional, see also optiset for a list): -0212 % solverOpts - specific SCIP options (list of pairs of parameter names and values) -0213 % maxiter - maximum LP solver iterations -0214 % maxnodes - maximum nodes to explore -0215 % maxtime - maximum execution time [s] -0216 % tolrfun - primal feasibility tolerance -0217 % display - solver display level [0-5] -0218 % probfile - write problem to given file -0219 % presolvedfile - write presolved problem to file -0220 % -0221 % Return Status: -0222 % 0 - Unknown -0223 % 1 - User Interrupted -0224 % 2 - Node Limit Reached -0225 % 3 - Total Node Limit Reached -0226 % 4 - Stall Node Limit Reached -0227 % 5 - Time Limit Reached -0228 % 6 - Memory Limit Reached -0229 % 7 - Gap Limit Reached -0230 % 8 - Solution Limit Reached -0231 % 9 - Solution Improvement Limit Reached -0232 % 10 - Restart Limit Reached -0233 % 11 - Problem Solved to Optimality -0234 % 12 - Problem is Infeasible -0235 % 13 - Problem is Unbounded -0236 % 14 - Problem is Either Infeasible or Unbounded -0237 -0238 res.origStat = exitflag; -0239 res.full = xopt; -0240 res.obj = fval; -0241 -0242 switch exitflag -0243 case 11 -0244 res.stat = 1; -0245 case {1, 5, 6, 7, 8, 9, 10, 13} % If user-interrupted, not sure if optimal -0246 res.stat = 2; -0247 otherwise -0248 res.stat = 0; -0249 end -0250 otherwise -0251 error('RAVEN solver not defined or unknown. Try using setRavenSolver(''solver'').'); -0252 end -0253 if res.stat>0 -0254 res.full=res.full(1:size(prob.a,2)); -0255 end -0256 end -0257 -0258 function s_merged=structUpdate(s_old,s_new) -0259 %Remove overlapping fields from first struct; -0260 %Obtain all unique names of remaining fields; -0261 %Merge both structs -0262 s_merged = rmfield(s_old, intersect(fieldnames(s_old), fieldnames(s_new))); -0263 names = [fieldnames(s_merged); fieldnames(s_new)]; -0264 s_merged = cell2struct([struct2cell(s_merged); struct2cell(s_new)], names, 1); +0093 % If run in parallel, then one thread per gurobi. getCurrentTask +0094 % requires the Parallel Computing Toolbox. +0095 inParallel = false; +0096 if license('test','Distrib_Computing_Toolbox') +0097 try +0098 inParallel = ~isempty(getCurrentTask); +0099 catch +0100 end +0101 end +0102 if inParallel +0103 solverparams.Threads=1; +0104 end +0105 solverparams = structUpdate(solverparams,params); +0106 +0107 % Restructering problem according to gurobi format +0108 if isfield(prob, 'csense') +0109 prob.sense = renameparams(prob.csense, {'L','G','E'}, {'<','>','='}); +0110 prob = rmfield(prob, {'csense'}); +0111 end +0112 if isfield(prob, 'osense') +0113 osense = prob.osense; +0114 prob.modelsense = renameparams(num2str(prob.osense), {'1','-1'}, {'min','max'}); +0115 prob = rmfield(prob, {'osense'}); +0116 end +0117 [prob.obj, prob.rhs, prob.vtype] = deal(prob.c, prob.b, prob.vartype); +0118 prob = rmfield(prob, {'c','b','vartype'}); +0119 +0120 resG = gurobi(prob,solverparams); +0121 +0122 try +0123 % Name output fields the same as COBRA does +0124 res.full = resG.x; +0125 res.obj = resG.objval; +0126 res.origStat = resG.status; +0127 if isfield(resG,{'pi','rc'}) +0128 res.dual = -resG.pi*osense; +0129 res.rcost = -resG.rc*osense; +0130 end +0131 if milp && strcmp(resG.status, 'TIME_LIMIT') +0132 % If res has the objval field, it succeeded, regardless of +0133 % time_limit status +0134 resG.status = 'OPTIMAL'; +0135 end +0136 switch resG.status +0137 case 'OPTIMAL' +0138 res.stat = 1; +0139 case 'UNBOUNDED' +0140 res.stat = 2; +0141 otherwise +0142 res.stat = 0; +0143 end +0144 if ~milp +0145 res.vbasis = resG.vbasis; +0146 res.cbasis = resG.cbasis; +0147 else +0148 res.mipgap = resG.mipgap; +0149 end +0150 catch +0151 res.stat = 0; +0152 res.origStat = resG.status; % useful information to have +0153 end +0154 %% Use GLPK using RAVEN-provided binary +0155 case 'glpk' +0156 solverparams.scale = 1; % Auto scaling +0157 %solverparams.tmlim = defaultparams.timeLimit; +0158 solverparams.tolbnd = defaultparams.feasTol; +0159 solverparams.toldj = defaultparams.optTol; +0160 solverparams.tolint = defaultparams.intTol; +0161 solverparams.tolobj = defaultparams.objTol; +0162 solverparams.msglev = 0; % Level of verbosity +0163 solverparams = structUpdate(solverparams,params); +0164 +0165 prob.csense = renameparams(prob.csense, {'L','G','E'}, {'U','L','S'}); +0166 +0167 if milp +0168 solverparams.tmlim = solverparams.tmlim*10; +0169 solverparams.msglev = 1; % Level of verbosity +0170 disp('Issues have been observed when using GLPK for MILP solving. Be advised to carefully observe the results, or us another solver.') +0171 end +0172 +0173 % Ensure that RAVEN glpk binary is used, return to original +0174 % directory afterwards +0175 [ravenDir,currDir]=findRAVENroot(); +0176 cd(fullfile(ravenDir,'software','GLPKmex')) +0177 [xopt, fmin, errnum, extra] = glpk(prob.c, prob.A, prob.b, prob.lb, prob.ub, prob.csense, prob.vartype, prob.osense, solverparams); +0178 cd(currDir) +0179 +0180 switch errnum % 1 = undefined; 2 = feasible; 3 = infeasible; 4 = no feasible solution; 5 = optimal; 6 = no unbounded solution +0181 case 5 +0182 res.stat = 1; % Optimal +0183 case 2 +0184 res.stat = 2; % Feasible, but not optimal +0185 otherwise +0186 res.stat = 0; +0187 end +0188 res.origStat = errnum; +0189 res.full = xopt; +0190 res.obj = fmin; +0191 res.dual = -extra.lambda*prob.osense; +0192 res.rcost = -extra.redcosts*prob.osense; +0193 %% Use scip +0194 case {'soplex','scip'} % Old 'soplex' option also allowed +0195 [xopt,fval,exitflag] = scip([], prob.c, prob.A,-prob.b, prob.b, prob.lb, prob.ub, prob.vartype); +0196 +0197 % [x,fval,exitflag,stats] = scip(H, f, A, rl, ru, lb, ub, xtype, sos, qc, nl, x0, opts) +0198 % +0199 % Input arguments*: +0200 % H - quadratic objective matrix (sparse, optional [NOT TRIL / TRIU]) +0201 % f - linear objective vector +0202 % A - linear constraint matrix (sparse) +0203 % rl - linear constraint lhs +0204 % ru - linear constraint rhs +0205 % lb - decision variable lower bounds +0206 % ub - decision variable upper bounds +0207 % xtype - string of variable integrality ('c' continuous, 'i' integer, 'b' binary) +0208 % sos - SOS structure with fields type, index and weight (see below) +0209 % qc - Quadratic Constraints structure with fields Q, l, qrl and qru (see below) +0210 % nl - Nonlinear Objective and Constraints structure (see below) +0211 % x0 - primal solution +0212 % opts - solver options (see below) +0213 % +0214 % Return arguments: +0215 % x - solution vector +0216 % fval - objective value at the solution +0217 % exitflag - exit status (see below) +0218 % stats - statistics structure +0219 % +0220 % Option Fields (all optional, see also optiset for a list): +0221 % solverOpts - specific SCIP options (list of pairs of parameter names and values) +0222 % maxiter - maximum LP solver iterations +0223 % maxnodes - maximum nodes to explore +0224 % maxtime - maximum execution time [s] +0225 % tolrfun - primal feasibility tolerance +0226 % display - solver display level [0-5] +0227 % probfile - write problem to given file +0228 % presolvedfile - write presolved problem to file +0229 % +0230 % Return Status: +0231 % 0 - Unknown +0232 % 1 - User Interrupted +0233 % 2 - Node Limit Reached +0234 % 3 - Total Node Limit Reached +0235 % 4 - Stall Node Limit Reached +0236 % 5 - Time Limit Reached +0237 % 6 - Memory Limit Reached +0238 % 7 - Gap Limit Reached +0239 % 8 - Solution Limit Reached +0240 % 9 - Solution Improvement Limit Reached +0241 % 10 - Restart Limit Reached +0242 % 11 - Problem Solved to Optimality +0243 % 12 - Problem is Infeasible +0244 % 13 - Problem is Unbounded +0245 % 14 - Problem is Either Infeasible or Unbounded +0246 +0247 res.origStat = exitflag; +0248 res.full = xopt; +0249 res.obj = fval; +0250 +0251 switch exitflag +0252 case 11 +0253 res.stat = 1; +0254 case {1, 5, 6, 7, 8, 9, 10, 13} % If user-interrupted, not sure if optimal +0255 res.stat = 2; +0256 otherwise +0257 res.stat = 0; +0258 end +0259 otherwise +0260 error('RAVEN solver not defined or unknown. Try using setRavenSolver(''solver'').'); +0261 end +0262 if res.stat>0 +0263 res.full=res.full(1:size(prob.a,2)); +0264 end 0265 end 0266 -0267 function paramlist = renameparams(paramlist,old,new) -0268 if ~iscell(paramlist) -0269 wasNoCell = true; -0270 paramlist={paramlist}; -0271 else -0272 wasNoCell = false; -0273 end -0274 for i=1:numel(old) -0275 paramlist = regexprep(paramlist,old{i},new{i}); -0276 end -0277 if wasNoCell -0278 paramlist=paramlist{1}; -0279 end -0280 end +0267 function s_merged=structUpdate(s_old,s_new) +0268 %Remove overlapping fields from first struct; +0269 %Obtain all unique names of remaining fields; +0270 %Merge both structs +0271 s_merged = rmfield(s_old, intersect(fieldnames(s_old), fieldnames(s_new))); +0272 names = [fieldnames(s_merged); fieldnames(s_new)]; +0273 s_merged = cell2struct([struct2cell(s_merged); struct2cell(s_new)], names, 1); +0274 end +0275 +0276 function paramlist = renameparams(paramlist,old,new) +0277 if ~iscell(paramlist) +0278 wasNoCell = true; +0279 paramlist={paramlist}; +0280 else +0281 wasNoCell = false; +0282 end +0283 for i=1:numel(old) +0284 paramlist = regexprep(paramlist,old{i},new{i}); +0285 end +0286 if wasNoCell +0287 paramlist=paramlist{1}; +0288 end +0289 end
Generated by m2html © 2005
\ No newline at end of file diff --git a/installation/checkInstallation.m b/installation/checkInstallation.m index de8c12363..40ef456fe 100755 --- a/installation/checkInstallation.m +++ b/installation/checkInstallation.m @@ -309,8 +309,37 @@ [~,res]=evalc("runtests('hmmerTests.m')"); res=interpretResults(res); if res==false - fprintf([' This is essential to run getKEGGModelFromHomology()\n'... - ' when using a FASTA file as input\n']) + fprintf([' HMMER is only required by getKEGGModelForOrganism() when\n'... + ' using a FASTA file as input. If you do not use this function,\n'... + ' this failure can be ignored.\n']) + % On Windows, HMMER runs through Windows Subsystem for Linux + if ispc + [status,~]=system('where wsl'); + if status~=0 + printOrange([' MATLAB cannot find wsl.exe. On Windows, HMMER runs through WSL.\n'... + ' Install it with "wsl --install" in a Command Prompt, then restart\n'... + ' MATLAB (it inherits the system PATH at startup) and rerun checkInstallation\n']) + else + [status,~]=system('wsl -e true'); + if status~=0 + printOrange([' MATLAB finds wsl.exe, but cannot start a Linux distribution\n'... + ' through it. Run "wsl -l -v" in a Command Prompt to check that a\n'... + ' distribution is installed and set up, then restart MATLAB and\n'... + ' rerun checkInstallation\n']) + else + try + [status,~]=system(['wsl "' getWSLpath(fullfile(ravenDir,'software','hmmer','hmmsearch')) '" -h']); + catch + status=1; + end + if status~=0 + printOrange([' WSL works, but the RAVEN HMMER binary cannot be run in it.\n'... + ' Check that the files in software/hmmer are present and that\n'... + ' the drive holding RAVEN is accessible from WSL\n']) + end + end + end + end end if developMode diff --git a/solver/optimizeProb.m b/solver/optimizeProb.m index 3430c45d3..bcb137e98 100755 --- a/solver/optimizeProb.m +++ b/solver/optimizeProb.m @@ -90,7 +90,16 @@ solverparams.FeasibilityTol = defaultparams.feasTol; solverparams.OptimalityTol = defaultparams.optTol; solverparams.Presolve = 2; - if ~isempty(getCurrentTask) % If run in parallel, then one thread per gurobi + % If run in parallel, then one thread per gurobi. getCurrentTask + % requires the Parallel Computing Toolbox. + inParallel = false; + if license('test','Distrib_Computing_Toolbox') + try + inParallel = ~isempty(getCurrentTask); + catch + end + end + if inParallel solverparams.Threads=1; end solverparams = structUpdate(solverparams,params); From 2156ff13ea86ac080e980c815e2b8a914e57df1c Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 30 Sep 2026 21:05:44 +0200 Subject: [PATCH 2/3] fix: outputs and metabolite merging in replaceMets (#756) --- core/replaceMets.m | 82 ++++---- doc/core/removeMets.html | 2 +- doc/core/replaceMets.html | 176 ++++++++---------- .../unit_tests/modelCurationTests.html | 92 ++++++++- testing/unit_tests/modelCurationTests.m | 88 +++++++++ 5 files changed, 296 insertions(+), 144 deletions(-) diff --git a/core/replaceMets.m b/core/replaceMets.m index 9301d8628..9cfa0804c 100755 --- a/core/replaceMets.m +++ b/core/replaceMets.m @@ -20,7 +20,7 @@ % Output: % model model structure with selected metabolites replaced % removedRxns identifiers of duplicate reactions that were removed -% idxDuplRxns index of removedRxns in original model +% idxDuplRxns indices of removedRxns in the input model % % Note: This function is useful when the model contains both 'oxygen' and % 'o2' as metabolite names. If 'oxygen' and 'o2' are identifiers instead, @@ -59,10 +59,11 @@ error('The to-be-replaced metabolite cannot be found in the model.'); end -rxnsWithMet = find(model.S(metIdx,:)); +[~, rxnsWithMet] = find(model.S(metIdx,:)); +rxnsWithMet = unique(rxnsWithMet); if verbose==true - fprintf('\n\nThe following reactions contain the to-be-replaced metabolite as reactant:\n') - fprintf(strjoin(model.rxns(rxnsWithMet),'\n')) + fprintf('\n\nThe following reactions involve the to-be-replaced metabolite:\n') + fprintf('%s\n',model.rxns{rxnsWithMet}) fprintf('\n') end @@ -88,8 +89,11 @@ idxDelete=[]; if identifiers + if any(ismember(metIdx,repIdx)) + error('The to-be-replaced metabolite and the replacement metabolite are the same.'); + end originalStoch = model.S(metIdx,rxnsWithMet); - model.S(repIdx,rxnsWithMet) = originalStoch; + model.S(repIdx,rxnsWithMet) = model.S(repIdx,rxnsWithMet) + originalStoch; model.S(metIdx,rxnsWithMet) = 0; idxDelete = metIdx; else @@ -99,54 +103,38 @@ % the to-be-replace metabolite ID deleted. % Build list of metaboliteName[compartment] - metCompsN =cellstr(num2str(model.metComps)); - map = containers.Map(cellstr(num2str(transpose(1:length(model.comps)))),model.comps); - metCompsN = map.values(metCompsN); - metCompsN = strcat(lower(model.metNames),'[',metCompsN,']'); + comps = model.comps(:); + metCompsN = strcat(lower(model.metNames(:)),'[',comps(model.metComps(:)),']'); - for i = 1:length(repIdx) - metCompsNidx=find(strcmp(metCompsN(repIdx(i)), metCompsN)); - if length(metCompsNidx)>1 - for j = 2:length(metCompsNidx) - model.S(metCompsNidx(1),:) = model.S(metCompsNidx(1),:) + model.S(metCompsNidx(j),:); - idxDelete=[idxDelete; metCompsNidx(j)]; % Make list of metabolite IDs to delete - end + % Within each metaboliteName[compartment] group, the first replacement + % metabolite is kept and the other members are merged into it. + groups = unique(metCompsN(repIdx),'stable'); + for i = 1:numel(groups) + metCompsNidx = find(strcmp(groups{i}, metCompsN)); + keepIdx = metCompsNidx(find(ismember(metCompsNidx,repIdx),1)); + others = metCompsNidx(metCompsNidx ~= keepIdx); + if ~isempty(others) + model.S(keepIdx,:) = model.S(keepIdx,:) + sum(model.S(others,:),1); + idxDelete = [idxDelete; others]; % Make list of metabolite IDs to delete end end end - if ~isempty(idxDelete) - model.S(idxDelete,:) =[]; - model.mets(idxDelete) = []; - model.metNames(idxDelete) = []; - model.metComps(idxDelete) = []; - model.b(idxDelete) = []; - if isfield(model,'metFormulas') - model.metFormulas(idxDelete) = []; - end - if isfield(model,'unconstrained') - model.unconstrained(idxDelete) = []; - end - if isfield(model,'metMiriams') - model.metMiriams(idxDelete) = []; - end - if isfield(model,'metCharges') - model.metCharges(idxDelete) = []; - end - if isfield(model,'metDeltaG') - model.metDeltaG(idxDelete) = []; - end - if isfield(model,'inchis') - model.inchis(idxDelete) = []; - end - if isfield(model,'metSmiles') - model.metSmiles(idxDelete) = []; - end - if isfield(model,'metFrom') - model.metFrom(idxDelete) = []; - end + repMets = model.mets(repIdx); + model = removeMets(model,idxDelete); + % Indices of the replacement metabolites shift when metabolites with a + % lower index are deleted. + [~, repIdx] = ismember(repMets,model.mets); + repIdx(repIdx==0) = []; end % This could now have created duplicate reactions. Contract model. -model=contractModel(model,[],repIdx); +if isempty(idxDelete) + removedRxns = cell(0,1); + idxDuplRxns = zeros(0,1); +else + originalRxns = model.rxns; + [model, removedRxns] = contractModel(model,[],repIdx); + [~, idxDuplRxns] = ismember(removedRxns,originalRxns); +end end diff --git a/doc/core/removeMets.html b/doc/core/removeMets.html index 4a3d2d269..2df79cdf2 100644 --- a/doc/core/removeMets.html +++ b/doc/core/removeMets.html @@ -58,7 +58,7 @@

CROSS-REFERENCE INFORMATION ^ <li><a href=convertCharArray convertCharArray
  • getIndexes getIndexes
  • removeReactions removeReactions
    • This function is called by: +
  • mergeCompartments mergeCompartments
  • removeReactions removeReactions
  • replaceMets replaceMets
  • simplifyModel simplifyModel
    • diff --git a/doc/core/replaceMets.html b/doc/core/replaceMets.html index 6f4319513..db6d7f9b3 100644 --- a/doc/core/replaceMets.html +++ b/doc/core/replaceMets.html @@ -48,7 +48,7 @@

    DESCRIPTION ^DESCRIPTION ^CROSS-REFERENCE INFORMATION ^

    This function calls: +
  • contractModel contractModel
  • removeMets removeMets
    • This function is called by:
    @@ -91,7 +91,7 @@

    SOURCE CODE ^% Output: 0021 % model model structure with selected metabolites replaced 0022 % removedRxns identifiers of duplicate reactions that were removed -0023 % idxDuplRxns index of removedRxns in original model +0023 % idxDuplRxns indices of removedRxns in the input model 0024 % 0025 % Note: This function is useful when the model contains both 'oxygen' and 0026 % 'o2' as metabolite names. If 'oxygen' and 'o2' are identifiers instead, @@ -130,97 +130,85 @@

    SOURCE CODE ^'The to-be-replaced metabolite cannot be found in the model.'); 0060 end 0061 -0062 rxnsWithMet = find(model.S(metIdx,:)); -0063 if verbose==true -0064 fprintf('\n\nThe following reactions contain the to-be-replaced metabolite as reactant:\n') -0065 fprintf(strjoin(model.rxns(rxnsWithMet),'\n')) -0066 fprintf('\n') -0067 end -0068 -0069 model.metNames(metIdx) = model.metNames(repIdx(1)); -0070 if isfield(model,'metFormulas') -0071 model.metFormulas(metIdx) = model.metFormulas(repIdx(1)); -0072 end -0073 if isfield(model,'metMiriams') -0074 model.metMiriams(metIdx) = model.metMiriams(repIdx(1)); -0075 end -0076 if isfield(model,'metCharges') -0077 model.metCharges(metIdx) = model.metCharges(repIdx(1)); -0078 end -0079 if isfield(model,'metDeltaG') -0080 model.metDeltaG(metIdx) = model.metDeltaG(repIdx(1)); -0081 end -0082 if isfield(model,'inchis') -0083 model.inchis(metIdx) = model.inchis(repIdx(1)); -0084 end -0085 if isfield(model,'metSmiles') -0086 model.metSmiles(metIdx) = model.metSmiles(repIdx(1)); -0087 end -0088 -0089 idxDelete=[]; -0090 if identifiers -0091 originalStoch = model.S(metIdx,rxnsWithMet); -0092 model.S(repIdx,rxnsWithMet) = originalStoch; -0093 model.S(metIdx,rxnsWithMet) = 0; -0094 idxDelete = metIdx; -0095 else -0096 % Run through replacement metabolites and their compartments. If any of the -0097 % to-be-replaced metabolites is already present (checked by -0098 % metaboliteName[compartment], then the replacement metabolite is kept and -0099 % the to-be-replace metabolite ID deleted. -0100 -0101 % Build list of metaboliteName[compartment] -0102 metCompsN =cellstr(num2str(model.metComps)); -0103 map = containers.Map(cellstr(num2str(transpose(1:length(model.comps)))),model.comps); -0104 metCompsN = map.values(metCompsN); -0105 metCompsN = strcat(lower(model.metNames),'[',metCompsN,']'); -0106 -0107 for i = 1:length(repIdx) -0108 metCompsNidx=find(strcmp(metCompsN(repIdx(i)), metCompsN)); -0109 if length(metCompsNidx)>1 -0110 for j = 2:length(metCompsNidx) -0111 model.S(metCompsNidx(1),:) = model.S(metCompsNidx(1),:) + model.S(metCompsNidx(j),:); -0112 idxDelete=[idxDelete; metCompsNidx(j)]; % Make list of metabolite IDs to delete -0113 end -0114 end -0115 end -0116 end -0117 -0118 if ~isempty(idxDelete) -0119 model.S(idxDelete,:) =[]; -0120 model.mets(idxDelete) = []; -0121 model.metNames(idxDelete) = []; -0122 model.metComps(idxDelete) = []; -0123 model.b(idxDelete) = []; -0124 if isfield(model,'metFormulas') -0125 model.metFormulas(idxDelete) = []; -0126 end -0127 if isfield(model,'unconstrained') -0128 model.unconstrained(idxDelete) = []; -0129 end -0130 if isfield(model,'metMiriams') -0131 model.metMiriams(idxDelete) = []; -0132 end -0133 if isfield(model,'metCharges') -0134 model.metCharges(idxDelete) = []; -0135 end -0136 if isfield(model,'metDeltaG') -0137 model.metDeltaG(idxDelete) = []; -0138 end -0139 if isfield(model,'inchis') -0140 model.inchis(idxDelete) = []; -0141 end -0142 if isfield(model,'metSmiles') -0143 model.metSmiles(idxDelete) = []; -0144 end -0145 if isfield(model,'metFrom') -0146 model.metFrom(idxDelete) = []; -0147 end -0148 end -0149 -0150 % This could now have created duplicate reactions. Contract model. -0151 model=contractModel(model,[],repIdx); -0152 end +0062 [~, rxnsWithMet] = find(model.S(metIdx,:)); +0063 rxnsWithMet = unique(rxnsWithMet); +0064 if verbose==true +0065 fprintf('\n\nThe following reactions involve the to-be-replaced metabolite:\n') +0066 fprintf('%s\n',model.rxns{rxnsWithMet}) +0067 fprintf('\n') +0068 end +0069 +0070 model.metNames(metIdx) = model.metNames(repIdx(1)); +0071 if isfield(model,'metFormulas') +0072 model.metFormulas(metIdx) = model.metFormulas(repIdx(1)); +0073 end +0074 if isfield(model,'metMiriams') +0075 model.metMiriams(metIdx) = model.metMiriams(repIdx(1)); +0076 end +0077 if isfield(model,'metCharges') +0078 model.metCharges(metIdx) = model.metCharges(repIdx(1)); +0079 end +0080 if isfield(model,'metDeltaG') +0081 model.metDeltaG(metIdx) = model.metDeltaG(repIdx(1)); +0082 end +0083 if isfield(model,'inchis') +0084 model.inchis(metIdx) = model.inchis(repIdx(1)); +0085 end +0086 if isfield(model,'metSmiles') +0087 model.metSmiles(metIdx) = model.metSmiles(repIdx(1)); +0088 end +0089 +0090 idxDelete=[]; +0091 if identifiers +0092 if any(ismember(metIdx,repIdx)) +0093 error('The to-be-replaced metabolite and the replacement metabolite are the same.'); +0094 end +0095 originalStoch = model.S(metIdx,rxnsWithMet); +0096 model.S(repIdx,rxnsWithMet) = model.S(repIdx,rxnsWithMet) + originalStoch; +0097 model.S(metIdx,rxnsWithMet) = 0; +0098 idxDelete = metIdx; +0099 else +0100 % Run through replacement metabolites and their compartments. If any of the +0101 % to-be-replaced metabolites is already present (checked by +0102 % metaboliteName[compartment], then the replacement metabolite is kept and +0103 % the to-be-replace metabolite ID deleted. +0104 +0105 % Build list of metaboliteName[compartment] +0106 comps = model.comps(:); +0107 metCompsN = strcat(lower(model.metNames(:)),'[',comps(model.metComps(:)),']'); +0108 +0109 % Within each metaboliteName[compartment] group, the first replacement +0110 % metabolite is kept and the other members are merged into it. +0111 groups = unique(metCompsN(repIdx),'stable'); +0112 for i = 1:numel(groups) +0113 metCompsNidx = find(strcmp(groups{i}, metCompsN)); +0114 keepIdx = metCompsNidx(find(ismember(metCompsNidx,repIdx),1)); +0115 others = metCompsNidx(metCompsNidx ~= keepIdx); +0116 if ~isempty(others) +0117 model.S(keepIdx,:) = model.S(keepIdx,:) + sum(model.S(others,:),1); +0118 idxDelete = [idxDelete; others]; % Make list of metabolite IDs to delete +0119 end +0120 end +0121 end +0122 if ~isempty(idxDelete) +0123 repMets = model.mets(repIdx); +0124 model = removeMets(model,idxDelete); +0125 % Indices of the replacement metabolites shift when metabolites with a +0126 % lower index are deleted. +0127 [~, repIdx] = ismember(repMets,model.mets); +0128 repIdx(repIdx==0) = []; +0129 end +0130 +0131 % This could now have created duplicate reactions. Contract model. +0132 if isempty(idxDelete) +0133 removedRxns = cell(0,1); +0134 idxDuplRxns = zeros(0,1); +0135 else +0136 originalRxns = model.rxns; +0137 [model, removedRxns] = contractModel(model,[],repIdx); +0138 [~, idxDuplRxns] = ismember(removedRxns,originalRxns); +0139 end +0140 end
    Generated by m2html © 2005
    \ No newline at end of file diff --git a/doc/testing/unit_tests/modelCurationTests.html b/doc/testing/unit_tests/modelCurationTests.html index dde7fe80c..0846789f0 100644 --- a/doc/testing/unit_tests/modelCurationTests.html +++ b/doc/testing/unit_tests/modelCurationTests.html @@ -42,7 +42,7 @@

    CROSS-REFERENCE INFORMATION ^ <h2><a name=SUBFUNCTIONS ^

    +
  • function addExchangeRxnsTest(testCase)
  • function addTransportTest(testCase)
  • function addGenesRavenTest(testCase)
  • function removeGenesTest(testCase)
  • function addMetsTest(testCase)
  • function addMets_oneCompTest(testCase)
  • function removeMetsTest(testCase)
  • function addRxnsTest(testCase)
  • function removeReactionsTest(testCase)
  • function getMetsInCompTest(testCase)
  • function getRxnsInCompTest(testCase)
  • function getAllRxnsFromGenesTest(testCase)
  • function changeRxnsTest(testCase)
  • function removeBadRxnsTest(testCase)
  • function changeGrRulesTest(testCase)
  • function copyToCompsTest(testCase)
  • function setExchangeBoundsTest(testCase)
  • function addRxnsGenesMetsTest(testCase)
  • function replaceMetsOutputsTest(testCase)
  • function replaceMetsSharedReactionTest(testCase)
  • function replaceMetsNamesTest(testCase)
  • function replaceMetsNamesRepeatedReplacementTest(testCase)

    • SOURCE CODE ^

    0001 %run this test case with the command
@@ -537,7 +537,95 @@ 
    SOURCE CODE ^end
    +0493 end +0494 +0495 function replaceMetsOutputsTest(testCase) +0496 sourceDir = fileparts(which(mfilename)); +0497 load(fullfile(sourceDir,'test_data','ecoli_textbook.mat'), 'model'); +0498 +0499 %Add a duplicate of atp_c and a copy of ATPM that uses it +0500 metsToAdd.mets={'atpDup_c'}; +0501 metsToAdd.metNames={'atp duplicate'}; +0502 metsToAdd.compartments={'c'}; +0503 evalc('model=addMets(model,metsToAdd);'); +0504 rxnsToAdd.rxns={'ATPM_dup'}; +0505 rxnsToAdd.equations={'atpDup_c + h2o_c => adp_c + h_c + pi_c'}; +0506 evalc('model=addRxns(model,rxnsToAdd,1);'); +0507 +0508 [modelNew,removedRxns,idxDuplRxns]=replaceMets(model,'atpDup_c','atp_c',false,true); +0509 +0510 verifyEqual(testCase,numel(modelNew.rxns),numel(model.rxns)-1) +0511 verifyEmpty(testCase,find(strcmp(modelNew.mets,'atpDup_c'))) +0512 verifyEqual(testCase,removedRxns,{'ATPM_dup'}) +0513 verifyEqual(testCase,idxDuplRxns,numel(model.rxns)) +0514 end +0515 +0516 function replaceMetsSharedReactionTest(testCase) +0517 %The replaced metabolite occurs in a reaction that already contains the +0518 %replacement metabolite; the stoichiometries are summed. +0519 sourceDir = fileparts(which(mfilename)); +0520 load(fullfile(sourceDir,'test_data','ecoli_textbook.mat'), 'model'); +0521 +0522 metsToAdd.mets={'atpDup_c'}; +0523 metsToAdd.metNames={'atp duplicate'}; +0524 metsToAdd.compartments={'c'}; +0525 evalc('model=addMets(model,metsToAdd);'); +0526 rxnsToAdd.rxns={'ATPshared'}; +0527 rxnsToAdd.equations={'atp_c + atpDup_c + h2o_c => adp_c + h_c + pi_c'}; +0528 evalc('model=addRxns(model,rxnsToAdd,1);'); +0529 +0530 modelNew=replaceMets(model,'atpDup_c','atp_c',false,true); +0531 +0532 rxnIdx=strcmp(modelNew.rxns,'ATPshared'); +0533 verifyEqual(testCase,full(modelNew.S(strcmp(modelNew.mets,'atp_c'),rxnIdx)),-2) +0534 verifyEqual(testCase,numel(modelNew.mets),numel(model.mets)-1) +0535 verifyEmpty(testCase,find(strcmp(modelNew.mets,'atpDup_c'))) +0536 end +0537 +0538 function replaceMetsNamesTest(testCase) +0539 %The replacement metabolite has a higher index than the to-be-replaced +0540 %metabolite, and the removed duplicate reaction is not the last reaction. +0541 sourceDir = fileparts(which(mfilename)); +0542 load(fullfile(sourceDir,'test_data','ecoli_textbook.mat'), 'model'); +0543 +0544 metsToAdd.mets={'atpNew_c'}; +0545 metsToAdd.metNames={'ATPnew'}; +0546 metsToAdd.compartments={'c'}; +0547 evalc('model=addMets(model,metsToAdd);'); +0548 rxnsToAdd.rxns={'ATPM_dup','ATPM_other'}; +0549 rxnsToAdd.equations={'atpNew_c + h2o_c => adp_c + h_c + pi_c','h2o_c => h_c + pi_c'}; +0550 evalc('model=addRxns(model,rxnsToAdd,1);'); +0551 +0552 [modelNew,removedRxns,idxDuplRxns]=replaceMets(model,'ATP','ATPnew'); +0553 +0554 verifyEqual(testCase,numel(modelNew.rxns),numel(model.rxns)-1) +0555 verifyEqual(testCase,removedRxns,{'ATPM_dup'}) +0556 verifyEqual(testCase,idxDuplRxns,numel(model.rxns)-1) +0557 verifyEqual(testCase,model.rxns(idxDuplRxns),removedRxns) +0558 verifyEmpty(testCase,find(strcmp(modelNew.mets,'atp_c'))) +0559 end +0560 +0561 function replaceMetsNamesRepeatedReplacementTest(testCase) +0562 %The replacement name occurs twice in the same compartment; the +0563 %stoichiometry is added once. +0564 sourceDir = fileparts(which(mfilename)); +0565 load(fullfile(sourceDir,'test_data','ecoli_textbook.mat'), 'model'); +0566 +0567 metsToAdd.mets={'atp2_c','atpNew_c'}; +0568 metsToAdd.metNames={'ATP2','ATPnew'}; +0569 metsToAdd.compartments={'c','c'}; +0570 evalc('model=addMets(model,metsToAdd);'); +0571 model.metNames{strcmp(model.mets,'atp2_c')}='ATP'; +0572 rxnsToAdd.rxns={'rNew'}; +0573 rxnsToAdd.equations={'atpNew_c => adp_c'}; +0574 evalc('model=addRxns(model,rxnsToAdd,1);'); +0575 +0576 modelNew=replaceMets(model,'ATPnew','ATP'); +0577 +0578 verifyEqual(testCase,full(modelNew.S(strcmp(modelNew.mets,'atp_c'),strcmp(modelNew.rxns,'rNew'))),-1) +0579 verifyEmpty(testCase,find(strcmp(modelNew.mets,'atpNew_c'))) +0580 verifyEmpty(testCase,find(strcmp(modelNew.mets,'atp2_c'))) +0581 end
    Generated by m2html © 2005
    \ No newline at end of file diff --git a/testing/unit_tests/modelCurationTests.m b/testing/unit_tests/modelCurationTests.m index 62e3dc962..bd6f80900 100755 --- a/testing/unit_tests/modelCurationTests.m +++ b/testing/unit_tests/modelCurationTests.m @@ -491,3 +491,91 @@ function addRxnsGenesMetsTest(testCase) verifyEqual(testCase,modelNew,modelManual) end + +function replaceMetsOutputsTest(testCase) +sourceDir = fileparts(which(mfilename)); +load(fullfile(sourceDir,'test_data','ecoli_textbook.mat'), 'model'); + +%Add a duplicate of atp_c and a copy of ATPM that uses it +metsToAdd.mets={'atpDup_c'}; +metsToAdd.metNames={'atp duplicate'}; +metsToAdd.compartments={'c'}; +evalc('model=addMets(model,metsToAdd);'); +rxnsToAdd.rxns={'ATPM_dup'}; +rxnsToAdd.equations={'atpDup_c + h2o_c => adp_c + h_c + pi_c'}; +evalc('model=addRxns(model,rxnsToAdd,1);'); + +[modelNew,removedRxns,idxDuplRxns]=replaceMets(model,'atpDup_c','atp_c',false,true); + +verifyEqual(testCase,numel(modelNew.rxns),numel(model.rxns)-1) +verifyEmpty(testCase,find(strcmp(modelNew.mets,'atpDup_c'))) +verifyEqual(testCase,removedRxns,{'ATPM_dup'}) +verifyEqual(testCase,idxDuplRxns,numel(model.rxns)) +end + +function replaceMetsSharedReactionTest(testCase) +%The replaced metabolite occurs in a reaction that already contains the +%replacement metabolite; the stoichiometries are summed. +sourceDir = fileparts(which(mfilename)); +load(fullfile(sourceDir,'test_data','ecoli_textbook.mat'), 'model'); + +metsToAdd.mets={'atpDup_c'}; +metsToAdd.metNames={'atp duplicate'}; +metsToAdd.compartments={'c'}; +evalc('model=addMets(model,metsToAdd);'); +rxnsToAdd.rxns={'ATPshared'}; +rxnsToAdd.equations={'atp_c + atpDup_c + h2o_c => adp_c + h_c + pi_c'}; +evalc('model=addRxns(model,rxnsToAdd,1);'); + +modelNew=replaceMets(model,'atpDup_c','atp_c',false,true); + +rxnIdx=strcmp(modelNew.rxns,'ATPshared'); +verifyEqual(testCase,full(modelNew.S(strcmp(modelNew.mets,'atp_c'),rxnIdx)),-2) +verifyEqual(testCase,numel(modelNew.mets),numel(model.mets)-1) +verifyEmpty(testCase,find(strcmp(modelNew.mets,'atpDup_c'))) +end + +function replaceMetsNamesTest(testCase) +%The replacement metabolite has a higher index than the to-be-replaced +%metabolite, and the removed duplicate reaction is not the last reaction. +sourceDir = fileparts(which(mfilename)); +load(fullfile(sourceDir,'test_data','ecoli_textbook.mat'), 'model'); + +metsToAdd.mets={'atpNew_c'}; +metsToAdd.metNames={'ATPnew'}; +metsToAdd.compartments={'c'}; +evalc('model=addMets(model,metsToAdd);'); +rxnsToAdd.rxns={'ATPM_dup','ATPM_other'}; +rxnsToAdd.equations={'atpNew_c + h2o_c => adp_c + h_c + pi_c','h2o_c => h_c + pi_c'}; +evalc('model=addRxns(model,rxnsToAdd,1);'); + +[modelNew,removedRxns,idxDuplRxns]=replaceMets(model,'ATP','ATPnew'); + +verifyEqual(testCase,numel(modelNew.rxns),numel(model.rxns)-1) +verifyEqual(testCase,removedRxns,{'ATPM_dup'}) +verifyEqual(testCase,idxDuplRxns,numel(model.rxns)-1) +verifyEqual(testCase,model.rxns(idxDuplRxns),removedRxns) +verifyEmpty(testCase,find(strcmp(modelNew.mets,'atp_c'))) +end + +function replaceMetsNamesRepeatedReplacementTest(testCase) +%The replacement name occurs twice in the same compartment; the +%stoichiometry is added once. +sourceDir = fileparts(which(mfilename)); +load(fullfile(sourceDir,'test_data','ecoli_textbook.mat'), 'model'); + +metsToAdd.mets={'atp2_c','atpNew_c'}; +metsToAdd.metNames={'ATP2','ATPnew'}; +metsToAdd.compartments={'c','c'}; +evalc('model=addMets(model,metsToAdd);'); +model.metNames{strcmp(model.mets,'atp2_c')}='ATP'; +rxnsToAdd.rxns={'rNew'}; +rxnsToAdd.equations={'atpNew_c => adp_c'}; +evalc('model=addRxns(model,rxnsToAdd,1);'); + +modelNew=replaceMets(model,'ATPnew','ATP'); + +verifyEqual(testCase,full(modelNew.S(strcmp(modelNew.mets,'atp_c'),strcmp(modelNew.rxns,'rNew'))),-1) +verifyEmpty(testCase,find(strcmp(modelNew.mets,'atpNew_c'))) +verifyEmpty(testCase,find(strcmp(modelNew.mets,'atp2_c'))) +end From ee5f1bb89a055856b8e7cc43dec53806a6d76511 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 7 Oct 2026 21:38:24 +0200 Subject: [PATCH 3/3] fix: read RAVEN 3 model files, and tolerate a non-double rev (#759) * fix: read RAVEN 3 model files, and tolerate a non-double rev RAVEN 3 and raven-toolbox write YAML this version cannot read. The metaData map is tagged (`- metaData: !!omap`), which the exact-match case missed, so the whole block was skipped and the model came back with no id, name, version or date. A reaction's free text moved from `rxnNotes` to `notes`, which aborted the read outright with "Unknown entry in yaml file". EC codes and SMILES moved inside the per-entry `annotation` block, where they were gathered as MIRIAMs instead of filling eccodes and metSmiles. An empty gene list is written inline, and the GECKO flag is spelled `gecko_light`. All of these are accepted on read now; what this version writes is unchanged, so its own output stays valid for older RAVEN. importModel also only found a model note inside a , which RAVEN's exportModel writes but a writer serialising a notes dictionary does not; such a writer emits one "

    key: value

    " per entry, so the labelled entry is taken when there is one and the text otherwise. mergeLinear combines rev with a sparse logical, which MATLAB refuses for an integer class ("Sparse integer array arithmetic operations are not supported"). A model built outside RAVEN can carry an integer rev, so it is coerced once on entry rather than failing several steps into ftINIT. metComps is documented as an integer column: a reader that turns the 1-based index into a position needs an integer, and MATLAB stores an integer-valued double narrowed to an integer type in any case. * docs: regenerate HTML for the changed functions updateDocumentation regenerates the whole doc/ tree, but the committed files were generated on Windows and every other one differs only in the path separator of its index link and title (external\kegg vs external/kegg). Only the three functions this branch changes are updated, so the diff stays about the change rather than about the platform it was regenerated on. * fix: compare rev rather than coerce it, and cover ftINIT too The earlier fix coerced reducedModel.rev to double on entry to mergeLinear, which changed the class of rev in the model it returns -- a type the caller did not ask to have changed. Comparing instead matches what every other use of rev in RAVEN already does (rev == 1, rev ~= 0, find(rev)), works for a double, an integer or a logical alike, and leaves the model untouched. On the released Human-GEM 2.1.0 .mat, whose rev is int64, both forms return the same 8751 reactions with the same S and the same rev values, but the comparison leaves rev as int64 instead of silently making it double. ftINIT had the same pattern and the same failure: any(prepData.refModel.S(metSel,:) < 0,1) & prepData.refModel.rev.' a sparse logical combined with rev, which MATLAB refuses for an integer array. It only runs when metabolomics data is supplied, which is why #1143 did not reach it. Every other use of rev compares first and is already class-agnostic. --- INIT/ftINIT.m | 2 +- INIT/mergeLinear.m | 8 +- doc/INIT/ftINIT.html | 2 +- doc/INIT/mergeLinear.html | 300 +++++----- doc/io/importModel.html | 807 +++++++++++++-------------- doc/io/readYAMLmodel.html | 879 ++++++++++++++++-------------- io/importModel.m | 21 +- io/readYAMLmodel.m | 51 +- readme/RAVEN_structure_fields.csv | 2 +- 9 files changed, 1090 insertions(+), 982 deletions(-) diff --git a/INIT/ftINIT.m b/INIT/ftINIT.m index ab5f3e030..9dba76ad1 100644 --- a/INIT/ftINIT.m +++ b/INIT/ftINIT.m @@ -135,7 +135,7 @@ %Get the matching mets metSel = ismember(upper(prepData.refModel.metNames),upper(metabolomicsData{i})); prodRxnsSel = any(prepData.refModel.S(metSel,:) > 0,1) | ... %direct producers - (any(prepData.refModel.S(metSel,:) < 0,1) & prepData.refModel.rev.'); %reversible reactions that are consumers + (any(prepData.refModel.S(metSel,:) < 0,1) & prepData.refModel.rev.' ~= 0); %reversible reactions that are consumers %convert the production rxns from refModel to minModel prepData.groupIds; [~,ia,ib] = intersect(prepData.minModel.rxns,prepData.refModel.rxns); diff --git a/INIT/mergeLinear.m b/INIT/mergeLinear.m index 1378353b2..114d7228d 100644 --- a/INIT/mergeLinear.m +++ b/INIT/mergeLinear.m @@ -62,8 +62,12 @@ for i=1:numel(twoNonZero) involvedRxns=find(reducedModel.S(twoNonZero(i),:)); %Check that we can have one positive and one negative - pos = sum(reducedModel.S(twoNonZero(i),involvedRxns).' > 0 | reducedModel.rev(involvedRxns)); - neg = sum(reducedModel.S(twoNonZero(i),involvedRxns).' < 0 | reducedModel.rev(involvedRxns)); + %rev is compared rather than used as-is: the other operand is a + %sparse logical, and MATLAB refuses to combine one with an integer + %array ("Sparse integer array arithmetic operations are not + %supported"). A model built outside RAVEN can carry an integer rev. + pos = sum(reducedModel.S(twoNonZero(i),involvedRxns).' > 0 | reducedModel.rev(involvedRxns) ~= 0); + neg = sum(reducedModel.S(twoNonZero(i),involvedRxns).' < 0 | reducedModel.rev(involvedRxns) ~= 0); %Check so that one or both of the reactions haven't been merged diff --git a/doc/INIT/ftINIT.html b/doc/INIT/ftINIT.html index 6133e894d..8e8a54983 100644 --- a/doc/INIT/ftINIT.html +++ b/doc/INIT/ftINIT.html @@ -270,7 +270,7 @@

    SOURCE CODE ^%Get the matching mets 0136 metSel = ismember(upper(prepData.refModel.metNames),upper(metabolomicsData{i})); 0137 prodRxnsSel = any(prepData.refModel.S(metSel,:) > 0,1) | ... %direct producers -0138 (any(prepData.refModel.S(metSel,:) < 0,1) & prepData.refModel.rev.'); %reversible reactions that are consumers +0138 (any(prepData.refModel.S(metSel,:) < 0,1) & prepData.refModel.rev.' ~= 0); %reversible reactions that are consumers 0139 %convert the production rxns from refModel to minModel 0140 prepData.groupIds; 0141 [~,ia,ib] = intersect(prepData.minModel.rxns,prepData.refModel.rxns); diff --git a/doc/INIT/mergeLinear.html b/doc/INIT/mergeLinear.html index 16467e80e..5d90456af 100644 --- a/doc/INIT/mergeLinear.html +++ b/doc/INIT/mergeLinear.html @@ -125,157 +125,161 @@

    SOURCE CODE ^for i=1:numel(twoNonZero) 0063 involvedRxns=find(reducedModel.S(twoNonZero(i),:)); 0064 %Check that we can have one positive and one negative -0065 pos = sum(reducedModel.S(twoNonZero(i),involvedRxns).' > 0 | reducedModel.rev(involvedRxns)); -0066 neg = sum(reducedModel.S(twoNonZero(i),involvedRxns).' < 0 | reducedModel.rev(involvedRxns)); -0067 -0068 -0069 %Check so that one or both of the reactions haven't been merged -0070 %already -0071 if numel(involvedRxns)==2 && isempty(intersect(bannedIndexes,involvedRxns)) && pos >= 1 && neg >= 1 -0072 %first, take care of a special case: If the first reaction is producing the metabolite and if it is reversible, -0073 %and the second is also producing it and is not reversible, change the order - the code below will not work otherwise -0074 if reducedModel.rev(involvedRxns(1)) && (~reducedModel.rev(involvedRxns(2))) && ... -0075 (reducedModel.S(twoNonZero(i),involvedRxns(1)) > 0) && (reducedModel.S(twoNonZero(i),involvedRxns(2)) > 0) -0076 involvedRxns = flip(involvedRxns); -0077 end -0078 -0079 %first make sure the first reaction is producing the metabolite -0080 if reducedModel.S(twoNonZero(i),involvedRxns(1)) < 0 -0081 %it is not producing the metabolite - fix that -0082 %first choice: use the second reaction as producer if it is producing -0083 if reducedModel.S(twoNonZero(i),involvedRxns(2)) > 0 -0084 involvedRxns = flip(involvedRxns);%make the second the first -0085 else -0086 %now we know that the second reaction is not producing, so we can safely try to make the first a producer -0087 if reducedModel.rev(involvedRxns(1)) == 1 -0088 [reducedModel,reversedRxns] = flipRxn(reducedModel, involvedRxns(1), groupIds, reversedRxns); -0089 else %ok, finally try to flip the second reaction -0090 if reducedModel.rev(involvedRxns(2)) == 1 -0091 [reducedModel,reversedRxns] = flipRxn(reducedModel, involvedRxns(2), groupIds, reversedRxns); -0092 involvedRxns = flip(involvedRxns);%make the second the first -0093 else -0094 error('We should never end up here!'); -0095 end -0096 end -0097 end -0098 end -0099 %Now, make sure the second rxn is a consumer -0100 if reducedModel.S(twoNonZero(i),involvedRxns(2)) > 0 -0101 if reducedModel.rev(involvedRxns(2)) == 1 -0102 [reducedModel,reversedRxns] = flipRxn(reducedModel, involvedRxns(2), groupIds, reversedRxns); -0103 else -0104 error('We should never end up here!'); -0105 end -0106 end -0107 -0108 %Calculate how many times the second reaction has to be -0109 %multiplied before being merged with the first -0110 stoichRatio=abs(reducedModel.S(twoNonZero(i),involvedRxns(1))/reducedModel.S(twoNonZero(i),involvedRxns(2))); +0065 %rev is compared rather than used as-is: the other operand is a +0066 %sparse logical, and MATLAB refuses to combine one with an integer +0067 %array ("Sparse integer array arithmetic operations are not +0068 %supported"). A model built outside RAVEN can carry an integer rev. +0069 pos = sum(reducedModel.S(twoNonZero(i),involvedRxns).' > 0 | reducedModel.rev(involvedRxns) ~= 0); +0070 neg = sum(reducedModel.S(twoNonZero(i),involvedRxns).' < 0 | reducedModel.rev(involvedRxns) ~= 0); +0071 +0072 +0073 %Check so that one or both of the reactions haven't been merged +0074 %already +0075 if numel(involvedRxns)==2 && isempty(intersect(bannedIndexes,involvedRxns)) && pos >= 1 && neg >= 1 +0076 %first, take care of a special case: If the first reaction is producing the metabolite and if it is reversible, +0077 %and the second is also producing it and is not reversible, change the order - the code below will not work otherwise +0078 if reducedModel.rev(involvedRxns(1)) && (~reducedModel.rev(involvedRxns(2))) && ... +0079 (reducedModel.S(twoNonZero(i),involvedRxns(1)) > 0) && (reducedModel.S(twoNonZero(i),involvedRxns(2)) > 0) +0080 involvedRxns = flip(involvedRxns); +0081 end +0082 +0083 %first make sure the first reaction is producing the metabolite +0084 if reducedModel.S(twoNonZero(i),involvedRxns(1)) < 0 +0085 %it is not producing the metabolite - fix that +0086 %first choice: use the second reaction as producer if it is producing +0087 if reducedModel.S(twoNonZero(i),involvedRxns(2)) > 0 +0088 involvedRxns = flip(involvedRxns);%make the second the first +0089 else +0090 %now we know that the second reaction is not producing, so we can safely try to make the first a producer +0091 if reducedModel.rev(involvedRxns(1)) == 1 +0092 [reducedModel,reversedRxns] = flipRxn(reducedModel, involvedRxns(1), groupIds, reversedRxns); +0093 else %ok, finally try to flip the second reaction +0094 if reducedModel.rev(involvedRxns(2)) == 1 +0095 [reducedModel,reversedRxns] = flipRxn(reducedModel, involvedRxns(2), groupIds, reversedRxns); +0096 involvedRxns = flip(involvedRxns);%make the second the first +0097 else +0098 error('We should never end up here!'); +0099 end +0100 end +0101 end +0102 end +0103 %Now, make sure the second rxn is a consumer +0104 if reducedModel.S(twoNonZero(i),involvedRxns(2)) > 0 +0105 if reducedModel.rev(involvedRxns(2)) == 1 +0106 [reducedModel,reversedRxns] = flipRxn(reducedModel, involvedRxns(2), groupIds, reversedRxns); +0107 else +0108 error('We should never end up here!'); +0109 end +0110 end 0111 -0112 %Add the second to the first -0113 reducedModel.S(:,involvedRxns(1))=reducedModel.S(:,involvedRxns(1))+reducedModel.S(:,involvedRxns(2))*stoichRatio; -0114 -0115 %Clear the second reaction -0116 reducedModel.S(:,involvedRxns(2))=0; -0117 -0118 %This is to prevent numerical issues. It should be 0 -0119 %already -0120 reducedModel.S(twoNonZero(i),involvedRxns(1))=0; +0112 %Calculate how many times the second reaction has to be +0113 %multiplied before being merged with the first +0114 stoichRatio=abs(reducedModel.S(twoNonZero(i),involvedRxns(1))/reducedModel.S(twoNonZero(i),involvedRxns(2))); +0115 +0116 %Add the second to the first +0117 reducedModel.S(:,involvedRxns(1))=reducedModel.S(:,involvedRxns(1))+reducedModel.S(:,involvedRxns(2))*stoichRatio; +0118 +0119 %Clear the second reaction +0120 reducedModel.S(:,involvedRxns(2))=0; 0121 -0122 %At this point the second reaction is certain to be deleted -0123 %in a later step and can therefore be ignored -0124 -0125 %Recalculate the bounds for the new reaction. This can be -0126 %problematic since the scale of the bounds may change -0127 %dramatically. Let the most constraining reaction determine -0128 %the new bound -0129 lb1=reducedModel.lb(involvedRxns(1)); -0130 lb2=reducedModel.lb(involvedRxns(2)); -0131 ub1=reducedModel.ub(involvedRxns(1)); -0132 ub2=reducedModel.ub(involvedRxns(2)); -0133 -0134 if lb2~=-inf -0135 reducedModel.lb(involvedRxns(1))=max(lb1,lb2/stoichRatio); -0136 end -0137 if ub2~=inf -0138 reducedModel.ub(involvedRxns(1))=min(ub1,ub2/stoichRatio); -0139 end -0140 -0141 %take care of the .rev flag - it could be that the combined rxn changes from rev to irrev -0142 reducedModel.rev(involvedRxns(1)) = reducedModel.rev(involvedRxns(1))*reducedModel.rev(involvedRxns(2));%this is a way to do an "and" operation with 0 and 1 numbers -0143 -0144 %Then recalculate the objective coefficient. The resulting -0145 %coefficient is the weighted sum of the previous -0146 reducedModel.c(involvedRxns(1))=reducedModel.c(involvedRxns(1))+reducedModel.c(involvedRxns(2))*stoichRatio; +0122 %This is to prevent numerical issues. It should be 0 +0123 %already +0124 reducedModel.S(twoNonZero(i),involvedRxns(1))=0; +0125 +0126 %At this point the second reaction is certain to be deleted +0127 %in a later step and can therefore be ignored +0128 +0129 %Recalculate the bounds for the new reaction. This can be +0130 %problematic since the scale of the bounds may change +0131 %dramatically. Let the most constraining reaction determine +0132 %the new bound +0133 lb1=reducedModel.lb(involvedRxns(1)); +0134 lb2=reducedModel.lb(involvedRxns(2)); +0135 ub1=reducedModel.ub(involvedRxns(1)); +0136 ub2=reducedModel.ub(involvedRxns(2)); +0137 +0138 if lb2~=-inf +0139 reducedModel.lb(involvedRxns(1))=max(lb1,lb2/stoichRatio); +0140 end +0141 if ub2~=inf +0142 reducedModel.ub(involvedRxns(1))=min(ub1,ub2/stoichRatio); +0143 end +0144 +0145 %take care of the .rev flag - it could be that the combined rxn changes from rev to irrev +0146 reducedModel.rev(involvedRxns(1)) = reducedModel.rev(involvedRxns(1))*reducedModel.rev(involvedRxns(2));%this is a way to do an "and" operation with 0 and 1 numbers 0147 -0148 %store which reactions that have been merged -0149 rxnInd1 = find(strcmp(origRxnIds, reducedModel.rxns(involvedRxns(1)))); -0150 rxnInd2 = find(strcmp(origRxnIds, reducedModel.rxns(involvedRxns(2)))); -0151 grpId = max(groupIds(rxnInd1),groupIds(rxnInd2)); -0152 if grpId == 0 -0153 grpId = nextGroupId; -0154 nextGroupId = nextGroupId + 1; -0155 end -0156 -0157 if groupIds(rxnInd1) ~= grpId -0158 if groupIds(rxnInd1) == 0 -0159 %not merged before, just set the group id -0160 groupIds(rxnInd1) = grpId; -0161 else -0162 %merged before - all rxns with the same group id should be changed -0163 groupIds(groupIds == groupIds(rxnInd1)) = grpId; -0164 end -0165 end -0166 if groupIds(rxnInd2) ~= grpId -0167 if groupIds(rxnInd2) == 0 -0168 %not merged before, just set the group id -0169 groupIds(rxnInd2) = grpId; -0170 else -0171 %merged before - all rxns with the same group id should be changed -0172 groupIds(groupIds == groupIds(rxnInd2)) = grpId; -0173 end -0174 end -0175 -0176 %Iterate again -0177 mergedSome=true; -0178 end -0179 end -0180 -0181 %All possible reactions merged -0182 if mergedSome==false -0183 break; -0184 end -0185 -0186 %Now delete all reactions that involve no metabolites -0187 I=find(sum(reducedModel.S~=0,1)==0); -0188 -0189 %Remove reactions -0190 reducedModel=removeReactions(reducedModel,I); -0191 -0192 %Remove metabolites -0193 notInUse=sum(reducedModel.S~=0,2)==0; -0194 reducedModel=removeMets(reducedModel,notInUse); -0195 end -0196 -0197 function [model1,reversedRxns1] = flipRxn(model1, rxnInd, groupIds1, reversedRxns1) -0198 model1.S(:,rxnInd) = model1.S(:,rxnInd)*-1; -0199 %swap the bounds -0200 ub = model1.ub(rxnInd); -0201 model1.ub(rxnInd) = -model1.lb(rxnInd); -0202 model1.lb(rxnInd) = -ub; -0203 %flip the objective -0204 model1.c(rxnInd) = -model1.c(rxnInd); -0205 -0206 %now take care of the reversedRxns - if this is a group, reverse all of the -0207 %reactions in the group in the reversedRxns index - they will all be reversed at the -0208 %same time since they are the same rxn. -0209 rxnIndices = rxnInd; -0210 if groupIds1(rxnInd) > 0 -0211 rxnIndices = find(groupIds1 == groupIds1(rxnInd)); -0212 end -0213 reversedRxns1(rxnIndices) = ~reversedRxns1(rxnIndices); -0214 end -0215 end +0148 %Then recalculate the objective coefficient. The resulting +0149 %coefficient is the weighted sum of the previous +0150 reducedModel.c(involvedRxns(1))=reducedModel.c(involvedRxns(1))+reducedModel.c(involvedRxns(2))*stoichRatio; +0151 +0152 %store which reactions that have been merged +0153 rxnInd1 = find(strcmp(origRxnIds, reducedModel.rxns(involvedRxns(1)))); +0154 rxnInd2 = find(strcmp(origRxnIds, reducedModel.rxns(involvedRxns(2)))); +0155 grpId = max(groupIds(rxnInd1),groupIds(rxnInd2)); +0156 if grpId == 0 +0157 grpId = nextGroupId; +0158 nextGroupId = nextGroupId + 1; +0159 end +0160 +0161 if groupIds(rxnInd1) ~= grpId +0162 if groupIds(rxnInd1) == 0 +0163 %not merged before, just set the group id +0164 groupIds(rxnInd1) = grpId; +0165 else +0166 %merged before - all rxns with the same group id should be changed +0167 groupIds(groupIds == groupIds(rxnInd1)) = grpId; +0168 end +0169 end +0170 if groupIds(rxnInd2) ~= grpId +0171 if groupIds(rxnInd2) == 0 +0172 %not merged before, just set the group id +0173 groupIds(rxnInd2) = grpId; +0174 else +0175 %merged before - all rxns with the same group id should be changed +0176 groupIds(groupIds == groupIds(rxnInd2)) = grpId; +0177 end +0178 end +0179 +0180 %Iterate again +0181 mergedSome=true; +0182 end +0183 end +0184 +0185 %All possible reactions merged +0186 if mergedSome==false +0187 break; +0188 end +0189 +0190 %Now delete all reactions that involve no metabolites +0191 I=find(sum(reducedModel.S~=0,1)==0); +0192 +0193 %Remove reactions +0194 reducedModel=removeReactions(reducedModel,I); +0195 +0196 %Remove metabolites +0197 notInUse=sum(reducedModel.S~=0,2)==0; +0198 reducedModel=removeMets(reducedModel,notInUse); +0199 end +0200 +0201 function [model1,reversedRxns1] = flipRxn(model1, rxnInd, groupIds1, reversedRxns1) +0202 model1.S(:,rxnInd) = model1.S(:,rxnInd)*-1; +0203 %swap the bounds +0204 ub = model1.ub(rxnInd); +0205 model1.ub(rxnInd) = -model1.lb(rxnInd); +0206 model1.lb(rxnInd) = -ub; +0207 %flip the objective +0208 model1.c(rxnInd) = -model1.c(rxnInd); +0209 +0210 %now take care of the reversedRxns - if this is a group, reverse all of the +0211 %reactions in the group in the reversedRxns index - they will all be reversed at the +0212 %same time since they are the same rxn. +0213 rxnIndices = rxnInd; +0214 if groupIds1(rxnInd) > 0 +0215 rxnIndices = find(groupIds1 == groupIds1(rxnInd)); +0216 end +0217 reversedRxns1(rxnIndices) = ~reversedRxns1(rxnIndices); +0218 end +0219 end
    Generated by m2html © 2005
    \ No newline at end of file diff --git a/doc/io/importModel.html b/doc/io/importModel.html index dff37de8c..dc124f99e 100644 --- a/doc/io/importModel.html +++ b/doc/io/importModel.html @@ -849,415 +849,426 @@

    SOURCE CODE ^end 0740 end 0741 if isfield(modelSBML,'notes') -0742 startString=strfind(modelSBML.notes,'xhtml">'); -0743 endString=strfind(modelSBML.notes,'</body>'); -0744 if any(startString) && any(endString) -0745 model.annotation.note=modelSBML.notes(startString(1)+7:endString-1); -0746 model.annotation.note=regexprep(model.annotation.note,'<p.*?>|</p.*?>',''); -0747 model.annotation.note=strtrim(model.annotation.note); -0748 if regexp(model.annotation.note,'This file was generated using the exportModel function in RAVEN Toolbox \d\.\d and OutputSBML in libSBML') -0749 model.annotation=rmfield(model.annotation,'note'); % Default note added when running exportModel -0750 end -0751 end -0752 end -0753 -0754 if any(~cellfun(@isempty,compartmentOutside)) -0755 model.compOutside=compartmentOutside; -0756 end -0757 -0758 model.rxnNames=reactionNames; -0759 model.metNames=metaboliteNames; -0760 -0761 %Match the compartments for metabolites -0762 [~, J]=ismember(metaboliteCompartments,model.comps); -0763 model.metComps=J; +0742 %Two shapes reach here. RAVEN's own exportModel writes the note as bare +0743 %text inside a <body>; a writer that serialises a notes dictionary emits +0744 %one '<p>key: value</p>' per entry and no <body> at all, carrying the note +0745 %under its own key. Take the labelled entry when there is one, so the +0746 %other entries (version, and whatever else the dictionary held) stay out +0747 %of the note, and fall back to the <body> text otherwise. +0748 note=parseNote(modelSBML.notes,'note'); +0749 if isempty(note) +0750 startString=strfind(modelSBML.notes,'xhtml">'); +0751 endString=strfind(modelSBML.notes,'</body>'); +0752 if any(startString) && any(endString) +0753 note=modelSBML.notes(startString(1)+7:endString-1); +0754 end +0755 end +0756 if ~isempty(note) +0757 model.annotation.note=regexprep(note,'<p.*?>|</p.*?>',''); +0758 model.annotation.note=strtrim(model.annotation.note); +0759 if regexp(model.annotation.note,'This file was generated using the exportModel function in RAVEN Toolbox \d\.\d and OutputSBML in libSBML') +0760 model.annotation=rmfield(model.annotation,'note'); % Default note added when running exportModel +0761 end +0762 end +0763 end 0764 -0765 %If any genes have been loaded (only for the new format) -0766 if ~isempty(geneNames) -0767 %In some rare cases geneNames may not necessarily be used in grRules. -0768 %That is true for Yeast 7.6. It's therefore important to change gene -0769 %systematic names to geneIDs in sophisticated way. Gene systematic -0770 %names are not unique, since exactly the same name may be in different -0771 %compartments -0772 if all(cellfun(@isempty,strfind(grRules,geneNames{1}))) -0773 geneShortNames=geneNames; -0774 %geneShortNames contain compartments as well, so these are removed -0775 geneShortNames=regexprep(geneShortNames,' \[.+$',''); -0776 %grRules obtained from modifier fields contain geneNames. These are -0777 %changed into geneIDs. grRulesFromModifier is a good way to have -0778 %geneIDs and rxns association when it's important to resolve -0779 %systematic name ambiguities -0780 grRulesFromModifier=regexprep(regexprep(grRulesFromModifier,'\[|\]','_'),regexprep(geneNames,'\[|\]','_'),geneIDs); -0781 grRules=regexprep(regexprep(grRules,'\[|\]','_'),regexprep(geneNames,'\[|\]','_'),geneIDs); -0782 -0783 %Yeast 7.6 contains several metabolites, which were used in gene -0784 %associations. For that reason, the list of species ID is created -0785 %and we then check whether any of them have kegg.genes annotation -0786 %thereby obtaining systematic gene names -0787 geneShortNames=vertcat(geneShortNames,metaboliteNames); -0788 geneIDs=vertcat(geneIDs,metaboliteIDs); -0789 geneSystNames=extractMiriam(vertcat(geneMiriams,metaboliteMiriams),'kegg.genes'); -0790 geneCompartments=vertcat(geneCompartments,metaboliteCompartments); -0791 geneMiriams=vertcat(geneMiriams,metaboliteMiriams); -0792 -0793 %Now we retain information for only these entries, which have -0794 %kegg.genes annotation -0795 geneShortNames=geneShortNames(~cellfun('isempty',geneSystNames)); -0796 geneIDs=geneIDs(~cellfun('isempty',geneSystNames)); -0797 geneSystNames=geneSystNames(~cellfun('isempty',geneSystNames)); -0798 geneCompartments=geneCompartments(~cellfun('isempty',geneSystNames)); -0799 geneMiriams=geneMiriams(~cellfun('isempty',geneSystNames)); -0800 %Now we reorder geneIDs and geneSystNames by geneSystNames string -0801 %length -0802 geneNames=geneIDs;%Backuping geneIDs, since we need unsorted order for later -0803 [~, Indx] = sort(cellfun('size', geneSystNames, 2), 'descend'); -0804 geneIDs = geneIDs(Indx); -0805 geneSystNames = geneSystNames(Indx); -0806 for i=1:numel(geneSystNames) -0807 for j=1:numel(grRules) -0808 if strfind(grRules{j},geneSystNames{i}) -0809 if ~isempty(grRules{j}) -0810 if sum(ismember(geneSystNames,geneSystNames{i}))==1 -0811 grRules{j}=regexprep(grRules{j},geneSystNames{i},geneIDs{i}); -0812 elseif sum(ismember(geneSystNames,geneSystNames{i}))>1 -0813 counter=0; -0814 ovrlpIDs=geneIDs(ismember(geneSystNames,geneSystNames{i})); -0815 for k=1:numel(ovrlpIDs) -0816 if strfind(grRulesFromModifier{j},ovrlpIDs{k}) -0817 counter=counter+1; -0818 grRules{j}=regexprep(grRules{j},geneSystNames{i},ovrlpIDs{k}); -0819 end -0820 if counter>1 -0821 EM=['Gene association is ambiguous for reaction ' modelSBML.reaction(j).id]; -0822 dispEM(EM); -0823 end -0824 end -0825 end -0826 end -0827 end -0828 end -0829 end -0830 end -0831 model.genes=geneNames; -0832 model.grRules=grRules; -0833 [grRules,rxnGeneMat] = standardizeGrRules(model,true); -0834 model.grRules = grRules; -0835 model.rxnGeneMat = rxnGeneMat; -0836 -0837 %Match the compartments for genes -0838 [~, J]=ismember(geneCompartments,model.comps); -0839 model.geneComps=J; -0840 else -0841 if ~all(cellfun(@isempty,grRules)) -0842 %If fbc_geneProduct exists, follow the specified gene order, such -0843 %that matching geneShortNames in function below will work -0844 if isfield(modelSBML,'fbc_geneProduct') -0845 genes={modelSBML.fbc_geneProduct.fbc_id}; -0846 -0847 %Get gene Miriams if they were not retrieved above (this occurs -0848 %when genes are stored as fbc_geneProduct instead of species) -0849 if isempty(geneMiriams) -0850 geneMiriams = cell(numel(genes),1); -0851 if isfield(modelSBML.fbc_geneProduct,'sboTerm') && isscalar(unique([modelSBML.fbc_geneProduct.sboTerm])) -0852 %If all the SBO terms are identical, don't add them to geneMiriams -0853 modelSBML.fbc_geneProduct = rmfield(modelSBML.fbc_geneProduct,'sboTerm'); -0854 end -0855 for i = 1:numel(genes) -0856 geneMiriams{i}=parseMiriam(modelSBML.fbc_geneProduct(i).annotation); -0857 if isfield(modelSBML.fbc_geneProduct(i),'sboTerm') && ~(modelSBML.fbc_geneProduct(i).sboTerm==-1) -0858 geneMiriams{i} = addSBOtoMiriam(geneMiriams{i},modelSBML.fbc_geneProduct(i).sboTerm); -0859 end -0860 end -0861 end -0862 proteins={modelSBML.fbc_geneProduct.fbc_name}; -0863 else -0864 genes=getGenesFromGrRules(grRules); -0865 end -0866 model.genes=genes; -0867 model.grRules=grRules; -0868 [grRules,rxnGeneMat] = standardizeGrRules(model,true); -0869 model.grRules = grRules; -0870 model.rxnGeneMat = rxnGeneMat; -0871 end -0872 end -0873 -0874 if all(cellfun(@isempty,geneShortNames)) -0875 if isfield(modelSBML,'fbc_geneProduct') -0876 for i=1:numel(genes) -0877 if ~isempty(modelSBML.fbc_geneProduct(i).fbc_label) -0878 geneShortNames{i,1}=modelSBML.fbc_geneProduct(i).fbc_label; -0879 elseif ~isempty(modelSBML.fbc_geneProduct(i).fbc_name) -0880 geneShortNames{i,1}=modelSBML.fbc_geneProduct(i).fbc_name; -0881 else -0882 geneShortNames{i,1}=''; -0883 end -0884 end -0885 end -0886 end -0887 -0888 %If any InChIs have been loaded -0889 if any(~cellfun(@isempty,metaboliteInChI)) -0890 model.inchis=metaboliteInChI; -0891 end -0892 -0893 %If any formulas have been loaded -0894 if any(~cellfun(@isempty,metaboliteFormula)) -0895 model.metFormulas=metaboliteFormula; -0896 end -0897 -0898 %If any charges have been loaded -0899 if ~isempty(metaboliteCharges) -0900 model.metCharges=metaboliteCharges; -0901 end -0902 -0903 %If any gene short names have been loaded -0904 if any(~cellfun(@isempty,geneShortNames)) -0905 model.geneShortNames=geneShortNames; -0906 end -0907 -0908 %If any Miriam strings for compartments have been loaded -0909 if any(~cellfun(@isempty,compartmentMiriams)) -0910 model.compMiriams=compartmentMiriams; -0911 end -0912 -0913 %If any Miriam strings for metabolites have been loaded -0914 if any(~cellfun(@isempty,metaboliteMiriams)) -0915 model.metMiriams=metaboliteMiriams; -0916 end -0917 -0918 %If any subsystems have been loaded -0919 if any(~cellfun(@isempty,subsystems)) -0920 model.subSystems=subsystems; -0921 end -0922 if any(rxnComps) -0923 if all(rxnComps) -0924 model.rxnComps=rxnComps; -0925 else -0926 if supressWarnings==false -0927 EM='The compartments for the following reactions could not be matched. Ignoring reaction compartment information'; -0928 dispEM(EM,false,model.rxns(rxnComps==0)); -0929 end -0930 end -0931 end -0932 -0933 %If any ec-codes have been loaded -0934 if any(~cellfun(@isempty,eccodes)) -0935 model.eccodes=eccodes; -0936 end -0937 -0938 %If any Miriam strings for reactions have been loaded -0939 if any(~cellfun(@isempty,rxnMiriams)) -0940 model.rxnMiriams=rxnMiriams; -0941 end -0942 -0943 %If any Miriam strings for genes have been loaded -0944 if any(~cellfun(@isempty,geneMiriams)) -0945 model.geneMiriams=geneMiriams; -0946 end -0947 -0948 %If any protein strings have been loaded -0949 if any(~cellfun(@isempty,proteins)) -0950 proteins = reshape(proteins,[],1); -0951 model.proteins=proteins; +0765 if any(~cellfun(@isempty,compartmentOutside)) +0766 model.compOutside=compartmentOutside; +0767 end +0768 +0769 model.rxnNames=reactionNames; +0770 model.metNames=metaboliteNames; +0771 +0772 %Match the compartments for metabolites +0773 [~, J]=ismember(metaboliteCompartments,model.comps); +0774 model.metComps=J; +0775 +0776 %If any genes have been loaded (only for the new format) +0777 if ~isempty(geneNames) +0778 %In some rare cases geneNames may not necessarily be used in grRules. +0779 %That is true for Yeast 7.6. It's therefore important to change gene +0780 %systematic names to geneIDs in sophisticated way. Gene systematic +0781 %names are not unique, since exactly the same name may be in different +0782 %compartments +0783 if all(cellfun(@isempty,strfind(grRules,geneNames{1}))) +0784 geneShortNames=geneNames; +0785 %geneShortNames contain compartments as well, so these are removed +0786 geneShortNames=regexprep(geneShortNames,' \[.+$',''); +0787 %grRules obtained from modifier fields contain geneNames. These are +0788 %changed into geneIDs. grRulesFromModifier is a good way to have +0789 %geneIDs and rxns association when it's important to resolve +0790 %systematic name ambiguities +0791 grRulesFromModifier=regexprep(regexprep(grRulesFromModifier,'\[|\]','_'),regexprep(geneNames,'\[|\]','_'),geneIDs); +0792 grRules=regexprep(regexprep(grRules,'\[|\]','_'),regexprep(geneNames,'\[|\]','_'),geneIDs); +0793 +0794 %Yeast 7.6 contains several metabolites, which were used in gene +0795 %associations. For that reason, the list of species ID is created +0796 %and we then check whether any of them have kegg.genes annotation +0797 %thereby obtaining systematic gene names +0798 geneShortNames=vertcat(geneShortNames,metaboliteNames); +0799 geneIDs=vertcat(geneIDs,metaboliteIDs); +0800 geneSystNames=extractMiriam(vertcat(geneMiriams,metaboliteMiriams),'kegg.genes'); +0801 geneCompartments=vertcat(geneCompartments,metaboliteCompartments); +0802 geneMiriams=vertcat(geneMiriams,metaboliteMiriams); +0803 +0804 %Now we retain information for only these entries, which have +0805 %kegg.genes annotation +0806 geneShortNames=geneShortNames(~cellfun('isempty',geneSystNames)); +0807 geneIDs=geneIDs(~cellfun('isempty',geneSystNames)); +0808 geneSystNames=geneSystNames(~cellfun('isempty',geneSystNames)); +0809 geneCompartments=geneCompartments(~cellfun('isempty',geneSystNames)); +0810 geneMiriams=geneMiriams(~cellfun('isempty',geneSystNames)); +0811 %Now we reorder geneIDs and geneSystNames by geneSystNames string +0812 %length +0813 geneNames=geneIDs;%Backuping geneIDs, since we need unsorted order for later +0814 [~, Indx] = sort(cellfun('size', geneSystNames, 2), 'descend'); +0815 geneIDs = geneIDs(Indx); +0816 geneSystNames = geneSystNames(Indx); +0817 for i=1:numel(geneSystNames) +0818 for j=1:numel(grRules) +0819 if strfind(grRules{j},geneSystNames{i}) +0820 if ~isempty(grRules{j}) +0821 if sum(ismember(geneSystNames,geneSystNames{i}))==1 +0822 grRules{j}=regexprep(grRules{j},geneSystNames{i},geneIDs{i}); +0823 elseif sum(ismember(geneSystNames,geneSystNames{i}))>1 +0824 counter=0; +0825 ovrlpIDs=geneIDs(ismember(geneSystNames,geneSystNames{i})); +0826 for k=1:numel(ovrlpIDs) +0827 if strfind(grRulesFromModifier{j},ovrlpIDs{k}) +0828 counter=counter+1; +0829 grRules{j}=regexprep(grRules{j},geneSystNames{i},ovrlpIDs{k}); +0830 end +0831 if counter>1 +0832 EM=['Gene association is ambiguous for reaction ' modelSBML.reaction(j).id]; +0833 dispEM(EM); +0834 end +0835 end +0836 end +0837 end +0838 end +0839 end +0840 end +0841 end +0842 model.genes=geneNames; +0843 model.grRules=grRules; +0844 [grRules,rxnGeneMat] = standardizeGrRules(model,true); +0845 model.grRules = grRules; +0846 model.rxnGeneMat = rxnGeneMat; +0847 +0848 %Match the compartments for genes +0849 [~, J]=ismember(geneCompartments,model.comps); +0850 model.geneComps=J; +0851 else +0852 if ~all(cellfun(@isempty,grRules)) +0853 %If fbc_geneProduct exists, follow the specified gene order, such +0854 %that matching geneShortNames in function below will work +0855 if isfield(modelSBML,'fbc_geneProduct') +0856 genes={modelSBML.fbc_geneProduct.fbc_id}; +0857 +0858 %Get gene Miriams if they were not retrieved above (this occurs +0859 %when genes are stored as fbc_geneProduct instead of species) +0860 if isempty(geneMiriams) +0861 geneMiriams = cell(numel(genes),1); +0862 if isfield(modelSBML.fbc_geneProduct,'sboTerm') && isscalar(unique([modelSBML.fbc_geneProduct.sboTerm])) +0863 %If all the SBO terms are identical, don't add them to geneMiriams +0864 modelSBML.fbc_geneProduct = rmfield(modelSBML.fbc_geneProduct,'sboTerm'); +0865 end +0866 for i = 1:numel(genes) +0867 geneMiriams{i}=parseMiriam(modelSBML.fbc_geneProduct(i).annotation); +0868 if isfield(modelSBML.fbc_geneProduct(i),'sboTerm') && ~(modelSBML.fbc_geneProduct(i).sboTerm==-1) +0869 geneMiriams{i} = addSBOtoMiriam(geneMiriams{i},modelSBML.fbc_geneProduct(i).sboTerm); +0870 end +0871 end +0872 end +0873 proteins={modelSBML.fbc_geneProduct.fbc_name}; +0874 else +0875 genes=getGenesFromGrRules(grRules); +0876 end +0877 model.genes=genes; +0878 model.grRules=grRules; +0879 [grRules,rxnGeneMat] = standardizeGrRules(model,true); +0880 model.grRules = grRules; +0881 model.rxnGeneMat = rxnGeneMat; +0882 end +0883 end +0884 +0885 if all(cellfun(@isempty,geneShortNames)) +0886 if isfield(modelSBML,'fbc_geneProduct') +0887 for i=1:numel(genes) +0888 if ~isempty(modelSBML.fbc_geneProduct(i).fbc_label) +0889 geneShortNames{i,1}=modelSBML.fbc_geneProduct(i).fbc_label; +0890 elseif ~isempty(modelSBML.fbc_geneProduct(i).fbc_name) +0891 geneShortNames{i,1}=modelSBML.fbc_geneProduct(i).fbc_name; +0892 else +0893 geneShortNames{i,1}=''; +0894 end +0895 end +0896 end +0897 end +0898 +0899 %If any InChIs have been loaded +0900 if any(~cellfun(@isempty,metaboliteInChI)) +0901 model.inchis=metaboliteInChI; +0902 end +0903 +0904 %If any formulas have been loaded +0905 if any(~cellfun(@isempty,metaboliteFormula)) +0906 model.metFormulas=metaboliteFormula; +0907 end +0908 +0909 %If any charges have been loaded +0910 if ~isempty(metaboliteCharges) +0911 model.metCharges=metaboliteCharges; +0912 end +0913 +0914 %If any gene short names have been loaded +0915 if any(~cellfun(@isempty,geneShortNames)) +0916 model.geneShortNames=geneShortNames; +0917 end +0918 +0919 %If any Miriam strings for compartments have been loaded +0920 if any(~cellfun(@isempty,compartmentMiriams)) +0921 model.compMiriams=compartmentMiriams; +0922 end +0923 +0924 %If any Miriam strings for metabolites have been loaded +0925 if any(~cellfun(@isempty,metaboliteMiriams)) +0926 model.metMiriams=metaboliteMiriams; +0927 end +0928 +0929 %If any subsystems have been loaded +0930 if any(~cellfun(@isempty,subsystems)) +0931 model.subSystems=subsystems; +0932 end +0933 if any(rxnComps) +0934 if all(rxnComps) +0935 model.rxnComps=rxnComps; +0936 else +0937 if supressWarnings==false +0938 EM='The compartments for the following reactions could not be matched. Ignoring reaction compartment information'; +0939 dispEM(EM,false,model.rxns(rxnComps==0)); +0940 end +0941 end +0942 end +0943 +0944 %If any ec-codes have been loaded +0945 if any(~cellfun(@isempty,eccodes)) +0946 model.eccodes=eccodes; +0947 end +0948 +0949 %If any Miriam strings for reactions have been loaded +0950 if any(~cellfun(@isempty,rxnMiriams)) +0951 model.rxnMiriams=rxnMiriams; 0952 end 0953 -0954 model.unconstrained=metaboliteUnconstrained; -0955 -0956 %Convert SBML IDs back into their original strings. Here we are using part -0957 %from convertSBMLID, originating from the COBRA Toolbox -0958 model.rxns=regexprep(model.rxns,'__([0-9]+)__','${char(str2num($1))}'); -0959 model.mets=regexprep(model.mets,'__([0-9]+)__','${char(str2num($1))}'); -0960 model.comps=regexprep(model.comps,'__([0-9]+)__','${char(str2num($1))}'); -0961 model.grRules=regexprep(model.grRules,'__([0-9]+)__','${char(str2num($1))}'); -0962 model.genes=regexprep(model.genes,'__([0-9]+)__','${char(str2num($1))}'); -0963 model.id=regexprep(model.id,'__([0-9]+)__','${char(str2num($1))}'); +0954 %If any Miriam strings for genes have been loaded +0955 if any(~cellfun(@isempty,geneMiriams)) +0956 model.geneMiriams=geneMiriams; +0957 end +0958 +0959 %If any protein strings have been loaded +0960 if any(~cellfun(@isempty,proteins)) +0961 proteins = reshape(proteins,[],1); +0962 model.proteins=proteins; +0963 end 0964 -0965 if removePrefix -0966 [model, hasChanged]=removeIdentifierPrefix(model); -0967 dispEM(['The following fields have prefixes removed from all entries. '... -0968 'If this is undesired, run importModel with removePrefix as false. Example: '... -0969 'importModel(''filename.xml'',[],false);'],false,hasChanged) -0970 end -0971 -0972 %Remove unused fields -0973 if isempty(model.annotation) -0974 model=rmfield(model,'annotation'); -0975 end -0976 if isempty(model.compOutside) -0977 model=rmfield(model,'compOutside'); -0978 end -0979 if isempty(model.compMiriams) -0980 model=rmfield(model,'compMiriams'); +0965 model.unconstrained=metaboliteUnconstrained; +0966 +0967 %Convert SBML IDs back into their original strings. Here we are using part +0968 %from convertSBMLID, originating from the COBRA Toolbox +0969 model.rxns=regexprep(model.rxns,'__([0-9]+)__','${char(str2num($1))}'); +0970 model.mets=regexprep(model.mets,'__([0-9]+)__','${char(str2num($1))}'); +0971 model.comps=regexprep(model.comps,'__([0-9]+)__','${char(str2num($1))}'); +0972 model.grRules=regexprep(model.grRules,'__([0-9]+)__','${char(str2num($1))}'); +0973 model.genes=regexprep(model.genes,'__([0-9]+)__','${char(str2num($1))}'); +0974 model.id=regexprep(model.id,'__([0-9]+)__','${char(str2num($1))}'); +0975 +0976 if removePrefix +0977 [model, hasChanged]=removeIdentifierPrefix(model); +0978 dispEM(['The following fields have prefixes removed from all entries. '... +0979 'If this is undesired, run importModel with removePrefix as false. Example: '... +0980 'importModel(''filename.xml'',[],false);'],false,hasChanged) 0981 end -0982 if isempty(model.rxnComps) -0983 model=rmfield(model,'rxnComps'); -0984 end -0985 if isempty(model.grRules) -0986 model=rmfield(model,'grRules'); -0987 end -0988 if isempty(model.rxnGeneMat) -0989 model=rmfield(model,'rxnGeneMat'); -0990 end -0991 if isempty(model.subSystems) -0992 model=rmfield(model,'subSystems'); -0993 else -0994 model.subSystems(cellfun(@isempty,subsystems))={{''}}; -0995 % If all subSystems have single entries, then unnest them -0996 if all(cellfun(@(x) iscell(x) && isscalar(x), model.subSystems)) -0997 model.subSystems = transpose([model.subSystems{:}]); -0998 end -0999 end -1000 if isempty(model.eccodes) -1001 model=rmfield(model,'eccodes'); -1002 end -1003 if isempty(model.rxnMiriams) -1004 model=rmfield(model,'rxnMiriams'); -1005 end -1006 if cellfun(@isempty,model.rxnNotes) -1007 model=rmfield(model,'rxnNotes'); -1008 end -1009 if cellfun(@isempty,model.rxnReferences) -1010 model=rmfield(model,'rxnReferences'); -1011 end -1012 if isempty(model.rxnConfidenceScores) || all(isnan(model.rxnConfidenceScores)) -1013 model=rmfield(model,'rxnConfidenceScores'); -1014 end -1015 if isempty(model.genes) -1016 model=rmfield(model,'genes'); -1017 elseif isrow(model.genes) -1018 model.genes=transpose(model.genes); +0982 +0983 %Remove unused fields +0984 if isempty(model.annotation) +0985 model=rmfield(model,'annotation'); +0986 end +0987 if isempty(model.compOutside) +0988 model=rmfield(model,'compOutside'); +0989 end +0990 if isempty(model.compMiriams) +0991 model=rmfield(model,'compMiriams'); +0992 end +0993 if isempty(model.rxnComps) +0994 model=rmfield(model,'rxnComps'); +0995 end +0996 if isempty(model.grRules) +0997 model=rmfield(model,'grRules'); +0998 end +0999 if isempty(model.rxnGeneMat) +1000 model=rmfield(model,'rxnGeneMat'); +1001 end +1002 if isempty(model.subSystems) +1003 model=rmfield(model,'subSystems'); +1004 else +1005 model.subSystems(cellfun(@isempty,subsystems))={{''}}; +1006 % If all subSystems have single entries, then unnest them +1007 if all(cellfun(@(x) iscell(x) && isscalar(x), model.subSystems)) +1008 model.subSystems = transpose([model.subSystems{:}]); +1009 end +1010 end +1011 if isempty(model.eccodes) +1012 model=rmfield(model,'eccodes'); +1013 end +1014 if isempty(model.rxnMiriams) +1015 model=rmfield(model,'rxnMiriams'); +1016 end +1017 if cellfun(@isempty,model.rxnNotes) +1018 model=rmfield(model,'rxnNotes'); 1019 end -1020 if isempty(model.geneComps) -1021 model=rmfield(model,'geneComps'); +1020 if cellfun(@isempty,model.rxnReferences) +1021 model=rmfield(model,'rxnReferences'); 1022 end -1023 if isempty(model.geneMiriams) -1024 model=rmfield(model,'geneMiriams'); +1023 if isempty(model.rxnConfidenceScores) || all(isnan(model.rxnConfidenceScores)) +1024 model=rmfield(model,'rxnConfidenceScores'); 1025 end -1026 if isempty(model.geneShortNames) -1027 model=rmfield(model,'geneShortNames'); -1028 end -1029 if isempty(model.proteins) -1030 model=rmfield(model,'proteins'); -1031 end -1032 if isempty(model.inchis) -1033 model=rmfield(model,'inchis'); -1034 end -1035 if isempty(model.metFormulas) -1036 model=rmfield(model,'metFormulas'); -1037 end -1038 if isempty(model.metMiriams) -1039 model=rmfield(model,'metMiriams'); -1040 end -1041 if ~any(model.metCharges) -1042 model=rmfield(model,'metCharges'); -1043 end -1044 -1045 %This just removes the grRules if no genes have been loaded -1046 if ~isfield(model,'genes') && isfield(model,'grRules') -1047 model=rmfield(model,'grRules'); +1026 if isempty(model.genes) +1027 model=rmfield(model,'genes'); +1028 elseif isrow(model.genes) +1029 model.genes=transpose(model.genes); +1030 end +1031 if isempty(model.geneComps) +1032 model=rmfield(model,'geneComps'); +1033 end +1034 if isempty(model.geneMiriams) +1035 model=rmfield(model,'geneMiriams'); +1036 end +1037 if isempty(model.geneShortNames) +1038 model=rmfield(model,'geneShortNames'); +1039 end +1040 if isempty(model.proteins) +1041 model=rmfield(model,'proteins'); +1042 end +1043 if isempty(model.inchis) +1044 model=rmfield(model,'inchis'); +1045 end +1046 if isempty(model.metFormulas) +1047 model=rmfield(model,'metFormulas'); 1048 end -1049 -1050 %Print warnings about bad structure -1051 if supressWarnings==false -1052 checkModelStruct(model,false); -1053 end -1054 -1055 if removeExcMets==true -1056 model=simplifyModel(model); -1057 end -1058 end -1059 -1060 function fieldContent=parseNote(searchString,fieldName) -1061 %The function obtains the particular information from 'notes' field, using -1062 %fieldName as the dummy string -1063 -1064 fieldContent=''; +1049 if isempty(model.metMiriams) +1050 model=rmfield(model,'metMiriams'); +1051 end +1052 if ~any(model.metCharges) +1053 model=rmfield(model,'metCharges'); +1054 end +1055 +1056 %This just removes the grRules if no genes have been loaded +1057 if ~isfield(model,'genes') && isfield(model,'grRules') +1058 model=rmfield(model,'grRules'); +1059 end +1060 +1061 %Print warnings about bad structure +1062 if supressWarnings==false +1063 checkModelStruct(model,false); +1064 end 1065 -1066 if strfind(searchString,fieldName) -1067 [~,targetString] = regexp(searchString,['<p>' fieldName '.*?</p>'],'tokens','match'); -1068 targetString=regexprep(targetString,'<p>|</p>',''); -1069 targetString=regexprep(targetString,[fieldName, ':'],''); -1070 for i=1:numel(targetString) -1071 fieldContent=[fieldContent ';' strtrim(targetString{1,i})]; -1072 end -1073 fieldContent=regexprep(fieldContent,'^;|;$',''); -1074 else -1075 fieldContent=''; -1076 end -1077 end -1078 -1079 function fieldContent=parseAnnotation(searchString,startString,midString,fieldName) -1080 -1081 fieldContent=''; -1082 -1083 %Removing whitespace characters from the ending strings, which may occur in -1084 %several cases -1085 searchString=regexprep(searchString,'" />','"/>'); -1086 [~,targetString] = regexp(searchString,['<rdf:li rdf:resource="' startString fieldName midString '.*?"/>'],'tokens','match'); -1087 targetString=regexprep(targetString,'<rdf:li rdf:resource="|"/>',''); -1088 targetString=regexprep(targetString,startString,''); -1089 targetString=regexprep(targetString,[fieldName midString],''); -1090 -1091 for i=1:numel(targetString) -1092 fieldContent=[fieldContent ';' strtrim(targetString{1,i})]; -1093 end -1094 -1095 fieldContent=regexprep(fieldContent,'^;|;$',''); -1096 end -1097 -1098 function miriamStruct=parseMiriam(searchString) -1099 %Generates miriam structure from annotation field -1100 -1101 %Finding whether miriams are written in the old or the new way -1102 if strfind(searchString,'urn:miriam:') -1103 startString='urn:miriam:'; -1104 elseif strfind(searchString,'http://identifiers.org/') -1105 startString='http://identifiers.org/'; -1106 elseif strfind(searchString,'https://identifiers.org/') -1107 startString='https://identifiers.org/'; -1108 else -1109 miriamStruct=[]; -1110 return; -1111 end -1112 -1113 miriamStruct=[]; -1114 -1115 searchString=regexprep(searchString,'" />','"/>'); -1116 [~,targetString] = regexp(searchString,'<rdf:li rdf:resource=".*?"/>','tokens','match'); -1117 targetString=regexprep(targetString,'<rdf:li rdf:resource="|"/>',''); -1118 targetString=regexprep(targetString,startString,''); -1119 -1120 fwdslash = contains(targetString,'/'); -1121 midString = cell(numel(targetString),1); -1122 midString(fwdslash) = {'/'}; -1123 midString(~fwdslash) = {':'}; -1124 -1125 counter=0; -1126 for i=1:numel(targetString) -1127 if isempty(regexp(targetString{1,i},'inchi|ec-code', 'once')) -1128 counter=counter+1; -1129 miriamStruct.name{counter,1} = regexprep(targetString{1,i},[midString{i} '.+'],'','once'); -1130 miriamStruct.value{counter,1} = regexprep(targetString{1,i},[miriamStruct.name{counter,1} midString{i}],'','once'); -1131 miriamStruct.name{counter,1} = regexprep(miriamStruct.name{counter,1},'^obo\.',''); -1132 end -1133 end -1134 end +1066 if removeExcMets==true +1067 model=simplifyModel(model); +1068 end +1069 end +1070 +1071 function fieldContent=parseNote(searchString,fieldName) +1072 %The function obtains the particular information from 'notes' field, using +1073 %fieldName as the dummy string +1074 +1075 fieldContent=''; +1076 +1077 if strfind(searchString,fieldName) +1078 [~,targetString] = regexp(searchString,['<p>' fieldName '.*?</p>'],'tokens','match'); +1079 targetString=regexprep(targetString,'<p>|</p>',''); +1080 targetString=regexprep(targetString,[fieldName, ':'],''); +1081 for i=1:numel(targetString) +1082 fieldContent=[fieldContent ';' strtrim(targetString{1,i})]; +1083 end +1084 fieldContent=regexprep(fieldContent,'^;|;$',''); +1085 else +1086 fieldContent=''; +1087 end +1088 end +1089 +1090 function fieldContent=parseAnnotation(searchString,startString,midString,fieldName) +1091 +1092 fieldContent=''; +1093 +1094 %Removing whitespace characters from the ending strings, which may occur in +1095 %several cases +1096 searchString=regexprep(searchString,'" />','"/>'); +1097 [~,targetString] = regexp(searchString,['<rdf:li rdf:resource="' startString fieldName midString '.*?"/>'],'tokens','match'); +1098 targetString=regexprep(targetString,'<rdf:li rdf:resource="|"/>',''); +1099 targetString=regexprep(targetString,startString,''); +1100 targetString=regexprep(targetString,[fieldName midString],''); +1101 +1102 for i=1:numel(targetString) +1103 fieldContent=[fieldContent ';' strtrim(targetString{1,i})]; +1104 end +1105 +1106 fieldContent=regexprep(fieldContent,'^;|;$',''); +1107 end +1108 +1109 function miriamStruct=parseMiriam(searchString) +1110 %Generates miriam structure from annotation field +1111 +1112 %Finding whether miriams are written in the old or the new way +1113 if strfind(searchString,'urn:miriam:') +1114 startString='urn:miriam:'; +1115 elseif strfind(searchString,'http://identifiers.org/') +1116 startString='http://identifiers.org/'; +1117 elseif strfind(searchString,'https://identifiers.org/') +1118 startString='https://identifiers.org/'; +1119 else +1120 miriamStruct=[]; +1121 return; +1122 end +1123 +1124 miriamStruct=[]; +1125 +1126 searchString=regexprep(searchString,'" />','"/>'); +1127 [~,targetString] = regexp(searchString,'<rdf:li rdf:resource=".*?"/>','tokens','match'); +1128 targetString=regexprep(targetString,'<rdf:li rdf:resource="|"/>',''); +1129 targetString=regexprep(targetString,startString,''); +1130 +1131 fwdslash = contains(targetString,'/'); +1132 midString = cell(numel(targetString),1); +1133 midString(fwdslash) = {'/'}; +1134 midString(~fwdslash) = {':'}; 1135 -1136 function miriam = addSBOtoMiriam(miriam,sboTerm) -1137 %Appends SBO term to miriam structure -1138 -1139 sboTerm = {['SBO:' sprintf('%07u',sboTerm)]}; % convert to proper format -1140 if isempty(miriam) -1141 miriam.name = {'sbo'}; -1142 miriam.value = sboTerm; -1143 elseif any(strcmp('sbo',miriam.name)) -1144 currSbo = strcmp('sbo',miriam.name); -1145 miriam.value(currSbo) = sboTerm; -1146 else -1147 miriam.name(end+1) = {'sbo'}; -1148 miriam.value(end+1) = sboTerm; -1149 end -1150 end +1136 counter=0; +1137 for i=1:numel(targetString) +1138 if isempty(regexp(targetString{1,i},'inchi|ec-code', 'once')) +1139 counter=counter+1; +1140 miriamStruct.name{counter,1} = regexprep(targetString{1,i},[midString{i} '.+'],'','once'); +1141 miriamStruct.value{counter,1} = regexprep(targetString{1,i},[miriamStruct.name{counter,1} midString{i}],'','once'); +1142 miriamStruct.name{counter,1} = regexprep(miriamStruct.name{counter,1},'^obo\.',''); +1143 end +1144 end +1145 end +1146 +1147 function miriam = addSBOtoMiriam(miriam,sboTerm) +1148 %Appends SBO term to miriam structure +1149 +1150 sboTerm = {['SBO:' sprintf('%07u',sboTerm)]}; % convert to proper format +1151 if isempty(miriam) +1152 miriam.name = {'sbo'}; +1153 miriam.value = sboTerm; +1154 elseif any(strcmp('sbo',miriam.name)) +1155 currSbo = strcmp('sbo',miriam.name); +1156 miriam.value(currSbo) = sboTerm; +1157 else +1158 miriam.name(end+1) = {'sbo'}; +1159 miriam.value(end+1) = sboTerm; +1160 end +1161 end
    Generated by m2html © 2005
    \ No newline at end of file diff --git a/doc/io/readYAMLmodel.html b/doc/io/readYAMLmodel.html index d1e03747e..c295b9e19 100644 --- a/doc/io/readYAMLmodel.html +++ b/doc/io/readYAMLmodel.html @@ -205,7 +205,7 @@

    SOURCE CODE ^% import different sections 0150 switch tline_raw -0151 case '- metaData:' +0151 case {'- metaData:', '- metaData: !!omap'} 0152 section = 1; 0153 if verbose 0154 fprintf('\t%d\n', section); @@ -225,7 +225,7 @@

    SOURCE CODE ^end 0169 pos=0; 0170 continue -0171 case '- genes:' +0171 case {'- genes:', '- genes: []'} 0172 section = 4; 0173 if verbose 0174 fprintf('\t%d\n', section); @@ -291,7 +291,7 @@

    SOURCE CODE ^case 'organization' 0236 model.annotation.organization = tline_value; -0237 case 'geckoLight' +0237 case {'geckoLight', 'gecko_light'} 0238 if strcmp(tline_value,'true') 0239 model.ec.geckoLight = true; 0240 end @@ -324,427 +324,466 @@

    SOURCE CODE ^readFieldValue(model, 'inchis', tline_value, pos); 0268 readList=''; miriamKey=''; 0269 case 'smiles' -0270 model = readFieldValue(model, 'metSmiles', tline_value, pos); -0271 readList=''; miriamKey=''; -0272 case 'deltaG' -0273 model = readFieldValue(model, 'metDeltaG', tline_value, pos); -0274 readList=''; miriamKey=''; -0275 case 'metFrom' -0276 model = readFieldValue(model, 'metFrom', tline_value, pos); -0277 readList=''; miriamKey=''; -0278 case 'annotation' -0279 readList = 'annotation'; -0280 otherwise -0281 switch readList -0282 case 'annotation' -0283 [metMiriams, miriamKey, metMirNo] = gatherAnnotation(pos,metMiriams,tline_key,tline_value,miriamKey,metMirNo); -0284 otherwise -0285 error(['Unknown entry in yaml file: ' tline_raw]) -0286 end -0287 end; continue -0288 end -0289 -0290 % import reactions: -0291 if section == 3 -0292 switch tline_key -0293 case 'id' -0294 pos = pos + 1; -0295 model = readFieldValue(model, 'rxns', tline_value,pos); -0296 readList=''; miriamKey=''; -0297 case 'name' -0298 model = readFieldValue(model, 'rxnNames', tline_value, pos); -0299 readList=''; miriamKey=''; -0300 case 'lower_bound' -0301 model.lb(pos,1) = {tline_value}; -0302 readList=''; miriamKey=''; -0303 case 'upper_bound' -0304 model.ub(pos,1) = {tline_value}; -0305 readList=''; miriamKey=''; -0306 case 'rev' -0307 model.rev(pos,1) = {tline_value}; -0308 readList=''; miriamKey=''; -0309 case 'gene_reaction_rule' -0310 model = readFieldValue(model, 'grRules', tline_value, pos); -0311 readList=''; miriamKey=''; -0312 case 'rxnNotes' -0313 model = readFieldValue(model, 'rxnNotes', tline_value, pos); -0314 readList=''; miriamKey=''; -0315 case 'rxnFrom' -0316 model = readFieldValue(model, 'rxnFrom', tline_value, pos); -0317 readList=''; miriamKey=''; -0318 case 'deltaG' -0319 model = readFieldValue(model, 'rxnDeltaG', tline_value, pos); -0320 readList=''; miriamKey=''; -0321 case 'objective_coefficient' -0322 model.c(pos,1) = 1; -0323 readList=''; miriamKey=''; -0324 case 'references' -0325 model = readFieldValue(model, 'rxnReferences', tline_value, pos); -0326 readList=''; miriamKey=''; -0327 case 'confidence_score' -0328 model = readFieldValue(model, 'rxnConfidenceScores', tline_value, pos); -0329 readList=''; miriamKey=''; -0330 case 'eccodes' -0331 if isempty(tline_value) -0332 readList = 'eccodes'; -0333 else -0334 eccodes(ecCodeNo,1:2)={pos,tline_value}; -0335 ecCodeNo=ecCodeNo+1; -0336 end -0337 case 'subsystem' -0338 if isempty(tline_value) -0339 readList = 'subsystem'; -0340 else -0341 subSystems(subSysNo,1:2)={pos,tline_value}; -0342 subSysNo=subSysNo+1; -0343 end -0344 case 'metabolites' -0345 readList = 'equation'; -0346 case 'annotation' -0347 readList = 'annotation'; -0348 -0349 otherwise -0350 switch readList -0351 case 'eccodes' -0352 eccodes(ecCodeNo,1:2)={pos,regexprep(tline_value,'^ +- "?(.*)"?$','$1')}; -0353 ecCodeNo=ecCodeNo+1; -0354 case 'subsystem' -0355 subSystems(subSysNo,1:2)={pos,regexprep(tline_value,'^ +- "?(.*)"?$','$1')}; -0356 subSysNo=subSysNo+1; -0357 case 'annotation' -0358 [rxnMiriams, miriamKey,rxnMirNo] = gatherAnnotation(pos,rxnMiriams,tline_key,tline_value,miriamKey,rxnMirNo); -0359 case 'equation' -0360 coeff = sscanf(tline_value,'%f'); -0361 equations(equatiNo,1:3)={pos,tline_key,coeff}; -0362 equatiNo=equatiNo+1; -0363 otherwise -0364 error(['Unknown entry in yaml file: ' tline_raw]) -0365 end -0366 end; continue -0367 end -0368 -0369 % import genes: -0370 if section == 4 -0371 switch tline_key -0372 case 'id' -0373 pos = pos + 1; -0374 model = readFieldValue(model, 'genes', tline_value, pos); -0375 readList = ''; -0376 miriamKey = ''; -0377 case 'name' -0378 model = readFieldValue(model, 'geneShortNames', tline_value, pos); -0379 case 'protein' -0380 model = readFieldValue(model, 'proteins', tline_value, pos); -0381 case 'annotation' -0382 readList = 'annotation'; -0383 otherwise -0384 switch readList -0385 case 'annotation' -0386 [geneMiriams, miriamKey,genMirNo] = gatherAnnotation(pos,geneMiriams,tline_key,tline_value,miriamKey,genMirNo); -0387 otherwise -0388 error(['Unknown entry in yaml file: ' tline_raw]) -0389 end -0390 end; continue -0391 end -0392 -0393 % import compartments: -0394 if section == 5 -0395 model.comps(end+1,1) = {tline_key}; -0396 model.compNames(end+1,1) = {tline_value}; -0397 end -0398 -0399 % import ec reaction info -0400 if section == 6 -0401 switch tline_key -0402 case 'id' -0403 pos = pos + 1; -0404 model.ec = readFieldValue(model.ec, 'rxns', tline_value, pos); -0405 readList=''; -0406 case 'kcat' -0407 model.ec = readFieldValue(model.ec, 'kcat', tline_value, pos); -0408 readList=''; -0409 case 'source' -0410 model.ec = readFieldValue(model.ec, 'source', tline_value, pos); -0411 readList=''; -0412 case 'notes' -0413 model.ec = readFieldValue(model.ec, 'notes', tline_value, pos); -0414 readList=''; -0415 case 'eccodes' -0416 if isempty(tline_value) -0417 readList = 'eccodes'; -0418 else -0419 ecGecko(ecGeckoNo,1:2)={pos,tline_value}; -0420 ecGeckoNo=ecGeckoNo+1; -0421 end -0422 case 'enzymes' -0423 readList = 'enzStoich'; -0424 otherwise -0425 switch readList -0426 case 'eccodes' -0427 ecGecko(ecGeckoNo,1:2)={pos,regexprep(tline_value,'^ +- "?(.*)"?$','$1')}; -0428 ecGeckoNo=ecGeckoNo+1; -0429 case 'enzStoich' -0430 coeff = sscanf(tline_value,'%f'); -0431 enzStoich(enzStoichNo,1:3)={pos,tline_key,coeff}; -0432 enzStoichNo=enzStoichNo+1; -0433 otherwise -0434 error(['Unknown entry in yaml file: ' tline_raw]) -0435 end -0436 end; continue -0437 end -0438 -0439 % import ec enzyme info -0440 if section == 7 -0441 switch tline_key -0442 case 'genes' -0443 pos = pos + 1; -0444 model.ec = readFieldValue(model.ec, 'genes', tline_value, pos); -0445 case 'enzymes' -0446 model.ec = readFieldValue(model.ec, 'enzymes', tline_value, pos); -0447 case 'mw' -0448 model.ec = readFieldValue(model.ec, 'mw', tline_value, pos); -0449 case 'sequence' -0450 model.ec = readFieldValue(model.ec, 'sequence', tline_value, pos); -0451 case 'concs' -0452 model.ec = readFieldValue(model.ec, 'concs', tline_value, pos); -0453 otherwise -0454 error(['Unknown entry in yaml file: ' tline_raw]) -0455 end; continue -0456 end -0457 end -0458 -0459 %Parse annotations -0460 if ~isempty(metMiriams) -0461 locs = cell2mat(metMiriams(:,1)); -0462 for i=unique(locs)' -0463 model.metMiriams{i,1}.name=metMiriams(locs==i,2); -0464 model.metMiriams{i,1}.value=metMiriams(locs==i,3); -0465 end -0466 end -0467 if ~isempty(rxnMiriams) -0468 locs = cell2mat(rxnMiriams(:,1)); -0469 for i=unique(locs)' -0470 model.rxnMiriams{i,1}.name=rxnMiriams(locs==i,2); -0471 model.rxnMiriams{i,1}.value=rxnMiriams(locs==i,3); -0472 end -0473 end -0474 if ~isempty(geneMiriams) -0475 locs = cell2mat(geneMiriams(:,1)); -0476 for i=unique(locs)' -0477 model.geneMiriams{i,1}.name=geneMiriams(locs==i,2); -0478 model.geneMiriams{i,1}.value=geneMiriams(locs==i,3); -0479 end -0480 end -0481 -0482 %Parse subSystems -0483 if ~isempty(subSystems) -0484 locs = cell2mat(subSystems(:,1)); -0485 for i=unique(locs)' -0486 model.subSystems{i,1}=subSystems(locs==i,2); -0487 end -0488 end -0489 -0490 %Parse ec-codes -0491 if ~isempty(eccodes) -0492 locs = cell2mat(eccodes(:,1)); -0493 for i=unique(locs)' -0494 eccodesCat=strjoin(eccodes(locs==i,2),';'); -0495 model.eccodes{i,1}=eccodesCat; -0496 end -0497 emptyEc=cellfun('isempty',model.eccodes); -0498 model.eccodes(emptyEc)={''}; -0499 end -0500 -0501 % follow-up data processing -0502 if verbose -0503 fprintf('\nimporting completed\nfollow-up processing...'); -0504 end -0505 [~, model.metComps] = ismember(model.metComps, model.comps); -0506 [~, model.geneComps] = ismember(model.geneComps, model.comps); -0507 [~, model.rxnComps] = ismember(model.rxnComps, model.comps); -0508 -0509 % Fill S-matrix -0510 rxnIdx = cellfun('isempty', equations(:,1)); -0511 equations(rxnIdx,:) = ''; -0512 rxnIdx = cell2mat(equations(:,1)); -0513 [~,metIdx] = ismember(equations(:,2),model.mets); -0514 coeffs = cell2mat(equations(:,3)); -0515 model.S=sparse(max(metIdx),max(rxnIdx)); -0516 linearIndices = sub2ind([max(metIdx), max(rxnIdx)],metIdx,rxnIdx); -0517 model.S(linearIndices) = coeffs; -0518 -0519 % Convert strings to numeric -0520 model.metCharges = str2double(model.metCharges); -0521 model.lb = str2double(model.lb); -0522 model.ub = str2double(model.ub); -0523 model.rxnConfidenceScores = str2double(model.rxnConfidenceScores); -0524 model.b = zeros(length(model.mets),1); -0525 model.metDeltaG = str2double(model.metDeltaG); -0526 model.rxnDeltaG = str2double(model.rxnDeltaG); -0527 -0528 % Fill some other fields -0529 model.annotation.defaultLB = min(model.lb); -0530 model.annotation.defaultUB = max(model.ub); -0531 if numel(model.lb)<numel(model.rxns) %No LB reported = min -0532 model.lb(end+1:numel(model.rxns)-numel(model.lb),1) = double(model.annotation.defaultLB); -0533 end -0534 if numel(model.ub)<numel(model.rxns) %No UB reported = max -0535 model.ub(end+1:numel(model.rxns)-numel(model.ub),1) = double(model.annotation.defaultUB); -0536 end -0537 if ~all(cellfun('isempty',model.rev)) -0538 model.rev = str2double(model.rev); -0539 else -0540 model.rev = []; -0541 end -0542 if numel(model.rev)<numel(model.rxns) %No rev reported, assume from LB and UB -0543 model.rev(end+1:numel(model.rxns)-numel(model.rev),1) = double(model.lb<0 & model.ub>0); -0544 end -0545 -0546 % Remove empty fields, otherwise fill to correct length -0547 % Reactions -0548 for i={'rxnNames','grRules','eccodes','rxnNotes','rxnReferences',... -0549 'rxnFrom','subSystems','rxnMiriams'} % Empty strings -0550 model = emptyOrFill(model,i{1},{''},'rxns'); -0551 end -0552 for i={'c'} % Zeros -0553 model = emptyOrFill(model,i{1},0,'rxns',true); -0554 end -0555 for i={'rxnConfidenceScores','rxnDeltaG'} % NaNs -0556 model = emptyOrFill(model,i{1},NaN,'rxns'); -0557 end -0558 for i={'rxnComps'} % Ones, assume first compartment -0559 model = emptyOrFill(model,i{1},1,'rxns'); -0560 end -0561 % Metabolites -0562 for i={'metNames','inchis','metFormulas','metMiriams','metFrom','metSmiles','metNotes'} % Empty strings -0563 model = emptyOrFill(model,i{1},{''},'mets'); -0564 end -0565 for i={'metCharges','unconstrained'} % Zeros -0566 model = emptyOrFill(model,i{1},0,'mets'); -0567 end -0568 for i={'metDeltaG'} % % NaNs -0569 model = emptyOrFill(model,i{1},NaN,'mets'); -0570 end -0571 for i={'metComps'} % Ones, assume first compartment -0572 model = emptyOrFill(model,i{1},1,'mets'); -0573 end -0574 % Genes -0575 for i={'geneMiriams','geneShortNames','proteins'} % Empty strings -0576 model = emptyOrFill(model,i{1},{''},'genes'); -0577 end -0578 for i={'geneComps'} % Ones, assume first compartment -0579 model = emptyOrFill(model,i{1},1,'genes'); +0270 if isempty(tline_value) +0271 readList = 'smilesInAnnotation'; +0272 else +0273 model = readFieldValue(model, 'metSmiles', tline_value, pos); +0274 readList=''; miriamKey=''; +0275 end +0276 case 'deltaG' +0277 model = readFieldValue(model, 'metDeltaG', tline_value, pos); +0278 readList=''; miriamKey=''; +0279 case 'metFrom' +0280 model = readFieldValue(model, 'metFrom', tline_value, pos); +0281 readList=''; miriamKey=''; +0282 case 'annotation' +0283 readList = 'annotation'; +0284 otherwise +0285 switch readList +0286 case 'smilesInAnnotation' +0287 %As for ec-code on a reaction: a bare list item +0288 %continues the list, anything else ends it. +0289 if isempty(tline_key) +0290 model = readFieldValue(model, 'metSmiles', ... +0291 regexprep(tline_value,'^ +- "?(.*)"?$','$1'), pos); +0292 else +0293 readList = 'annotation'; +0294 [metMiriams, miriamKey, metMirNo] = gatherAnnotation(pos,metMiriams,tline_key,tline_value,miriamKey,metMirNo); +0295 end +0296 case 'annotation' +0297 [metMiriams, miriamKey, metMirNo] = gatherAnnotation(pos,metMiriams,tline_key,tline_value,miriamKey,metMirNo); +0298 otherwise +0299 error(['Unknown entry in yaml file: ' tline_raw]) +0300 end +0301 end; continue +0302 end +0303 +0304 % import reactions: +0305 if section == 3 +0306 switch tline_key +0307 case 'id' +0308 pos = pos + 1; +0309 model = readFieldValue(model, 'rxns', tline_value,pos); +0310 readList=''; miriamKey=''; +0311 case 'name' +0312 model = readFieldValue(model, 'rxnNames', tline_value, pos); +0313 readList=''; miriamKey=''; +0314 case 'lower_bound' +0315 model.lb(pos,1) = {tline_value}; +0316 readList=''; miriamKey=''; +0317 case 'upper_bound' +0318 model.ub(pos,1) = {tline_value}; +0319 readList=''; miriamKey=''; +0320 case 'rev' +0321 model.rev(pos,1) = {tline_value}; +0322 readList=''; miriamKey=''; +0323 case 'gene_reaction_rule' +0324 model = readFieldValue(model, 'grRules', tline_value, pos); +0325 readList=''; miriamKey=''; +0326 case {'rxnNotes', 'notes'} +0327 model = readFieldValue(model, 'rxnNotes', tline_value, pos); +0328 readList=''; miriamKey=''; +0329 case 'rxnFrom' +0330 model = readFieldValue(model, 'rxnFrom', tline_value, pos); +0331 readList=''; miriamKey=''; +0332 case 'deltaG' +0333 model = readFieldValue(model, 'rxnDeltaG', tline_value, pos); +0334 readList=''; miriamKey=''; +0335 case 'objective_coefficient' +0336 model.c(pos,1) = 1; +0337 readList=''; miriamKey=''; +0338 case 'references' +0339 model = readFieldValue(model, 'rxnReferences', tline_value, pos); +0340 readList=''; miriamKey=''; +0341 case 'confidence_score' +0342 model = readFieldValue(model, 'rxnConfidenceScores', tline_value, pos); +0343 readList=''; miriamKey=''; +0344 case 'eccodes' +0345 if isempty(tline_value) +0346 readList = 'eccodes'; +0347 else +0348 eccodes(ecCodeNo,1:2)={pos,tline_value}; +0349 ecCodeNo=ecCodeNo+1; +0350 end +0351 case 'subsystem' +0352 if isempty(tline_value) +0353 readList = 'subsystem'; +0354 else +0355 subSystems(subSysNo,1:2)={pos,tline_value}; +0356 subSysNo=subSysNo+1; +0357 end +0358 case 'metabolites' +0359 readList = 'equation'; +0360 case 'ec-code' +0361 %EC numbers emitted inside the annotation block rather than +0362 %under the legacy top-level 'eccodes' key, which is still +0363 %accepted above. A block list follows when the value is +0364 %empty; it is read under its own state so that gathering can +0365 %resume for the rest of the annotation block once the list +0366 %ends, instead of swallowing every later entry as an EC code. +0367 if isempty(tline_value) +0368 readList = 'eccodesInAnnotation'; +0369 else +0370 eccodes(ecCodeNo,1:2)={pos,tline_value}; +0371 ecCodeNo=ecCodeNo+1; +0372 end +0373 case 'annotation' +0374 readList = 'annotation'; +0375 +0376 otherwise +0377 switch readList +0378 case 'eccodes' +0379 eccodes(ecCodeNo,1:2)={pos,regexprep(tline_value,'^ +- "?(.*)"?$','$1')}; +0380 ecCodeNo=ecCodeNo+1; +0381 case 'eccodesInAnnotation' +0382 %A bare list item ("- 1.1.2.4") has no colon, so +0383 %tline_key is empty and the list continues. Anything +0384 %else is the annotation block's next entry, so the +0385 %list has ended. +0386 if isempty(tline_key) +0387 eccodes(ecCodeNo,1:2)={pos,regexprep(tline_value,'^ +- "?(.*)"?$','$1')}; +0388 ecCodeNo=ecCodeNo+1; +0389 else +0390 readList = 'annotation'; +0391 [rxnMiriams, miriamKey,rxnMirNo] = gatherAnnotation(pos,rxnMiriams,tline_key,tline_value,miriamKey,rxnMirNo); +0392 end +0393 case 'subsystem' +0394 subSystems(subSysNo,1:2)={pos,regexprep(tline_value,'^ +- "?(.*)"?$','$1')}; +0395 subSysNo=subSysNo+1; +0396 case 'annotation' +0397 [rxnMiriams, miriamKey,rxnMirNo] = gatherAnnotation(pos,rxnMiriams,tline_key,tline_value,miriamKey,rxnMirNo); +0398 case 'equation' +0399 coeff = sscanf(tline_value,'%f'); +0400 equations(equatiNo,1:3)={pos,tline_key,coeff}; +0401 equatiNo=equatiNo+1; +0402 otherwise +0403 error(['Unknown entry in yaml file: ' tline_raw]) +0404 end +0405 end; continue +0406 end +0407 +0408 % import genes: +0409 if section == 4 +0410 switch tline_key +0411 case 'id' +0412 pos = pos + 1; +0413 model = readFieldValue(model, 'genes', tline_value, pos); +0414 readList = ''; +0415 miriamKey = ''; +0416 case 'name' +0417 model = readFieldValue(model, 'geneShortNames', tline_value, pos); +0418 case 'protein' +0419 model = readFieldValue(model, 'proteins', tline_value, pos); +0420 case 'annotation' +0421 readList = 'annotation'; +0422 otherwise +0423 switch readList +0424 case 'annotation' +0425 [geneMiriams, miriamKey,genMirNo] = gatherAnnotation(pos,geneMiriams,tline_key,tline_value,miriamKey,genMirNo); +0426 otherwise +0427 error(['Unknown entry in yaml file: ' tline_raw]) +0428 end +0429 end; continue +0430 end +0431 +0432 % import compartments: +0433 if section == 5 +0434 model.comps(end+1,1) = {tline_key}; +0435 model.compNames(end+1,1) = {tline_value}; +0436 end +0437 +0438 % import ec reaction info +0439 if section == 6 +0440 switch tline_key +0441 case 'id' +0442 pos = pos + 1; +0443 model.ec = readFieldValue(model.ec, 'rxns', tline_value, pos); +0444 readList=''; +0445 case 'kcat' +0446 model.ec = readFieldValue(model.ec, 'kcat', tline_value, pos); +0447 readList=''; +0448 case 'source' +0449 model.ec = readFieldValue(model.ec, 'source', tline_value, pos); +0450 readList=''; +0451 case 'notes' +0452 model.ec = readFieldValue(model.ec, 'notes', tline_value, pos); +0453 readList=''; +0454 case 'eccodes' +0455 if isempty(tline_value) +0456 readList = 'eccodes'; +0457 else +0458 ecGecko(ecGeckoNo,1:2)={pos,tline_value}; +0459 ecGeckoNo=ecGeckoNo+1; +0460 end +0461 case 'enzymes' +0462 readList = 'enzStoich'; +0463 otherwise +0464 switch readList +0465 case 'eccodes' +0466 ecGecko(ecGeckoNo,1:2)={pos,regexprep(tline_value,'^ +- "?(.*)"?$','$1')}; +0467 ecGeckoNo=ecGeckoNo+1; +0468 case 'enzStoich' +0469 coeff = sscanf(tline_value,'%f'); +0470 enzStoich(enzStoichNo,1:3)={pos,tline_key,coeff}; +0471 enzStoichNo=enzStoichNo+1; +0472 otherwise +0473 error(['Unknown entry in yaml file: ' tline_raw]) +0474 end +0475 end; continue +0476 end +0477 +0478 % import ec enzyme info +0479 if section == 7 +0480 switch tline_key +0481 case 'genes' +0482 pos = pos + 1; +0483 model.ec = readFieldValue(model.ec, 'genes', tline_value, pos); +0484 case 'enzymes' +0485 model.ec = readFieldValue(model.ec, 'enzymes', tline_value, pos); +0486 case 'mw' +0487 model.ec = readFieldValue(model.ec, 'mw', tline_value, pos); +0488 case 'sequence' +0489 model.ec = readFieldValue(model.ec, 'sequence', tline_value, pos); +0490 case 'concs' +0491 model.ec = readFieldValue(model.ec, 'concs', tline_value, pos); +0492 otherwise +0493 error(['Unknown entry in yaml file: ' tline_raw]) +0494 end; continue +0495 end +0496 end +0497 +0498 %Parse annotations +0499 if ~isempty(metMiriams) +0500 locs = cell2mat(metMiriams(:,1)); +0501 for i=unique(locs)' +0502 model.metMiriams{i,1}.name=metMiriams(locs==i,2); +0503 model.metMiriams{i,1}.value=metMiriams(locs==i,3); +0504 end +0505 end +0506 if ~isempty(rxnMiriams) +0507 locs = cell2mat(rxnMiriams(:,1)); +0508 for i=unique(locs)' +0509 model.rxnMiriams{i,1}.name=rxnMiriams(locs==i,2); +0510 model.rxnMiriams{i,1}.value=rxnMiriams(locs==i,3); +0511 end +0512 end +0513 if ~isempty(geneMiriams) +0514 locs = cell2mat(geneMiriams(:,1)); +0515 for i=unique(locs)' +0516 model.geneMiriams{i,1}.name=geneMiriams(locs==i,2); +0517 model.geneMiriams{i,1}.value=geneMiriams(locs==i,3); +0518 end +0519 end +0520 +0521 %Parse subSystems +0522 if ~isempty(subSystems) +0523 locs = cell2mat(subSystems(:,1)); +0524 for i=unique(locs)' +0525 model.subSystems{i,1}=subSystems(locs==i,2); +0526 end +0527 end +0528 +0529 %Parse ec-codes +0530 if ~isempty(eccodes) +0531 locs = cell2mat(eccodes(:,1)); +0532 for i=unique(locs)' +0533 eccodesCat=strjoin(eccodes(locs==i,2),';'); +0534 model.eccodes{i,1}=eccodesCat; +0535 end +0536 emptyEc=cellfun('isempty',model.eccodes); +0537 model.eccodes(emptyEc)={''}; +0538 end +0539 +0540 % follow-up data processing +0541 if verbose +0542 fprintf('\nimporting completed\nfollow-up processing...'); +0543 end +0544 [~, model.metComps] = ismember(model.metComps, model.comps); +0545 [~, model.geneComps] = ismember(model.geneComps, model.comps); +0546 [~, model.rxnComps] = ismember(model.rxnComps, model.comps); +0547 +0548 % Fill S-matrix +0549 rxnIdx = cellfun('isempty', equations(:,1)); +0550 equations(rxnIdx,:) = ''; +0551 rxnIdx = cell2mat(equations(:,1)); +0552 [~,metIdx] = ismember(equations(:,2),model.mets); +0553 coeffs = cell2mat(equations(:,3)); +0554 model.S=sparse(max(metIdx),max(rxnIdx)); +0555 linearIndices = sub2ind([max(metIdx), max(rxnIdx)],metIdx,rxnIdx); +0556 model.S(linearIndices) = coeffs; +0557 +0558 % Convert strings to numeric +0559 model.metCharges = str2double(model.metCharges); +0560 model.lb = str2double(model.lb); +0561 model.ub = str2double(model.ub); +0562 model.rxnConfidenceScores = str2double(model.rxnConfidenceScores); +0563 model.b = zeros(length(model.mets),1); +0564 model.metDeltaG = str2double(model.metDeltaG); +0565 model.rxnDeltaG = str2double(model.rxnDeltaG); +0566 +0567 % Fill some other fields +0568 model.annotation.defaultLB = min(model.lb); +0569 model.annotation.defaultUB = max(model.ub); +0570 if numel(model.lb)<numel(model.rxns) %No LB reported = min +0571 model.lb(end+1:numel(model.rxns)-numel(model.lb),1) = double(model.annotation.defaultLB); +0572 end +0573 if numel(model.ub)<numel(model.rxns) %No UB reported = max +0574 model.ub(end+1:numel(model.rxns)-numel(model.ub),1) = double(model.annotation.defaultUB); +0575 end +0576 if ~all(cellfun('isempty',model.rev)) +0577 model.rev = str2double(model.rev); +0578 else +0579 model.rev = []; 0580 end -0581 % Comps -0582 for i={'compNames'} % Empty strings -0583 model = emptyOrFill(model,i{1},{''},'comps'); -0584 end -0585 for i={'compOutside'} % First comp -0586 model = emptyOrFill(model,i{1},model.comps{1},'comps'); -0587 end -0588 % Single fields are kept, even if empty -0589 % for i={'description','name','version','date','annotation'} -0590 % if isempty(model.(i{1})) -0591 % model = rmfield(model,i{1}); -0592 % end -0593 % end -0594 -0595 if isfield(model,'subSystems') -0596 % If all entries are 1x1, then flatten -0597 if all(cellfun(@(x) numel(x) <= 1, model.subSystems)) -0598 model.subSystems = transpose([model.subSystems{:}]); -0599 end -0600 end -0601 -0602 % Make rxnGeneMat fields and map to the existing model.genes field -0603 if isfield(model,'grRules') -0604 [genes, rxnGeneMat] = getGenesFromGrRules(model.grRules); -0605 model.rxnGeneMat = sparse(numel(model.rxns),numel(model.genes)); -0606 [~,geneOrder] = ismember(genes,model.genes); -0607 if any(geneOrder == 0) -0608 error(['The grRules includes the following gene(s), that are not in '... -0609 'the list of model genes: ', genes{~geneOrder}]) -0610 end -0611 model.rxnGeneMat(:,geneOrder) = rxnGeneMat; +0581 if numel(model.rev)<numel(model.rxns) %No rev reported, assume from LB and UB +0582 model.rev(end+1:numel(model.rxns)-numel(model.rev),1) = double(model.lb<0 & model.ub>0); +0583 end +0584 +0585 % Remove empty fields, otherwise fill to correct length +0586 % Reactions +0587 for i={'rxnNames','grRules','eccodes','rxnNotes','rxnReferences',... +0588 'rxnFrom','subSystems','rxnMiriams'} % Empty strings +0589 model = emptyOrFill(model,i{1},{''},'rxns'); +0590 end +0591 for i={'c'} % Zeros +0592 model = emptyOrFill(model,i{1},0,'rxns',true); +0593 end +0594 for i={'rxnConfidenceScores','rxnDeltaG'} % NaNs +0595 model = emptyOrFill(model,i{1},NaN,'rxns'); +0596 end +0597 for i={'rxnComps'} % Ones, assume first compartment +0598 model = emptyOrFill(model,i{1},1,'rxns'); +0599 end +0600 % Metabolites +0601 for i={'metNames','inchis','metFormulas','metMiriams','metFrom','metSmiles','metNotes'} % Empty strings +0602 model = emptyOrFill(model,i{1},{''},'mets'); +0603 end +0604 for i={'metCharges','unconstrained'} % Zeros +0605 model = emptyOrFill(model,i{1},0,'mets'); +0606 end +0607 for i={'metDeltaG'} % % NaNs +0608 model = emptyOrFill(model,i{1},NaN,'mets'); +0609 end +0610 for i={'metComps'} % Ones, assume first compartment +0611 model = emptyOrFill(model,i{1},1,'mets'); 0612 end -0613 -0614 % Finalize GECKO model -0615 if isGECKO -0616 % Fill in empty fields and empty entries -0617 for i={'kcat','source','notes','eccodes'} % Even keep empty -0618 model.ec = emptyOrFill(model.ec,i{1},{''},'rxns',true); -0619 end -0620 for i={'enzymes','mw','sequence'} -0621 model.ec = emptyOrFill(model.ec,i{1},{''},'genes',true); -0622 end -0623 model.ec = emptyOrFill(model.ec,'concs',{'NaN'},'genes',true); -0624 model.ec = emptyOrFill(model.ec,'kcat',{'0'},'genes',true); -0625 % Change string to double -0626 for i={'kcat','mw','concs'} -0627 if isfield(model.ec,i{1}) -0628 model.ec.(i{1}) = str2double(model.ec.(i{1})); -0629 end -0630 end -0631 % Fill rxnEnzMat -0632 rxnIdx = cellfun('isempty', enzStoich(:,1)); -0633 enzStoich(rxnIdx,:) = ''; -0634 rxnIdx = cell2mat(enzStoich(:,1)); -0635 [~,enzIdx] = ismember(enzStoich(:,2),model.ec.enzymes); -0636 coeffs = cell2mat(enzStoich(:,3)); -0637 model.ec.rxnEnzMat = zeros(numel(model.ec.rxns), numel(model.ec.genes)); -0638 linearIndices = sub2ind([numel(model.ec.rxns), numel(model.ec.genes)], rxnIdx, enzIdx); -0639 model.ec.rxnEnzMat(linearIndices) = coeffs; -0640 %Parse ec-codes -0641 if ~isempty(ecGecko) -0642 locs = cell2mat(ecGecko(:,1)); -0643 for i=unique(locs)' -0644 ecGeckoCat=strjoin(ecGecko(locs==i,2),';'); -0645 model.ec.eccodes{i,1}=ecGeckoCat; -0646 end -0647 emptyEc=cellfun('isempty',model.ec.eccodes); -0648 model.ec.eccodes(emptyEc)={''}; +0613 % Genes +0614 for i={'geneMiriams','geneShortNames','proteins'} % Empty strings +0615 model = emptyOrFill(model,i{1},{''},'genes'); +0616 end +0617 for i={'geneComps'} % Ones, assume first compartment +0618 model = emptyOrFill(model,i{1},1,'genes'); +0619 end +0620 % Comps +0621 for i={'compNames'} % Empty strings +0622 model = emptyOrFill(model,i{1},{''},'comps'); +0623 end +0624 for i={'compOutside'} % First comp +0625 model = emptyOrFill(model,i{1},model.comps{1},'comps'); +0626 end +0627 % Single fields are kept, even if empty +0628 % for i={'description','name','version','date','annotation'} +0629 % if isempty(model.(i{1})) +0630 % model = rmfield(model,i{1}); +0631 % end +0632 % end +0633 +0634 if isfield(model,'subSystems') +0635 % If all entries are 1x1, then flatten +0636 if all(cellfun(@(x) numel(x) <= 1, model.subSystems)) +0637 model.subSystems = transpose([model.subSystems{:}]); +0638 end +0639 end +0640 +0641 % Make rxnGeneMat fields and map to the existing model.genes field +0642 if isfield(model,'grRules') +0643 [genes, rxnGeneMat] = getGenesFromGrRules(model.grRules); +0644 model.rxnGeneMat = sparse(numel(model.rxns),numel(model.genes)); +0645 [~,geneOrder] = ismember(genes,model.genes); +0646 if any(geneOrder == 0) +0647 error(['The grRules includes the following gene(s), that are not in '... +0648 'the list of model genes: ', genes{~geneOrder}]) 0649 end -0650 end -0651 -0652 if verbose -0653 fprintf(' Done!\n'); -0654 end -0655 end -0656 -0657 function model = emptyOrFill(model,field,emptyEntry,type,keepEmpty) -0658 if nargin<5 -0659 keepEmpty=false; -0660 end -0661 if isnumeric(emptyEntry) -0662 emptyCells=isempty(model.(field)); -0663 else -0664 emptyCells=cellfun('isempty',model.(field)); -0665 end -0666 if all(emptyCells) && ~keepEmpty -0667 model = rmfield(model, field); -0668 elseif numel(model.(field))<numel(model.(type)) -0669 model.(field)(end+1:numel(model.(type)),1)=emptyEntry; -0670 end -0671 end -0672 -0673 function model = readFieldValue(model, fieldName, value, pos) -0674 if numel(model.(fieldName))<pos-1 -0675 model.(fieldName)(end+1:pos,1) = {''}; -0676 end -0677 model.(fieldName)(pos,1) = {value}; -0678 end -0679 -0680 function [miriams,miriamKey,entryNumber] = gatherAnnotation(pos,miriams,key,value,miriamKey,entryNumber) -0681 if isempty(key) -0682 key=miriamKey; -0683 else -0684 miriamKey=key; -0685 end -0686 if ~isempty(value) -0687 miriams(entryNumber,1:3) = {pos, key, strip(value)}; -0688 entryNumber = entryNumber + 1; +0650 model.rxnGeneMat(:,geneOrder) = rxnGeneMat; +0651 end +0652 +0653 % Finalize GECKO model +0654 if isGECKO +0655 % Fill in empty fields and empty entries +0656 for i={'kcat','source','notes','eccodes'} % Even keep empty +0657 model.ec = emptyOrFill(model.ec,i{1},{''},'rxns',true); +0658 end +0659 for i={'enzymes','mw','sequence'} +0660 model.ec = emptyOrFill(model.ec,i{1},{''},'genes',true); +0661 end +0662 model.ec = emptyOrFill(model.ec,'concs',{'NaN'},'genes',true); +0663 model.ec = emptyOrFill(model.ec,'kcat',{'0'},'genes',true); +0664 % Change string to double +0665 for i={'kcat','mw','concs'} +0666 if isfield(model.ec,i{1}) +0667 model.ec.(i{1}) = str2double(model.ec.(i{1})); +0668 end +0669 end +0670 % Fill rxnEnzMat +0671 rxnIdx = cellfun('isempty', enzStoich(:,1)); +0672 enzStoich(rxnIdx,:) = ''; +0673 rxnIdx = cell2mat(enzStoich(:,1)); +0674 [~,enzIdx] = ismember(enzStoich(:,2),model.ec.enzymes); +0675 coeffs = cell2mat(enzStoich(:,3)); +0676 model.ec.rxnEnzMat = zeros(numel(model.ec.rxns), numel(model.ec.genes)); +0677 linearIndices = sub2ind([numel(model.ec.rxns), numel(model.ec.genes)], rxnIdx, enzIdx); +0678 model.ec.rxnEnzMat(linearIndices) = coeffs; +0679 %Parse ec-codes +0680 if ~isempty(ecGecko) +0681 locs = cell2mat(ecGecko(:,1)); +0682 for i=unique(locs)' +0683 ecGeckoCat=strjoin(ecGecko(locs==i,2),';'); +0684 model.ec.eccodes{i,1}=ecGeckoCat; +0685 end +0686 emptyEc=cellfun('isempty',model.ec.eccodes); +0687 model.ec.eccodes(emptyEc)={''}; +0688 end 0689 end -0690 end +0690 +0691 if verbose +0692 fprintf(' Done!\n'); +0693 end +0694 end +0695 +0696 function model = emptyOrFill(model,field,emptyEntry,type,keepEmpty) +0697 if nargin<5 +0698 keepEmpty=false; +0699 end +0700 if isnumeric(emptyEntry) +0701 emptyCells=isempty(model.(field)); +0702 else +0703 emptyCells=cellfun('isempty',model.(field)); +0704 end +0705 if all(emptyCells) && ~keepEmpty +0706 model = rmfield(model, field); +0707 elseif numel(model.(field))<numel(model.(type)) +0708 model.(field)(end+1:numel(model.(type)),1)=emptyEntry; +0709 end +0710 end +0711 +0712 function model = readFieldValue(model, fieldName, value, pos) +0713 if numel(model.(fieldName))<pos-1 +0714 model.(fieldName)(end+1:pos,1) = {''}; +0715 end +0716 model.(fieldName)(pos,1) = {value}; +0717 end +0718 +0719 function [miriams,miriamKey,entryNumber] = gatherAnnotation(pos,miriams,key,value,miriamKey,entryNumber) +0720 if isempty(key) +0721 key=miriamKey; +0722 else +0723 miriamKey=key; +0724 end +0725 if ~isempty(value) +0726 miriams(entryNumber,1:3) = {pos, key, strip(value)}; +0727 entryNumber = entryNumber + 1; +0728 end +0729 end
    Generated by m2html © 2005
    \ No newline at end of file diff --git a/io/importModel.m b/io/importModel.m index dae96063c..b8fdade73 100755 --- a/io/importModel.m +++ b/io/importModel.m @@ -739,11 +739,22 @@ end end if isfield(modelSBML,'notes') - startString=strfind(modelSBML.notes,'xhtml">'); - endString=strfind(modelSBML.notes,''); - if any(startString) && any(endString) - model.annotation.note=modelSBML.notes(startString(1)+7:endString-1); - model.annotation.note=regexprep(model.annotation.note,'|',''); + %Two shapes reach here. RAVEN's own exportModel writes the note as bare + %text inside a ; a writer that serialises a notes dictionary emits + %one '

    key: value

    ' per entry and no at all, carrying the note + %under its own key. Take the labelled entry when there is one, so the + %other entries (version, and whatever else the dictionary held) stay out + %of the note, and fall back to the text otherwise. + note=parseNote(modelSBML.notes,'note'); + if isempty(note) + startString=strfind(modelSBML.notes,'xhtml">'); + endString=strfind(modelSBML.notes,''); + if any(startString) && any(endString) + note=modelSBML.notes(startString(1)+7:endString-1); + end + end + if ~isempty(note) + model.annotation.note=regexprep(note,'|',''); model.annotation.note=strtrim(model.annotation.note); if regexp(model.annotation.note,'This file was generated using the exportModel function in RAVEN Toolbox \d\.\d and OutputSBML in libSBML') model.annotation=rmfield(model.annotation,'note'); % Default note added when running exportModel diff --git a/io/readYAMLmodel.m b/io/readYAMLmodel.m index 1eb65f3b6..ea353f39f 100755 --- a/io/readYAMLmodel.m +++ b/io/readYAMLmodel.m @@ -148,7 +148,7 @@ tline_value = line_value{i}; % import different sections switch tline_raw - case '- metaData:' + case {'- metaData:', '- metaData: !!omap'} section = 1; if verbose fprintf('\t%d\n', section); @@ -168,7 +168,7 @@ end pos=0; continue - case '- genes:' + case {'- genes:', '- genes: []'} section = 4; if verbose fprintf('\t%d\n', section); @@ -234,7 +234,7 @@ model.annotation.email = tline_value; case 'organization' model.annotation.organization = tline_value; - case 'geckoLight' + case {'geckoLight', 'gecko_light'} if strcmp(tline_value,'true') model.ec.geckoLight = true; end @@ -267,8 +267,12 @@ model = readFieldValue(model, 'inchis', tline_value, pos); readList=''; miriamKey=''; case 'smiles' - model = readFieldValue(model, 'metSmiles', tline_value, pos); - readList=''; miriamKey=''; + if isempty(tline_value) + readList = 'smilesInAnnotation'; + else + model = readFieldValue(model, 'metSmiles', tline_value, pos); + readList=''; miriamKey=''; + end case 'deltaG' model = readFieldValue(model, 'metDeltaG', tline_value, pos); readList=''; miriamKey=''; @@ -279,6 +283,16 @@ readList = 'annotation'; otherwise switch readList + case 'smilesInAnnotation' + %As for ec-code on a reaction: a bare list item + %continues the list, anything else ends it. + if isempty(tline_key) + model = readFieldValue(model, 'metSmiles', ... + regexprep(tline_value,'^ +- "?(.*)"?$','$1'), pos); + else + readList = 'annotation'; + [metMiriams, miriamKey, metMirNo] = gatherAnnotation(pos,metMiriams,tline_key,tline_value,miriamKey,metMirNo); + end case 'annotation' [metMiriams, miriamKey, metMirNo] = gatherAnnotation(pos,metMiriams,tline_key,tline_value,miriamKey,metMirNo); otherwise @@ -309,7 +323,7 @@ case 'gene_reaction_rule' model = readFieldValue(model, 'grRules', tline_value, pos); readList=''; miriamKey=''; - case 'rxnNotes' + case {'rxnNotes', 'notes'} model = readFieldValue(model, 'rxnNotes', tline_value, pos); readList=''; miriamKey=''; case 'rxnFrom' @@ -343,6 +357,19 @@ end case 'metabolites' readList = 'equation'; + case 'ec-code' + %EC numbers emitted inside the annotation block rather than + %under the legacy top-level 'eccodes' key, which is still + %accepted above. A block list follows when the value is + %empty; it is read under its own state so that gathering can + %resume for the rest of the annotation block once the list + %ends, instead of swallowing every later entry as an EC code. + if isempty(tline_value) + readList = 'eccodesInAnnotation'; + else + eccodes(ecCodeNo,1:2)={pos,tline_value}; + ecCodeNo=ecCodeNo+1; + end case 'annotation' readList = 'annotation'; @@ -351,6 +378,18 @@ case 'eccodes' eccodes(ecCodeNo,1:2)={pos,regexprep(tline_value,'^ +- "?(.*)"?$','$1')}; ecCodeNo=ecCodeNo+1; + case 'eccodesInAnnotation' + %A bare list item ("- 1.1.2.4") has no colon, so + %tline_key is empty and the list continues. Anything + %else is the annotation block's next entry, so the + %list has ended. + if isempty(tline_key) + eccodes(ecCodeNo,1:2)={pos,regexprep(tline_value,'^ +- "?(.*)"?$','$1')}; + ecCodeNo=ecCodeNo+1; + else + readList = 'annotation'; + [rxnMiriams, miriamKey,rxnMirNo] = gatherAnnotation(pos,rxnMiriams,tline_key,tline_value,miriamKey,rxnMirNo); + end case 'subsystem' subSystems(subSysNo,1:2)={pos,regexprep(tline_value,'^ +- "?(.*)"?$','$1')}; subSysNo=subSysNo+1; diff --git a/readme/RAVEN_structure_fields.csv b/readme/RAVEN_structure_fields.csv index f471abf7b..4c95e19cb 100644 --- a/readme/RAVEN_structure_fields.csv +++ b/readme/RAVEN_structure_fields.csv @@ -49,7 +49,7 @@ geneMiriams;;Struct;Structure with MIRIAM information about the genes;Non-mandat ;value;Column Cell Array of Strings;The ID for particular MIRIAM entry;Non-mandatory field;P12345; geneShortNames;;Column Cell Array of Strings;Gene alternative names (e.g. ERG10);Non-mandatory field;ERG10; metNames;;Column Cell Array of Strings;Metabolite names;Non-mandatory field;alpha-D-glucose; -metComps;;Column Vector of Doubles;Compartments for metabolites;Mandatory field;1;Shows the index of affiliated compartment in comps, not the compartment abbreviation +metComps;;Column Vector of Integers;Compartments for metabolites;Mandatory field;1;Written as an integer so a reader that converts the 1-based index into a position gets one; MATLAB stores an integer-valued double narrowed to an integer type anyway, and both classes index comps identically. inchis;;Column Cell Array of Strings;InChI-codes for metabolites;Non-mandatory field;1/C6H12O6/c7-1-2-3(8)4(9)5(10)6(11)12-2/h2-11H,1H2/t2-,3-,4+,5-,6+/m1/s1; metFormulas;;Column Cell Array of Strings;Chemical formulas for metabolites;Non-mandatory field;C6H12O6; metMiriams;;Struct;Structure with MIRIAM information about the metabolites;Non-mandatory field;;If there are multiple miriams, 'name' indexes must match 'value' indexes