Skip to content

Commit

Permalink
Merge pull request #24 from SysBioChalmers/devel
Browse files Browse the repository at this point in the history 
Devel
  • Loading branch information
@IVANDOMENZAIN
IVANDOMENZAIN authored Apr 11, 2018
2 parents 1c23f50 + 5e70249 commit 5b1ef84
Show file tree
Hide file tree
Showing 46 changed files with 85,796 additions and 76,363 deletions.
Binary file modified Databases/PhylDist.mat
View file
Binary file not shown.
Binary file modified Databases/ProtDatabase.mat
100644 → 100755
View file
Binary file not shown.
44,030 changes: 22,871 additions & 21,159 deletions Databases/max_KCAT.txt
View file

Large diffs are not rendered by default.

38,748 changes: 22,619 additions & 16,129 deletions Databases/max_MW.txt
View file

Large diffs are not rendered by default.

27,780 changes: 14,044 additions & 13,736 deletions Databases/max_SA.txt
View file

Large diffs are not rendered by default.

12,866 changes: 6,433 additions & 6,433 deletions Databases/uniprot-organism%3Ayeast.tab
View file

Large diffs are not rendered by default.

104 changes: 104 additions & 0 deletions Matlab_Module/Kcat_sensitivity_analysis/changeMedia_batch.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,104 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% model = changeMedia(model,media,flux)
%
% function that modifies the ecModel and makes it suitable for batch growth
% simulations on different carbon sources.
%
% INPUT:
% - model: An enzyme constrained model
% - meadia: Media type ('YEP' for complex, 'MAA' minimal with Aminoacids,
% 'Min' for minimal media)
% - flux: (Optional) A cell array with measured uptake fluxes in mmol/gDwh
%
% OUTPUT:
% - model: The ECmodel with
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [model,pos] = changeMedia_batch(model,c_source,media,flux)

% Give the carbon source (c_source) input variable with the following
% format: c_source = 'D-glucose exchange (reversible)'

% first block any uptake
pos = getUptakeIndexes(model,c_source);
model.ub(pos) = 0;

% Block O2 and glucose production (for avoiding multiple solutions):
model.ub(strcmp(model.rxnNames,'oxygen exchange')) = 0;
model.ub(strcmp(model.rxnNames,'D-glucose exchange')) = 0;

%Find substrate production rxn and block it:
pos_rev = strcmpi(model.rxnNames,c_source(1:strfind(c_source,...
' (reversible)')-1));
model.ub(pos_rev) = 0;
%For growth on fructose and mannose the transport takes place in a passive
%way. [Boles & Hollenberg, 2006]
if strcmp(c_source,'D-fructose exchange (reversible)')
model.S(strcmp(model.mets,'s_0796'),strcmp(model.rxns,'r_1134')) = 0;
model.S(strcmp(model.mets,'s_0794'),strcmp(model.rxns,'r_1134')) = 0;
elseif strcmp(c_source,'D-mannose exchange (reversible)')
model.S(strcmp(model.mets,'s_0796'),strcmp(model.rxns,'r_1139')) = 0;
model.S(strcmp(model.mets,'s_0794'),strcmp(model.rxns,'r_1139')) = 0;
end
%The media will define which rxns to fix:
if strcmpi(media,'YEP')
N = 25; %Aminoacids + Nucleotides
elseif strcmpi(media,'MAA')
N = 21; %Aminoacids
elseif strcmpi(media,'Min')
N = 1; %Only the carbon source
end

%UB parameter (manually optimized for glucose on Min+AA):
b = 0.08;
%UB parameter (manually optimized for glucose complex media):
c = 2;


%Define fluxes in case of ec model:
if nargin < 5 %Limited protein
if N>1
flux = b*ones(1,N);
if N>21
flux(22:25) = c;
end
end
flux(1) = 1000;
end

%Fix values as UBs:
for i = 1:N
model.ub(pos(i)) = flux(i);
end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function pos = getUptakeIndexes(model,c_source)
pos(1) = find(strcmpi(model.rxnNames,c_source));
pos(2) = find(strcmpi(model.rxnNames,'L-alanine exchange (reversible)'));
pos(3) = find(strcmpi(model.rxnNames,'L-arginine exchange (reversible)'));
pos(4) = find(strcmpi(model.rxnNames,'L-asparagine exchange (reversible)'));
pos(5) = find(strcmpi(model.rxnNames,'L-aspartate exchange (reversible)'));
pos(6) = find(strcmpi(model.rxnNames,'L-cysteine exchange (reversible)'));
pos(7) = find(strcmpi(model.rxnNames,'L-glutamine exchange (reversible)'));
pos(8) = find(strcmpi(model.rxnNames,'L-glutamate exchange (reversible)'));
pos(9) = find(strcmpi(model.rxnNames,'glycine exchange (reversible)'));
pos(10) = find(strcmpi(model.rxnNames,'L-histidine exchange (reversible)'));
pos(11) = find(strcmpi(model.rxnNames,'L-isoleucine exchange (reversible)'));
pos(12) = find(strcmpi(model.rxnNames,'L-leucine exchange (reversible)'));
pos(13) = find(strcmpi(model.rxnNames,'L-lysine exchange (reversible)'));
pos(14) = find(strcmpi(model.rxnNames,'L-methionine exchange (reversible)'));
pos(15) = find(strcmpi(model.rxnNames,'L-phenylalanine exchange (reversible)'));
pos(16) = find(strcmpi(model.rxnNames,'L-proline exchange (reversible)'));
pos(17) = find(strcmpi(model.rxnNames,'L-serine exchange (reversible)'));
pos(18) = find(strcmpi(model.rxnNames,'L-threonine exchange (reversible)'));
pos(19) = find(strcmpi(model.rxnNames,'L-tryptophan exchange (reversible)'));
pos(20) = find(strcmpi(model.rxnNames,'L-tyrosine exchange (reversible)'));
pos(21) = find(strcmpi(model.rxnNames,'L-valine exchange (reversible)'));
pos(22) = find(strcmpi(model.rxnNames,'2''-deoxyadenosine exchange (reversible)'));
pos(23) = find(strcmpi(model.rxnNames,'2''-deoxyguanosine exchange (reversible)'));
pos(24) = find(strcmpi(model.rxnNames,'thymidine exchange (reversible)'));
pos(25) = find(strcmpi(model.rxnNames,'deoxycytidine exchange (reversible)'));
pos(26) = find(strcmpi(model.rxnNames,'D-glucose exchange (reversible)'));
end
67 changes: 67 additions & 0 deletions Matlab_Module/Kcat_sensitivity_analysis/findMaxValue.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,67 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%function [value,organism,parameter] = findMaxValue(EC_cell,BRENDA, SA_cell)
%
% Function that gets the maximum kinetic parameter (Kcat or S.A.*Mw) from
% the BRENDA files for the specified set of EC numbers. The algorithm also
% returns the organism and the parameter type (Kcat or S.A.) of the query.
%
% Ivan Domenzain Last edited. 2018-02-06
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [value,organism,parameter] = findMaxValue(EC_cell,BRENDA,SA_cell)
%Looks for the maximum turnover number available for the EC# associated
%with the uniprot code
EC_cell = strsplit(EC_cell,' ');
value = [];
organism = [];
parameter = [];
for i=1:length(EC_cell)
find_flag = false;
%In case that wild cards are present in the EC number the search on
%the BRENDA file will be relaxed.
if ~isempty(strfind(EC_cell{i},'-'))
EC_cell{i} = EC_cell{i}(strfind(EC_cell{i},'-')-1:end);
find_flag = true;
end
ECnumber = ['EC' EC_cell{i}];
Kcat = 0; orgK = '';
if find_flag == true
matching = find(~cellfun(@isempty,strfind(BRENDA{1},ECnumber)));
else
% If no wild cards are present the EC number search in the
% BRENDA file will look for an exact match
matching = find(strcmpi(ECnumber,BRENDA{1}));
end
%Gets the maximum Kcat value for the queried EC#
if ~isempty(matching)
[Kcat, maxIndx] = max(BRENDA{4}(matching));
orgK = BRENDA{3}(matching(maxIndx));
end
% Looks for the maximum SA*Mw value available for the EC number
SA_Mw = 0; orgS = '';
if find_flag == true
matching = find(~cellfun(@isempty,strfind(SA_cell{1},ECnumber)));
else
matching = find(strcmpi(ECnumber,SA_cell{1}));
end
%Gets the maximum SA*Mw value for the queried EC#
if ~isempty(matching)
[SA_Mw, maxIndx] = max(SA_cell{3}(matching));
SA_Mw = SA_Mw; %[1/hr]
orgS = SA_cell{2}(matching(maxIndx));
end
%Choose the maximal available value as a turnover number for the EC
value = [value; max(Kcat,SA_Mw)];

if Kcat>SA_Mw
organism = [organism; {orgK}];
parameter = [parameter; {'K_cat'}];
else
organism = [organism; {orgS}];
parameter = [parameter; {'SA*Mw'}];
end
end
[value, maxIndx] = max(value);
organism = organism(maxIndx);
parameter = parameter(maxIndx);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
203 changes: 203 additions & 0 deletions Matlab_Module/Kcat_sensitivity_analysis/findTopLimitations.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,203 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function [limKcats,limRxns,breakFlag] = TopLimitingKcat(model,prevUnicodes,...
% prevRxns)
%
% Function that gets an EC model and returns the top growth limiting Kcat
% value in the network based on a sensitivity analysis on each of the
% assigned parameters (growthRate control coefficients)
%
% If the model is not able to grow then it performs the analysis on reaction
% basis (relaxing all of its matched kcats) and gets the top growth
% limiting one, if any.
%
% Ivan Domenzain Last edited. 2018-02-07
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [limitations,breakFlag] = findTopLimitations(model,prevUnicodes,option)

%Get an initial simulation optimizing for biomass production
base_sol = solveLP(model);
%Extract the metabolite indexes corresponding to the model enzymes
%(avoiding the protein pool)
enz_indxs = find(~cellfun(@isempty,strfind(model.metNames,'prot_')));
enz_indxs = enz_indxs(1:end-1);
%Extract the enzyme usage pseudoreactions indexes
enzUsageIndxs = find(~cellfun(@isempty,strfind(model.rxns,'prot_')));
enzUsageIndxs = enzUsageIndxs(1:end-1);
limitations = [];
breakFlag = false;
%If the model is growing, analyse the individual kcat values of the flux
%carrier enzymes
if ~isempty(base_sol.f)% ~= 0
%First, find the growth limiting Kcats in the Batch model
limKcats = findLimitingKcat(model,base_sol,enzUsageIndxs);
if ~isempty(limKcats{1})
%Sort the limiting Kcats cell according to their control coefficient
limKcats = sortByCCoeff(limKcats,limKcats{5});
j = 1;
pos = 0;
% Choose the most controlling Kcat that is not included in the cell
% prevUnicodes for avoiding infinite loops
while j<=length(limKcats{1})
%The top limiting coefficient is identified with its corresponding
%uniprot code and the rxn index in the model because there might
%be promiscous enzymes catalysing different reactions at different rates
Protname = [limKcats{1}{j}(6:end) '_' num2str(limKcats{3}(j))];
if ~ismember(Protname,prevUnicodes)
pos = j;
break
else
j = j+1;
end
end

for i=1:length(limKcats)
if pos~=0
limKcats{i} = limKcats{1,i}(pos);
else
% If all the significant Limiting Kcats have been previously
% modified then the iterations should be stopped
limKcats{i} = [];
breakFlag = true;
end
end
end

if ~isempty(limKcats{1})
limitations = limKcats;
else
breakFlag = true;
end
%If the model is not growing, analyse all of the kinetic parameters
%associated to each of the metabolic reactions.
else
breakFlag = true;
if option == 1
limRxns = findLimitingRxn(model,enz_indxs,enzUsageIndxs);
limitations = limRxns;
elseif option == 2
limEnz = findLimitingEnz(model,enz_indxs,enzUsageIndxs);
limitations = limEnz;
end
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function output_cell = sortByCCoeff(input_cell,values)
[~,indexes] = sort(values,'descend');
for i=1:length(input_cell)
output_cell{i} = input_cell{i}(indexes);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function kcatIndxs = findLimitingKcat(model,base_sol,enzUsageIndxs)


kcatIndxs = cell(1,6);
%Get the original rxn indexes of the flux carrier enzymes (enzyme usages)
enzUsageFlux = base_sol.x(enzUsageIndxs);
fluxCarriers = find(enzUsageFlux);
fluxCarriers = enzUsageIndxs(fluxCarriers);
%Extract the metabolite index for the flux carrying enzymes
usedEnzIndxs = zeros(size(fluxCarriers));
for i=1:length(usedEnzIndxs)
usedEnzIndxs(i) = find(model.S(:,fluxCarriers(i))==1);
end
%For each flux carrying enzyme calculate the growthRate control
%coefficient (gRCC) on each of its kinetic coefficients
for i = 1:length(usedEnzIndxs)
enzPos = usedEnzIndxs(i);
enzName = model.mets(enzPos);
%Get the nonzero kinetic coefficients
enz_rxnsCoeffs = find(model.S(enzPos,:));
%For each enzymatic reaction (avoid enzyme usage pseudoRxn)
if ~isempty(enz_rxnsCoeffs)
for j=1:length(enz_rxnsCoeffs)-1
temp_model = model;
coeffPos = enz_rxnsCoeffs(j);
temp_model.S(enzPos,coeffPos) = model.S(enzPos,coeffPos)/1000;
new_sol = solveLP(temp_model);
if ~isempty(new_sol.f)
gRCC = (new_sol.f - base_sol.f)/base_sol.f;
if abs(gRCC) > 1e-2
Kcat = (-1/model.S(enzPos,coeffPos))/3600;
kcatIndxs{1} = [kcatIndxs{1}; enzName];
kcatIndxs{2} = [kcatIndxs{2}; enzPos];
kcatIndxs{3} = [kcatIndxs{3}; coeffPos];
kcatIndxs{4} = [kcatIndxs{4}; Kcat];
kcatIndxs{5} = [kcatIndxs{5}; gRCC];
kcatIndxs{6} = [kcatIndxs{6}; model.rxnNames(coeffPos)];
end
end
end
end
end

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [limRxns] = findLimitingRxn(model,enz_indxs,enzUsageIndxs)

limRxns = cell(1,3);
for i=1:(enzUsageIndxs(1)-1)
%disp(['Relaxing Kcat coefficients for Rxn#' num2str(i)])
temp_model = model;
%Search if there are nonZero Kcat coefficients matched to the
%i-th reaction
nonZero = find(temp_model.S(enz_indxs,i));
if ~isempty(nonZero)
disp(['Analyzing rxn: #' num2str(i)])
%Flexibilize all non-zero coefficients in the i-th
%metabolic reaction and check if any growth can be obtained
pos = enz_indxs(nonZero);
temp_model.S(pos,i) = temp_model.S(pos,i)/1e6;
new_sol = solveLP(temp_model);
deltaGR = abs(new_sol.f);
if deltaGR>0
disp(['limiting rxn found: #' num2str(i)])
disp(deltaGR)
limRxns{1} = [limRxns{1}; model.rxnNames(i)];
limRxns{2} = [limRxns{2}; i];
limRxns{3} = [limRxns{3}; deltaGR];
end
end
end
%If limiting rxns were found, sort them according to the growthRate
%difference
if ~isempty(limRxns{3})
limRxns = sortByCCoeff(limRxns,limRxns{3});
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [limEnz] = findLimitingEnz(model,enz_indxs,enzUsageIndxs)

limEnz = cell(1,3);
for i=1:length(enz_indxs)
enzPos = enz_indxs(i);
%disp(['Relaxing Kcat coefficients for Enzyme #' num2str(i)])
temp_model = model;
%Search if there are nonZero Kcat coefficients matched to the
%i-th enzyme
nonZero = find(temp_model.S(enzPos,1:enzUsageIndxs(1)-1));
if ~isempty(nonZero)
disp(['Analyzing enzyme: #' num2str(i)])
%Flexibilize all non-zero coefficients and check if any growth
%can be obtained
temp_model.S(enzPos,nonZero) = temp_model.S(enzPos,nonZero)/1e6;

new_sol = solveLP(temp_model);
deltaGR = abs(new_sol.f);

if deltaGR>0
disp(['limiting enzyme found: #' num2str(i)])
disp(deltaGR)
limEnz{1} = [limEnz{1}; model.metNames(i)];
limEnz{2} = [limEnz{2}; i];
limEnz{3} = [limEnz{3}; deltaGR];
end
end
end
%If limiting rxns were found, sort them according to the growthRate
%difference
if ~isempty(limEnz{3})
limEnz = sortByCCoeff(limEnz,limEnz{3});
end
end
Loading

0 comments on commit 5b1ef84

Please sign in to comment.