Merge pull request #71 from SysBioChalmers/devel

Devel

ComplementaryScripts/dependencies.txt

@@ -1,4 +1,5 @@
-SBML_toolbox	v4.1.0
+RAVEN_toolbox	v2.0.0-rc.1
+COBRA_toolbox	commit b80de23
 SBML_level	3
 SBML_version	1
 fbc_version	2

ComplementaryScripts/increaseVersion.m

@@ -8,7 +8,7 @@ function increaseVersion(version)
 fid = fopen('../history.md','r');
 history = fscanf(fid,'%s');
 fclose(fid);
-if isempty(strfind(history,['yeast' version ':']))
+if contains(history,['yeast' version ':'])
     error('ERROR: update history.md first')
 end
 
@@ -21,7 +21,6 @@ function increaseVersion(version)
 saveYeastModel(model)
 
 %Store model as .mat (only for releases):
-mkdir('../ModelFiles/mat');
 save('../ModelFiles/mat/yeastGEM.mat','model');
 
 %Update version file:

ComplementaryScripts/modelCuration/checkMetBalance.m

@@ -0,0 +1,84 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% checkMetBalance(model,metName,flux,show_all)
+% Shows for a given metabolite all involved reactions on all compartments.
+%
+% INPUT:    model       A GEM as a structure object
+%           metName     The name of the desired metabolite (e.g. 'NADH')
+%           flux        (OPTIONAL) A solution of the model (nx1 vector)
+%           show_all    (OPTIONAL) If false, only non-zero fluxes will be
+%                       displayed (default = false; if no flux is given as
+%                       input then default = true)
+%
+% OUTPUT:   On the command window displays all involved reactions,
+%           including the flux value, the reaction ID and the formula. 
+%           Results are displayed by compartment, and separated depending
+%           if the metabolite is being consumed, produced or transported 
+%           in each rxn. If no flux was given as input, a zero will appear
+%           instead.
+%
+% Benjamín J. Sánchez
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function checkMetBalance(model,metName,flux,show_all)
+
+lb = model.lb;
+ub = model.ub;
+
+if nargin < 4
+    show_all = false;
+    if nargin < 3
+        flux     = zeros(size(model.rxns));
+        show_all = true;
+    end
+end
+
+pos = zeros(size(model.mets));
+for i = 1:length(model.metNames)
+    if strcmp(model.metNames{i},metName)
+        pos(i) = 1;
+    end
+end
+pos = find(pos);
+
+for i = 1:length(pos)
+    rxns_cons = {};
+    rxns_prod = {};
+    for j = 1:length(model.rxns)
+        if model.S(pos(i),j) ~= 0 && (show_all || abs(flux(j)) > 1e-10)
+            rxn         = model.rxns{j};
+            rxn_formula = printRxnFormula(model,model.rxns(j),false,false,true);
+            rxn_formula = rxn_formula{1};
+            if model.S(pos(i),j) < 0
+                rxns_cons = [rxns_cons;{lb(j) flux(j) ub(j) rxn rxn_formula}];
+            else
+                rxns_prod = [rxns_prod;{lb(j) flux(j) ub(j) rxn rxn_formula}];
+            end
+        end
+    end
+    if ~isempty(rxns_cons) || ~isempty(rxns_prod)
+        print_group(rxns_cons,'Consumption')
+        print_group(rxns_prod,'Production')
+    end
+end
+
+fprintf('\n')
+
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function print_group(rxns,name)
+
+if ~isempty(rxns)
+    fprintf([name ':\n'])
+    [~,order] = sort(cell2mat(rxns(:,1)),'descend');
+    rxns      = rxns(order,:);
+    for i = 1:length(order)
+        fprintf([sprintf('%6.2e',rxns{i,1}) '\t' sprintf('%6.2e',rxns{i,2}) ...
+            '\t' sprintf('%6.2e',rxns{i,3}) '\t' rxns{i,4} '\t\t' rxns{i,5} '\n'])
+    end
+end
+
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ComplementaryScripts/modelCuration/makeFormulasCompliant.m

@@ -0,0 +1,71 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% makeFormulasCompliant.m
+% Corrects SBML-incompatible metabolite formulas.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%aa codes:
+aas = {'s_0955[c]'	'ala'   	% A     Alanine
+       's_0981[c]'	'cys'       % C     Cysteine
+       's_0973[c]'	'asp'       % D     Aspartic acid
+       's_0991[c]'	'glu'       % E     Glutamic acid
+       's_1032[c]'	'phe'   	% F     Phenylalanine
+       's_1003[c]'	'gly'       % G     Glycine
+       's_1006[c]'	'his'   	% H     Histidine
+       's_1016[c]'	'ile'       % I     Isoleucine
+       's_1025[c]'	'lys'       % K     Lysine
+       's_1021[c]'	'leu'   	% L     Leucine
+       's_1029[c]'	'met'       % M     Methionine
+       's_0969[c]'	'asn'   	% N     Asparagine
+       's_1035[c]'	'pro'       % P     Proline
+       's_0999[c]'	'gln'   	% Q     Glutamine
+       's_0965[c]'	'arg'       % R     Arginine
+       's_1039[c]'	'ser'   	% S     Serine
+       's_1045[c]'	'thr'       % T     Threonine
+       's_1056[c]'	'val'       % V     Valine
+       's_1048[c]'	'trp'       % W     Tryptophan
+       's_1051[c]'	'tyr'};     % Y     Tyrosine
+
+%tRNA's are just produced and consumed in single loops, so we can just
+%replace the radical part with an "R" and still maintain mass balances,
+%i.e. aa-tRNA(aa) will have "R" + the formula of the corresponding aa and
+%tRNA(aa) will just have an "R". Example:
+%     tRNA(Ala): C10H17 O10PR2(C5H8O6PR)n -> R
+% Ala-tRNA(Ala): C13H22NO11PR2(C5H8O6PR)n -> C3H6NOR (ala-R without an -OH)
+% Cycle in which the 2 are involved:
+% r_4047: 0.4193 Ala-tRNA(Ala) + ...  -> 0.4193 tRNA(Ala) + ... + protein (growth)
+% r_0157: ATP           + L-alanine + tRNA(Ala) -> Ala-tRNA(Ala) + AMP         + diphosphate
+%         C10H12N5O13P3   C3H7NO2     R      	   C3H6NOR		   C10H12N5O7P	 HO7P2
+
+for i = 1:length(model.mets)
+    name    = model.metNames{i};
+    formula = model.metFormulas{i};
+    if contains(model.metNames{i},'-tRNA(')
+        %Correct metabolite:
+        try
+            aaName = lower(name(1:strfind(name,'-')-1));
+            aaID   = aas{strcmp(aas(:,2),aaName),1};
+            aaPos  = strcmp(model.mets,aaID);
+            model.metFormulas{i} = [model.metFormulas{aaPos} 'R'];
+            model.metFormulas{i} = takeOutFromFormula(model.metFormulas{i},'OH');
+        catch
+            %formyl-met: C6H11NO3S
+            model.metFormulas{i} = 'C6H10NO2SR';
+        end
+
+        %Correct associated tRNA pair:
+        pairName = name(strfind(name,'-')+1:end);
+        pairPos  = strcmp(model.metNames,pairName);
+        if sum(pairPos) > 0
+            model.metFormulas{pairPos} = 'R';
+        end
+    end
+end
+
+%Display the remaining problems:
+for i = 1:length(model.mets)
+    if contains(model.metFormulas{i},')n')
+        disp([model.metNames{i} ': ' model.metFormulas{i}]);
+    end
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ComplementaryScripts/modelCuration/takeOutFromFormula.m

@@ -0,0 +1,50 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% formula = takeOutFromFormula(formula,elements)
+% Takes away from formula each of the elements specified.
+%
+% Benjamín J. Sánchez
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function formula = takeOutFromFormula(formula,elements)
+
+for i = 1:length(elements)
+    posInFormula = strfind(formula,elements(i));
+    if posInFormula == length(formula)
+        formula = strrep(formula,elements(i),'');               %only 1 element at the end of the formula
+    else
+        for j = posInFormula+1:length(formula)
+            if isletter(formula(j))
+                if j == posInFormula+1
+                    formula = strrep(formula,elements(i),'');   %only 1 element in the middle of the formula
+                else
+                    formula = goDown(formula,posInFormula,j);	%more than 1 element in the formula
+                end
+                break
+            elseif j == length(formula)
+                formula = goDown(formula,posInFormula,j+1);     %more than 1 element at the end of the formula
+            end
+        end
+    end
+end
+
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function formula = goDown(formula,startPos,endPos)
+
+numbElem = str2double(formula(startPos+1:endPos-1));
+if numbElem == 2
+    numbElem = '';
+else
+    numbElem = num2str(numbElem - 1);
+end
+if endPos == length(formula)+1
+    formula = [formula(1:startPos) numbElem];
+else
+   formula = [formula(1:startPos) numbElem formula(endPos:end)];
+end
+
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ComplementaryScripts/anaerobicModel.m → ComplementaryScripts/otherChanges/anaerobicModel.m

@@ -31,6 +31,8 @@
 GAM   = 30.49;  %Data from Nissen et al. 1997
 P     = 0.461;  %Data from Nissen et al. 1997
 NGAM  = 0;      %Refit done in Jouthen et al. 2012
+cd ../modelCuration
 model = changeBiomass(model,P,GAM,NGAM);
+cd ../otherChanges
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ComplementaryScripts/saveYeastModel.m

@@ -22,24 +22,22 @@ function saveYeastModel(model)
 end
 fclose(fid);
 
-%Retrieve SBML toolbox version:
-SBMLpath = which('SBMLToolbox.m');
-slashPos = strfind(SBMLpath,'\');
-if isempty(slashPos)
-    slashPos = strfind(SBMLpath,'/');
-end
-try
-    SBMLpath = SBMLpath(1:slashPos(end-1));
-    fid      = fopen([SBMLpath 'VERSION.txt'],'r');
-    SBMLTver = fscanf(fid,'%s');
-    fclose(fid);
-catch
-    SBMLTver = '?';
-end
+%Retrieve RAVEN version:
+RAVENver = getVersion('checkInstallation.m','version.txt');
+
+%Retrieve latest COBRA commit:
+COBRApath   = which('initCobraToolbox.m');
+slashPos    = getSlashPos(COBRApath);
+COBRApath   = COBRApath(1:slashPos(end)-1);
+currentPath = pwd;
+cd(COBRApath)
+COBRAcommit = git('log -n 1 --format=%H');
+cd(currentPath)
 
 %Save file with versions:
 fid = fopen('dependencies.txt','wt');
-fprintf(fid,['SBML_toolbox\tv' SBMLTver '\n']);
+fprintf(fid,['RAVEN_toolbox\tv' RAVENver '\n']);
+fprintf(fid,['COBRA_toolbox\tcommit ' COBRAcommit(1:7) '\n']);
 fields = fieldnames(model.modelVersion);
 for i = 1:length(fields)
     value = model.modelVersion.(fields{i});
@@ -49,4 +47,32 @@ function saveYeastModel(model)
 
 end
 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function version = getVersion(IDfileName,VERfileName)
+
+try
+    path     = which(IDfileName);
+    slashPos = getSlashPos(path);
+    path     = path(1:slashPos(end-1));
+    fid      = fopen([path VERfileName],'r');
+    version  = fscanf(fid,'%s');
+    fclose(fid);
+catch
+    version = '?';
+end
+
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function slashPos = getSlashPos(path)
+
+slashPos = strfind(path,'\');       %Windows
+if isempty(slashPos)
+    slashPos = strfind(path,'/');   %MAC/Linux
+end
+
+end
+
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%