Merge pull request #179 from SysBioChalmers/removeBoundaryComp

feat: remove the boundary compartment

code/addBoundaryMets.m

@@ -4,17 +4,18 @@
 %   Boundary metabolites are pseudometabolites that exist at the system
 %   boundary, and are used for the import/export of mass into/out of the
 %   system. 
-%   For example: "glc[boundary] <==> glc[extracellular]"
+%   For example: "glc[extracellular] <==> glc[boundary]"
 %
-%   Some model formats (e.g., Cobra) do not use boundary metabolites, but
-%   instead formulate these exchange reactions (or "demand/DM" or "sink"
-%   reactions) without an explicit reactant. 
+%   Some models do not use boundary metabolites, but instead formulate
+%   these exchange reactions (or "demand/DM" or "sink" reactions) without
+%   an explicit product (or reactant).
 %   For example: "glc[extracellular] <==> "
 %
-%   This function identifies all reactions that are formulated with an 
-%   exchange of an extracellular metabolite into nothing as a reaction with
-%   a boundary metabolite, and adds any new boundary metabolites to the 
-%   model, if they did not yet exist.
+%   This function identifies all reactions that involve the conversion of 
+%   a single metabolite into nothing (or vice versa), and balances the
+%   reaction with the same metabolite in a different (boundary) comparment.
+%   The function will also add any new boundary metabolites to the model if
+%   they did not yet exist.
 %
 % USAGE:
 %
@@ -25,11 +26,12 @@
 %
 %   model       Model structure.
 %
-%   exch_only   (Optional, default = TRUE) 
+%   exch_only   (Optional, default = FALSE) 
 %               If TRUE, only exchange rxns (i.e., into/out of the
 %               extracellular compartment) will be considered.
 %               If FALSE, reactions involving mets in other compartments
 %               will also be considered. For example:
+%
 %               "met[nucleus] -->" becomes "met[nucleus] --> met[boundary]"
 %
 %
@@ -38,10 +40,9 @@
 %   new_model   New model structure, with added boundary metabolites.
 %
 
-
 % handle input arguments
 if nargin < 2
-    exch_only = true;
+    exch_only = false;
 end
 
 % add a boundary compartment to the model if it does not yet exist

code/io/importHumanYaml.m

@@ -52,7 +52,7 @@
 model.metComps={};
 model.inchis={};
 model.metFormulas={};
-model.unconstrained=[];
+%model.unconstrained=[]; %abandoned
 model.rxnReferences={};
 model.rxnFrom={};
 model.metFrom={};
@@ -198,8 +198,8 @@
 
             case 'metabolites'
                 readEquation = true;
-                leftEquation  = '';
-                rightEquation = '';
+                leftEquation  = {''};
+                rightEquation = {''};
 
             otherwise
                 if readSubsystems
@@ -277,11 +277,11 @@
 [~, metIdx] = ismember(model.mets, newMets);
 model.S = S(metIdx, :);
 
-% Although this works with HumanGEM, but it is NOT a generic solution of
-% dealing with the `unconstrained` field for other models!
-model.unconstrained = double(endsWith(model.mets, 'x'));
-        if ~silentMode
-            fprintf(' Done!\n');
+% Now the `unconstrained` field is abandoned in HumanGEM
+%model.unconstrained = double(endsWith(model.mets, 'x'));
+
+if ~silentMode
+    fprintf(' Done!\n');
 end
 end
 

code/modelCuration/removeBoundaryCompartment.m

@@ -0,0 +1,52 @@
+%
+% FILE NAME:    removeBoundaryCompartment.m
+%
+% PURPOSE: This script removes the boundary compartment [x] and all
+%          metabolites within this compartment from the model, as discussed
+%          in #172
+
+% load model
+ihuman = importHumanYaml('../../modelFiles/yml/HumanGEM.yml');
+
+% delete unconstrained (boundary) metabolites
+model_del = simplifyModel(ihuman);
+del_mets = setdiff(ihuman.mets, model_del.mets);
+fprintf('\nRemoved %u boundary metabolites from the model.\n\n', numel(del_mets));
+
+% remove the boundary compartment itself and update the metComps field
+metCompSymbols = model_del.comps(model_del.metComps);
+[~,comp_ind] = ismember('boundary', lower(model_del.compNames));
+model_del.comps(comp_ind) = [];
+model_del.compNames(comp_ind) = [];
+[~,model_del.metComps] = ismember(metCompSymbols, model_del.comps);
+
+
+% write new model to yml
+writeHumanYaml(model_del, '../../modelFiles/yml/HumanGEM.yml');
+
+
+% import metabolite annotations
+metAssoc = jsondecode(fileread('../../data/annotation/humanGEMMetAssoc.JSON'));
+
+% remove boundary metabolites
+del_ind = ismember(metAssoc.mets, del_mets);
+f = fieldnames(metAssoc);
+for i = 1:numel(f)
+    metAssoc.(f{i})(del_ind) = [];
+end
+
+% verify that annotation structure is aligned with model
+if ~isequal(model_del.mets, metAssoc.mets)
+    error('Model and metabolite annotation structures not synced!');
+end
+
+% export metabolite annotations
+jsonStr = jsonencode(metAssoc);
+fid = fopen('../../data/annotation/humanGEMMetAssoc.JSON', 'w');
+fwrite(fid, prettyJson(jsonStr));
+fclose(fid);
+
+
+
+
+