Merge pull request #8 from SysBioChalmers/fix/yeastGEMcompliance

Fix/yeast GEM compliance

models/SLIMEr.m

@@ -1,7 +1,7 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % model = SLIMEr(model,data,includeTails)
 %
-% Benjamín J. Sánchez. Last update: 2018-03-30
+% Benjamín J. Sánchez. Last update: 2018-05-20
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 function model = SLIMEr(model,data,includeTails)
@@ -12,13 +12,19 @@
 for i = 1:length(metIDs)
     backName = getBackboneName(metNames{i});
     if ~isempty(backName)
-        if ~ismember(backName,model.metNames)
+        if ismember(backName,model.metNames)
+            if ~startsWith(backName,'ergosterol [')
+                model.metFormulas{strcmp(model.metNames,backName)} = '';
+            end
+        else
             model = addLipidSpecies(model,backName,'',false);
         end
         %Add transport rxn to cytoplasm for non cytoplasmic backbones:
         if ~contains(backName,'cytoplasm')
             cytoName = [backName(1:strfind(backName,'[')) 'cytoplasm]'];
-            if ~ismember(cytoName,model.metNames)
+            if ismember(cytoName,model.metNames)
+                model.metFormulas{strcmp(model.metNames,cytoName)} = '';
+            else
                 model = addLipidSpecies(model,cytoName,'',false);
             end
             backID    = model.mets(strcmp(model.metNames,backName));
@@ -53,12 +59,12 @@
 end
 
 %Create lipid pseudo-rxn for backbones:
-model = addLipidSpecies(model,'lipid - backbones [cytoplasm]','NA',includeTails);
+model = addLipidSpecies(model,'lipid - backbones [cytoplasm]','',includeTails);
 model = changeLipidComp(model,data.lipidData);
 
 %Create lipid pseudo-rxn for tails:
 if includeTails
-    model = addLipidSpecies(model,'lipid - tails [cytoplasm]','NA',includeTails);
+    model = addLipidSpecies(model,'lipid - tails [cytoplasm]','',includeTails);
     model = changeChainComp(model,data.chainData);
 end
 
@@ -121,8 +127,9 @@
 model.S(Hpos,bioRxn)   = +GAM;
 model.S(Ppos,bioRxn)   = +GAM;
 
-%Remove wrongly created field:
-model = rmfield(model,'grRules');
+%Correct some wrongly created fields:
+model   = rmfield(model,'grRules');
+model.S = full(model.S);
 
 end
 

models/addSLIMErxn.m

@@ -1,7 +1,7 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % [model,toDelete] = addSLIMErxn(model,rxnID,specID)
 %
-% Benjamín J. Sánchez. Last update: 2018-03-24
+% Benjamín J. Sánchez. Last update: 2018-05-20
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 function [model,toDelete] = addSLIMErxn(model,rxnID,specID)
@@ -39,7 +39,8 @@
 specName = specName(1:strfind(specName,'[')-2);
 
 %Stoich. coeff. for backbone: molecular weight of specific species
-specMW = getMWfromFormula(model.metFormulas(specPos));
+specFormula = model.metFormulas(specPos);
+specMW      = getMWfromFormula(specFormula);
 
 %Define number of tails to add:
 switch backName
@@ -67,9 +68,9 @@
           'inositol phosphomannosylinositol phosphoceramide'    %MIPC
           'mannosylinositol phosphorylceramide'}                %M(IP)2C
         if contains(specName,'(C24)')
-            tailsRxn = {'C18:0';'C24:0'};
+            tailsRxn = {'C18:0','C24:0'};
         elseif contains(specName,'(C26)')
-            tailsRxn = {'C18:0';'C26:0'};
+            tailsRxn = {'C18:0','C26:0'};
         end
 
     %Cases with specific names:
@@ -100,19 +101,38 @@
 end
 
 %Find tail metabolites in model:
-tailPos    = find(~cellfun(@isempty,strfind(model.metNames,' chain [cytoplasm]')));
-tailIDs    = model.mets(tailPos)';
-tailsModel = model.metNames(tailPos)';
-tailsMWs   = getMWfromFormula(model.metFormulas(tailPos)');
-tailCoeffs = zeros(size(tailIDs));
+tailPos       = contains(model.metNames,' chain [cytoplasm]');
+tailIDs       = model.mets(tailPos)';
+tailsModel    = model.metNames(tailPos)';
+tailsFormulas = model.metFormulas(tailPos)';
+tailsMWs      = getMWfromFormula(tailsFormulas);
+tailCoeffs    = zeros(size(tailIDs));
 
 %Match to corresponding tail:
+prodFormulas = cell(size(tailsRxn));
 for i = 1:length(tailsRxn)
-    tailName   = [tailsRxn{i} ' chain [cytoplasm]'];
-    tailMatch  = strcmp(tailsModel,tailName);
-    tailCoeffs = tailCoeffs + tailMatch.*tailsMWs;
+    tailName        = [tailsRxn{i} ' chain [cytoplasm]'];
+    tailMatch       = strcmp(tailsModel,tailName);
+    tailCoeffs      = tailCoeffs + tailMatch.*tailsMWs;
+    prodFormulas{i} = tailsFormulas{tailMatch};
 end
 
+%Assign molecular formula to backbone by balancing out the SLIME rxns:
+backName     = [backName ' ['];
+formula      = '';
+elements     = {'C','H','N','O','P','S'};
+for j = 1:length(elements)
+    Nin  = getStoichFromFormula(specFormula, elements{j});
+    Nout = getStoichFromFormula(prodFormulas, elements{j});
+    diff = Nin - sum(Nout);
+    if diff == 1
+        formula = [formula elements{j}];
+    elseif diff > 1
+        formula = [formula elements{j} num2str(diff)];
+    end
+end
+model.metFormulas(startsWith(model.metNames,backName)) = {formula};
+
 %Create SLIME rxn (with same ID as previous ISA rxn but different name):
 model = addReaction(model,rxnID, ...
                     'reactionName', rxnName, ...

models/createNewModels.m

@@ -1,7 +1,7 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % createNewModels
 %
-% Benjamín J. Sánchez. Last update: 2018-03-30
+% Benjamín J. Sánchez. Last update: 2018-05-18
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 clear variables
@@ -12,7 +12,7 @@
 initCobraToolbox
 
 %Original model:
-model = load('yeastGEM.mat');
+model = load('yeast780.mat');
 model = model.model;
 
 %Create 2 models for each of the 10 conditions in the stress dataset:

models/getMWfromFormula.m

@@ -2,34 +2,21 @@
 % MWs = getMWfromFormula(metFormulas)
 % Returns the MW (in g/mmol) of a chemical formula
 %
-% Benjamín J. Sánchez. Last update: 2018-03-26
+% Benjamín J. Sánchez. Last update: 2018-05-20
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 function MWs = getMWfromFormula(metFormulas)
 
 %Atomic weights:
-AWs = {'H' 1.00794
-       'C' 12.011
-       'N' 14.00674
-       'O' 15.9994
-       'P' 30.973762
-       'S' 32.066};
-
+elements = {'C',   'H',      'O',     'N',      'P',       'S'};
+AWs      = [12.011, 1.00794, 15.9994, 14.00674, 30.973762, 32.066];
+
 %Calculate the MW for each formula:
 MWs = zeros(size(metFormulas));
-for i = 1:length(MWs)
-    formula = metFormulas{i};
-    elePos  = [regexp(formula,'[A-Z]') length(formula)+1];
+for i = 1:length(elements)
     %Find stochiometry of each element and multiply by atomic weight:
-    for j = 1:length(elePos)-1
-        element = formula(elePos(j));
-        AW      = AWs{strcmp(AWs(:,1),element),2}/1000; %[g/mmol]
-        stoich  = str2double(formula(elePos(j)+1:elePos(j+1)-1));
-        if isnan(stoich)
-            stoich = 1;
-        end
-        MWs(i) = MWs(i) + stoich*AW;
-    end
+    stoich = getStoichFromFormula(metFormulas,elements{i});
+    MWs    = MWs + stoich*AWs(i)/1000;
 end
 
 end

models/getStoichFromFormula.m

@@ -0,0 +1,36 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% stoich = getStoichFromFormula(metFormulas,element)
+% Returns the stoichiometry of a given element in a set of chemical formulas
+%
+% Benjamín J. Sánchez. Last update: 2018-05-20
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function stoich = getStoichFromFormula(metFormulas,element)
+
+%Calculate the MW for each formula:
+stoich = zeros(size(metFormulas));
+for i = 1:length(stoich)
+    formula = metFormulas{i};
+    elePos  = strfind(formula,element);
+    if length(elePos) > 1
+        error(['Non-standard formula: ' formula])
+    elseif isempty(elePos)
+        stoich(i) = 0;
+    elseif elePos == length(formula)
+        stoich(i) = 1;
+    else
+        rest    = formula(elePos+1:end);
+        nextPos = regexp(rest,'[A-Z]');
+        if isempty(nextPos)
+            stoich(i) = str2double(rest);
+        elseif nextPos(1) == 1
+            stoich(i) = 1;
+        else
+            stoich(i) = str2double(rest(1:nextPos(1)-1));
+        end
+    end
+end
+
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%