Merge pull request #31 from SysBioChalmers/devel

GECKO 1.2.1

.gitignore

@@ -53,3 +53,4 @@ docs/_*
 #########
 geckopy/geckopy/data_files/*.xml
 *~
+Matlab_Module/*.mat

Matlab_Module/Kcat_sensitivity_analysis/modifyKcats.m

@@ -78,6 +78,7 @@
                          'error','gRControlCoeff'};
 
         modifications = cell2table(modifications,'VariableNames',varNamesTable);
+        modifications = truncateValues(modifications,4);
         writetable(modifications,['../../Models/' name '/data/' name '_kcatModifications.txt']);
 
     else
@@ -90,11 +91,13 @@
         if ~isempty(limRxns)
             varNamesTable = {'rxnNames','rxnPos','gRControlCoeff'};
             modifications = cell2table(limRxns,'VariableNames',varNamesTable);
+            modifications = truncateValues(modifications,4);
             writetable(modifications,['../../Models/' name '/data/' name '_limitingRxns.txt']);
         end
         if ~isempty(limEnz)
             varNamesTable = {'EnzNames','EnzPos','gRControlCoeff'};
             modifications = cell2table(limEnz,'VariableNames',varNamesTable);
+            modifications = truncateValues(modifications,4);
             writetable(modifications,['../../Models/' name '/data/' name '_limitingEnzymes.txt']);
         end
     end

Matlab_Module/Kcat_sensitivity_analysis/topUsedEnzymes.m

@@ -41,6 +41,7 @@
     end
     %Write the top-ten used enzyme names and their percentage usages for
     %every condition on the output file
+    outputFile = truncateValues(outputFile,1);
     writetable(outputFile,['../../Models/' name '/data/' name '_topUsedEnzymes.txt'])
 end
 

Matlab_Module/Kcat_sensitivity_analysis/truncateValues.m

@@ -0,0 +1,19 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% table = truncateValues(table,nCols)
+%
+%
+% Benjamin Sanchez    Last edited: 2018-05-25
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function table = truncateValues(table,nCols)
+
+[m,n] = size(table);
+for i = 1:m
+    for j = (n-nCols+1):n
+        orderMagn  = max([ceil(log10(abs(table{i,j}))),0]);
+        table{i,j} = round(table{i,j},6-orderMagn);
+    end
+end
+
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Matlab_Module/change_model/addArmReaction.m

@@ -9,7 +9,7 @@
 % OUTPUTS:
 % model             Modified GEM structure (1x1 struct)
 % 
-% Cheng Zhang. Last edited: 2016-12-21
+% Cheng Zhang & Ivan Domenzain. Last edited: 2018-05-28
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 function model = addArmReaction(model,rxn)
@@ -25,16 +25,32 @@
 obj      = model.c(rxnIndex);
 coeffsS  = model.S(sub_pos,rxnIndex)';
 coeffsP  = model.S(pro_pos,rxnIndex)';
+grRule   = model.grRules{rxnIndex};
+
+subSystem = '';
+if isfield(model,'subSystems')
+    if ~isempty(model.subSystems{rxnIndex}{1})
+        subSystem = model.subSystems{rxnIndex};
+    end
+end
 
 %Create new rxn:
-mets   = [metS,['pmet_' rxn]];
-coeffs = [coeffsS,1];
-name   = {['arm_' rxn],[model.rxnNames{rxnIndex} ' (arm)']};
-model  = addReaction(model,name,mets,coeffs,true,LB,UB,obj);
+rxnID = ['arm_' rxn];
+model = addReaction(model,rxnID, ...
+                    'reactionName', [model.rxnNames{rxnIndex} ' (arm)'], ...
+                    'metaboliteList', [metS,['pmet_' rxn]], ...
+                    'stoichCoeffList', [coeffsS,1], ...
+                    'lowerBound', LB, ...
+                    'upperBound', UB, ...
+                    'objectiveCoef', obj, ...
+                    'subSystem', subSystem);
+model.grRules{strcmp(model.rxns,rxnID)} = grRule;
 
 %Change old rxn:
-name   = {rxn,model.rxnNames{rxnIndex}};
-model  = addReaction(model,name,[['pmet_' rxn], metP],[-1,coeffsP]);
+model = addReaction(model,rxn, ...
+                    'reactionName', model.rxnNames{rxnIndex}, ...
+                    'metaboliteList', [['pmet_' rxn], metP], ...
+                    'stoichCoeffList', [-1,coeffsP]);
 
 %Update metComps:
 pos = strcmp(model.mets,['pmet_' rxn]);

Matlab_Module/change_model/addEnzymesToRxn.m

@@ -1,5 +1,5 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% model = addEnzymesToRxn(model,kvalues,rxn,newMets,newRxnName)
+% model = addEnzymesToRxn(model,kvalues,rxn,newMets,newRxnName,kegg,swissprot)
 % Adds new metabolite to the left side of a selected reaction in the model.
 % If the reaction does not exist it will create a new one.
 %
@@ -13,29 +13,52 @@
 % OUTPUTS:
 % model             Modified GEM structure (1x1 struct)
 % 
-% Cheng Zhang. Last edited: 2016-03-22
+% Cheng Zhang & Ivan Domenzain. Last edited: 2018-05-28
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-function model = addEnzymesToRxn(model,kvalues,rxn,newMets,newRxnName)
+function model = addEnzymesToRxn(model,kvalues,rxn,newMets,newRxnName,kegg,swissprot)
 
 %Define all neccesary parts for new (or changed) rxn:
-rxnIndex = find(ismember(model.rxns,rxn)); 
-metS     = model.mets(model.S(:,rxnIndex) < 0)';
-metP     = model.mets(model.S(:,rxnIndex) > 0)';
-LB       = model.lb(rxnIndex);
-UB       = model.ub(rxnIndex);    
-obj      = model.c(rxnIndex);
-coeffsS  = model.S(model.S(:,rxnIndex)<0,rxnIndex)';
-coeffsP  = model.S(model.S(:,rxnIndex)>0,rxnIndex)';
+rxnIndex  = strcmp(model.rxns,rxn); 
+metS      = model.mets(model.S(:,rxnIndex) < 0)';
+metP      = model.mets(model.S(:,rxnIndex) > 0)';
+LB        = model.lb(rxnIndex);
+UB        = model.ub(rxnIndex);    
+obj       = model.c(rxnIndex);
+coeffsS   = model.S(model.S(:,rxnIndex)<0,rxnIndex)';
+coeffsP   = model.S(model.S(:,rxnIndex)>0,rxnIndex)';
+genes     = cell(size(newMets));
 
-%Include enzyme in reaction:
+subSystem = '';
+if isfield(model,'subSystems')
+    if ~isempty(model.subSystems{rxnIndex}{1})
+        subSystem = model.subSystems{rxnIndex};
+    end
+end
+
+%Find genes either in swissprot or in kegg and with them construct the gene rule:
 for i = 1:length(newMets)
-    metS    = [metS,newMets{i}];
-    coeffsS = [coeffsS,-1/kvalues(i)];
+    protein = strrep(newMets{i},'prot_','');
+    try
+        geneIDs      = swissprot{strcmp(swissprot(:,1),protein),3};
+        geneIDs      = strsplit(geneIDs,' ');
+        [genes(i),~] = intersect(geneIDs,model.genes);
+    catch
+        genes{i} = kegg{strcmp(kegg(:,1),protein),3};
+    end
 end
-mets   = [metS,metP];
-coeffs = [coeffsS,coeffsP];
-model  = addReaction(model,newRxnName,mets,coeffs,true,LB,UB,obj,'','','','',false);
+grRule = strjoin(genes,' and ');
+
+%Include enzyme in reaction:
+model = addReaction(model,newRxnName{1}, ...
+                    'reactionName', newRxnName{2}, ...
+                    'metaboliteList', [metS,newMets,metP], ...
+                    'stoichCoeffList', [coeffsS,-kvalues.^-1,coeffsP], ...
+                    'lowerBound', LB, ...
+                    'upperBound', UB, ...
+                    'objectiveCoef', obj, ...
+                    'subSystem', subSystem);
+model.grRules{strcmp(model.rxns,newRxnName{1})} = grRule;
 
 end
 

Matlab_Module/change_model/addProtein.m

@@ -12,7 +12,7 @@
 % OUTPUTS:
 % model             Model with the added protein
 % 
-% Cheng Zhang & Benjamín Sánchez. Last edited: 2018-03-15
+% Cheng Zhang & Benjamin J. Sanchez. Last edited: 2018-05-28
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 function model = addProtein(model,P,kegg,swissprot)
@@ -23,19 +23,23 @@
 pos_e         = strcmp(model.enzymes,P);        %position in model.enzymes
 
 %Update model.MWs & model.sequences vectors:
-match_geneName = false;
+match_gen      = false;
+match_enzName = false;
 match_MW       = false;
 match_seq      = false;
+gene           = [];
 for i = 1:length(swissprot)
     %Gene name:
-    if strcmp(P,swissprot{i,1}) && ~isempty(swissprot{i,3}) && ~match_geneName
-        match_geneName  = true;
-        geneName        = swissprot{i,3};
-        pos_space       = strfind(geneName,' ');
-        if isempty(pos_space)
-            pos_space = length(geneName)+1;
+    if strcmp(P,swissprot{i,1}) && ~isempty(swissprot{i,3}) && ~match_enzName
+        match_enzName     = true;
+        geneNames         = swissprot{i,3};
+        geneIDs           = strsplit(geneNames,' ');
+        [geneSwissprot,~] = intersect(geneIDs,model.genes);
+        if ~isempty(geneSwissprot)
+            match_gen = true;
+            gene      = geneSwissprot{1};
         end
-        model.geneNames{pos_e,1} = geneName(1:pos_space(1)-1);
+        model.enzNames{pos_e,1} = geneIDs{1};
     end
     %Molecular Weight:
     if strcmp(P,swissprot{i,1}) && swissprot{i,5} ~= 0 && ~match_MW
@@ -48,8 +52,8 @@
         model.sequences{pos_e,1} = swissprot{i,6};
     end
 end
-if ~match_geneName
-    model.geneNames{pos_e,1} = '-';
+if ~match_enzName
+    model.enzNames{pos_e,1} = '-';
 end
 if ~match_MW
     model.MWs(pos_e,1) = mean(cell2mat(swissprot(:,5)))/1000;	%average g/mmol
@@ -59,14 +63,13 @@
 end
 
 %Update model.genes & model.pathways vectors:
-match_gen  = false;
 match_path = false;
 for i = 1:length(kegg)
     if strcmp(P,kegg{i,1})
         %Gene:
         if ~isempty(kegg{i,3}) && ~match_gen
-            match_gen            = true;
-            model.genes{pos_e,1} = kegg{i,3};
+            match_gen = true;
+            gene      = kegg{i,3};
         end
         %Pathway:
         if ~isempty(kegg{i,6}) && ~match_path
@@ -90,34 +93,38 @@
     if sum(unknowns) == 0
         idx = 0;
     else
-        unknowns  = model.genes(unknowns);
+        unknowns  = model.enzGenes(unknowns);
         pos_final = strfind(unknowns{end},'_')+1;
         idx       = str2double(unknowns{end}(pos_final:end));
     end
-    model.genes{pos_e,1} = ['unknown_' num2str(idx+1)];
+    gene = ['unknown_' num2str(idx+1)];
 end
 if ~match_path
     model.pathways{pos_e,1} = '-';
 end
 
+%Add gene to gene list if non-existing previously:
+if ~ismember(gene,model.genes)
+    model.enzNames(pos_e,1)
+    model = addGenes(model,{gene},'geneNames',model.enzNames(pos_e,1));
+end
+
 %Add exchange reaction of protein: -> P
-model = addReaction(model, ...                      %model
-                    ['prot_' P '_exchange'], ...    %rxn name
-                    {prot_name}, ...                %metabolites
-                    1, ...                          %stoichiometry
-                    true, ...                       %reversibility
-                    0, ...                          %LB
-                    Inf, ...                        %UB
-                    0, ...                          %c
-                    {''}, ...                       %subsystem
-                    model.genes{pos_e,1});          %gene rule
+rxnID = ['prot_' P '_exchange'];
+model = addReaction(model,rxnID, ...
+                    'metaboliteList', {prot_name}, ...
+                    'stoichCoeffList', 1, ...
+                    'lowerBound', 0, ...
+                    'upperBound', Inf);
+model.grRules{strcmp(model.rxns,rxnID)} = gene;
+model.enzGenes{pos_e,1} = gene;
 
 %Update metComps:
 pos_m = strcmp(model.mets,prot_name);   %position in model.mets
 if isfield(model,'compNames')
     cytIndex = find(strcmpi(model.compNames,'cytoplasm'),1);
     if ~isempty(cytIndex)
-        model.metComps(pos_m) = 2;              %For simplification all proteins are in cytosol
+        model.metComps(pos_m) = 2;	%For simplification all proteins are in cytosol
     else
         model.metComps(pos_m) = 1; 
     end

Matlab_Module/change_model/convertToEnzymeModel.m

@@ -1,7 +1,7 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % eModel = convertToEnzymeModel(model,uniprots,kcats)
 % Converts standard GEM to GEM accounting for enzymes as pseudo
-% metabolites, with -(MW/kcat) as the corresponding stoich. coeffs.
+% metabolites, with -(1/kcat) as the corresponding stoich. coeffs.
 %
 % INPUT:
 % model             The GEM structure (1x1 struct)
@@ -11,18 +11,26 @@
 % OUTPUT:
 % eModel            Modified GEM structure (1x1 struct)
 % 
-% Cheng Zhang. Last edited: 2016-10-30
+% Cheng Zhang. Last edited: 2018-05-24
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 function eModel = convertToEnzymeModel(irrevModel,uniprots,kcats)
 
+%Load databases:
+cd ../../Databases
+data      = load('ProtDatabase.mat');
+swissprot = data.swissprot;
+kegg      = data.kegg;
+cd ../Matlab_Module/change_model
+
 eModel  = irrevModel;
 enzymes = cell(5000,1);
 [m,n]   = size(uniprots);
 y       = 0;
+
 for i = 1:m
     rxnID = irrevModel.rxns{i};
-    x = 0;
+    x     = 0;
     for j = 1:n
         %Update vector enzymes and count the number of isozymes (x):
         if ~isempty(uniprots{i,j}) && kcats(i,j) > 0 % if kcat=0 no mets will be added
@@ -39,7 +47,7 @@
         %>1 enzyme: Will include an "arm reaction" for controlling the
         %total flux through the system of parallel rxns.
         if x > 1
-            eModel = addArmReaction(eModel,rxnID);
+            eModel  = addArmReaction(eModel,rxnID);
         end
         x = 0;
         %For each enzyme adds a new parallel rxn with that enzyme as pseudo metabolite.
@@ -52,8 +60,8 @@
                 newMets = cell(size(uniprots{i,j}));
                 for k = 1:length(newMets)
                     newMets{k} = ['prot_' uniprots{i,j}{k}];
-                end 
-                eModel = addEnzymesToRxn(eModel,kvalues,rxnID,newMets,{newID,newName}); 
+                end
+                eModel = addEnzymesToRxn(eModel,kvalues,rxnID,newMets,{newID,newName},kegg,swissprot); 
             end
         end
         eModel = removeRxns(eModel,{rxnID});  %Remove the original rxn
@@ -64,29 +72,23 @@
 enzymes(y+1:end) = [];
 enzymes          = unique(enzymes)';
 
-%Load databases:
-cd ../../Databases
-data      = load('ProtDatabase.mat');
-swissprot = data.swissprot;
-kegg      = data.kegg;
-cd ../Matlab_Module/change_model
-
-%Reset gene rules:
-eModel = rmfield(eModel,'grRules');
-eModel = rmfield(eModel,'rxnGeneMat');
-eModel.rules = cell(size(eModel.rxns));
+%Update rxnGeneMat:
+cd ../get_enzyme_data
+[~,rxnGeneMat]    = standardizeGrRules(eModel);
+eModel.rxnGeneMat = rxnGeneMat;
 
 %Create additional fields in model:
 eModel.enzymes   = cell(0,1);
-eModel.genes     = cell(0,1);
-eModel.geneNames = cell(0,1);
+eModel.enzGenes  = cell(0,1);
+eModel.enzNames  = cell(0,1);
 eModel.MWs       = zeros(0,1);
 eModel.sequences = cell(0,1);
 eModel.pathways  = cell(0,1);
+cd ../change_model
 for i = 1:length(enzymes)
     eModel = addProtein(eModel,enzymes{i},kegg,swissprot);
 end
-
+eModel.rules = eModel.grRules;
 end
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%