From 34052a04ac17b480df3037d92983bfeccef34ef9 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 15 Oct 2026 13:25:53 +0200 Subject: [PATCH] Read the model adapter from model_adapter.toml, the file geckopy reads A project's parameters can now live in model_adapter.toml, so one file configures GECKO and geckopy. ModelAdapterManager.getAdapter and setDefault take the project folder or the TOML file and return a TomlAdapter; adapter classes keep working as before. - readAdapterToml and writeAdapterToml read and write the subset of TOML that adapters use, with an error naming the line for anything else. - loadAdapterParams maps the snake_case keys to GECKO's field names and fills geckopy's defaults, both from adapterParameters.tsv, and rejects unknown keys. Evotune source groups become kcatSources and sigma0logSource. - TomlAdapter takes spontaneous reactions from [spontaneous] and the anaerobic and protein-content functions from [matlab]; parameters only GECKO uses go under [matlab.params]. A class that needs more code can subclass TomlAdapter. - convertAdapterToToml writes an existing adapter class to model_adapter.toml and lists the methods it cannot carry over. - createGECKOproject writes model_adapter.toml by default ('adapterFormat','class' for the old template). - The full_ecModel tutorial uses model_adapter.toml; HumanGEMAdapter in light_ecModel subclasses TomlAdapter. It now marks all 14 reactions that data/spontaneousReactions.tsv flags as spontaneous; the old class found one. Original-commit: https://github.com/SysBioChalmers/GECKO/commit/5b32cf61053c8cf0f611dc2311f4822ec290f820 --- src/model_adapter/ModelAdapterManager.m | 21 +- src/model_adapter/TomlAdapter.m | 113 +++++ src/model_adapter/adapterParameters.tsv | 36 ++ src/model_adapter/adapterTemplate.toml | 113 +++++ src/model_adapter/convertAdapterToToml.m | 228 +++++++++ src/model_adapter/loadAdapterParams.m | 249 +++++++++ src/model_adapter/readAdapterToml.m | 471 ++++++++++++++++++ src/model_adapter/writeAdapterToml.m | 143 ++++++ src/utilities/createGECKOproject.m | 26 +- test/VerificationMatrix.tsv | 11 +- .../manual_tests/TestYeastGEMSimpleWorkflow.m | 2 +- test/unit_tests/adapterTomlTests.m | 237 +++++++++ tutorials/full_ecModel/YeastGEMAdapter.m | 89 ---- tutorials/full_ecModel/code/plotlightVSfull.m | 2 +- tutorials/full_ecModel/model_adapter.toml | 85 ++++ tutorials/full_ecModel/protocol.m | 16 +- tutorials/light_ecModel/HumanGEMAdapter.m | 82 +-- tutorials/light_ecModel/model_adapter.toml | 58 +++ 18 files changed, 1804 insertions(+), 178 deletions(-) create mode 100644 src/model_adapter/TomlAdapter.m create mode 100644 src/model_adapter/adapterParameters.tsv create mode 100644 src/model_adapter/adapterTemplate.toml create mode 100644 src/model_adapter/convertAdapterToToml.m create mode 100644 src/model_adapter/loadAdapterParams.m create mode 100644 src/model_adapter/readAdapterToml.m create mode 100644 src/model_adapter/writeAdapterToml.m create mode 100644 test/unit_tests/adapterTomlTests.m delete mode 100644 tutorials/full_ecModel/YeastGEMAdapter.m create mode 100644 tutorials/full_ecModel/model_adapter.toml create mode 100644 tutorials/light_ecModel/model_adapter.toml diff --git a/src/model_adapter/ModelAdapterManager.m b/src/model_adapter/ModelAdapterManager.m index e50fcdc1..a75b0caa 100644 --- a/src/model_adapter/ModelAdapterManager.m +++ b/src/model_adapter/ModelAdapterManager.m @@ -4,6 +4,10 @@ % stores a process-wide default adapter, so that GECKO functions can % retrieve the active adapter without it being passed explicitly. % +% An adapter is loaded from model_adapter.toml (give the file or the +% project folder; see TomlAdapter) or from an adapter class (give its .m +% file). +% % Notes % ----- % All methods are static. The default adapter is held in a persistent @@ -20,9 +24,22 @@ end ModelAdapterManager.showGECKO4Notice(); + adapterPath = char(adapterPath); + if isfolder(adapterPath) + tomlFile = fullfile(adapterPath, 'model_adapter.toml'); + if ~isfile(tomlFile) + error('No model_adapter.toml found in %s.', adapterPath); + end + adapter = TomlAdapter(tomlFile); + return + end [adapterFolder, adapterClassName, extension] = fileparts(adapterPath); - if ~strcmp(extension, '.m') - error('Please provide the full path to the adapter file, including the file extension.'); + if strcmp(extension, '.toml') + adapter = TomlAdapter(adapterPath); + return + elseif ~strcmp(extension, '.m') + error(['Please provide the project folder, its model_adapter.toml, or the ' ... + 'full path to the adapter class file, including the .m extension.']); else s = pathsep; pathStr = [s, path, s]; diff --git a/src/model_adapter/TomlAdapter.m b/src/model_adapter/TomlAdapter.m new file mode 100644 index 00000000..28b7fef0 --- /dev/null +++ b/src/model_adapter/TomlAdapter.m @@ -0,0 +1,113 @@ +% TomlAdapter Model adapter configured by model_adapter.toml. +% +% The model adapter for a project whose parameters are in +% model_adapter.toml, the configuration file that GECKO and geckopy share. +% No adapter class has to be written: what used to be code in an adapter +% class is set in the TOML file. +% +% - Spontaneous reactions: [spontaneous] name_contains (reactions whose +% name contains the text) and/or file (a TSV of reaction identifiers, +% optionally with a second 0/1 column). +% - Anaerobic conditions and protein content during kcat tuning: +% [matlab] make_model_anaerobic and change_protein_biomass name MATLAB +% functions, for instance in the project's code folder, which is added +% to the MATLAB path. Each takes the ecModel (and, for +% change_protein_biomass, the protein content) and returns the ecModel. +% +% A project that needs other custom code can keep an adapter class and +% read its parameters from the TOML file in its constructor: +% +% obj.params = loadAdapterParams(fullfile(projectPath,'model_adapter.toml')); +% +% Examples +% -------- +% adapter = TomlAdapter(fullfile(projectPath,'model_adapter.toml')); +% % or, as the default adapter: +% ModelAdapterManager.setDefault(projectPath); +% +% See also +% -------- +% ModelAdapter, ModelAdapterManager, loadAdapterParams, convertAdapterToToml +classdef TomlAdapter < ModelAdapter + methods + function obj = TomlAdapter(tomlFile) + % TomlAdapter Load the adapter from model_adapter.toml. + % + % Parameters + % ---------- + % tomlFile : char + % path to model_adapter.toml, or to the project folder that + % holds it. + obj.params = loadAdapterParams(tomlFile); + codeFolder = fullfile(obj.params.path, 'code'); + if isfolder(codeFolder) + addpath(codeFolder); + end + end + + function [spont,spontRxnNames] = getSpontaneousReactions(obj,model) + % getSpontaneousReactions Reactions set as spontaneous in [spontaneous]. + % + % Returns + % ------- + % spont : logical + % true for each reaction of model that is spontaneous. + % spontRxnNames : cell + % identifiers of those reactions. + spont = false(numel(model.rxns),1); + if isfield(obj.params,'spontaneous') + s = obj.params.spontaneous; + if isfield(s,'name_contains') && ~isempty(s.name_contains) + spont = spont | contains(model.rxnNames, s.name_contains); + end + if isfield(s,'file') && ~isempty(s.file) + spont = spont | ismember(model.rxns, readSpontaneousFile(s.file)); + end + end + spontRxnNames = model.rxns(spont); + end + + function ecModel = makeModelAnaerobic(obj,ecModel) + % makeModelAnaerobic Run the function named in [matlab] make_model_anaerobic. + fun = hookName(obj, 'makeModelAnaerobic'); + if ~isempty(fun) + ecModel = feval(fun, ecModel); + end + end + + function ecModel = changeProteinBiomass(obj,ecModel,Ptot) + % changeProteinBiomass Run the function named in [matlab] change_protein_biomass. + fun = hookName(obj, 'changeProteinBiomass'); + if ~isempty(fun) + ecModel = feval(fun, ecModel, Ptot); + end + end + end +end + +function fun = hookName(obj, field) +fun = ''; +if isfield(obj.params,'matlab') && isfield(obj.params.matlab, field) + fun = obj.params.matlab.(field); +end +end + +function rxns = readSpontaneousFile(file) +% Reaction identifiers in the first column; with a second column, only +% rows where it is 1 or true. The first row is a header. +fid = fopen(file, 'r'); +if fid < 0 + error('GECKO:TomlAdapter', 'Cannot read the spontaneous reactions file %s', file); +end +header = fgetl(fid); +nCols = numel(strsplit(header, '\t')); +if nCols >= 2 + cols = textscan(fid, ['%q %q' repmat(' %*q', 1, nCols-2)], 'Delimiter', '\t'); + keep = ismember(lower(strtrim(cols{2})), {'1','true'}); + rxns = cols{1}(keep); +else + cols = textscan(fid, '%q', 'Delimiter', '\t'); + rxns = cols{1}; +end +fclose(fid); +end diff --git a/src/model_adapter/adapterParameters.tsv b/src/model_adapter/adapterParameters.tsv new file mode 100644 index 00000000..c4c45f80 --- /dev/null +++ b/src/model_adapter/adapterParameters.tsv @@ -0,0 +1,36 @@ +toml_key matlab_field default description +conv_gem convGEM Path to the conventional GEM, relative to the project folder +org_name org_name Scientific name of the organism +sigma sigma 0.5 Average enzyme saturation factor +p_tot Ptot 0.5 Total protein content [g protein/gDW] +f f 0.5 Fraction of enzymes in the model [g enzyme/g protein] +gr_exp gR_exp 0.41 Growth rate the model should reach when not limited by nutrient uptake [1/h] +c_source c_source "" Exchange reaction of the main carbon source +bio_rxn bioRxn "" Biomass pseudoreaction +enzyme_comp enzyme_comp "cytoplasm" Compartment of the protein pseudometabolites +kcat_aggregate_brenda kcatAggregation "max" BRENDA aggregate read by fuzzyKcatMatching: max or median +kcat_aggregate_candidates "max" Aggregation of several candidate kcats for one reaction (geckopy) +kcat_aggregate_isozymes "mean" Aggregation of isozyme kcats when filling a missing value (geckopy) +kegg.id kegg.ID "sce" KEGG organism code +kegg.gene_id kegg.geneID "kegg" Field in the KEGG entry that matches model gene identifiers +uniprot.type uniprot.type "taxonomy" Whether uniprot.id is a taxonomy or a proteome identifier +uniprot.id uniprot.ID "559292" Taxonomy or proteome identifier +uniprot.gene_id_field uniprot.geneIDfield "gene_oln" UniProt field that matches model gene identifiers +uniprot.reviewed uniprot.reviewed false Use only reviewed (Swiss-Prot) entries +complex.taxonomic_id complex.taxonomicID Taxonomy identifier for Complex Portal +spontaneous.name_contains spontaneous.name_contains Reactions whose name contains this text are spontaneous +spontaneous.file spontaneous.file TSV of spontaneous reactions, relative to the project folder: reaction identifiers, optionally a 0/1 column +evotune.sigma0_log_default evotune.sigma0logDefault 0.5 Prior standard deviation (log) for kcats whose source is in no group +evotune.source_groups evotune.kcatSources Trust tiers: group name to the ec.source values it covers +evotune.sigma0_log_source evotune.sigma0logSource Prior standard deviation (log) per group +evotune.rmse_threshold evotune.rmseThreshold 0.2 Stop when the best RMSE reaches this +evotune.max_generations evotune.maxGenerations 150 Maximum number of CMA-ES generations +evotune.max_growth_weight evotune.maxGrowthWeight 1.0 Weight of the maximum-growth RMSE against the flux RMSE +evotune.prior_penalty_weight evotune.priorPenaltyWeight 0.03 Weight of the prior term in the search objective +evotune.tie_isozymes evotune.tieIsozymes true Give isozyme copies of one reaction a single kcat +okp.method okp.method "CataPro" OpenKineticsPredictor method +okp.handle_long_sequences okp.handleLongSequences "truncate" How OpenKineticsPredictor handles sequences that are too long +okp.include_similarity_columns okp.includeSimilarityColumns true Return similarity to the method's training data +okp.canonicalize_substrates okp.canonicalizeSubstrates true Canonicalize substrate SMILES before prediction +matlab.make_model_anaerobic matlab.makeModelAnaerobic MATLAB function that constrains the ecModel to anaerobic growth +matlab.change_protein_biomass matlab.changeProteinBiomass MATLAB function that changes the protein content of the biomass diff --git a/src/model_adapter/adapterTemplate.toml b/src/model_adapter/adapterTemplate.toml new file mode 100644 index 00000000..f75c2143 --- /dev/null +++ b/src/model_adapter/adapterTemplate.toml @@ -0,0 +1,113 @@ +# Model adapter for KEY_NAME. +# Read by GECKO (ModelAdapterManager.setDefault on the project folder) and +# by geckopy (ModelAdapter.from_folder). Paths are relative to the project +# folder. Keys that are left out get the default listed, with their MATLAB +# names, in GECKO's src/model_adapter/adapterParameters.tsv. + +# Conventional GEM that this ecModel is built from. +conv_gem = "models/yourModel.xml" + +# Scientific name of the organism. +org_name = "genus species" + +# Average enzyme saturation factor. +sigma = 0.5 +# Total protein content in the cell [g protein/gDW]. +p_tot = 0.5 +# Fraction of enzymes in the model [g enzyme/g protein]. +f = 0.5 +# Growth rate the model should reach when not limited by nutrient uptake, +# for example the maximum growth rate [1/h]. +gr_exp = 0.41 + +# Exchange reaction of glucose, or another preferred carbon source. +c_source = "r_1714" +# Biomass pseudoreaction. +bio_rxn = "r_4041" +# Compartment of the protein pseudometabolites. They can all sit in one +# compartment, whichever compartments the enzymes act in; typically the +# cytoplasm. +enzyme_comp = "cytoplasm" + +# Which BRENDA aggregate fuzzyKcatMatching reads for each (EC, substrate, +# organism): "max" (the highest turnover reported) or "median". +kcat_aggregate_brenda = "max" + +[kegg] +# KEGG organism code, from https://www.genome.jp/kegg/catalog/org_list.html +id = "sce" +# Field for the KEGG gene identifier; it should match the gene identifiers +# of the model. "kegg" takes the KEGG entry identifier (for example +# YER023W: https://www.genome.jp/dbget-bin/www_bget?sce:YER023W). Gene +# identifiers from the "Other DBs" section of the KEGG page can also be +# used, for example "NCBI-GeneID", "UniProt" or "Ensembl"; not all are +# available for every organism and gene. +gene_id = "kegg" + +[uniprot] +# "proteome" (https://www.uniprot.org/proteomes/) or "taxonomy" +# (https://www.uniprot.org/taxonomy), and the matching identifier. +type = "taxonomy" +id = "559292" +# UniProt field that matches the gene identifiers of the model: one of the +# "Returned Field" entries under "Names & Taxonomy" at +# https://www.uniprot.org/help/return_fields +gene_id_field = "gene_oln" +# Use only reviewed (Swiss-Prot) entries. They have the highest confidence, +# but coverage can be (very) low for non-model organisms. +reviewed = false + +[complex] +# Taxonomy identifier for Complex Portal; leave out to skip it. +# taxonomic_id = 559292 + +# Reactions without an enzyme: those whose name contains name_contains, +# and those listed in file (a TSV with a header row, reaction identifiers +# in the first column and optionally a 0/1 column). +[spontaneous] +name_contains = "spontaneous" +# file = "data/spontaneousReactions.tsv" + +# Hyperparameters for evotune (CMA-ES) kcat tuning. +[evotune] +# Prior standard deviation (log) for kcats whose source is in no group. +sigma0_log_default = 0.5 +# Weight on the max-growth RMSE against the flux RMSE: +# (rmse_flux + w * rmse_max_growth) / (w + 1). At 1 both count equally. +max_growth_weight = 1.0 +# Weight on a prior term in the search objective, +# rmse + w * mean((log(k/k0) / sigma0_log)^2). Keeps large corrections +# reproducible across seeds; 0 scores on RMSE alone. +prior_penalty_weight = 0.03 +# Give isozyme copies of one reaction (sharing a prior value and a source) +# a single shared kcat, so the search cannot invent a distinction the kcat +# assignment never made. +tie_isozymes = true +# Stop when the RMSE falls below this. +rmse_threshold = 0.2 +# Maximum number of CMA-ES generations. +max_generations = 150 + +# Kcat sources, from least to most trusted, each with its prior width +# (lower is more trusted). +[evotune.source_groups.dlkcat] +sources = ["dlkcat"] + +[evotune.source_groups.brenda] +sources = ["brenda"] + +[evotune.source_groups.custom] +sources = ["custom"] + +[evotune.sigma0_log_source] +dlkcat = 0.4 +brenda = 0.2 +custom = 0.1 + +# Read by GECKO only. Names of MATLAB functions, for instance in the code +# folder, used during kcat tuning: one that constrains the ecModel to +# anaerobic growth, and one that changes the protein content of the +# biomass (it receives the ecModel and the protein content). +[matlab] +# make_model_anaerobic = "anaerobicModel" +# change_protein_biomass = "scaleBiomassProtein" diff --git a/src/model_adapter/convertAdapterToToml.m b/src/model_adapter/convertAdapterToToml.m new file mode 100644 index 00000000..fe24c636 --- /dev/null +++ b/src/model_adapter/convertAdapterToToml.m @@ -0,0 +1,228 @@ +function report = convertAdapterToToml(adapter, tomlFile) +% convertAdapterToToml Write the parameters of an adapter class to model_adapter.toml. +% +% Converts an existing ModelAdapter class (GECKO 3 or 4 style) to +% model_adapter.toml, the configuration file shared by GECKO and geckopy. +% Parameter names are translated with adapterParameters.tsv (for instance +% Ptot becomes p_tot), paths inside the project folder are written +% relative to it, and parameters without a geckopy name are kept under +% [matlab.params], so that loadAdapterParams gives back the same params. +% +% Methods that the adapter class overrides cannot be written to TOML. +% They are listed in the report, with how to carry them over: +% - getSpontaneousReactions: set [spontaneous] name_contains or file; +% - makeModelAnaerobic, changeProteinBiomass: set [matlab] +% make_model_anaerobic or change_protein_biomass to a function name; +% - other methods: keep the adapter class and let its constructor call +% loadAdapterParams. +% +% Parameters +% ---------- +% adapter : ModelAdapter or char +% an adapter object, or the path to its .m file. +% tomlFile : char, optional +% file to write (default: model_adapter.toml in the adapter's +% params.path). An existing file is not overwritten. +% +% Returns +% ------- +% report : struct +% - file: the file written; +% - unmapped: parameters written under [matlab.params]; +% - methods: overridden methods that need carrying over by hand. +% +% Examples +% -------- +% report = convertAdapterToToml(fullfile(projectPath,'MyGEMAdapter.m')); +% adapter = TomlAdapter(projectPath); +% +% See also +% -------- +% TomlAdapter, loadAdapterParams, writeAdapterToml + +if ischar(adapter) || isstring(adapter) + adapter = ModelAdapterManager.getAdapter(char(adapter)); +end +params = adapter.params; +if isfield(params,'path') && isstruct(params.path) + params.path = params.path.Name; % a struct path, as some old adapters keep it +end +projectPath = ''; +if isfield(params,'path') + projectPath = char(params.path); +end +if nargin < 2 || isempty(tomlFile) + tomlFile = fullfile(projectPath, 'model_adapter.toml'); +end +if isfile(tomlFile) + error('GECKO:convertAdapterToToml', '%s exists already; it is not overwritten.', tomlFile); +end + +map = readMap(); +data = struct(); +used = {}; +for k = 1:numel(map.tomlKey) + field = map.matlabField{k}; + if isempty(field) || any(strcmp(field, {'evotune.kcatSources','evotune.sigma0logSource'})) + continue + end + [found, value] = getDotted(params, field); + if ~found || (isempty(value) && ~ischar(value)) + if found + used{end+1} = field; %#ok + end + continue + end + used{end+1} = field; %#ok + if any(strcmp(map.tomlKey{k}, {'conv_gem','spontaneous.file'})) + value = relativePath(char(value), projectPath); + end + data = setDotted(data, map.tomlKey{k}, value); +end +used{end+1} = 'path'; + +% Trust tiers: one group per ec.source value. +if isfield(params,'evotune') && isfield(params.evotune,'kcatSources') + used = [used, {'evotune.kcatSources','evotune.sigma0logSource'}]; + sources = params.evotune.kcatSources; + widths = params.evotune.sigma0logSource; + groups = struct(); + sigma0 = struct(); + for k = 1:numel(sources) + name = matlab.lang.makeValidName(sources{k}); + groups.(name).sources = sources(k); + sigma0.(name) = widths(k); + end + data.evotune.source_groups = groups; + data.evotune.sigma0_log_source = sigma0; +end + +% Everything else is kept for GECKO under [matlab.params]. +extra = unusedFields(params, '', used); +report.unmapped = {}; +if ~isempty(fieldnames(extra)) + data.matlab.params = extra; + report.unmapped = flattenNames(extra, ''); +end + +report.methods = overriddenMethods(adapter); +report.file = tomlFile; + +header = {sprintf('Model adapter for %s, converted from %s.', ... + getOr(params,'org_name','this project'), class(adapter)), ... + 'Shared by GECKO (MATLAB) and geckopy; see adapterParameters.tsv in GECKO for every key.'}; +writeAdapterToml(data, tomlFile, header); + +fprintf('Wrote %s.\n', tomlFile); +if ~isempty(report.unmapped) + fprintf('Kept for GECKO only, under [matlab.params]: %s\n', strjoin(report.unmapped, ', ')); +end +if ~isempty(report.methods) + printOrange(sprintf(['%s overrides methods that a TOML file cannot hold: %s.\n' ... + 'Set [spontaneous] for getSpontaneousReactions, [matlab] function names for\n' ... + 'makeModelAnaerobic and changeProteinBiomass, or keep the class and let its\n' ... + 'constructor call loadAdapterParams.\n'], class(adapter), strjoin(report.methods, ', '))); +end +end + +function map = readMap() +file = fullfile(fileparts(mfilename('fullpath')), 'adapterParameters.tsv'); +fid = fopen(file, 'r'); +cols = textscan(fid, '%s %s %s %s', 'Delimiter', '\t', 'HeaderLines', 1, 'Whitespace', ''); +fclose(fid); +map.tomlKey = cols{1}; +map.matlabField = cols{2}; +end + +function names = overriddenMethods(adapter) +% Methods defined by the adapter class itself rather than ModelAdapter. +names = {}; +mc = metaclass(adapter); +for k = 1:numel(mc.MethodList) + m = mc.MethodList(k); + if strcmp(m.DefiningClass.Name, mc.Name) && ~strcmp(m.Name, mc.Name) && ~m.Static + names{end+1} = m.Name; %#ok + end +end +names = sort(names); +end + +function extra = unusedFields(s, prefix, used) +extra = struct(); +fields = fieldnames(s); +for k = 1:numel(fields) + name = fields{k}; + if ~isempty(prefix) + name = [prefix '.' name]; + end + if any(strcmp(name, used)) + continue + end + value = s.(fields{k}); + if isstruct(value) && isscalar(value) && any(startsWith(used, [name '.'])) + sub = unusedFields(value, name, used); + if ~isempty(fieldnames(sub)) + extra.(fields{k}) = sub; + end + elseif ~(isempty(value) && ~ischar(value) && ~iscell(value)) + extra.(fields{k}) = value; + end +end +end + +function names = flattenNames(s, prefix) +names = {}; +fields = fieldnames(s); +for k = 1:numel(fields) + name = fields{k}; + if ~isempty(prefix) + name = [prefix '.' name]; + end + if isstruct(s.(fields{k})) + names = [names, flattenNames(s.(fields{k}), name)]; %#ok + else + names{end+1} = name; %#ok + end +end +end + +function p = relativePath(p, projectPath) +p = strrep(p, '\', '/'); +base = strrep(projectPath, '\', '/'); +if ~isempty(base) && startsWith(p, [base '/']) + p = p(numel(base)+2:end); +end +end + +function [found, value] = getDotted(s, key) +parts = strsplit(key, '.'); +value = s; +found = true; +for k = 1:numel(parts) + if ~isstruct(value) || ~isfield(value, parts{k}) + found = false; + value = []; + return + end + value = value.(parts{k}); +end +end + +function s = setDotted(s, key, value) +parts = strsplit(key, '.'); +if numel(parts) == 1 + s.(parts{1}) = value; + return +end +if ~isfield(s, parts{1}) || ~isstruct(s.(parts{1})) + s.(parts{1}) = struct(); +end +s.(parts{1}) = setDotted(s.(parts{1}), strjoin(parts(2:end), '.'), value); +end + +function v = getOr(s, field, default) +v = default; +if isfield(s, field) && ~isempty(s.(field)) + v = s.(field); +end +end diff --git a/src/model_adapter/loadAdapterParams.m b/src/model_adapter/loadAdapterParams.m new file mode 100644 index 00000000..94cd4bbd --- /dev/null +++ b/src/model_adapter/loadAdapterParams.m @@ -0,0 +1,249 @@ +function params = loadAdapterParams(tomlFile) +% loadAdapterParams Read model adapter parameters from a TOML file. +% +% Reads model_adapter.toml, the configuration file shared by GECKO and +% geckopy, and returns the parameters with GECKO's field names, so that +% they can be used as ModelAdapter.params. Keys are translated with +% adapterParameters.tsv (for instance p_tot becomes Ptot); keys that are +% not given get the default listed there, which is geckopy's default. +% +% Parameters +% ---------- +% tomlFile : char +% path to model_adapter.toml, or to the project folder that holds it. +% +% Returns +% ------- +% params : struct +% the adapter parameters, as in ModelAdapter.params: +% - path: the project folder (the folder of the TOML file); +% - convGEM and spontaneous.file: absolute paths; +% - evotune.kcatSources and evotune.sigma0logSource: built from the +% evotune.source_groups and evotune.sigma0_log_source tables, one +% entry per ec.source value; +% - fields under [matlab.params] are copied as they are, for +% parameters that only GECKO uses. +% +% Examples +% -------- +% params = loadAdapterParams(fullfile(projectPath,'model_adapter.toml')); +% +% See also +% -------- +% TomlAdapter, readAdapterToml, convertAdapterToToml + +if isfolder(tomlFile) + tomlFile = fullfile(tomlFile, 'model_adapter.toml'); +end +if ~isfile(tomlFile) + error('GECKO:loadAdapterParams', 'No model adapter file found at %s', tomlFile); +end +tomlFile = char(tomlFile); +[projectPath, ~, ~] = fileparts(getAbsolutePath(tomlFile)); + +data = readAdapterToml(tomlFile); +map = adapterParameterMap(); + +% Unknown keys are errors, as in geckopy, so that a typo is not ignored. +known = [map.tomlKey; {'evotune.source_groups'; 'evotune.sigma0_log_source'}]; +checkKnownKeys(data, '', known, tomlFile); + +params = struct(); +params.path = projectPath; +for k = 1:numel(map.tomlKey) + key = map.tomlKey{k}; + if any(strcmp(key, {'evotune.source_groups', 'evotune.sigma0_log_source'})) + continue + end + [found, value] = getDotted(data, key); + if ~found + if isempty(map.default{k}) + continue + end + value = map.default{k}; + end + if isempty(map.matlabField{k}) + field = key; % a key that only geckopy reads keeps its name + else + field = map.matlabField{k}; + end + params = setDotted(params, field, value); +end + +if ~isfield(params, 'complex') + params.complex.taxonomicID = []; % geckopy's default: no Complex Portal query +end + +for f = {'convGEM'} + if isfield(params, f{1}) + params.(f{1}) = projectFile(projectPath, params.(f{1})); + end +end +if isfield(params, 'spontaneous') && isfield(params.spontaneous, 'file') + params.spontaneous.file = projectFile(projectPath, params.spontaneous.file); +end + +params = sourceGroups(params, data, tomlFile); + +if isfield(data, 'matlab') && isfield(data.matlab, 'params') + params = mergeStruct(params, data.matlab.params); +end +end + +function map = adapterParameterMap() +% Reads adapterParameters.tsv into toml keys, MATLAB fields and defaults. +file = fullfile(fileparts(mfilename('fullpath')), 'adapterParameters.tsv'); +fid = fopen(file, 'r'); +cols = textscan(fid, '%s %s %s %s', 'Delimiter', '\t', 'HeaderLines', 1, ... + 'Whitespace', '', 'ReturnOnError', false); +fclose(fid); +map.tomlKey = cols{1}; +map.matlabField = cols{2}; +map.default = cell(size(cols{3})); +for k = 1:numel(cols{3}) + if ~isempty(cols{3}{k}) + map.default{k} = parseDefault(cols{3}{k}); + end +end +end + +function value = parseDefault(text) +% A default is written as a TOML value. +tmp = [tempname '.toml']; +fid = fopen(tmp, 'w'); +fprintf(fid, 'x = %s\n', text); +fclose(fid); +cleanup = onCleanup(@() delete(tmp)); +data = readAdapterToml(tmp); +value = data.x; +if ischar(value) && isempty(value) + value = ''; +end +end + +function checkKnownKeys(s, prefix, known, tomlFile) +fields = fieldnames(s); +for k = 1:numel(fields) + key = fields{k}; + if ~isempty(prefix) + key = [prefix '.' key]; + end + if strcmp(key, 'matlab.params') + continue % free-form parameters for GECKO only + end + if any(strcmp(key, known)) || any(startsWith(known, [key '.'])) + if isstruct(s.(fields{k})) && ~any(strcmp(key, known)) + checkKnownKeys(s.(fields{k}), key, known, tomlFile); + end + continue + end + if any(startsWith(key, {'evotune.source_groups.', 'evotune.sigma0_log_source.'})) + continue + end + error('GECKO:loadAdapterParams', ['Unknown key "%s" in %s. Parameters that only ' ... + 'GECKO uses go under [matlab.params].'], key, tomlFile); +end +end + +function params = sourceGroups(params, data, tomlFile) +% evotune.source_groups and evotune.sigma0_log_source become the +% kcatSources and sigma0logSource lists that evotune reads: one entry per +% ec.source value, with its group's prior width. +groups = struct(); +widths = struct(); +if isfield(data, 'evotune') && isfield(data.evotune, 'source_groups') + groups = data.evotune.source_groups; + if isfield(data.evotune, 'sigma0_log_source') + widths = data.evotune.sigma0_log_source; + end +else + % geckopy's defaults (EvotuneParams.source_groups, sigma0_log_source). + groups.dlkcat.sources = {'dlkcat'}; + groups.brenda.sources = {'brenda'}; + groups.custom.sources = {'custom'}; + widths.dlkcat = 0.4; + widths.brenda = 0.2; + widths.custom = 0.1; +end +names = fieldnames(groups); +if ~isempty(setxor(names, fieldnames(widths))) + error('GECKO:loadAdapterParams', ['%s: evotune.sigma0_log_source must have ' ... + 'the same group names as evotune.source_groups'], tomlFile); +end +kcatSources = {}; +sigma0logSource = []; +for k = 1:numel(names) + group = groups.(names{k}); + sources = {}; + if isfield(group, 'sources') + sources = group.sources; + if ischar(sources) + sources = {sources}; + end + end + if isfield(group, 'match_okp') && group.match_okp && isfield(params, 'okp') + sources{end+1} = params.okp.method; %#ok + end + kcatSources = [kcatSources, sources]; %#ok + sigma0logSource = [sigma0logSource; repmat(widths.(names{k}), numel(sources), 1)]; %#ok +end +params.evotune.kcatSources = kcatSources; +params.evotune.sigma0logSource = sigma0logSource; +end + +function [found, value] = getDotted(s, key) +parts = strsplit(key, '.'); +value = s; +found = true; +for k = 1:numel(parts) + if ~isstruct(value) || ~isfield(value, parts{k}) + found = false; + value = []; + return + end + value = value.(parts{k}); +end +end + +function s = setDotted(s, key, value) +parts = strsplit(key, '.'); +if numel(parts) == 1 + s.(parts{1}) = value; + return +end +if ~isfield(s, parts{1}) || ~isstruct(s.(parts{1})) + s.(parts{1}) = struct(); +end +s.(parts{1}) = setDotted(s.(parts{1}), strjoin(parts(2:end), '.'), value); +end + +function s = mergeStruct(s, extra) +fields = fieldnames(extra); +for k = 1:numel(fields) + value = extra.(fields{k}); + if isstruct(value) && isscalar(value) && isfield(s, fields{k}) && isstruct(s.(fields{k})) + s.(fields{k}) = mergeStruct(s.(fields{k}), value); + else + s.(fields{k}) = value; + end +end +end + +function p = projectFile(projectPath, p) +if isempty(p) + return +end +if ~isAbsolute(p) + p = fullfile(projectPath, p); +end +end + +function tf = isAbsolute(p) +tf = startsWith(p, '/') || startsWith(p, '\') || ~isempty(regexp(p, '^[A-Za-z]:', 'once')); +end + +function p = getAbsolutePath(p) +if ~isAbsolute(p) + p = fullfile(pwd, p); +end +end diff --git a/src/model_adapter/readAdapterToml.m b/src/model_adapter/readAdapterToml.m new file mode 100644 index 00000000..eca857f2 --- /dev/null +++ b/src/model_adapter/readAdapterToml.m @@ -0,0 +1,471 @@ +function data = readAdapterToml(filename) +% readAdapterToml Read a model adapter TOML file into a struct. +% +% Reads the subset of TOML that model adapter files use: comments, bare +% keys, basic and literal strings, integers, floats, booleans, arrays of +% these values (also over several lines), tables ([kegg], [a.b.c]) and +% arrays of tables ([[a]]). Anything else, such as dates, inline tables, +% multi-line strings or quoted keys, stops with an error that names the +% line. +% +% Parameters +% ---------- +% filename : char +% path to the TOML file. +% +% Returns +% ------- +% data : struct +% one field per key, with tables as nested structs and arrays of +% tables as struct arrays. Strings are char, numbers double, booleans +% logical, arrays of strings a 1xn cell and arrays of numbers or +% booleans a 1xn vector. An empty array is an empty cell. +% +% Examples +% -------- +% data = readAdapterToml(fullfile(projectPath,'model_adapter.toml')); +% +% See also +% -------- +% writeAdapterToml, loadAdapterParams + +text = fileread(filename); +text = strrep(text, sprintf('\r\n'), newline); +lines = strsplit(text, newline, 'CollapseDelimiters', false); + +data = struct(); +tablePath = {}; % keys of the current table, each element {name, arrayIndex} +arrayNames = {}; % dotted names declared as arrays of tables +i = 0; +while i < numel(lines) + i = i + 1; + lineNo = i; + line = stripComment(lines{i}, filename, lineNo); + if isempty(strtrim(line)) + continue + end + trimmed = strtrim(line); + if startsWith(trimmed, '[[') + if ~endsWith(trimmed, ']]') + tomlError(filename, lineNo, 'array of tables header must end with ]]'); + end + keys = splitKey(trimmed(3:end-2), filename, lineNo); + [data, tablePath] = appendTableArray(data, keys, arrayNames, filename, lineNo); + arrayNames = union(arrayNames, {strjoin(keys, '.')}); + elseif startsWith(trimmed, '[') + if ~endsWith(trimmed, ']') + tomlError(filename, lineNo, 'table header must end with ]'); + end + keys = splitKey(trimmed(2:end-1), filename, lineNo); + tablePath = resolvePath(data, keys, arrayNames); + data = ensureTable(data, tablePath, filename, lineNo); + else + eq = findUnquoted(line, '='); + if isempty(eq) + tomlError(filename, lineNo, 'expected key = value'); + end + keys = splitKey(line(1:eq-1), filename, lineNo); + valueText = strtrim(line(eq+1:end)); + % Arrays may span several lines: read on until the brackets close. + while startsWith(valueText, '[') && bracketDepth(valueText) > 0 + if i >= numel(lines) + tomlError(filename, lineNo, 'array is not closed'); + end + i = i + 1; + valueText = [valueText ' ' strtrim(stripComment(lines{i}, filename, i))]; %#ok + end + value = parseValue(valueText, filename, lineNo); + data = setKey(data, [tablePath, keyPath(keys)], value, filename, lineNo); + end +end +end + +function line = stripComment(line, filename, lineNo) +% Removes a # comment that is not inside a string. +inBasic = false; inLiteral = false; k = 1; +while k <= numel(line) + c = line(k); + if inBasic + if c == '\' + k = k + 1; + elseif c == '"' + inBasic = false; + end + elseif inLiteral + if c == '''' + inLiteral = false; + end + elseif c == '"' + if k + 2 <= numel(line) && strcmp(line(k:k+2), '"""') + tomlError(filename, lineNo, 'multi-line strings are not supported'); + end + inBasic = true; + elseif c == '''' + if k + 2 <= numel(line) && strcmp(line(k:k+2), '''''''') + tomlError(filename, lineNo, 'multi-line strings are not supported'); + end + inLiteral = true; + elseif c == '#' + line = line(1:k-1); + return + end + k = k + 1; +end +end + +function idx = findUnquoted(line, ch) +% Position of the first ch outside a string, or empty. +idx = []; +inBasic = false; inLiteral = false; k = 1; +while k <= numel(line) + c = line(k); + if inBasic + if c == '\' + k = k + 1; + elseif c == '"' + inBasic = false; + end + elseif inLiteral + if c == '''' + inLiteral = false; + end + elseif c == '"' + inBasic = true; + elseif c == '''' + inLiteral = true; + elseif c == ch + idx = k; + return + end + k = k + 1; +end +end + +function depth = bracketDepth(text) +% Open [ minus closed ] outside strings. +depth = 0; inBasic = false; inLiteral = false; k = 1; +while k <= numel(text) + c = text(k); + if inBasic + if c == '\' + k = k + 1; + elseif c == '"' + inBasic = false; + end + elseif inLiteral + if c == '''' + inLiteral = false; + end + elseif c == '"' + inBasic = true; + elseif c == '''' + inLiteral = true; + elseif c == '[' + depth = depth + 1; + elseif c == ']' + depth = depth - 1; + end + k = k + 1; +end +end + +function keys = splitKey(text, filename, lineNo) +% Dotted bare key into its parts. +keys = strtrim(strsplit(strtrim(text), '.')); +for k = 1:numel(keys) + if isempty(regexp(keys{k}, '^[A-Za-z][A-Za-z0-9_]*$', 'once')) + tomlError(filename, lineNo, sprintf(['key "%s" is not supported: use letters, ' ... + 'digits and underscores, starting with a letter'], keys{k})); + end +end +end + +function p = keyPath(keys) +p = cellfun(@(k) {k, 0}, keys, 'UniformOutput', false); +end + +function path = resolvePath(data, keys, arrayNames) +% Path of {name, index} pairs; a key that names an array of tables refers +% to its last element. +path = cell(1, numel(keys)); +for k = 1:numel(keys) + idx = 0; + if ismember(strjoin(keys(1:k), '.'), arrayNames) + idx = numel(getKey(data, [path(1:k-1), {{keys{k}, 0}}])); + end + path{k} = {keys{k}, idx}; +end +end + +function [data, tablePath] = appendTableArray(data, keys, arrayNames, filename, lineNo) +% Appends a new element to the array of tables named by keys. +parent = resolvePath(data, keys(1:end-1), arrayNames); +data = ensureTable(data, parent, filename, lineNo); +name = strjoin(keys, '.'); +existing = getKey(data, [parent, {{keys{end}, 0}}]); +if isempty(existing) || ~ismember(name, arrayNames) + if ~isempty(existing) + tomlError(filename, lineNo, sprintf('"%s" is already a table or value', name)); + end + n = 1; +else + n = numel(existing) + 1; +end +tablePath = [parent, {{keys{end}, n}}]; +data = setTableElement(data, parent, keys{end}, n); +end + +function data = ensureTable(data, path, filename, lineNo) +if isempty(path) + return +end +name = path{1}{1}; idx = path{1}{2}; +if ~isfield(data, name) + data.(name) = struct(); +elseif ~isstruct(data.(name)) + tomlError(filename, lineNo, sprintf('"%s" is already a value, not a table', name)); +end +child = data.(name); +if idx == 0 + child = ensureTable(child, path(2:end), filename, lineNo); + data.(name) = child; +else + element = child(idx); + element = ensureTable(element, path(2:end), filename, lineNo); + data.(name) = assignElement(child, idx, element); +end +end + +function data = setTableElement(data, parent, name, n) +% Adds element n (an empty struct) to the table array name under parent. +if isempty(parent) + if ~isfield(data, name) || isempty(data.(name)) || n == 1 + data.(name) = struct(); + else + data.(name) = assignElement(data.(name), n, struct()); + end + return +end +pname = parent{1}{1}; pidx = parent{1}{2}; +child = data.(pname); +if pidx == 0 + data.(pname) = setTableElement(child, parent(2:end), name, n); +else + element = setTableElement(child(pidx), parent(2:end), name, n); + data.(pname) = assignElement(child, pidx, element); +end +end + +function value = getKey(data, path) +value = []; +node = data; +for k = 1:numel(path) + name = path{k}{1}; idx = path{k}{2}; + if ~isstruct(node) || ~isfield(node, name) + value = []; + return + end + node = node.(name); + if idx > 0 + node = node(idx); + end +end +value = node; +end + +function data = setKey(data, path, value, filename, lineNo) +name = path{1}{1}; idx = path{1}{2}; +if numel(path) == 1 + % [] only fills a field that another element of a table array has; + % TOML itself has no empty value. + if isfield(data, name) && ~isequal(data.(name), []) + tomlError(filename, lineNo, sprintf('key "%s" is defined twice', name)); + end + data.(name) = value; + return +end +if ~isfield(data, name) + data.(name) = struct(); +elseif ~isstruct(data.(name)) + tomlError(filename, lineNo, sprintf('"%s" is already a value, not a table', name)); +end +child = data.(name); +if idx == 0 + data.(name) = setKey(child, path(2:end), value, filename, lineNo); +else + element = setKey(child(idx), path(2:end), value, filename, lineNo); + data.(name) = assignElement(child, idx, element); +end +end + +function s = assignElement(s, idx, element) +% Assigns element to s(idx), adding fields so the struct array stays +% consistent. +if idx > numel(s) || isempty(fieldnames(s)) + if numel(s) == 1 && isempty(fieldnames(s)) && idx == 1 + s = element; + return + end +end +fs = fieldnames(element); +for f = 1:numel(fs) + if ~isfield(s, fs{f}) + [s.(fs{f})] = deal([]); + end +end +fs = fieldnames(s); +for f = 1:numel(fs) + if ~isfield(element, fs{f}) + element.(fs{f}) = []; + end +end +element = orderfields(element, s(1)); +s(idx) = element; +end + +function value = parseValue(text, filename, lineNo) +text = strtrim(text); +if isempty(text) + tomlError(filename, lineNo, 'missing value'); +end +switch text(1) + case '"' + [value, rest] = parseBasicString(text, filename, lineNo); + case '''' + close = find(text(2:end) == '''', 1) + 1; + if isempty(close) + tomlError(filename, lineNo, 'string is not closed'); + end + value = text(2:close-1); + rest = text(close+1:end); + case '[' + [value, rest] = parseArray(text, filename, lineNo); + case '{' + tomlError(filename, lineNo, 'inline tables are not supported; use a [table] instead'); + otherwise + [value, rest] = parseScalar(text, filename, lineNo); +end +if ~isempty(strtrim(rest)) + tomlError(filename, lineNo, sprintf('unexpected text after the value: %s', strtrim(rest))); +end +end + +function [value, rest] = parseBasicString(text, filename, lineNo) +value = ''; +k = 2; +while k <= numel(text) + c = text(k); + if c == '"' + rest = text(k+1:end); + return + elseif c == '\' + if k == numel(text) + break + end + e = text(k+1); + switch e + case 'n', value(end+1) = newline; + case 't', value(end+1) = sprintf('\t'); + case '"', value(end+1) = '"'; + case '\', value(end+1) = '\'; + case 'u' + value(end+1) = char(hex2dec(text(k+2:k+5))); + k = k + 4; + otherwise + tomlError(filename, lineNo, sprintf('escape \\%s is not supported', e)); + end + k = k + 2; + continue + else + value(end+1) = c; %#ok + end + k = k + 1; +end +tomlError(filename, lineNo, 'string is not closed'); +end + +function [value, rest] = parseArray(text, filename, lineNo) +% Parses [a, b, ...]; returns a cell for strings, a vector otherwise. +items = {}; +k = 2; +while true + while k <= numel(text) && (isspace(text(k)) || text(k) == ',') + k = k + 1; + end + if k > numel(text) + tomlError(filename, lineNo, 'array is not closed'); + end + if text(k) == ']' + rest = text(k+1:end); + break + end + if text(k) == '[' + tomlError(filename, lineNo, 'nested arrays are not supported'); + end + % Find the end of this element. + if text(k) == '"' + [item, after] = parseBasicString(text(k:end), filename, lineNo); + k = numel(text) - numel(after) + 1; + elseif text(k) == '''' + close = find(text(k+1:end) == '''', 1) + k; + if isempty(close) + tomlError(filename, lineNo, 'string is not closed'); + end + item = text(k+1:close-1); + k = close + 1; + else + stop = k; + while stop <= numel(text) && text(stop) ~= ',' && text(stop) ~= ']' + stop = stop + 1; + end + [item, after] = parseScalar(strtrim(text(k:stop-1)), filename, lineNo); + if ~isempty(strtrim(after)) + tomlError(filename, lineNo, 'unexpected text in array'); + end + k = stop; + end + items{end+1} = item; %#ok +end +if isempty(items) + value = {}; +elseif all(cellfun(@ischar, items)) + value = items; +elseif all(cellfun(@(x) isnumeric(x) && isscalar(x), items)) + value = cell2mat(items); +elseif all(cellfun(@islogical, items)) + value = cell2mat(items); +else + tomlError(filename, lineNo, 'arrays must hold values of one type'); +end +end + +function [value, rest] = parseScalar(text, filename, lineNo) +token = regexp(text, '^[^\s,\]]+', 'match', 'once'); +rest = text(numel(token)+1:end); +switch token + case 'true' + value = true; + case 'false' + value = false; + case {'inf', '+inf'} + value = Inf; + case '-inf' + value = -Inf; + case {'nan', '+nan', '-nan'} + value = NaN; + otherwise + clean = strrep(token, '_', ''); + if ~isempty(regexp(clean, '^[+-]?(0|[1-9][0-9]*)$', 'once')) || ... + ~isempty(regexp(clean, '^[+-]?(0|[1-9][0-9]*)(\.[0-9]+)?([eE][+-]?[0-9]+)?$', 'once')) + value = str2double(clean); + elseif ~isempty(regexp(token, '^\d{4}-\d{2}-\d{2}', 'once')) + tomlError(filename, lineNo, 'dates are not supported'); + else + tomlError(filename, lineNo, sprintf('cannot read value "%s"; strings need quotes', token)); + end +end +end + +function tomlError(filename, lineNo, message) +error('GECKO:readAdapterToml', '%s, line %d: %s', filename, lineNo, message); +end diff --git a/src/model_adapter/writeAdapterToml.m b/src/model_adapter/writeAdapterToml.m new file mode 100644 index 00000000..68991a4f --- /dev/null +++ b/src/model_adapter/writeAdapterToml.m @@ -0,0 +1,143 @@ +function writeAdapterToml(data, filename, header) +% writeAdapterToml Write a struct as a model adapter TOML file. +% +% Writes the subset of TOML that readAdapterToml reads: values first, +% then each nested struct as a [table] and each struct array as an array +% of tables ([[name]]). +% +% Parameters +% ---------- +% data : struct +% values may be char or string scalars, numeric or logical scalars or +% vectors, cell arrays of char, structs and struct arrays. +% filename : char +% path of the TOML file to write; an existing file is overwritten. +% header : char or cell, optional +% comment lines written at the top of the file (without '#'). +% +% Examples +% -------- +% writeAdapterToml(struct('sigma',0.5), 'model_adapter.toml', ... +% {'Adapter for ecYeastGEM'}); +% +% See also +% -------- +% readAdapterToml, convertAdapterToToml + +if nargin < 3 + header = {}; +elseif ischar(header) || isstring(header) + header = cellstr(header); +end + +out = {}; +for k = 1:numel(header) + out{end+1} = ['# ' header{k}]; %#ok +end +if ~isempty(header) + out{end+1} = ''; +end +out = writeTable(out, data, {}); + +fid = fopen(filename, 'w'); +if fid < 0 + error('GECKO:writeAdapterToml', 'Cannot write %s', filename); +end +fprintf(fid, '%s\n', out{:}); +fclose(fid); +end + +function out = writeTable(out, s, prefix) +fields = fieldnames(s); +tables = {}; +for k = 1:numel(fields) + value = s.(fields{k}); + if isstruct(value) + tables{end+1} = fields{k}; %#ok + elseif isempty(value) && ~iscell(value) && ~ischar(value) + continue % unset value; TOML has no null + else + out{end+1} = [fields{k} ' = ' formatValue(value, strjoin([prefix fields(k)], '.'))]; %#ok + end +end +for k = 1:numel(tables) + value = s.(tables{k}); + name = strjoin([prefix tables(k)], '.'); + if numel(value) == 1 + out{end+1} = ''; %#ok + out{end+1} = ['[' name ']']; %#ok + out = writeTable(out, value, [prefix tables(k)]); + else + for e = 1:numel(value) + out{end+1} = ''; %#ok + out{end+1} = ['[[' name ']]']; %#ok + out = writeTable(out, value(e), [prefix tables(k)]); + end + end +end +end + +function text = formatValue(value, name) +if isstring(value) && isscalar(value) + value = char(value); +end +if ischar(value) + text = quote(value); +elseif iscell(value) + if ~all(cellfun(@(x) ischar(x) || (isstring(x) && isscalar(x)), value)) + error('GECKO:writeAdapterToml', '%s: only cell arrays of text can be written', name); + end + parts = cellfun(@(x) quote(char(x)), value, 'UniformOutput', false); + text = ['[' strjoin(parts(:)', ', ') ']']; +elseif islogical(value) || isnumeric(value) + if isscalar(value) + text = formatScalar(value); + else + parts = arrayfun(@formatScalar, value(:)', 'UniformOutput', false); + text = ['[' strjoin(parts, ', ') ']']; + end +else + error('GECKO:writeAdapterToml', '%s: values of class %s cannot be written', name, class(value)); +end +end + +function text = formatScalar(x) +if islogical(x) + if x + text = 'true'; + else + text = 'false'; + end +elseif isinf(x) + if x > 0 + text = 'inf'; + else + text = '-inf'; + end +elseif isnan(x) + text = 'nan'; +elseif x == round(x) && abs(x) < 1e15 + text = sprintf('%d', x); +else + text = sprintf('%.17g', x); + % Use the shortest form that reads back to the same double. + for p = 1:17 + candidate = sprintf('%.*g', p, x); + if str2double(candidate) == x + text = candidate; + break + end + end + if isempty(regexp(text, '[.eE]', 'once')) + text = [text '.0']; + end +end +end + +function text = quote(s) +s = strrep(s, '\', '\\'); +s = strrep(s, '"', '\"'); +s = strrep(s, sprintf('\t'), '\t'); +s = strrep(s, newline, '\n'); +text = ['"' s '"']; +end diff --git a/src/utilities/createGECKOproject.m b/src/utilities/createGECKOproject.m index e1d127fb..5204c1be 100644 --- a/src/utilities/createGECKOproject.m +++ b/src/utilities/createGECKOproject.m @@ -8,23 +8,30 @@ function createGECKOproject(varargin) % -------------------- % name : char % a name for the folder structure used in GECKO. Also creates a basic -% adapter class, which must be manually adjusted. If not defined, a +% model adapter, which must be manually adjusted. If not defined, a % dialog box will appear. % path : char % a path where to create the folder. If not defined, a dialog box will % appear. +% adapterFormat : char +% 'toml' (default) writes model_adapter.toml, which GECKO and geckopy +% both read (see TomlAdapter); 'class' writes an adapter class +% Adapter.m instead. % % Examples % -------- % % optional arguments may be given positionally or as name-value pairs: % createGECKOproject(name, path); % createGECKOproject('name', name); +% createGECKOproject('name', name, 'adapterFormat', 'class'); p = parseGECKOargs(varargin, { ... 'name', []; ... - 'path', []}); + 'path', []; ... + 'adapterFormat', 'toml'}); name = p.name; path = p.path; +adapterFormat = validatestring(p.adapterFormat, {'toml','class'}); if isempty(name) prompt = {'Please provide a project name (e.g. ecYeastGEM)'}; @@ -51,9 +58,15 @@ function createGECKOproject(varargin) fclose(fid); end - % Read the template adapter class - fid = fopen(fullfile(findGECKOroot, 'src', ... - 'model_adapter', 'adapterTemplate.m')); + % Read the template adapter + if strcmp(adapterFormat, 'toml') + template = 'adapterTemplate.toml'; + filename = fullfile(path, name, 'model_adapter.toml'); + else + template = 'adapterTemplate.m'; + filename = fullfile(path, name, [name 'Adapter.m']); + end + fid = fopen(fullfile(findGECKOroot, 'src', 'model_adapter', template)); f = fread(fid, '*char')'; fclose(fid); @@ -62,8 +75,7 @@ function createGECKOproject(varargin) f = strrep(f, 'KEY_PATH', path); f = strrep(f, 'KEY_NAME', name); - % Save the class file - filename = fullfile(path, name, [name 'Adapter.m']); + % Save the adapter fid = fopen(filename, 'w'); fwrite(fid, f); fclose(fid); diff --git a/test/VerificationMatrix.tsv b/test/VerificationMatrix.tsv index fad92a5a..0ae3c5de 100644 --- a/test/VerificationMatrix.tsv +++ b/test/VerificationMatrix.tsv @@ -40,13 +40,20 @@ limit_proteins/flexibilizeProtConcs.m tc0013 limit_proteins/getConcControlCoeffs.m tc0013 limit_proteins/updateProtPool.m model_adapter/ModelAdapter.m Base class for adapters, not explicitly tested. -model_adapter/ModelAdapterManager.m tc0008 +model_adapter/ModelAdapterManager.m tc0008, tc0067 +model_adapter/readAdapterToml.m tc0061, tc0062 +model_adapter/writeAdapterToml.m tc0063 +model_adapter/loadAdapterParams.m tc0064, tc0065 +model_adapter/convertAdapterToToml.m tc0066 +model_adapter/TomlAdapter.m tc0066, tc0067 +model_adapter/adapterParameters.tsv tc0064, tc0065 +model_adapter/adapterTemplate.toml tc0068 utilities/* Not explicitly tested - The code has existed for a long time. There are a lot of functions in here, unclear which of them are still useful and functional. io/loadConventionalGEM.m tc0009 io/loadEcModel.m tc0009 io/saveEcModel.m tc0009 tutorials/tutorial_Human-GEM/HumanGEMAdapter.m tcm0002 - does not test all aspects of the adapter, but some. -tutorials/tutorial_yeast-GEM/YeastGEMAdapter.m tcm0001 - does not test all aspects of the adapter, but some. +tutorials/full_ecModel/model_adapter.toml tcm0001 - does not test all aspects of the adapter, but some. tutorials/tutorial_yeast-GEM/code/changeMedia_batch.m Not explicitly tested - The code has existed for a long time tutorials/tutorial_yeast-GEM/code/getModelParameters_obsolete.m Deprecated tutorials/tutorial_yeast-GEM/code/manualModifications.m Deprecated diff --git a/test/manual_tests/TestYeastGEMSimpleWorkflow.m b/test/manual_tests/TestYeastGEMSimpleWorkflow.m index f641ad2d..70c312cb 100644 --- a/test/manual_tests/TestYeastGEMSimpleWorkflow.m +++ b/test/manual_tests/TestYeastGEMSimpleWorkflow.m @@ -2,7 +2,7 @@ GECKORoot = findGECKOroot(); -yeastAdapter = ModelAdapterManager.getAdapter(fullfile(GECKORoot, 'tutorials', 'full_ecModel', 'YeastGEMAdapter.m')); +yeastAdapter = ModelAdapterManager.getAdapter(fullfile(GECKORoot, 'tutorials', 'full_ecModel')); yeastGEM = readYAMLmodel(fullfile(yeastAdapter.getParameters().path,'models','yeast-GEM.yml')); %Full model diff --git a/test/unit_tests/adapterTomlTests.m b/test/unit_tests/adapterTomlTests.m new file mode 100644 index 00000000..2bfb1cf3 --- /dev/null +++ b/test/unit_tests/adapterTomlTests.m @@ -0,0 +1,237 @@ +%run this test case with the command +%results = runtests('adapterTomlTests.m') + +function tests = adapterTomlTests + tests = functiontests(localfunctions); +end + +function testReadAdapterToml_tc0061(testCase) + % Every TOML feature that model adapter files use. + file = writeTemp({ ... + '# comment', ... + 'name = "a \"quoted\" # not a comment" # comment', ... + 'literal = ''C:\path''', ... + 'int = 559_292', ... + 'float = 2.5e-1', ... + 'yes = true', ... + 'list = ["a", "b",', ... + ' "c"] # over two lines', ... + 'nums = [0.4, 0.2, 0.1]', ... + 'empty = []', ... + '', ... + '[kegg]', ... + 'id = "sce"', ... + '', ... + '[evotune.source_groups.okp]', ... + 'sources = ["OKP"]', ... + 'match_okp = true', ... + '', ... + '[[stage]]', ... + 'action = "select"', ... + '[[stage]]', ... + 'action = "fill"', ... + 'sources = ["isozymes"]'}); + d = readAdapterToml(file); + verifyEqual(testCase, d.name, 'a "quoted" # not a comment') + verifyEqual(testCase, d.literal, 'C:\path') + verifyEqual(testCase, d.int, 559292) + verifyEqual(testCase, d.float, 0.25) + verifyEqual(testCase, d.yes, true) + verifyEqual(testCase, d.list, {'a','b','c'}) + verifyEqual(testCase, d.nums, [0.4 0.2 0.1]) + verifyEqual(testCase, d.empty, {}) + verifyEqual(testCase, d.kegg.id, 'sce') + verifyEqual(testCase, d.evotune.source_groups.okp.sources, {'OKP'}) + verifyTrue(testCase, d.evotune.source_groups.okp.match_okp) + verifyEqual(testCase, numel(d.stage), 2) + verifyEqual(testCase, d.stage(2).action, 'fill') + verifyEqual(testCase, d.stage(2).sources, {'isozymes'}) +end + +function testReadAdapterTomlErrors_tc0062(testCase) + % Unsupported TOML stops with the line number. + verifyError(testCase, @() readAdapterToml(writeTemp({'x = { a = 1 }'})), 'GECKO:readAdapterToml') + verifyError(testCase, @() readAdapterToml(writeTemp({'x = 2026-10-15'})), 'GECKO:readAdapterToml') + verifyError(testCase, @() readAdapterToml(writeTemp({'x = unquoted'})), 'GECKO:readAdapterToml') + verifyError(testCase, @() readAdapterToml(writeTemp({'x = 1', 'x = 2'})), 'GECKO:readAdapterToml') + try + readAdapterToml(writeTemp({'a = 1', 'b = """x"""'})); + catch err + verifySubstring(testCase, err.message, 'line 2') + end +end + +function testWriteAdapterToml_tc0063(testCase) + % What writeAdapterToml writes, readAdapterToml reads back unchanged. + s.text = 'with "quotes" and \ backslash'; + s.number = 0.1; + s.integer = 150; + s.flag = false; + s.list = {'dlkcat','brenda'}; + s.vector = [0.4 0.2]; + s.table.inner.value = 1e-8; + file = [tempname '.toml']; + writeAdapterToml(s, file, 'header'); + verifyEqual(testCase, readAdapterToml(file), s) +end + +function testLoadAdapterParamsDefaults_tc0064(testCase) + % A minimal file gets geckopy's defaults, under GECKO's field names. + folder = tempname; mkdir(folder); + writeTemp({'conv_gem = "models/x.xml"', 'org_name = "testus testus"', 'p_tot = 0.4'}, ... + fullfile(folder, 'model_adapter.toml')); + params = loadAdapterParams(folder); + verifyEqual(testCase, params.path, folder) + verifyEqual(testCase, params.convGEM, fullfile(folder, 'models/x.xml')) + verifyEqual(testCase, params.Ptot, 0.4) + verifyEqual(testCase, params.sigma, 0.5) + verifyEqual(testCase, params.kegg.ID, 'sce') + verifyEqual(testCase, params.kcatAggregation, 'max') + verifyEqual(testCase, params.evotune.kcatSources, {'dlkcat','brenda','custom'}) + verifyEqual(testCase, params.evotune.sigma0logSource, [0.4; 0.2; 0.1]) + verifyEqual(testCase, params.evotune.maxGenerations, 150) + verifyEqual(testCase, params.okp.handleLongSequences, 'truncate') + verifyEmpty(testCase, params.complex.taxonomicID) + % A misspelled key is an error, as in geckopy. + writeTemp({'conv_gem = "x.xml"', 'org_name = "x"', 'sigmaa = 1'}, fullfile(folder, 'model_adapter.toml')); + verifyError(testCase, @() loadAdapterParams(folder), 'GECKO:loadAdapterParams') +end + +function testLoadAdapterParamsSourceGroups_tc0065(testCase) + % Trust tiers with several sources and the OKP method. + folder = tempname; mkdir(folder); + writeTemp({'conv_gem = "x.xml"', 'org_name = "x"', ... + '[okp]', 'method = "EITLEM"', ... + '[evotune.source_groups.okp]', 'sources = ["OpenKineticsPredictor"]', 'match_okp = true', ... + '[evotune.source_groups.database]', 'sources = ["brenda", "okp_brenda"]', ... + '[evotune.sigma0_log_source]', 'okp = 0.25', 'database = 0.2', ... + '[matlab]', 'make_model_anaerobic = "anaerobicModel_GECKO"', ... + '[matlab.params]', 'c_UptakeExp = 0.64'}, fullfile(folder, 'model_adapter.toml')); + params = loadAdapterParams(fullfile(folder, 'model_adapter.toml')); + verifyEqual(testCase, params.evotune.kcatSources, {'OpenKineticsPredictor','EITLEM','brenda','okp_brenda'}) + verifyEqual(testCase, params.evotune.sigma0logSource, [0.25; 0.25; 0.2; 0.2]) + verifyEqual(testCase, params.matlab.makeModelAnaerobic, 'anaerobicModel_GECKO') + verifyEqual(testCase, params.c_UptakeExp, 0.64) +end + +function testConvertAdapterToToml_tc0066(testCase) + % Converting an adapter class and loading the TOML gives the same + % parameters back. + geckoPath = findGECKOroot; + legacyFolder = tempname; mkdir(legacyFolder); + writeTemp(legacyAdapterClass(), fullfile(legacyFolder, 'LegacyFixtureAdapter.m')); + adapters = {fullfile(geckoPath,'test','unit_tests','ecTestGEM','TestGEMAdapter.m'), ... + fullfile(legacyFolder,'LegacyFixtureAdapter.m')}; + for k = 1:numel(adapters) + original = ModelAdapterManager.getAdapter(adapters{k}); + folder = tempname; mkdir(folder); + file = fullfile(folder, 'model_adapter.toml'); + report = convertAdapterToToml(original, file); + verifyTrue(testCase, any(strcmp(report.methods, 'getSpontaneousReactions'))) + converted = TomlAdapter(file); + expected = original.params; + actual = converted.params; + % The project folder is where the TOML file is; the model path is + % relative to it. + [~, convName] = fileparts(expected.convGEM); + [~, convName2] = fileparts(actual.convGEM); + verifyEqual(testCase, convName2, convName) + expected = rmfield(expected, {'path','convGEM'}); + verifyEqual(testCase, keepFields(actual, expected), expected) + end + verifyEqual(testCase, report.unmapped, {'c_UptakeExp'}) +end + +function lines = legacyAdapterClass() +% An adapter class as written for GECKO 3 and 4: evotune trust tiers and a +% parameter that geckopy does not know. +lines = { ... + 'classdef LegacyFixtureAdapter < ModelAdapter', ... + ' methods', ... + ' function obj = LegacyFixtureAdapter()', ... + ' obj.params.path = fileparts(mfilename(''fullpath''));', ... + ' obj.params.convGEM = fullfile(obj.params.path,''models'',''human-GEM.yml'');', ... + ' obj.params.sigma = 0.1;', ... + ' obj.params.Ptot = 0.505717472;', ... + ' obj.params.f = 0.412;', ... + ' obj.params.gR_exp = 0.020663429;', ... + ' obj.params.c_UptakeExp = 0.641339301;', ... + ' obj.params.org_name = ''homo sapiens'';', ... + ' obj.params.kcatAggregation = ''median'';', ... + ' obj.params.complex.taxonomicID = 9606;', ... + ' obj.params.kegg.ID = ''hsa'';', ... + ' obj.params.kegg.geneID = ''Ensembl'';', ... + ' obj.params.uniprot.type = ''proteome'';', ... + ' obj.params.uniprot.ID = ''UP000005640'';', ... + ' obj.params.uniprot.geneIDfield = ''gene_primary'';', ... + ' obj.params.uniprot.reviewed = true;', ... + ' obj.params.c_source = ''MAR09034'';', ... + ' obj.params.bioRxn = ''MAR13082'';', ... + ' obj.params.enzyme_comp = ''Cytosol'';', ... + ' obj.params.evotune.sigma0logDefault = 0.5;', ... + ' obj.params.evotune.kcatSources = {''dlkcat'',''brenda'',''custom''};', ... + ' obj.params.evotune.sigma0logSource = [0.4; 0.2; 0.1];', ... + ' obj.params.evotune.maxGrowthWeight = 2.0;', ... + ' obj.params.evotune.priorPenaltyWeight = 0.03;', ... + ' obj.params.evotune.tieIsozymes = true;', ... + ' obj.params.evotune.rmseThreshold = 0.2;', ... + ' obj.params.evotune.maxGenerations = 212;', ... + ' end', ... + ' function [spont,spontRxnNames] = getSpontaneousReactions(obj,model)', ... + ' spont = contains(model.rxnNames,''spontaneous'');', ... + ' spontRxnNames = model.rxnNames(spont);', ... + ' end', ... + ' end', ... + 'end'}; +end + +function testTomlAdapterSpontaneous_tc0067(testCase) + folder = tempname; mkdir(fullfile(folder,'data')); + writeTemp({'rxns spontaneous', 'R1 0', 'R2 1', 'R3 1'}, fullfile(folder,'data','spont.tsv')); + writeTemp({'conv_gem = "x.xml"', 'org_name = "x"', ... + '[spontaneous]', 'name_contains = "spontaneous"', 'file = "data/spont.tsv"'}, ... + fullfile(folder, 'model_adapter.toml')); + adapter = ModelAdapterManager.getAdapter(folder); + verifyClass(testCase, adapter, 'TomlAdapter') + model.rxns = {'R1';'R2';'R4';'R5'}; + model.rxnNames = {'a';'b';'spontaneous c';'d'}; + [spont, names] = adapter.getSpontaneousReactions(model); + verifyEqual(testCase, spont, [false; true; true; false]) + verifyEqual(testCase, names, {'R2';'R4'}) +end + +function testCreateGECKOprojectToml_tc0068(testCase) + % A new project gets model_adapter.toml, which loads as it is. + here = pwd; + restore = onCleanup(@() cd(here)); + root = tempname; mkdir(root); + createGECKOproject('name', 'ecTestProject', 'path', root); + folder = fullfile(root, 'ecTestProject'); + verifyTrue(testCase, isfile(fullfile(folder, 'model_adapter.toml'))) + adapter = ModelAdapterManager.getAdapter(folder); + verifyEqual(testCase, adapter.params.org_name, 'genus species') + verifyEqual(testCase, adapter.params.evotune.kcatSources, {'dlkcat','brenda','custom'}) + verifyEqual(testCase, adapter.params.spontaneous.name_contains, 'spontaneous') +end + +function s = keepFields(s, template) +% s reduced to the fields (recursively) that template has. +fields = fieldnames(s); +for k = 1:numel(fields) + if ~isfield(template, fields{k}) + s = rmfield(s, fields{k}); + elseif isstruct(s.(fields{k})) && isstruct(template.(fields{k})) + s.(fields{k}) = keepFields(s.(fields{k}), template.(fields{k})); + end +end +s = orderfields(s, template); +end + +function file = writeTemp(lines, file) +if nargin < 2 + file = [tempname '.toml']; +end +fid = fopen(file, 'w'); +fprintf(fid, '%s\n', lines{:}); +fclose(fid); +end diff --git a/tutorials/full_ecModel/YeastGEMAdapter.m b/tutorials/full_ecModel/YeastGEMAdapter.m deleted file mode 100644 index b9ad7414..00000000 --- a/tutorials/full_ecModel/YeastGEMAdapter.m +++ /dev/null @@ -1,89 +0,0 @@ -classdef YeastGEMAdapter < ModelAdapter - methods - function obj = YeastGEMAdapter() - obj.params.path = fullfile(findGECKOroot,'tutorials','full_ecModel'); - addpath(fullfile(obj.params.path,'code')); - - obj.params.convGEM = fullfile(obj.params.path,'models','yeast-GEM.yml'); - - obj.params.sigma = 0.5; - - obj.params.Ptot = 0.5; - - obj.params.f = 0.5; - - obj.params.gR_exp = 0.41; - - obj.params.org_name = 'saccharomyces cerevisiae'; - - obj.params.complex.taxonomicID = 559292; - - obj.params.kegg.ID = 'sce'; - - obj.params.kegg.geneID = 'kegg'; - - obj.params.uniprot.type = 'proteome'; - - obj.params.uniprot.ID = 'UP000002311'; - - obj.params.uniprot.geneIDfield = 'gene_oln'; - - obj.params.uniprot.reviewed = true; - - obj.params.c_source = 'r_1714'; - - obj.params.bioRxn = 'r_4041'; - - % Name of the compartment where the protein pseudometabolites - % should be located (all be located in the same compartment, - % this does not interfere with them catalyzing reactions in - % different compartments). Typically, cytoplasm is chosen. - obj.params.enzyme_comp = 'cytoplasm'; - - %% Hyperparameters for evotune (CMA-ES) kcat fitting - % Default initial uncertainty (standard deviation in log-space) for kcat values - obj.params.evotune.sigma0logDefault = 0.5; - % Data sources for kcat values, ordered from least to most trusted - obj.params.evotune.kcatSources = {'dlkcat','brenda','custom'}; - % Initial uncertainty for each source (lower = more trusted data) - obj.params.evotune.sigma0logSource = [0.4; 0.2; 0.1]; - - % Weight on the max-growth RMSE against the flux RMSE: - % (rmseFlux + w*rmseMaxGrowth) / (w+1). At 2 the max-growth - % conditions count double against the flux conditions. - obj.params.evotune.maxGrowthWeight = 2.0; - % Weight on a prior term in the search objective, - % rmse + w*mean((log(k/k0)/sigma0log)^2). Keeps large corrections - % reproducible across seeds; 0 scores on RMSE alone. - obj.params.evotune.priorPenaltyWeight = 0.03; - % Give isozyme copies of one reaction (sharing a prior value and - % a source) a single shared kcat, so the search cannot invent a - % distinction the kcat assignment never made. - obj.params.evotune.tieIsozymes = true; - - % Stop optimization when RMSE falls below this threshold - obj.params.evotune.rmseThreshold = 0.2; - % Maximum number of CMA-ES generations before termination - obj.params.evotune.maxGenerations = 212; - end - - function ecModel = makeModelAnaerobic(obj,ecModel) - % Constrains ecModel to anaerobic growth conditions (blocked - % oxygen uptake, enabled sterol/fatty acid exchanges, and other - % curations), as implemented in anaerobicModel_GECKO. - ecModel = anaerobicModel_GECKO(ecModel); - end - function ecModel = changeProteinBiomass(obj,ecModel,Ptot) - % Currently a no-op: returns ecModel unchanged. Intended to - % rescale the biomass reaction's protein content to Ptot (as - % scaleBioMass_GECKO would), but this step is disabled for now. - ecModel = ecModel; - % Taken from yeast-GEM 9.0.2 - %ecModel = scaleBioMass_GECKO(ecModel,'protein',Ptot,'carbohydrate',false); - end - function [spont,spontRxnNames] = getSpontaneousReactions(obj,model) - spont = contains(model.rxnNames,'spontaneous'); - spontRxnNames = model.rxnNames(spont); - end - end -end diff --git a/tutorials/full_ecModel/code/plotlightVSfull.m b/tutorials/full_ecModel/code/plotlightVSfull.m index 216fed1e..f497351d 100644 --- a/tutorials/full_ecModel/code/plotlightVSfull.m +++ b/tutorials/full_ecModel/code/plotlightVSfull.m @@ -18,7 +18,7 @@ % mapped onto the conventional (non-ec) GEM reactions % -adapterLocation = fullfile(findGECKOroot,'tutorials','full_ecModel','YeastGEMAdapter.m'); +adapterLocation = fullfile(findGECKOroot,'tutorials','full_ecModel'); ModelAdapterManager.setDefault(adapterLocation); model = loadConventionalGEM(); model = setParam(model,'lb','r_1714',-1000); diff --git a/tutorials/full_ecModel/model_adapter.toml b/tutorials/full_ecModel/model_adapter.toml new file mode 100644 index 00000000..ea7d466a --- /dev/null +++ b/tutorials/full_ecModel/model_adapter.toml @@ -0,0 +1,85 @@ +# Model adapter for the full_ecModel tutorial (yeast-GEM). +# Read by GECKO (ModelAdapterManager.setDefault on this folder) and by +# geckopy (ModelAdapter.from_folder). Paths are relative to this folder. +# Every key, its MATLAB name and its default are listed in GECKO's +# src/model_adapter/adapterParameters.tsv. + +# Conventional GEM that the ecModel is built from. +conv_gem = "models/yeast-GEM.yml" + +org_name = "saccharomyces cerevisiae" + +# Average enzyme saturation factor. +sigma = 0.5 +# Total protein content [g protein/gDW]. +p_tot = 0.5 +# Fraction of enzymes in the model [g enzyme/g protein]. +f = 0.5 +# Growth rate the model should reach when not limited by nutrient uptake [1/h]. +gr_exp = 0.41 + +c_source = "r_1714" +bio_rxn = "r_4041" +# Compartment of the protein pseudometabolites. They can all sit in one +# compartment, whichever compartments the enzymes act in; typically the +# cytoplasm. +enzyme_comp = "cytoplasm" + +[kegg] +id = "sce" +gene_id = "kegg" + +[uniprot] +type = "proteome" +id = "UP000002311" +gene_id_field = "gene_oln" +reviewed = true + +[complex] +taxonomic_id = 559292 + +# Reactions with "spontaneous" in their name have no enzyme. +[spontaneous] +name_contains = "spontaneous" + +# Hyperparameters for evotune (CMA-ES) kcat tuning. +[evotune] +# Prior standard deviation (log) for kcats whose source is in no group. +sigma0_log_default = 0.5 +# Weight on the max-growth RMSE against the flux RMSE: +# (rmse_flux + w * rmse_max_growth) / (w + 1). At 2 the max-growth +# conditions count double against the flux conditions. +max_growth_weight = 2.0 +# Weight on a prior term in the search objective, +# rmse + w * mean((log(k/k0) / sigma0_log)^2). Keeps large corrections +# reproducible across seeds; 0 scores on RMSE alone. +prior_penalty_weight = 0.03 +# Give isozyme copies of one reaction (sharing a prior value and a source) +# a single shared kcat, so the search cannot invent a distinction the kcat +# assignment never made. +tie_isozymes = true +# Stop when the RMSE falls below this. +rmse_threshold = 0.2 +# Maximum number of CMA-ES generations. +max_generations = 212 + +# Kcat sources, from least to most trusted, each with its prior width +# (lower is more trusted). +[evotune.source_groups.dlkcat] +sources = ["dlkcat"] + +[evotune.source_groups.brenda] +sources = ["brenda"] + +[evotune.source_groups.custom] +sources = ["custom"] + +[evotune.sigma0_log_source] +dlkcat = 0.4 +brenda = 0.2 +custom = 0.1 + +# Read by GECKO only. During kcat tuning, anaerobic conditions are set by +# anaerobicModel_GECKO in the code folder, which GECKO adds to the path. +[matlab] +make_model_anaerobic = "anaerobicModel_GECKO" diff --git a/tutorials/full_ecModel/protocol.m b/tutorials/full_ecModel/protocol.m index a6379254..84395d3a 100644 --- a/tutorials/full_ecModel/protocol.m +++ b/tutorials/full_ecModel/protocol.m @@ -54,17 +54,19 @@ % STEP 3-7 Modify the model adapter % In the Nature Protocols paper is explained how to decide on the organism- % and model-specific parameters in the model adapter, which for this -% tutorial is located at tutorials/full_ecModel/YeastGEMAdapter.m. +% tutorial is located at tutorials/full_ecModel/model_adapter.toml. The +% same file configures geckopy, the Python version of GECKO. %% STAGE 1: Expansion from a starting metabolic model to an ecModel structure % STEP 8 Set modelAdapter -adapterLocation = fullfile(findGECKOroot,'tutorials','full_ecModel','YeastGEMAdapter.m'); +adapterLocation = fullfile(findGECKOroot,'tutorials','full_ecModel'); ModelAdapter = ModelAdapterManager.setDefault(adapterLocation); -% With the above line, we set the YeastGEMAdapter as the default adapter -% from here onwards, which means that any GECKO function that requires a -% model adapter as input will use YeastGEMAdapter. This is until you have -% set a new default adapter, or closed MATLAB. +% With the above line, we set the model adapter of this project folder, +% read from its model_adapter.toml, as the default adapter from here +% onwards, which means that any GECKO function that requires a model +% adapter as input will use it. This is until you have set a new default +% adapter, or closed MATLAB. % We can also load the ModelAdapter into the MATLAB Workspace, so that you % can quickly see its content, for instance the parameters. @@ -457,7 +459,7 @@ % STEP 67-68 Simulate Crabtree effect with protein pool % We will below run a custom plotCrabtree function that is kept in the code -% subfolder, which YeastGEMAdapter adds to the MATLAB path. +% subfolder, which the model adapter adds to the MATLAB path. [fluxes, gRate] = plotCrabtree(ecModel); % fluxes has all the predicted fluxes, while gRate is a vector with the % corresponding growth rates that were simulated, as visualized on the diff --git a/tutorials/light_ecModel/HumanGEMAdapter.m b/tutorials/light_ecModel/HumanGEMAdapter.m index 13756b80..43d94503 100644 --- a/tutorials/light_ecModel/HumanGEMAdapter.m +++ b/tutorials/light_ecModel/HumanGEMAdapter.m @@ -1,76 +1,20 @@ -classdef HumanGEMAdapter < ModelAdapter +% HumanGEMAdapter Model adapter for the light_ecModel tutorial. +% +% The parameters are in model_adapter.toml in this folder, which geckopy +% reads as well. This class only adds getHumanGEMRootPath, so it shows how +% a project keeps custom code next to a TOML adapter: subclass TomlAdapter +% and pass the project folder to its constructor. +% +% See also +% -------- +% TomlAdapter, loadAdapterParams +classdef HumanGEMAdapter < TomlAdapter methods function obj = HumanGEMAdapter() - geckoPath = findGECKOroot; - obj.params.path = fullfile(geckoPath,'tutorials','light_ecModel'); - - % The model distributed with the light_ecModel tutorial is Human-GEM - % is version 1.15.0, available from - % https://github.com/SysBioChalmers/Human-GEM/releases/tag/v1.15.0 - % In addition, the following lines were run to reduce its size - % before storing it in this GECKO tutorial: - % ihuman = simplifyModel(ihuman,false,false,true,true); - % [ihuman.grRules,skipped] = simplifyGrRules(ihuman.grRules,true); - % ihuman = deleteUnusedGenes(ihuman); - % ihuman = rmfield(ihuman,{'unconstrained','rxnReferences','rxnFrom','metFrom','rxnConfidenceScores'}); - % ihuman.name = ihuman.id; - % writeYAMLmodel(ihuman,'human-GEM.yml') - obj.params.convGEM = fullfile(obj.params.path,'models','human-GEM.yml'); - - obj.params.sigma = 0.1; %This was changed to a low number to give a reasonable growth rate - this should be investigated more - - obj.params.Ptot = 0.505717472; %Average across NCI60 cell lines - - obj.params.f = 0.412; %estimated as TPM of model genes / all genes in GTEx - - obj.params.gR_exp = 0.020663429; %[g/gDw h]/Average across NCI60 cell lines - - obj.params.c_UptakeExp = 0.641339301; %[mmol/gDw h]/Average across NCI60 cell lines - - obj.params.org_name = 'homo sapiens'; - - obj.params.complex.taxonomicID = 9606; - - obj.params.kegg.ID = 'hsa'; - - obj.params.kegg.geneID = 'Ensembl'; - - obj.params.uniprot.type = 'proteome'; - - obj.params.uniprot.ID = 'UP000005640'; - - obj.params.uniprot.geneIDfield = 'gene_primary'; - - obj.params.uniprot.reviewed = true; - - obj.params.c_source = 'MAR09034'; - - obj.params.bioRxn = 'MAR13082'; - - % Name of the compartment where the protein pseudometabolites - % should be located (all be located in the same compartment, - % this does not interfere with them catalyzing reactions in - % different compartments). Typically, cytoplasm is chosen. - obj.params.enzyme_comp = 'Cytosol'; - - end - - function [spont,spontRxnNames] = getSpontaneousReactions(obj,model) - fID=fopen(fullfile(obj.params.path,'data','spontaneousReactions.tsv'),'r'); - rxns_tsv = textscan(fID,'%q %q','Delimiter','\t','HeaderLines',1); - fclose(fID); - % The above file is derived from the reactions.tsv that is available - % in the Human-GEM repository, matching version 1.15.0 of the model. - % https://github.com/SysBioChalmers/Human-GEM/blob/194ebe5431c83e25f78df61caacad2fa485b5cb4/model/reactions.tsv - % Only the columns with reaction identifiers and assignment of - % spontaneous reactions were kept in spontaneousReactions.tsv. - spont = logical(str2double(cell2mat(rxns_tsv{2}))); - spontRxnNames = rxns_tsv{1}; - [spont,~] = ismember(model.rxns,spontRxnNames(spont)); - spontRxnNames = model.rxns(spont); + obj@TomlAdapter(fullfile(findGECKOroot,'tutorials','light_ecModel')); end end - + methods(Static) function path = getHumanGEMRootPath() path = fileparts(which('Human-GEM.mat'));% This file should be on the path diff --git a/tutorials/light_ecModel/model_adapter.toml b/tutorials/light_ecModel/model_adapter.toml new file mode 100644 index 00000000..9f1a45a0 --- /dev/null +++ b/tutorials/light_ecModel/model_adapter.toml @@ -0,0 +1,58 @@ +# Model adapter for the light_ecModel tutorial (Human-GEM). +# Read by GECKO through HumanGEMAdapter, which adds custom code, and by +# geckopy (ModelAdapter.from_folder). Paths are relative to this folder. +# Every key, its MATLAB name and its default are listed in GECKO's +# src/model_adapter/adapterParameters.tsv. + +# Human-GEM 1.15.0 +# (https://github.com/SysBioChalmers/Human-GEM/releases/tag/v1.15.0), +# reduced in size before it was stored in this tutorial: +# ihuman = simplifyModel(ihuman,false,false,true,true); +# [ihuman.grRules,skipped] = simplifyGrRules(ihuman.grRules,true); +# ihuman = deleteUnusedGenes(ihuman); +# ihuman = rmfield(ihuman,{'unconstrained','rxnReferences','rxnFrom','metFrom','rxnConfidenceScores'}); +# ihuman.name = ihuman.id; +# writeYAMLmodel(ihuman,'human-GEM.yml') +conv_gem = "models/human-GEM.yml" + +org_name = "homo sapiens" + +# Set low to give a reasonable growth rate; this should be investigated +# further. +sigma = 0.1 +# Average across NCI60 cell lines [g protein/gDW]. +p_tot = 0.505717472 +# Estimated as TPM of model genes / all genes in GTEx. +f = 0.412 +# Average across NCI60 cell lines [1/h]. +gr_exp = 0.020663429 + +c_source = "MAR09034" +bio_rxn = "MAR13082" +# Compartment of the protein pseudometabolites. They can all sit in one +# compartment, whichever compartments the enzymes act in; typically the +# cytoplasm. +enzyme_comp = "Cytosol" + +[kegg] +id = "hsa" +gene_id = "Ensembl" + +[uniprot] +type = "proteome" +id = "UP000005640" +gene_id_field = "gene_primary" +reviewed = true + +[complex] +taxonomic_id = 9606 + +# Derived from reactions.tsv of Human-GEM 1.15.0 +# (https://github.com/SysBioChalmers/Human-GEM/blob/194ebe5431c83e25f78df61caacad2fa485b5cb4/model/reactions.tsv), +# keeping the reaction identifiers and the spontaneous column. +[spontaneous] +file = "data/spontaneousReactions.tsv" + +[matlab.params] +# Glucose uptake, average across NCI60 cell lines [mmol/gDW/h]. +c_UptakeExp = 0.641339301