Read the model adapter from model_adapter.toml, shared with geckopy

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];

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

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

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"