From 0ddabeecfc5091b8a9cc4619b22eedce9efa39e2 Mon Sep 17 00:00:00 2001
From: BenjaSanchez <bensan@chalmers.se>
Date: Fri, 18 May 2018 16:20:38 +0200
Subject: [PATCH 1/4] style: model file name

---
 models/createNewModels.m              |   4 ++--
 models/{yeastGEM.mat => yeast780.mat} | Bin
 2 files changed, 2 insertions(+), 2 deletions(-)
 rename models/{yeastGEM.mat => yeast780.mat} (100%)

diff --git a/models/createNewModels.m b/models/createNewModels.m
index 6201737..fd03795 100644
--- a/models/createNewModels.m
+++ b/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:
diff --git a/models/yeastGEM.mat b/models/yeast780.mat
similarity index 100%
rename from models/yeastGEM.mat
rename to models/yeast780.mat

From 2f0df9f0bacd99b5ba9b77abe106e744474e25c4 Mon Sep 17 00:00:00 2001
From: BenjaSanchez <bensan@chalmers.se>
Date: Sun, 20 May 2018 13:42:09 +0200
Subject: [PATCH 2/4] refactor: getMWfromFormula

included getStoichFromFormula
---
 models/getMWfromFormula.m     | 27 +++++++-------------------
 models/getStoichFromFormula.m | 36 +++++++++++++++++++++++++++++++++++
 2 files changed, 43 insertions(+), 20 deletions(-)
 create mode 100644 models/getStoichFromFormula.m

diff --git a/models/getMWfromFormula.m b/models/getMWfromFormula.m
index 2ab99f4..6033847 100644
--- a/models/getMWfromFormula.m
+++ b/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
diff --git a/models/getStoichFromFormula.m b/models/getStoichFromFormula.m
new file mode 100644
index 0000000..b5c8514
--- /dev/null
+++ b/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
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\ No newline at end of file

From 2a8dc319ad6e4224701be703f28d0a2531182a37 Mon Sep 17 00:00:00 2001
From: BenjaSanchez <bensan@chalmers.se>
Date: Sun, 20 May 2018 14:06:50 +0200
Subject: [PATCH 3/4] refactor: addSLIMErxns

---
 models/addSLIMErxn.m | 20 +++++++++++---------
 1 file changed, 11 insertions(+), 9 deletions(-)

diff --git a/models/addSLIMErxn.m b/models/addSLIMErxn.m
index 6288a65..4ea3c43 100644
--- a/models/addSLIMErxn.m
+++ b/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,11 +101,12 @@
 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:
 for i = 1:length(tailsRxn)

From 9368d42f8357402ba89b467b1a95b12b9fa0ac0f Mon Sep 17 00:00:00 2001
From: BenjaSanchez <bensan@chalmers.se>
Date: Sun, 20 May 2018 18:34:31 +0200
Subject: [PATCH 4/4] fix: compliance with yeast-GEM

molecular formulas assigned to all backbones based on mass balances. A couple of other changes in the model format
---
 models/SLIMEr.m      | 21 ++++++++++++++-------
 models/addSLIMErxn.m | 24 +++++++++++++++++++++---
 2 files changed, 35 insertions(+), 10 deletions(-)

diff --git a/models/SLIMEr.m b/models/SLIMEr.m
index a5ce0e7..666d760 100644
--- a/models/SLIMEr.m
+++ b/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
 
diff --git a/models/addSLIMErxn.m b/models/addSLIMErxn.m
index 4ea3c43..5f3973b 100644
--- a/models/addSLIMErxn.m
+++ b/models/addSLIMErxn.m
@@ -109,12 +109,30 @@
 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, ...