diff --git a/core/applyCondition.m b/core/applyCondition.m new file mode 100644 index 00000000..8494485f --- /dev/null +++ b/core/applyCondition.m @@ -0,0 +1,130 @@ +function model = applyCondition(model, condition) +% applyCondition +% Apply a deterministic "condition" to a model: a prelude that resets +% exchange bounds, optional metabolite removals + automatic charge +% rebalancing of a pseudoreaction, optional biomass-stoichiometry +% delta, and a per-reaction bounds diff. The schema is intentionally +% narrow so a condition can be reviewed as data. +% +% Yeast-GEM was the first consumer; the same schema works for any +% GEM that keeps its condition presets as data rather than as code. +% Project-specific extensions (e.g. yeast-GEM's amino_acid_ratio +% step that rewrites a protein pseudoreaction's stoichiometry from a +% side-car TSV) are handled by the *caller* before / after this +% function — kept upstream-narrow on purpose. +% +% Inputs: +% model RAVEN model struct. +% condition Either a path to a YAML condition file or a struct +% already produced by parseYAML. The expected schema +% (all keys optional): +% +% prelude: +% reset_exchanges: out % truthy -> reset all +% +% cofactor_pseudoreaction: +% rxn_id: r_4598 +% remove_mets: +% - { met: s_3714 } +% charge_balance_met: s_0794 +% +% biomass_stoichiometry_delta: +% rxn_id: r_4041 +% add: +% - { met: s_0689, coef: 0.08 } +% - { met: s_0687, coef: -0.08 } +% - { met: s_0794, coef: -0.16 } +% +% bounds: +% - { rxn: r_1654, lb: -1000 } +% - { rxn: r_1992, lb: 0 } +% - { rxn: r_1663, lb: 0, ub: 0 } +% +% expected_uptake_count: 15 +% +% Output: +% model Modified model. +% +% Usage: model = applyCondition(model, 'data/conditions/anaerobic.yml') +% model = applyCondition(model, parseYAML('data/conditions/anaerobic.yml')) + +if ischar(condition) || isstring(condition) + cond = parseYAML(char(condition)); +elseif isstruct(condition) + cond = condition; +else + error('applyCondition:invalidCondition', ... + 'condition must be a YAML file path or a struct.'); +end + +% --- Step 1: prelude --------------------------------------------------- +if isfield(cond, 'prelude') && isfield(cond.prelude, 'reset_exchanges') + [~, exchangeRxns] = getExchangeRxns(model, cond.prelude.reset_exchanges); + model.lb(exchangeRxns) = 0; + model.ub(exchangeRxns) = 1000; +end + +% --- Step 2: cofactor pseudoreaction edits ---------------------------- +if isfield(cond, 'cofactor_pseudoreaction') + cp = cond.cofactor_pseudoreaction; + cofacIdx = getIndexes(model, cp.rxn_id, 'rxns'); + if isfield(cp, 'remove_mets') + for i = 1:numel(cp.remove_mets) + metIdx = getIndexes(model, cp.remove_mets{i}.met, 'mets'); + model.S(metIdx, cofacIdx) = 0; + end + end + if isfield(cp, 'charge_balance_met') + balanceIdx = find(strcmp(model.mets, cp.charge_balance_met)); + model.S(balanceIdx, cofacIdx) = 0; + model.S(balanceIdx, cofacIdx) = ... + -sum(model.S(:, cofacIdx) .* model.metCharges, 'omitnan'); + end +end + +% --- Step 3: biomass stoichiometry delta ------------------------------ +if isfield(cond, 'biomass_stoichiometry_delta') + delta = cond.biomass_stoichiometry_delta; + bioIdx = getIndexes(model, delta.rxn_id, 'rxns'); + if isfield(delta, 'add') + for i = 1:numel(delta.add) + entry = delta.add{i}; + metIdx = getIndexes(model, entry.met, 'mets'); + model.S(metIdx, bioIdx) = full(model.S(metIdx, bioIdx)) + entry.coef; + end + end +end + +% --- Step 4: bounds --------------------------------------------------- +nUptake = 0; +if isfield(cond, 'bounds') + for i = 1:numel(cond.bounds) + b = cond.bounds{i}; + rxnIdx = find(strcmp(model.rxns, b.rxn)); + if isempty(rxnIdx) + warning('applyCondition:missingRxn', ... + 'Reaction %s not found in model; skipping.', b.rxn); + continue; + end + if isfield(b, 'lb') + model.lb(rxnIdx) = b.lb; + if b.lb == -1000 + nUptake = nUptake + 1; + end + end + if isfield(b, 'ub') + model.ub(rxnIdx) = b.ub; + end + end +end + +% --- Step 5: uptake sanity check -------------------------------------- +if isfield(cond, 'expected_uptake_count') + if nUptake ~= cond.expected_uptake_count + warning('applyCondition:uptakeCountMismatch', ... + 'Expected %d uptake reactions, applied %d. Some may be missing from the model.', ... + cond.expected_uptake_count, nUptake); + end +end + +end diff --git a/core/assignSBOterms.m b/core/assignSBOterms.m new file mode 100644 index 00000000..0c0e1755 --- /dev/null +++ b/core/assignSBOterms.m @@ -0,0 +1,159 @@ +function model = assignSBOterms(model, opts) +% assignSBOterms +% Assign SBO terms to metabolites and reactions following a generic +% rule set. Mirrors raven_python.annotation.add_sbo_terms; +% organism-agnostic, parameterised entirely by `opts`. The +% yeast-GEM port of this function is the legacy addSBOterms.m, +% which becomes a thin shim here. +% +% Rules +% ----- +% Metabolites: +% SBO:0000649 (Biomass) when met.name is in opts.biomassMetNames, +% or ends with any of opts.biomassMetSuffixes. Otherwise +% SBO:0000247 (Simple chemical). +% +% Reactions (default → override → pseudoreaction override): +% SBO:0000176 (Metabolic reaction) default. +% Single-reactant reactions become: +% SBO:0000627 (exchange) if the lone metabolite is +% extracellular (compartment 'e' or compartment name +% containing 'extracellular'), +% SBO:0000632 (sink) if coef < 0, +% SBO:0000628 (demand) otherwise. +% Transport reactions (detected by opts.transportDetector or +% the default heuristic: same metName in ≥ 2 compartments +% in a single reaction) → SBO:0000655. +% Reactions whose name matches opts.biomassRxnName → SBO:0000629. +% Reactions whose name matches opts.ngamRxnName → SBO:0000630. +% Reactions whose name contains any of +% opts.pseudoreactionSubstrings → SBO:0000395. +% +% "fill" semantic — SBO is written via editMiriam(..., 'fill') so +% pre-existing SBO annotations are preserved. +% +% Inputs: +% model RAVEN model struct. +% opts (opt) struct with any of the following fields. Missing +% fields take the defaults shown: +% biomassMetNames {'biomass','DNA','RNA','protein', +% 'carbohydrate','lipid','cofactor','ion'} +% biomassMetSuffixes {' backbone',' chain'} +% biomassRxnName 'biomass pseudoreaction' +% ngamRxnName 'non-growth associated maintenance reaction' +% pseudoreactionSubstrings {'pseudoreaction','SLIME rxn'} +% onlyLastReactionForPseudo false. yeast-GEM bug-compat +% flag — replicates the +% legacy `for i=numel(...)` +% typo (pseudoreaction +% SBOs applied only to the +% last reaction). Off by +% default; turn ON for +% byte-equivalent yeast-GEM +% output. +% +% Output: +% model Modified model. +% +% Usage: model = assignSBOterms(model) +% model = assignSBOterms(model, struct('onlyLastReactionForPseudo', true)) + +if nargin < 2 || isempty(opts) + opts = struct(); +end +opts = applyDefaults(opts); + +% Metabolite SBO ------------------------------------------------------ +metsSBO = cell(size(model.mets)); +for i = 1:length(model.mets) + metName = model.metNames{i}; + if any(strcmp(opts.biomassMetNames, metName)) || endsWithAny(metName, opts.biomassMetSuffixes) + metsSBO{i} = 'SBO:0000649'; + else + metsSBO{i} = 'SBO:0000247'; + end +end + +% Reaction SBO -------------------------------------------------------- +rxnSBO = cell(size(model.rxns)); +rxnSBO(:) = {'SBO:0000176'}; + +% Single-reactant reactions +reactantNumber = sum(model.S ~= 0, 1); +singleRxns = find(reactantNumber == 1); +for k = 1:numel(singleRxns) + idx = singleRxns(k); + metRow = find(model.S(:, idx)); + compName = model.compNames{model.metComps(metRow)}; + compShort = model.comps{model.metComps(metRow)}; + if strcmp(compShort, 'e') || strcmp(compName, 'extracellular') + rxnSBO{idx} = 'SBO:0000627'; + elseif sum(model.S(:, idx)) < 0 + rxnSBO{idx} = 'SBO:0000632'; + else + rxnSBO{idx} = 'SBO:0000628'; + end +end + +% Transport reactions +if isfield(opts, 'transportRxnIdxs') && ~isempty(opts.transportRxnIdxs) + transportIdxs = opts.transportRxnIdxs; +else + transportIdxs = getTransportRxns(model); +end +rxnSBO(transportIdxs) = {'SBO:0000655'}; + +% Pseudoreaction overrides +if opts.onlyLastReactionForPseudo + pseudoTargets = numel(model.rxns); +else + pseudoTargets = 1:numel(model.rxns); +end +for ii = pseudoTargets + name = model.rxnNames{ii}; + if strcmp(name, opts.biomassRxnName) + rxnSBO{ii} = 'SBO:0000629'; + elseif strcmp(name, opts.ngamRxnName) + rxnSBO{ii} = 'SBO:0000630'; + else + for k = 1:numel(opts.pseudoreactionSubstrings) + if contains(name, opts.pseudoreactionSubstrings{k}) + rxnSBO{ii} = 'SBO:0000395'; + break; + end + end + end +end + +model = editMiriam(model, 'met', 'all', 'sbo', metsSBO, 'fill'); +model = editMiriam(model, 'rxn', 'all', 'sbo', rxnSBO, 'fill'); +end + + +function opts = applyDefaults(opts) +defaults = struct( ... + 'biomassMetNames', {{'biomass','DNA','RNA','protein','carbohydrate','lipid','cofactor','ion'}}, ... + 'biomassMetSuffixes', {{' backbone',' chain'}}, ... + 'biomassRxnName', 'biomass pseudoreaction', ... + 'ngamRxnName', 'non-growth associated maintenance reaction', ... + 'pseudoreactionSubstrings', {{'pseudoreaction','SLIME rxn'}}, ... + 'onlyLastReactionForPseudo', false); +fields = fieldnames(defaults); +for k = 1:numel(fields) + f = fields{k}; + if ~isfield(opts, f) || isempty(opts.(f)) + opts.(f) = defaults.(f); + end +end +end + + +function tf = endsWithAny(s, suffixes) +tf = false; +for i = 1:numel(suffixes) + if endsWith(s, suffixes{i}) + tf = true; + return; + end +end +end diff --git a/core/curateModelFromTables.m b/core/curateModelFromTables.m new file mode 100644 index 00000000..742b50bd --- /dev/null +++ b/core/curateModelFromTables.m @@ -0,0 +1,345 @@ +function newModel=curateModelFromTables(model,metsInfo,genesInfo,rxnsCoeffs,rxnsInfo,metPrefix,rxnPrefix) +% curateModelFromTables +% Curate existing and/or add new metabolites, reactions and genes +% from tabular data files. Originally extracted from yeast-GEM's +% curateMetsRxnsGenes; generalised here so any GEM project can drive +% batch curation from the same set of *.tsv files. +% +% If the *.tsv files contain metabolites, reactions and/or genes that are +% already present in the model, then information in the model will be +% overwritten. Note that this includes empty annotations in the *.tsv +% files! Metabolites are matched by metaboliteName[comp]; reactions by +% the stoichiometry of its reactants and products; genes by their gene +% name. This function can therefore be used to add new entities in the +% model, or curate those already existing in the model. +% +% Input: +% model RAVEN model structure to be curated. +% metsInfo path to a *.tsv file with metabolite information, or +% 'none' to skip metabolite curation. Columns: +% metNames, comps, formula, charge, inchi, metNotes, +% then any number of MIRIAM-namespace columns. +% genesInfo path to a *.tsv file with gene information, or +% 'none'. Columns: genes, geneShortNames, then MIRIAM. +% rxnsCoeffs path to a *.tsv file with reaction stoichiometric +% coefficients, or 'none'. Columns: rxnIdx, rxnNames, +% metNames, comps, coefficient. One row per +% (reaction, metabolite) pair. +% rxnsInfo path to a *.tsv file with reaction information, or +% 'none'. Columns: rxnIdx, rxnNames, grRules, lb, ub, +% rev, subSystems, eccodes, rxnNotes, rxnReferences, +% rxnConfidenceScores, then MIRIAM. +% metPrefix prefix used to mint fresh metabolite ids (e.g. 's_' +% for yeast-GEM, 'M_' for the cobrapy/BiGG default). +% Default: 'M_'. +% rxnPrefix prefix used to mint fresh reaction ids. Default: 'R_'. +% +% Output: +% newModel curated RAVEN model structure. +% +% The 'everything after the core columns is MIRIAM' convention applies +% to all three info tables: any column whose header is not one of the +% listed core fields is treated as a MIRIAM annotation namespace and +% stored on the matching entity. +% +% Usage: newModel = curateModelFromTables(model, metsInfo, genesInfo, ... +% rxnsCoeffs, rxnsInfo, metPrefix, rxnPrefix) + +if nargin==4 + error('Provide both a ''rxnsInfo'' and a ''rxnsCoeffs'' file') +end +if nargin<4 + rxnsInfo='none'; + rxnsCoeffs='none'; +end +if nargin<3 + genesInfo='none'; +end +if nargin<6 || isempty(metPrefix) + metPrefix = 'M_'; +end +if nargin<7 || isempty(rxnPrefix) + rxnPrefix = 'R_'; +end +newModel=model; + +%% Metabolites +if ~strcmp(metsInfo,'none') + fid = fopen(metsInfo); + raw = textscan(fid,['%q' repmat(' %q',1,17)],'Delimiter','\t'); + fclose(fid); + metsToAdd.metNames = raw{1}(2:end); + metsToAdd.compartments = raw{2}(2:end); + metsToAdd.metFormulas = raw{3}(2:end); + metsToAdd.metCharges = cellfun(@str2num,raw{4}(2:end),'UniformOutput',false); + emptyEntry = cellfun(@isempty,metsToAdd.metCharges); + if all(emptyEntry) + metsToAdd = rmfield(metsToAdd,'metCharges'); + elseif any(emptyEntry) % If some charges are given, assume 0 for those without charges specified + metsToAdd.metCharges(emptyEntry) = {0}; + metsToAdd.metCharges = cell2mat(metsToAdd.metCharges); + else + metsToAdd.metCharges = cell2mat(metsToAdd.metCharges); + end + metsToAdd.inchis = raw{5}(2:end); + metsToAdd.metNotes = raw{6}(2:end); + + % Check if metabolite already exists (check by metName[comp]) + existingMets = []; + existingMetsIdx = []; + newMets = []; + for i=1:numel(metsToAdd.metNames) + metIdx = find(strcmp(newModel.metNames,metsToAdd.metNames{i})); + existMet = strcmp(metsToAdd.compartments{i},newModel.comps(newModel.metComps(metIdx))); + if any(existMet) + existingMets = [existingMets, i]; + existingMetsIdx = [existingMetsIdx, metIdx(existMet)]; + else + newMets = [newMets, i]; + end + end + + % Overwrite annotation of existing entries + if any(existingMets) + if isfield(newModel,'metFormulas') + newModel.metFormulas(existingMetsIdx) = metsToAdd.metFormulas(existingMets); + end + if isfield(newModel,'metCharges') && isfield(metsToAdd,'metCharges') + newModel.metCharges(existingMetsIdx) = metsToAdd.metCharges(existingMets); + end + if isfield(newModel,'inchis') + newModel.inchis(existingMetsIdx) = metsToAdd.inchis(existingMets); + end + newModel = extractAndAddMiriam(newModel,raw(7:end),existingMets,existingMetsIdx,'met'); + warning(['The following metabolites are already present in the model, '... + 'their annotation will be overwritten to match the metsInfo file. '... + 'If you do not particularly want to curate their annotations, it '... + 'would be better to removes these metabolites from metsInfo:\n\t\t%s'],... + strjoin(strcat(metsToAdd.metNames(existingMets),'[',metsToAdd.compartments(existingMets),']'),'\n\t\t')); + end + + % Continue with new metabolites + if numel(newMets)>0 + metsToAdd.metNames(existingMets) = []; + metsToAdd.compartments(existingMets) = []; + metsToAdd.metFormulas(existingMets) = []; + if all(cellfun(@isempty,metsToAdd.metFormulas)) + metsToAdd = rmfield(metsToAdd,'metFormulas'); + end + if isfield(metsToAdd,'metCharges') + metsToAdd.metCharges(existingMets) = []; + end + metsToAdd.inchis(existingMets) = []; + if all(cellfun(@isempty,metsToAdd.inchis)) + metsToAdd = rmfield(metsToAdd,'inchis'); + end + metsToAdd.metNotes(existingMets) = []; + if all(cellfun(@isempty,metsToAdd.metNotes)) + metsToAdd = rmfield(metsToAdd,'metNotes'); + end + % Add metabolites + newModel = addMets(newModel,metsToAdd,true,metPrefix); + addedIdx = numel(newModel.mets)-numel(newMets)+1:numel(newModel.mets); + newModel = extractAndAddMiriam(newModel,raw(7:end),newMets,addedIdx,'met'); + end +end +%% Genes +if ~strcmp(genesInfo,'none') + % Gather all data, first about reaction stoichiometries... + fid = fopen(genesInfo); + raw = textscan(fid,['%q' repmat(' %q',1,10)],'Delimiter','\t'); + fclose(fid); + genesToAdd.genes = raw{1}(2:end); + genesToAdd.geneShortNames = raw{2}(2:end); + + [~,notNewGene,existingGene] = intersect(genesToAdd.genes,newModel.genes); + if ~isempty(notNewGene) + warning(['The following genes are already present in the model, their annotation '... + 'will be overwritten to match the genesInfo file: \n\t\t%s'],... + strjoin(genesToAdd.genes(notNewGene),'\n\t\t')); + if isfield(newModel,'geneShortNames') + newModel.geneShortNames(existingGene) = genesToAdd.geneShortNames(notNewGene); + end + newModel = extractAndAddMiriam(newModel,raw(3:end),notNewGene,existingGene,'gene'); + end + + % Continue with new genes + toAdd = 1:numel(genesToAdd.genes); + toAdd(notNewGene) = []; + if ~isempty(toAdd) + genesToAdd.genes = genesToAdd.genes(toAdd); + genesToAdd.geneShortNames = genesToAdd.geneShortNames(toAdd); + if all(cellfun(@isempty,genesToAdd.geneShortNames)) + genesToAdd = rmfield(genesToAdd,'geneShortNames'); + end + % Add genes + newModel = addGenesRaven(newModel,genesToAdd); + addedIdx = numel(newModel.genes)-numel(toAdd)+1:numel(newModel.genes); + newModel = extractAndAddMiriam(newModel,raw(3:end),toAdd,addedIdx,'gene'); + end +end +%% Reactions +if ~any(strcmp({rxnsCoeffs,rxnsInfo},'none')) + % Gather all data, first about reaction stoichiometries... + fid = fopen(rxnsCoeffs); + raw = textscan(fid,'%q %q %q %q %f','Delimiter','\t','HeaderLines',1); + fclose(fid); + rxnCheck.coeffsIdx = str2double(raw{1}); + rxns.rxnNames = raw{2}; + rxns.metNames = raw{3}; + rxns.comps = raw{4}; + rxns.coeffs = raw{5}; + rxnCheck.coeffs = strcat(raw{1},'***',rxns.rxnNames); + + % ... and then additional reaction-specific data. + fid = fopen(rxnsInfo); + raw = textscan(fid,['%q' repmat(' %q',1,20)],'Delimiter','\t'); + fclose(fid); + rxnCheck.rxnsIdx = str2double(raw{1}(2:end)); + rxnCheck.rxns = strcat(raw{1}(2:end),'***',raw{2}(2:end)); + + notMatching = setxor(rxnCheck.coeffs,rxnCheck.rxns); + if numel(notMatching)>1 + error(['The following reactions ánd/or their indices are not matched '... + 'between the rxnsInfo and rxnsCoeffs files:\n\t\t%s'],... + strjoin(regexprep(notMatching,'^\d+***',''),'\n\t\t')) + end + + rxnsToAdd.rxnNames = raw{2}(2:end); + rxnsToAdd.grRules = raw{3}(2:end); + rxnsToAdd.lb = cellfun(@str2num,raw{4}(2:end)); + rxnsToAdd.ub = cellfun(@str2num,raw{5}(2:end)); + rxnsToAdd.rev = cellfun(@str2num,raw{6}(2:end)); + rxnsToAdd.subSystems = raw{7}(2:end); + rxnsToAdd.eccodes = raw{8}(2:end); + rxnsToAdd.rxnNotes = raw{9}(2:end); + rxnsToAdd.rxnReferences = raw{10}(2:end); + rxnsToAdd.rxnConfidenceScores = cellfun(@str2num,raw{11}(2:end),'UniformOutput',false); + emptyEntry = cellfun(@isempty,rxnsToAdd.rxnConfidenceScores); + if all(emptyEntry) + rxnsToAdd = rmfield(rxnsToAdd,'rxnConfidenceScores'); + elseif any(emptyEntry) % If some rxnConfidenceScores are given, assume 0 for those without rxnConfidenceScores specified + rxnsToAdd.rxnConfidenceScores(emptyEntry) = {0}; + rxnsToAdd.rxnConfidenceScores = cell2mat(rxnsToAdd.rxnConfidenceScores); + else + rxnsToAdd.rxnConfidenceScores = cell2mat(rxnsToAdd.rxnConfidenceScores); + end + + existingRxn=[]; + notNewRxn=[]; + for i=1:numel(rxnCheck.rxnsIdx) + rxnRows = find(rxnCheck.rxnsIdx(i)==rxnCheck.coeffsIdx); + rxnsToAdd.mets{i,1} = cell(1,1); + rxnsToAdd.stoichCoeffs{i,1} = cell(1,1); + for j=1:numel(rxnRows) + newMetComp = [rxns.metNames{rxnRows(j)}, '[', rxns.comps{rxnRows(j)}, ']']; + try + newMetComp = getIndexes(newModel,newMetComp,'metcomps'); + catch + error(['Not all metabolites in reaction "%s" are present in the '... + 'model or in the provided table with new metabolites.'],rxnsToAdd.rxnNames{i}) + end + rxnsToAdd.mets{i}(j) = newModel.mets(newMetComp); + rxnsToAdd.stoichCoeffs{i}{j} = rxns.coeffs(rxnRows(j)); + end + rxnsToAdd.stoichCoeffs{i}=cell2mat(rxnsToAdd.stoichCoeffs{i}); + %Check if the reaction not already exists, by stoichiometry of its + %products and reactants + modelCoeffs = transpose(newModel.S(getIndexes(newModel,rxnsToAdd.mets{i},'mets'),:)); + modelCoeffs2 = find(~all(modelCoeffs==0,2)); + modelCoeffs = modelCoeffs(modelCoeffs2,:); + [~, duplicateRxn] = intersect(modelCoeffs,rxnsToAdd.stoichCoeffs{i},'rows'); + if ~isempty(duplicateRxn) + existingRxn = [existingRxn,modelCoeffs2(duplicateRxn)]; + notNewRxn = [notNewRxn,i]; + end + end + + % If some reactions already existed, then replace its information + if ~isempty(existingRxn) + warning(['The following reactions are the same as existing model '... + 'reactions, their annotation will be overwritten to match the '... + 'rxnCoeffs file: \n\t\t(Existing reaction ID): New reaction name\n\t\t%s'],... + strjoin(strcat('(',newModel.rxns(existingRxn),{'): '},rxnsToAdd.rxnNames(notNewRxn)),'\n\t\t')); + newModel.rxnNames(existingRxn) = rxnsToAdd.rxnNames(notNewRxn); + newModel.lb(existingRxn) = rxnsToAdd.lb(notNewRxn); + newModel.ub(existingRxn) = rxnsToAdd.ub(notNewRxn); + newModel.rev(existingRxn) = rxnsToAdd.rev(notNewRxn); + if isfield(newModel,'subSystems') + newModel.subSystems(existingRxn) = rxnsToAdd.subSystems(notNewRxn); + end + if isfield(newModel,'eccodes') + newModel.eccodes(existingRxn) = rxnsToAdd.eccodes(notNewRxn); + end + if isfield(newModel,'rxnNotes') + newModel.rxnNotes(existingRxn) = rxnsToAdd.rxnNotes(notNewRxn); + end + if isfield(newModel,'rxnConfidenceScores') && isfield(rxnsToAdd,'rxnConfidenceScores') + newModel.rxnConfidenceScores(existingRxn) = rxnsToAdd.rxnConfidenceScores(notNewRxn); + end + emptyEntries = cellfun(@isempty,rxnsToAdd.grRules); + if ~all(emptyEntries) + newModel = changeGrRules(newModel,newModel.rxns(existingRxn(~emptyEntries)),rxnsToAdd.grRules(notNewRxn(~emptyEntries)),true); + end + newModel = extractAndAddMiriam(newModel,raw(11:end),notNewRxn,existingRxn,'rxn'); + end + + % Continue with new reactions + toAdd = 1:numel(rxnsToAdd.rxnNames); + toAdd(notNewRxn) = []; + if ~isempty(toAdd) + rxnsToAdd.rxnNames = rxnsToAdd.rxnNames(toAdd); + rxnsToAdd.lb = rxnsToAdd.lb(toAdd); + rxnsToAdd.ub = rxnsToAdd.ub(toAdd); + rxnsToAdd.rev = rxnsToAdd.rev(toAdd); + rxnsToAdd.subSystems = rxnsToAdd.subSystems(toAdd); + if all(cellfun(@isempty,rxnsToAdd.subSystems)) + rxnsToAdd = rmfield(rxnsToAdd,'subSystems'); + end + rxnsToAdd.eccodes = rxnsToAdd.eccodes(toAdd); + if all(cellfun(@isempty,rxnsToAdd.eccodes)) + rxnsToAdd = rmfield(rxnsToAdd,'eccodes'); + end + rxnsToAdd.rxnNotes = rxnsToAdd.rxnNotes(toAdd); + if all(cellfun(@isempty,rxnsToAdd.rxnNotes)) + rxnsToAdd = rmfield(rxnsToAdd,'rxnNotes'); + end + if isfield(rxnsToAdd,'rxnConfidenceScores') + rxnsToAdd.rxnConfidenceScores = rxnsToAdd.rxnConfidenceScores(toAdd); + end + rxnsToAdd.rxnReferences = rxnsToAdd.rxnReferences(toAdd); + if all(cellfun(@isempty,rxnsToAdd.rxnReferences)) + rxnsToAdd = rmfield(rxnsToAdd,'rxnReferences'); + end + rxnsToAdd.grRules = rxnsToAdd.grRules(toAdd); + rxnsToAdd.mets = rxnsToAdd.mets(toAdd); + rxnsToAdd.stoichCoeffs = rxnsToAdd.stoichCoeffs(toAdd); + rxnsToAdd.rxns = generateNewIds(model,'rxns',rxnPrefix,numel(rxnsToAdd.rxnNames)); + + newModel = addRxns(newModel,rxnsToAdd,1,[],false,false); + rxnsModelIdx = numel(newModel.rxns)-numel(toAdd)+1:numel(newModel.rxns); + newModel = extractAndAddMiriam(newModel,raw(11:end),toAdd,rxnsModelIdx,'rxn'); + end +end +end + +function newModel = extractAndAddMiriam(model,raw,inputIndex,modelIndex,type) +newModel=model; +miriamName = cell(numel(inputIndex),1); +miriamValues = cell(1,numel(inputIndex)); +for i=1:numel(raw) + miriamName{i} = raw{i}{1}; + miriamValues(1:numel(inputIndex),i) = raw{i}(inputIndex+1); +end +emptyMiriam = all(cellfun(@isempty,miriamValues),1); +miriamName(emptyMiriam) = []; +miriamValues(:,emptyMiriam) = []; +if ~isfield(newModel,[type 'Miriams']); + newModel.([type 'Miriams'])=cell(numel(newModel.([type 's'])),1); +end +if ~isempty(miriamName) + for i=1:numel(miriamName) + newModel = editMiriam(newModel,type,modelIndex,miriamName{i},miriamValues(:,i),'replace'); + end +end +end diff --git a/core/findDuplicateRxns.m b/core/findDuplicateRxns.m new file mode 100644 index 00000000..44d3011b --- /dev/null +++ b/core/findDuplicateRxns.m @@ -0,0 +1,40 @@ +function pairs = findDuplicateRxns(model, ignoreDirection) +% findDuplicateRxns +% Find reactions that share identical stoichiometry. Counterpart of +% raven_python.manipulation.find_duplicate_reactions, and the +% upstream version of yeast-GEM's findDuplicatedRxns. +% +% Only stoichiometry is compared — bounds, GPRs, and annotations +% are ignored. The default treats A→B and B→A as duplicates +% (typical curation use case: "find reactions that could be +% merged"). +% +% Inputs: +% model RAVEN model struct. +% ignoreDirection (opt, default true) Treat A→B and B→A as +% duplicates. +% +% Output: +% pairs Nx2 numeric array of reaction-index pairs +% (i, j) where reactions i and j share the +% same (possibly negated) stoichiometry, with +% i < j. Empty if the model has no duplicates. +% +% Usage: pairs = findDuplicateRxns(model) +% pairs = findDuplicateRxns(model, false) + +if nargin < 2 + ignoreDirection = true; +end + +pairs = zeros(0, 2); +n = numel(model.rxns); +for i = 1:n-1 + for j = i+1:n + if isequal(model.S(:, i), model.S(:, j)) || ... + (ignoreDirection && isequal(model.S(:, i), -model.S(:, j))) + pairs(end+1, :) = [i, j]; %#ok + end + end +end +end diff --git a/core/getBiomassFractions.m b/core/getBiomassFractions.m new file mode 100644 index 00000000..592e7830 --- /dev/null +++ b/core/getBiomassFractions.m @@ -0,0 +1,125 @@ +function fractions = getBiomassFractions(model, biomassConfig) +% getBiomassFractions +% Compute the mass fraction (g/gDW) per biomass component plus the +% total. Mirrors raven_python.biomass.sum_biomass; the MATLAB +% counterpart of yeast-GEM's legacy sumBioMass. +% +% The biomassConfig struct describes the per-organism biomass +% layout — see "Inputs" below. Components whose pseudoreaction is +% missing from the model contribute 0. +% +% Inputs: +% model RAVEN model struct. +% biomassConfig struct with fields: +% biomass_rxn rxn id of the top-level +% biomass pseudoreaction. +% proton_met met id of cytosolic H+ (used +% only by rescalePseudoreaction; +% may be unused here). +% components cell array of component +% structs with fields: +% .name component name +% (e.g. 'protein'). +% .pseudoreaction_name model.rxnNames +% entry to identify +% the pseudoreaction. +% .mass_strategy 'mw' | 'mw_minus_2h' +% | 'mw_minus_water' +% | 'grams' — see +% NOTES below. +% +% Output: +% fractions struct keyed by component name plus 'total': +% fractions.protein, fractions.RNA, ... etc. +% All values are in g/gDW. +% +% NOTES on mass_strategy: +% 'mw' MW from chemical formula +% 'mw_minus_2h' MW − 2.016 g/mol (two protons released per +% charged tRNA — protein-pseudoreaction substrates) +% 'mw_minus_water' MW − 18.015 g/mol (water released per +% polymerisation step — RNA / DNA) +% 'grams' stoichiometry already in g/gDW (lipid backbone) +% +% Usage: fractions = getBiomassFractions(model, biomassConfig) + +fractions = struct(); +total = 0; +for i = 1:numel(biomassConfig.components) + comp = biomassConfig.components{i}; + f = computeComponentFraction(model, comp); + fractions.(comp.name) = f; + total = total + f; +end +fractions.total = total; +end + +function f = computeComponentFraction(model, comp) +rxnPos = strcmp(model.rxnNames, comp.pseudoreaction_name); +if ~any(rxnPos) + f = 0; + return; +end +S_col = model.S(:, rxnPos); +isSub = find(S_col < 0); +if isempty(isSub) + f = 0; + return; +end +if strcmp(comp.mass_strategy, 'grams') + f = full(-sum(S_col(isSub))); + return; +end +offset = mwOffset(comp.mass_strategy); +formulas = model.metFormulas(isSub); +MWs = zeros(numel(formulas), 1); +for i = 1:numel(formulas) + MWs(i) = computeFormulaMW(formulas{i}); +end +zeroMW = MWs == 0; +if any(zeroMW) + error('getBiomassFractions:emptyFormula', ... + 'Biomass metabolite %s has an empty metFormula field.', ... + model.mets{isSub(find(zeroMW, 1))}); +end +MWs = MWs + offset; +f = full(-sum(S_col(isSub) .* MWs) / 1000); +end + +function offset = mwOffset(strategy) +switch strategy + case 'mw' + offset = 0; + case 'mw_minus_2h' + offset = -2.016; + case 'mw_minus_water' + offset = -18.015; + otherwise + error('getBiomassFractions:unknownStrategy', ... + 'Unknown mass_strategy: %s', strategy); +end +end + +function mw = computeFormulaMW(formula) +% Molecular weight in g/mol from a Hill-style chemical formula. Same +% element table as the legacy yeast-GEM sumBioMass. +tokens = regexp(formula, '([A-Z][a-z]*)(\d*)', 'tokens'); +if isempty(tokens) + mw = 0; + return; +end +tokensMatrix = vertcat(tokens{:}); +tokensMatrix(cellfun(@isempty, tokensMatrix(:,2)), 2) = {'1'}; +elements = tokensMatrix(:, 1); +counts = str2double(tokensMatrix(:, 2)); +elem = {'C', 12.01; 'H', 1.008; 'N', 14.007; 'O', 15.999; ... + 'P', 30.974; 'S', 32.06; 'R', 0; ... + 'Fe', 55.845; 'K', 39.098; 'Na', 22.99; 'Cl', 35.45; ... + 'Mn', 54.938; 'Zn', 65.38; 'Ca', 40.078; 'Mg', 24.305; 'Cu', 63.546}; +[~, elemMatch] = ismember(elements, elem(:,1)); +if any(elemMatch == 0) + error('getBiomassFractions:unknownElement', ... + 'Unknown element in formula %s', formula); +end +mw = sum(counts .* transpose([elem{elemMatch, 2}]), 'all'); +end diff --git a/core/scaleBiomassFraction.m b/core/scaleBiomassFraction.m new file mode 100644 index 00000000..0e4165f2 --- /dev/null +++ b/core/scaleBiomassFraction.m @@ -0,0 +1,50 @@ +function model = scaleBiomassFraction(model, biomassConfig, componentName, newValue, balanceOut) +% scaleBiomassFraction +% Rescale a biomass component to a target g/gDW value, optionally +% balancing a second component so the total biomass mass stays at +% 1 g/gDW. Mirrors raven_python.biomass.scale_biomass and yeast-GEM's +% legacy scaleBioMass. +% +% Inputs: +% model RAVEN model struct. +% biomassConfig struct (see getBiomassFractions). +% componentName Component to rescale. +% newValue Target fraction in g/gDW. +% balanceOut (opt) Second component name to adjust so the +% biomass total remains 1 g/gDW. Empty / omit +% to skip balancing. +% +% Output: +% model Modified model. +% +% Usage: model = scaleBiomassFraction(model, biomassConfig, 'protein', 0.5, 'carbohydrate') + +if nargin < 5 + balanceOut = ''; +end + +fractions = getBiomassFractions(model, biomassConfig); +if ~isfield(fractions, componentName) + error('scaleBiomassFraction:unknownComponent', ... + 'biomassConfig has no component named %s', componentName); +end +current = fractions.(componentName); +if current == 0 + error('scaleBiomassFraction:zeroCurrent', ... + ['Cannot scale %s to %g: current fraction is 0 ' ... + '(pseudoreaction missing or empty).'], componentName, newValue); +end +factor = newValue / current; +model = scaleBiomassPseudoreaction(model, biomassConfig, componentName, factor); + +if ~isempty(balanceOut) + fractions = getBiomassFractions(model, biomassConfig); + total = fractions.total; + balanceCurrent = fractions.(balanceOut); + if balanceCurrent == 0 + return; + end + balanceFactor = (balanceCurrent + (1 - total)) / balanceCurrent; + model = scaleBiomassPseudoreaction(model, biomassConfig, balanceOut, balanceFactor); +end +end diff --git a/core/scaleBiomassPseudoreaction.m b/core/scaleBiomassPseudoreaction.m new file mode 100644 index 00000000..4636297c --- /dev/null +++ b/core/scaleBiomassPseudoreaction.m @@ -0,0 +1,58 @@ +function model = scaleBiomassPseudoreaction(model, biomassConfig, componentName, factor) +% scaleBiomassPseudoreaction +% Multiply the substrate coefficients of one biomass component +% pseudoreaction by `factor` and rebalance H+ to preserve charge +% neutrality. Mirrors raven_python.biomass.rescale_pseudoreaction +% and yeast-GEM's legacy rescalePseudoReaction. +% +% "Substrate" means every metabolite in the pseudoreaction whose +% metabolite name does NOT match the component name (the +% component's product is left untouched). After rescaling, the +% coefficient of biomassConfig.proton_met is recomputed so the +% pseudoreaction's total ionic charge sums to zero. +% +% Inputs: +% model RAVEN model struct. +% biomassConfig struct (see getBiomassFractions). +% componentName Name of the component to rescale (must match +% biomassConfig.components{i}.name for some i, +% AND be the model.metNames of the produced +% metabolite in the matching pseudoreaction). +% factor Multiplicative factor. +% +% Output: +% model Modified model. +% +% Usage: model = scaleBiomassPseudoreaction(model, biomassConfig, 'protein', 0.9) + +comp = findComponent(biomassConfig, componentName); +rxnPos = find(strcmp(model.rxnNames, comp.pseudoreaction_name)); +if isempty(rxnPos) + error('scaleBiomassPseudoreaction:missingPseudoreaction', ... + 'No reaction named %s in model.', comp.pseudoreaction_name); +end + +for i = 1:length(model.mets) + S_ir = model.S(i, rxnPos); + isProd = strcmp(model.metNames{i}, componentName); + if S_ir ~= 0 && ~isProd + model.S(i, rxnPos) = factor * S_ir; + end +end + +% Rebalance H+ to keep charge neutrality. +Hc = find(strcmp(model.mets, biomassConfig.proton_met)); +model.S(Hc, rxnPos) = 0; +model.S(Hc, rxnPos) = -sum(model.S(:, rxnPos) .* model.metCharges, 'omitnan'); +end + +function comp = findComponent(cfg, name) +for i = 1:numel(cfg.components) + if strcmp(cfg.components{i}.name, name) + comp = cfg.components{i}; + return; + end +end +error('scaleBiomassPseudoreaction:unknownComponent', ... + 'biomassConfig has no component named %s', name); +end diff --git a/core/setGAM.m b/core/setGAM.m new file mode 100644 index 00000000..319c04df --- /dev/null +++ b/core/setGAM.m @@ -0,0 +1,54 @@ +function model = setGAM(model, value, biomassRxn, cofactorMetNames, ngamRxn, ngamValue) +% setGAM +% Set the growth-associated maintenance (GAM) coefficient in the +% biomass pseudoreaction, and optionally fix the non-growth +% maintenance (NGAM) reaction's bounds. Mirrors +% raven_python.biomass.set_gam and yeast-GEM's legacy changeGAM. +% +% For every metabolite in the biomass pseudoreaction whose +% `model.metNames` entry is in `cofactorMetNames`, the +% stoichiometric coefficient is set to ±`value` preserving the sign +% of the current coefficient. Yeast-GEM scales ATP, ADP, H2O, H+ +% and phosphate (with ATP and H2O on the substrate side, ADP / H+ / +% phosphate on the product side). +% +% Inputs: +% model RAVEN model struct. +% value New GAM value (mmol ATP / gDW per growth unit). +% biomassRxn Reaction id of the biomass pseudoreaction. +% cofactorMetNames Cell array of metabolite NAMES (not IDs) +% to rescale, e.g. {'ATP','ADP','H2O','H+', +% 'phosphate'}. +% ngamRxn (opt) NGAM reaction id. Required when +% ngamValue is supplied. +% ngamValue (opt) NGAM flux to fix. Sets the NGAM +% reaction's bounds to (ngamValue, ngamValue). +% +% Output: +% model Modified model. +% +% Usage: model = setGAM(model, 80, 'r_4041', {'ATP','ADP','H2O','H+','phosphate'}) + +if nargin < 4 + error('setGAM:missingArgs', ... + 'biomassRxn and cofactorMetNames are required.'); +end + +bioPos = strcmp(model.rxns, biomassRxn); +if ~any(bioPos) + error('setGAM:missingBiomassRxn', ... + 'Reaction %s not found in model.', biomassRxn); +end + +for i = 1:length(model.mets) + S_ix = model.S(i, bioPos); + isCofactor = any(strcmp(cofactorMetNames, model.metNames{i})); + if S_ix ~= 0 && isCofactor + model.S(i, bioPos) = sign(S_ix) * value; + end +end + +if nargin >= 6 && ~isempty(ngamRxn) && ~isempty(ngamValue) + model = setParam(model, 'eq', ngamRxn, ngamValue); +end +end diff --git a/io/loadDeltaGfromCSV.m b/io/loadDeltaGfromCSV.m new file mode 100644 index 00000000..da56e766 --- /dev/null +++ b/io/loadDeltaGfromCSV.m @@ -0,0 +1,61 @@ +function model = loadDeltaGfromCSV(model, metCsv, rxnCsv) +% loadDeltaGfromCSV +% Populate model.metDeltaG and model.rxnDeltaG from project CSV +% files. Mirrors raven_python.annotation.load_delta_g_csv and is +% the upstream version of yeast-GEM's loadDeltaG. +% +% Each CSV is a two-column table: identifier, deltaG. Rows whose +% identifier doesn't appear in the model are silently skipped. +% Pass an empty string ('') for either argument to skip that side. +% +% Inputs: +% model RAVEN model struct. +% metCsv Path to metabolite ΔG CSV (id, ΔG), or '' to skip. +% rxnCsv Path to reaction ΔG CSV (id, ΔG), or '' to skip. +% +% Output: +% model Model with metDeltaG and/or rxnDeltaG fields added. +% +% Usage: model = loadDeltaGfromCSV(model, ... +% 'data/databases/model_metDeltaG.csv', ... +% 'data/databases/model_rxnDeltaG.csv') + +if nargin < 3 + rxnCsv = ''; +end +if nargin < 2 + metCsv = ''; +end + +if ~isempty(metCsv) + if isfield(model, 'metDeltaG') + disp('Existing metDeltaG field will be overwritten.') + else + model.metDeltaG = nan(numel(model.mets), 1); + end + metG = readtable(metCsv); + [a, b] = ismember(model.mets, metG.(metG.Properties.VariableNames{1})); + model.metDeltaG(a) = metG.(metG.Properties.VariableNames{2})(b(a)); + if any(~a) + fprintf(['Not all metabolite identifiers are matched to %s; the latter\n' ... + 'file might have to be supplemented with deltaG values for new metabolites.\n'], ... + metCsv); + end +end + +if ~isempty(rxnCsv) + if isfield(model, 'rxnDeltaG') + disp('Existing rxnDeltaG field will be overwritten.') + else + model.rxnDeltaG = nan(numel(model.rxns), 1); + end + rxnG = readtable(rxnCsv); + [a, b] = ismember(model.rxns, rxnG.(rxnG.Properties.VariableNames{1})); + model.rxnDeltaG(a) = rxnG.(rxnG.Properties.VariableNames{2})(b(a)); + if any(~a) + fprintf(['Not all reaction identifiers are matched to %s; the latter\n' ... + 'file might have to be supplemented with deltaG values for new reactions.\n'], ... + rxnCsv); + end +end +end diff --git a/io/parseYAML.m b/io/parseYAML.m new file mode 100644 index 00000000..5c93b76a --- /dev/null +++ b/io/parseYAML.m @@ -0,0 +1,91 @@ +function out = parseYAML(filename) +% parseYAML +% Read an arbitrary YAML file into a MATLAB struct / cell tree. +% +% Use this for parsing arbitrary YAML configuration / data files +% (e.g. yeast-GEM's data/conditions/*.yml). For loading a cobra-format +% model YAML, use readYAMLmodel instead — that function knows the +% model schema and returns a populated RAVEN model struct. +% +% Implementation: delegates to Python's yaml.safe_load, then +% recursively converts the py.dict / py.list tree to native MATLAB +% struct / cell. Requires a working MATLAB-Python bridge and the +% pyyaml package in the linked Python environment: +% +% pip install pyyaml % from the MATLAB-linked Python env +% +% Input: +% filename path to the YAML file. +% +% Output: +% out MATLAB representation of the document: +% py.dict -> struct +% py.list -> cell column vector +% py.str -> char +% py.int -> double +% py.float -> double +% py.bool -> logical +% py.None -> [] +% +% Usage: cfg = parseYAML('data/conditions/anaerobic.yml') + +if ~isfile(filename) + error('parseYAML:fileNotFound', 'File not found: %s', filename); +end + +try + py.importlib.import_module('yaml'); +catch ME + error('parseYAML:pyyamlMissing', ... + ['pyyaml is required to read arbitrary YAML files. Install it ' ... + 'in your MATLAB-linked Python environment (`pip install pyyaml`).' ... + '\nUnderlying error: %s'], ME.message); +end + +f = py.builtins.open(filename, 'r'); +cleanup = onCleanup(@() f.close()); +data = py.yaml.safe_load(f); + +out = pyToMatlab(data); +end + + +function v = pyToMatlab(obj) +% Recursively convert pyyaml-loaded Python objects into MATLAB types. +if isa(obj, 'py.NoneType') + v = []; +elseif isa(obj, 'py.bool') + v = logical(obj); +elseif isa(obj, 'py.int') || isa(obj, 'py.float') + v = double(obj); +elseif isa(obj, 'py.str') + v = char(obj); +elseif isa(obj, 'py.dict') + v = struct(); + keys = cell(py.list(obj.keys())); + vals = cell(py.list(obj.values())); + for i = 1:numel(keys) + v.(matlabFieldName(char(keys{i}))) = pyToMatlab(vals{i}); + end +elseif isa(obj, 'py.list') || isa(obj, 'py.tuple') + cells = cell(obj); + v = cell(numel(cells), 1); + for i = 1:numel(cells) + v{i} = pyToMatlab(cells{i}); + end +else + % Fallback: best-effort + v = obj; +end +end + + +function name = matlabFieldName(key) +% Sanitise a YAML key into a valid MATLAB field name. Replaces non- +% alphanumeric characters with underscores; prefixes a digit-starting +% key with 'f_'. +name = regexprep(key, '[^A-Za-z0-9_]', '_'); +if isempty(name) || ~isstrprop(name(1), 'alpha') + name = ['f_' name]; +end +end diff --git a/io/readYAMLmodel.m b/io/readYAMLmodel.m index 1eb65f3b..f0027f0e 100755 --- a/io/readYAMLmodel.m +++ b/io/readYAMLmodel.m @@ -110,8 +110,9 @@ model.(modelFields{i,1})=modelFields{i,2}; end -% If GECKO model -if any(contains(line_key,'geckoLight')) +% If GECKO model — accept both the legacy `geckoLight` (inside metaData) +% and the cobrapy / raven_python style top-level `gecko_light` key. +if any(contains(line_key,'geckoLight')) || any(contains(line_key,'gecko_light')) isGECKO=true; ecFields = {'geckoLight', false;... 'rxns', {};... @@ -146,9 +147,10 @@ tline_raw = line_raw{i}; tline_key = line_key{i}; tline_value = line_value{i}; - % import different sections + % import different sections — accept the !!omap-tagged variant as + % well so cobrapy / raven_python output is recognized. switch tline_raw - case '- metaData:' + case {'- metaData:', '- metaData: !!omap'} section = 1; if verbose fprintf('\t%d\n', section); @@ -198,6 +200,31 @@ continue end + % cobrapy-style root-level keys (id, name, version, gecko_light). + % Cobra writes these at the top level (no metaData section); RAVEN's + % own writer normally puts them inside metaData, but we accept both + % so a cobra-written YAML can be ingested directly. + if isempty(regexp(tline_raw, '^ {2,}', 'once')) + switch tline_key + case 'id' + if isempty(model.id), model.id = tline_value; end + continue + case 'name' + if isempty(model.name), model.name = tline_value; end + continue + case 'version' + if ~isfield(model,'version') || isempty(model.version) + model.version = tline_value; + end + continue + case 'gecko_light' + if isGECKO && strcmp(tline_value,'true') + model.ec.geckoLight = true; + end + continue + end + end + % skip over empty keys if isempty(tline_raw) || (isempty(tline_key) && contains(tline_raw,'!!omap')) continue; @@ -234,7 +261,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 +294,12 @@ model = readFieldValue(model, 'inchis', tline_value, pos); readList=''; miriamKey=''; case 'smiles' + % Top-level (legacy MATLAB) and inside-annotation + % (cobrapy / current writer) layouts both land here. + % Don't reset readList — preserves the annotation + % gathering state if SMILES is nested inside the + % annotation block alongside other entries. model = readFieldValue(model, 'metSmiles', tline_value, pos); - readList=''; miriamKey=''; case 'deltaG' model = readFieldValue(model, 'metDeltaG', tline_value, pos); readList=''; miriamKey=''; @@ -309,7 +340,11 @@ case 'gene_reaction_rule' model = readFieldValue(model, 'grRules', tline_value, pos); readList=''; miriamKey=''; - case 'rxnNotes' + case {'notes','rxnNotes'} + % `notes` is the canonical reaction-side key (matches + % cobrapy and raven_python); `rxnNotes` is the legacy + % RAVEN MATLAB key, kept here for backward-compatible + % reads. model = readFieldValue(model, 'rxnNotes', tline_value, pos); readList=''; miriamKey=''; case 'rxnFrom' diff --git a/io/saveDeltaGtoCSV.m b/io/saveDeltaGtoCSV.m new file mode 100644 index 00000000..ee682940 --- /dev/null +++ b/io/saveDeltaGtoCSV.m @@ -0,0 +1,59 @@ +function saveDeltaGtoCSV(model, metCsv, rxnCsv, verbose) +% saveDeltaGtoCSV +% Persist model.metDeltaG and model.rxnDeltaG to project CSV files. +% Counterpart of loadDeltaGfromCSV and the upstream version of +% yeast-GEM's saveDeltaG. Mirrors raven_python.annotation.save_delta_g_csv. +% +% Each CSV gets two columns: identifier, deltaG. Rows are written +% in model order (one row per entity); identifiers without a +% matching field get NaN. Pass an empty string for metCsv or rxnCsv +% to skip that side. +% +% Inputs: +% model RAVEN model struct. +% metCsv Output path for the metabolite ΔG CSV, or '' to skip. +% rxnCsv Output path for the reaction ΔG CSV, or '' to skip. +% verbose (opt, default false) Print "wrote ..." per file. +% +% Usage: saveDeltaGtoCSV(model, ... +% 'data/databases/model_metDeltaG.csv', ... +% 'data/databases/model_rxnDeltaG.csv') + +if nargin < 4 + verbose = false; +end +if nargin < 3 + rxnCsv = ''; +end +if nargin < 2 + metCsv = ''; +end + +if ~isempty(metCsv) + if ~isfield(model, 'metDeltaG') + if verbose + fprintf('No metDeltaG field found, %s will not be changed.\n', metCsv); + end + else + metG = array2table([model.mets, num2cell(model.metDeltaG)]); + writetable(metG, metCsv); + if verbose + fprintf('Wrote %s\n', metCsv); + end + end +end + +if ~isempty(rxnCsv) + if ~isfield(model, 'rxnDeltaG') + if verbose + fprintf('No rxnDeltaG field found, %s will not be changed.\n', rxnCsv); + end + else + rxnG = array2table([model.rxns, num2cell(model.rxnDeltaG)]); + writetable(rxnG, rxnCsv); + if verbose + fprintf('Wrote %s\n', rxnCsv); + end + end +end +end diff --git a/io/writeYAMLmodel.m b/io/writeYAMLmodel.m index 2ada7074..ae76afdd 100755 --- a/io/writeYAMLmodel.m +++ b/io/writeYAMLmodel.m @@ -1,16 +1,24 @@ function writeYAMLmodel(model,fileName,preserveQuotes,sortIds) % writeYAMLmodel -% Writes a yaml file matching (roughly) the cobrapy yaml structure +% Writes a yaml file matching cobrapy's YAML structure. The format is +% cobrapy's native !!omap layout, extended with RAVEN-only top-level +% per-entry keys (inchis, deltaG, metFrom, eccodes, rxnFrom, +% references, confidence_score, protein) and the GECKO ec-rxns / +% ec-enzymes sections. Output is byte-stable with raven_python's +% io.yaml.write_yaml_model when called with the same model. % % model a model structure -% fileName name that the file will have. A dialog window will +% fileName name that the file will have. A dialog window will % open if no file name is specified. -% preserveQuotes if quotes should be preserved for strings -% (logical, default=true) +% preserveQuotes if all string values should be wrapped in double +% quotes. cobrapy emits quotes only where YAML +% requires them, so the default is false (matches +% cobrapy / raven-python). +% (logical, default=false) % sortIds if metabolites, reactions, genes and compartments % should be sorted alphabetically by their identifier, % otherwise they are kept in their original order -% (logical, default=false) +% (logical, default=false) % % Usage: writeYAMLmodel(model,fileName,preserveQuotes,sortIds) if nargin<2|| isempty(fileName) @@ -24,7 +32,7 @@ function writeYAMLmodel(model,fileName,preserveQuotes,sortIds) fileName=char(fileName); if nargin < 3 - preserveQuotes = true; + preserveQuotes = false; end if nargin < 4 sortIds = false; @@ -65,29 +73,42 @@ function writeYAMLmodel(model,fileName,preserveQuotes,sortIds) if fid == -1 error(['Cannot write to ' fileName ', does the directory exist?']) end -fprintf(fid,'---\n!!omap\n'); +% cobrapy emits a bare `!!omap` root with no document-start marker; +% match that for byte-stable round-tripping. +fprintf(fid,'!!omap\n'); %Insert file header (metadata) -writeMetadata(model,fid); +writeMetadata(model,fid,preserveQuotes); %Metabolites: +% Field order matches cobrapy + raven_python.io.yaml: +% id, name, compartment, charge, formula, notes, annotation, +% then RAVEN-only extras (inchis, deltaG, metFrom). +% SMILES goes inside the annotation block (cobrapy convention), not at +% metabolite top level — the reader still accepts top-level `smiles:` +% for backward compatibility with older yeast-GEM files. fprintf(fid,'- metabolites:\n'); for i = 1:length(model.mets) fprintf(fid,' - !!omap\n'); writeField(model, fid, 'mets', 'txt', i, ' - id', preserveQuotes) writeField(model, fid, 'metNames', 'txt', i, ' - name', preserveQuotes) writeField(model, fid, 'metComps', 'txt', i, ' - compartment', preserveQuotes) - writeField(model, fid, 'metFormulas', 'txt', i, ' - formula', preserveQuotes) writeField(model, fid, 'metCharges', 'num', i, ' - charge', preserveQuotes) + writeField(model, fid, 'metFormulas', 'txt', i, ' - formula', preserveQuotes) + writeField(model, fid, 'metNotes', 'txt', i, ' - notes', preserveQuotes) + writeAnnotation(model, fid, 'met', i, preserveQuotes) writeField(model, fid, 'inchis', 'txt', i, ' - inchis', preserveQuotes) - writeField(model, fid, 'metSmiles', 'txt', i, ' - smiles', preserveQuotes) - writeField(model, fid, 'metMiriams', 'txt', i, ' - annotation', preserveQuotes) writeField(model, fid, 'metDeltaG', 'num', i, ' - deltaG', preserveQuotes) - writeField(model, fid, 'metNotes', 'txt', i, ' - notes', preserveQuotes) writeField(model, fid, 'metFrom', 'txt', i, ' - metFrom', preserveQuotes) end %Reactions: +% Field order matches cobrapy + raven_python.io.yaml: +% id, name, metabolites, lower_bound, upper_bound, gene_reaction_rule, +% objective_coefficient, subsystem, notes, annotation, +% then RAVEN-only extras (eccodes, references, rxnFrom, deltaG, +% confidence_score). The notes key is the canonical `notes` (no +% longer `rxnNotes`); the reader still accepts the legacy key. fprintf(fid,'- reactions:\n'); for i = 1:length(model.rxns) fprintf(fid,' - !!omap\n'); @@ -97,17 +118,17 @@ function writeYAMLmodel(model,fileName,preserveQuotes,sortIds) writeField(model, fid, 'lb', 'num', i, ' - lower_bound', preserveQuotes) writeField(model, fid, 'ub', 'num', i, ' - upper_bound', preserveQuotes) writeField(model, fid, 'grRules', 'txt', i, ' - gene_reaction_rule', preserveQuotes) - writeField(model, fid, 'rxnFrom', 'txt', i, ' - rxnFrom', preserveQuotes) if model.c(i)~=0 - writeField(model, fid, 'c', 'num', i, ' - objective_coefficient', preserveQuotes) + writeField(model, fid, 'c', 'num', i, ' - objective_coefficient', preserveQuotes) end - writeField(model, fid, 'eccodes', 'txt', i, ' - eccodes', preserveQuotes) - writeField(model, fid, 'rxnReferences', 'txt', i, ' - references', preserveQuotes) writeField(model, fid, 'subSystems', 'txt', i, ' - subsystem', preserveQuotes) + writeField(model, fid, 'rxnNotes', 'txt', i, ' - notes', preserveQuotes) writeField(model, fid, 'rxnMiriams', 'txt', i, ' - annotation', preserveQuotes) + writeField(model, fid, 'eccodes', 'txt', i, ' - eccodes', preserveQuotes) + writeField(model, fid, 'rxnReferences', 'txt', i, ' - references', preserveQuotes) + writeField(model, fid, 'rxnFrom', 'txt', i, ' - rxnFrom', preserveQuotes) writeField(model, fid, 'rxnDeltaG', 'num', i, ' - deltaG', preserveQuotes) writeField(model, fid, 'rxnConfidenceScores', 'num', i, ' - confidence_score', preserveQuotes) - writeField(model, fid, 'rxnNotes', 'txt', i, ' - rxnNotes', preserveQuotes) end %Genes: @@ -132,6 +153,12 @@ function writeYAMLmodel(model,fileName,preserveQuotes,sortIds) %EC-model: if isfield(model,'ec') + % gecko_light flag at the top level (matches + % raven_python.io.yaml — keeps the metaData block a pure provenance + % container). The reader accepts both this key and the legacy + % geckoLight key inside metaData. + if model.ec.geckoLight; geckoLightStr = 'true'; else; geckoLightStr = 'false'; end + fprintf(fid,'- gecko_light: %s\n', geckoLightStr); fprintf(fid,'- ec-rxns:\n'); for i = 1:length(model.ec.rxns) fprintf(fid,' - !!omap\n'); @@ -219,9 +246,12 @@ function writeField(model,fid,fieldName,type,pos,name,preserveQuotes) end elseif strcmp(fieldName,'S') - %S: create header & write each metabolite in a new line - fprintf(fid,' %s: !!omap\n',name); + %S: create header & write each metabolite in a new line. Reactions + %with no metabolites emit `metabolites: !!omap []` (the flow-style + %empty omap cobrapy uses) so the file remains a valid YAML 1.2 + %document. if sum(field(:,pos) ~= 0) > 0 + fprintf(fid,' %s: !!omap\n',name); model.mets = model.mets(field(:,pos) ~= 0); model.coeffs = field(field(:,pos) ~= 0,pos); %Sort metabolites: @@ -230,6 +260,8 @@ function writeField(model,fid,fieldName,type,pos,name,preserveQuotes) for i = 1:length(model.mets) writeField(model, fid, 'coeffs', 'num', i, [' - ' model.mets{i}], preserveQuotes) end + else + fprintf(fid,' %s: !!omap []\n',name); end elseif strcmp(fieldName,'rxnEnzMat') @@ -249,7 +281,17 @@ function writeField(model,fid,fieldName,type,pos,name,preserveQuotes) elseif sum(strcmp({'subSystems','newMetMiriams','newRxnMiriams','newGeneMiriams','newCompMiriams','eccodes'},fieldName)) > 0 %eccodes/rxnNotes: if 1 write in 1 line, if more create header and list if strcmp(fieldName,'subSystems') - list = field{pos}; %subSystems already comes in a cell array + % The reader collapses an all-singleton subSystems field to + % a char column; defend against that (and against length + % mismatches caused by partial subsystem coverage) so the + % writer doesn't crash on shorter-than-rxns subsystem lists. + if iscell(field) + if pos > numel(field); return; end + list = field{pos}; + else + if pos > size(field, 1); return; end + list = field(pos, :); + end if isempty(list) return end @@ -306,14 +348,16 @@ function writeField(model,fid,fieldName,type,pos,name,preserveQuotes) %All other fields: if strcmp(type,'txt') value = field{pos}; - if preserveQuotes && ~isempty(value) - value = ['"',value,'"']; + if ~isempty(value) + if preserveQuotes || needsYamlQuoting(value) + value = ['"', escapeForDoubleQuoted(value), '"']; + end end elseif strcmp(type,'num') if isnan(field(pos)) value = []; else - value = sprintf('%.15g',full(field(pos))); + value = formatNumber(full(field(pos))); end end if ~isempty(value) @@ -323,64 +367,185 @@ function writeField(model,fid,fieldName,type,pos,name,preserveQuotes) end end -function writeMetadata(model,fid) -% Writes model metadata to the yaml file. This information will eventually -% be extracted entirely from the model, but for now, many of the entries -% are hard-coded defaults for HumanGEM. +function writeMetadata(model, fid, preserveQuotes) +% Writes the metaData block. Honors preserveQuotes so the rest of the +% file (which defaults to no surrounding quotes for cobra parity) stays +% consistent. The `date` field is preserved when the model carries one +% (model.date), so round-trips don't churn on every write; if absent +% it's filled with the current date. -fprintf(fid, '- metaData:\n'); -if isfield(model,'id') - fprintf(fid, ' id: "%s"\n', model.id); -else - fprintf(fid, ' id: "blankID"\n'); +fprintf(fid, '- metaData: !!omap\n'); +emitMetaField(fid, 'id', valueOrDefault(model,'id','blankID'), preserveQuotes); +emitMetaField(fid, 'name', valueOrDefault(model,'name','blankName'),preserveQuotes); +if isfield(model,'version') + emitMetaField(fid, 'version', model.version, preserveQuotes); end -if isfield(model,'name') - fprintf(fid, ' name: "%s"\n',model.name); +if isfield(model,'date') && ~isempty(model.date) + dateValue = model.date; else - fprintf(fid, ' name: "blankName"\n'); + dateValue = datestr(now, 29); %#ok % 29 = yyyy-mm-dd end -if isfield(model,'version') - fprintf(fid, ' version: "%s"\n',model.version); -end -fprintf(fid, ' date: "%s"\n',datestr(now,29)); % 29=YYYY-MM-DD +emitMetaField(fid, 'date', dateValue, preserveQuotes); if isfield(model,'annotation') - if isfield(model.annotation,'defaultLB') - fprintf(fid, ' defaultLB: "%g"\n', model.annotation.defaultLB); - end - if isfield(model.annotation,'defaultUB') - fprintf(fid, ' defaultUB: "%g"\n', model.annotation.defaultUB); - end - if isfield(model.annotation,'givenName') - fprintf(fid, ' givenName: "%s"\n', model.annotation.givenName); - end - if isfield(model.annotation,'familyName') - fprintf(fid, ' familyName: "%s"\n', model.annotation.familyName); - end - if isfield(model.annotation,'authors') - fprintf(fid, ' authors: "%s"\n', model.annotation.authors); - end - if isfield(model.annotation,'email') - fprintf(fid, ' email: "%s"\n', model.annotation.email); - end - if isfield(model.annotation,'organization') - fprintf(fid, ' organization: "%s"\n',model.annotation.organization); - end - if isfield(model.annotation,'taxonomy') - fprintf(fid, ' taxonomy: "%s"\n', model.annotation.taxonomy); - end - if isfield(model.annotation,'note') - fprintf(fid, ' note: "%s"\n', model.annotation.note); - end - if isfield(model.annotation,'sourceUrl') - fprintf(fid, ' sourceUrl: "%s"\n', model.annotation.sourceUrl); + annoFields = {'defaultLB','defaultUB','givenName','familyName', ... + 'authors','email','organization','taxonomy','note','sourceUrl'}; + for k = 1:numel(annoFields) + f = annoFields{k}; + if isfield(model.annotation, f) + emitMetaField(fid, f, model.annotation.(f), preserveQuotes); + end end end -if isfield(model,'ec') - if model.ec.geckoLight - geckoLight = 'true'; +% gecko_light is emitted at the top level (see geckoLight emission near +% the GECKO ec-* sections) to match raven_python.io.yaml; keeping it +% out of metaData lets cobrapy/ruamel keep the section a pure +% provenance block. +end + +function v = valueOrDefault(model, field, defaultVal) +if isfield(model, field) && ~isempty(model.(field)) + v = model.(field); +else + v = defaultVal; +end +end + +function emitMetaField(fid, key, value, preserveQuotes) +% Emit one ` - key: value` line inside the metaData omap block. +if islogical(value) + if value; value = 'true'; else; value = 'false'; end +end +if isnumeric(value) + value = formatNumber(double(value)); +elseif ~ischar(value) && ~isstring(value) + value = char(value); +end +if preserveQuotes + fprintf(fid, ' - %s: "%s"\n', key, value); +else + fprintf(fid, ' - %s: %s\n', key, value); +end +end + +function tf = needsYamlQuoting(s) +% A defensive subset of when a plain (unquoted) YAML scalar would be +% misparsed: leading YAML indicator chars, or any character that +% triggers flow-style collection / mapping parsing. Matches what +% ruamel.yaml's "round-trip" emitter quotes automatically. +if isempty(s) + tf = false; + return; +end +% Leading-character cases that turn into a flow indicator / tag / etc. +first = s(1); +if any(first == '[]{},&*!|>%@`#') + tf = true; return; +end +if first == '-' || first == '?' || first == ':' + tf = true; return; +end +% In-string cases: ': ' (key/value confusion), ' #' (comment), any flow +% bracket, leading or trailing whitespace, or anything outside ASCII +% printable. +if contains(s, ': ') || contains(s, ' #') || any(ismember(s, '[]{},')) + tf = true; return; +end +if ~strcmp(strip(string(s)), string(s)) + tf = true; return; +end +% YAML reserves certain tokens (true, false, null, …) as bare scalars +% reading as booleans / nulls. Quote when the entire value matches. +if any(strcmpi(s, {'true','false','null','yes','no','on','off','~'})) + tf = true; return; +end +tf = false; +end + +function out = escapeForDoubleQuoted(s) +% Escape backslashes and double quotes for emission inside YAML's +% double-quoted style. Conservative — full YAML escapes (Unicode etc.) +% are out of scope for the strings model curators normally use. +out = strrep(s, '\', '\\'); +out = strrep(out, '"', '\"'); +end + +function s = formatNumber(x) +% Format a finite number the way cobrapy / Python's float repr would — +% so whole-number floats round-trip as "1000.0", not "1000". Matches the +% ruamel.yaml output used by raven_python.io.yaml.write_yaml_model. +if isinf(x) + if x > 0 + s = '.inf'; else - geckoLight = 'false'; + s = '-.inf'; + end + return; +end +if x == floor(x) && abs(x) < 1e16 + s = sprintf('%.1f', x); % e.g. 1000 -> "1000.0" +else + s = sprintf('%.15g', x); +end +end + +function writeAnnotation(model, fid, kind, pos, preserveQuotes) +% Emit the per-entry `annotation` block, fusing MIRIAM cross-references +% with non-MIRIAM cobrapy-style annotation keys (currently: SMILES for +% metabolites). cobrapy expects SMILES inside `annotation.smiles`, not +% as a top-level metabolite key; this helper keeps the YAML aligned. +switch kind + case 'met' + miriamsField = 'metMiriams'; + extraName = 'smiles'; + extraField = 'metSmiles'; + otherwise + error('writeAnnotation:unsupportedKind', 'Unsupported kind: %s', kind); +end + +hasMiriams = isfield(model, miriamsField) && ~isempty(model.(miriamsField){pos}); +hasExtra = isfield(model, extraField) && ~isempty(model.(extraField){pos}); + +if ~hasMiriams && ~hasExtra + return; +end + +fprintf(fid, ' - annotation: !!omap\n'); +if hasMiriams + % Re-use the writeField MIRIAM path but suppress the block header + % it would emit (we already wrote it above). Tap into the same + % extractMiriam intermediate via a flat fprintf loop. + miriamNames = model.newMetMiriamNames; + miriamValues = model.newMetMiriams; + for j = 1:size(miriamValues, 2) + v = miriamValues{pos, j}; + if isempty(v); continue; end + list = strsplit(strrep(v, ' ', ''), ';'); + list = strip(list); + if numel(list) == 1 + valueOut = quoteIfNeeded(list{1}, preserveQuotes); + fprintf(fid, ' - %s: %s\n', miriamNames{j}, valueOut); + else + fprintf(fid, ' - %s:\n', miriamNames{j}); + for k = 1:numel(list) + fprintf(fid, ' - %s\n', ... + quoteIfNeeded(list{k}, preserveQuotes)); + end + end end - fprintf(fid,' geckoLight: "%s"\n',geckoLight); +end +if hasExtra + extraVal = quoteIfNeeded(model.(extraField){pos}, preserveQuotes); + fprintf(fid, ' - %s: %s\n', extraName, extraVal); +end +end + +function out = quoteIfNeeded(value, preserveQuotes) +% Quote a YAML scalar exactly when the surrounding writer would have +% lost it as a flow sequence / boolean / null otherwise. Conservative +% wrapper around needsYamlQuoting that respects preserveQuotes=true. +if preserveQuotes || needsYamlQuoting(value) + out = ['"', escapeForDoubleQuoted(value), '"']; +else + out = value; end end diff --git a/version.txt b/version.txt new file mode 100644 index 00000000..c200bec4 --- /dev/null +++ b/version.txt @@ -0,0 +1 @@ +2.11.1 \ No newline at end of file