diff --git a/src/design/TrimGdel/GRPRchecker.m b/src/design/TrimGdel/GRPRchecker.m new file mode 100644 index 0000000000..5a99f540e3 --- /dev/null +++ b/src/design/TrimGdel/GRPRchecker.m @@ -0,0 +1,92 @@ +function [GR ,PR] = GRPRchecker(model, targetMet, gvalue) +% GRPRchecker calculates the maximum GR and the minimu PR +% under the GR maximization when a constraint-based model, a target +% metabolite, and a gene deletion stratety are given. +% +% USAGE: +% +% function [GR, PR] +% = GRPRchecker(model, targetMet, givenGvalue) +% +% INPUTS: +% model: COBRA model structure containing the following required fields to perform gDel_minRN. +% +% *.rxns: Rxns in the model +% *.mets: Metabolites in the model +% *.genes: Genes in the model +% *.grRules: Gene-protein-reaction relations in the model +% *.S: Stoichiometric matrix (sparse) +% *.b: RHS of Sv = b (usually zeros) +% *.c: Objective coefficients +% *.lb: Lower bounds for fluxes +% *.ub: Upper bounds for fluxes +% *.rev: Reversibility of fluxes +% +% targetMet: target metabolites (e.g., 'btn_c') +% gvalue: The first column is the list of genes in the original model. +% The second column contains a 0/1 vector indicating which genes should be deleted. +% 0: indicates genes to be deleted. +% 1: indecates genes to be remained. +% +% OUTPUTS: +% GR: the growth rate obained when the gene deletion strategy is +% applied and the growth rate is maximized. +% PR: the minimum target metabolite production rate obained +% when the gene deletion strategy is applied and the growth rate is maximized. +% +% .. Author: - Takeyuki Tamura, Mar 06, 2025 +% + + +[model, targetRID, extype] = modelSetting(model, targetMet); + +m = size(model.mets, 1); +n = size(model.rxns, 1); +g = size(model.genes, 1); +gid = find(model.c); +pid = targetRID; + + +model2 = model; +[grRules0] = calculateGR(model, gvalue); +lb2 = model.lb; +ub2 = model.ub; + +for i=1:n + if grRules0{i, 4} == 0 + lb2(i) = 0; + ub2(i) = 0; + end +end + +gm.A = sparse(model.S); +gm.obj = -model.c; +gm.modelsense = 'Min'; +gm.sense = repmat('=', 1, size(model.S, 1)); +gm.lb = lb2; +gm.ub = ub2; +opt0 = gurobi(gm); + +[opt0.x(gid) opt0.x(pid)] + +GR0 = -opt0.objval; +lb2(gid) = GR0; +ub2(gid) = GR0; +model2.c(gid) = 0; +model2.c(pid) = 1; + +gm2.A = sparse(model.S); +gm2.obj = model2.c; +gm2.modelsense = 'Min'; +gm2.sense = repmat('=', 1, size(model.S, 1)); +gm2.lb = lb2; +gm2.ub = ub2; +opt1 = gurobi(gm2); + +GR = GR0 +PR = opt1.x(pid) +[GR PR] + +return; +end + diff --git a/src/design/TrimGdel/TrimGdel.m b/src/design/TrimGdel/TrimGdel.m new file mode 100644 index 0000000000..0e57b41e10 --- /dev/null +++ b/src/design/TrimGdel/TrimGdel.m @@ -0,0 +1,89 @@ +function [gvalue, GR, PR, size1, size2, size3, success] = TrimGdel(model, targetMet, maxLoop, PRLB, GRLB) +% TrimGdel appropriately considers GPR rules and determines +% a minimal gene deletion strategies to achieve growth-coupled production +% for a given target metabolite and a genome-scale model. +% even in the worst-case analysis (ensures the weak-growth-coupled production). +% +% Gurobi is required for this version. +% The CPLEX version is available on https://github.com/MetNetComp/TrimGdel +% +% USAGE: +% +% function [gvalue, GR, PR, size1, size2, size3, success] +% = TrimGdel(model, targetMet, maxLoop, PRLB, GRLB) +% +% INPUTS: +% model: COBRA model structure containing the following required fields to perform gDel_minRN. +% +% *.rxns: Rxns in the model +% *.mets: Metabolites in the model +% *.genes: Genes in the model +% *.grRules: Gene-protein-reaction relations in the model +% *.S: Stoichiometric matrix (sparse) +% *.b: RHS of Sv = b (usually zeros) +% *.c: Objective coefficients +% *.lb: Lower bounds for fluxes +% *.ub: Upper bounds for fluxes +% *.rev: Reversibility of fluxes +% +% targetMet: target metabolites (e.g., 'btn_c') +% maxLoop: the maximum number of iterations in gDel_minRN +% PRLB: the minimum required production rates of the target metabolites +% when gDel-minRN searches the gene deletion +% strategy candidates. +% (But it is not ensured to achieve this minimum required value +% when GR is maximized withoug PRLB.) +% GRLB: the minimum required growth rate +% when gDel-minRN searches the gene deletion +% strategy candidates. +% +% OUTPUTS: +% gvalue: a small gene deletion strategy (obtained by TrimGdel). +% The first column is the list of genes. +% The second column is a 0/1 vector indicating which genes should be deleted. +% 0: indicates genes to be deleted. +% 1: indecates genes to be remained. +% GR: the maximum growth rate when the obtained gene deletion +% strategy represented by gvalue is applied. +% PR: the minimum production rate of the target metabolite under +% the maximization of the growth rate when the obtained gene deletion +% strategy represented by gvalue is applied. +% size1: the number of gene deletions after Step1. +% size2: the number of gene deletions after Step2. +% size3: the number of gene deletions after Step3. +% success: indicates whether TrimGdel obained an appropriate gene +% deletion strategy. (1:success, 0:failure) +% +% NOTE: +% +% T. Tamura, "Trimming Gene Deletion Strategies for Growth-Coupled +% Production in Constraint-Based Metabolic Networks: TrimGdel," +% in IEEE/ACM Transactions on Computational Biology and Bioinformatics, +% vol. 20, no. 2, pp. 1540-1549, 2023. +% +% Comprehensive computational results are accumulated in MetNetComp +% database. +% https://metnetcomp.github.io/database1/indexFiles/index.html +% +% T. Tamura, "MetNetComp: Database for Minimal and Maximal Gene-Deletion Strategies +% for Growth-Coupled Production of Genome-Scale Metabolic Networks," +% in IEEE/ACM Transactions on Computational Biology and Bioinformatics, +% vol. 20, no. 6, pp. 3748-3758, 2023, +% +% .. Author: - Takeyuki Tamura, Mar 06, 2025 +% + +[gvalue gr pr it success] = gDel_minRN(model, targetMet, maxLoop, PRLB, GRLB) % Step 1 +if success + [gvalue, GR, PR, size1, size2, size3] = step2and3(model, targetMet, gvalue) % Step 2 and 3 +else + gvalue = []; + GR = 0; + PR = 0; + size1 = 0; + size2 = 0; + size3 = 0; +end + +end + diff --git a/src/design/TrimGdel/calculateGR.m b/src/design/TrimGdel/calculateGR.m new file mode 100644 index 0000000000..c81216a08c --- /dev/null +++ b/src/design/TrimGdel/calculateGR.m @@ -0,0 +1,77 @@ +function [grRules] = calculateGR(model, gvalue) +% calculateGR is a function of gDel_minRN that reads +% a COBRA model and a 0-1 assignment for genes and outputs whether +% each reaction is repressed or not. +% +% USAGE: +% +% function [grRules] = calculateGR(model, gvalue) +% +% INPUTS: +% model: COBRA model structure containing the following required fields to perform gDel_minRN. +% +% *.rxns: Rxns in the model +% *.mets: Metabolites in the model +% *.genes: Genes in the model +% *.grRules: Gene-protein-reaction relations in the model +% *.S: Stoichiometric matrix (sparse) +% *.b: RHS of Sv = b (usually zeros) +% *.c: Objective coefficients +% *.lb: Lower bounds for fluxes +% *.ub: Upper bounds for fluxes +% *.rev: Reversibility of fluxes +% +% gvalue: The first column is the list of genes in the original model. +% The second column contains a 0/1 vector indicating which genes should be deleted. +% 0: indicates genes to be deleted. +% 1: indecates genes to be remained. +% +% OUTPUT: +% grRules: The first column is the list of GPR-rules. If a reaction does +% not have a GPR-rule, it is represented as 1. +% In the second columun, each gene is converted to 0 or 1 based +% on the given 0-1 assignment, with AND converted to * and OR +% converted to +. +% The third comumn contains the calculation results from the +% second column. +% The fourth column is 0 if the third column is 0 and 1 if it is +% greater. If it is 0, the reaction is repressed; if it is 1, it +% is not repressed. +% +% .. Author: - Takeyuki Tamura, Mar 06, 2025 +% + +grRules = cell(size(model.rxns)); +for i=1:size(model.grRules, 1) + grRules{i, 1} = model.grRules{i,1}; +end +for i = 1:size(model.rxns, 1) + if isempty(grRules{i, 1})==1 + grRules{i,1} = '1'; + end +end +grRules(:, 2) = strrep(grRules, 'or', '+'); +grRules(:,2) = strrep(grRules(:,2), 'and', '*'); + +[xname2, index] = sortrows(gvalue(:,1), 'descend'); +for i=1:size(index, 1) + sorted_gvalue(i, 1) = gvalue{index(i, 1), 2}; +end +for i = 1:size(model.genes, 1) + grRules(:, 2) = strrep(grRules(:, 2), xname2{i, 1},num2str(sorted_gvalue(i, 1))); +end +for i = 1:size(grRules, 1) + %i + if isempty(grRules{i, 2}) == 0 + grRules{i, 3} = eval(grRules{i, 2}); + if grRules{i, 3} > 0.9 + grRules{i, 4} = 1; + else + grRules{i, 4} = 0; + end + else + grRules{i, 4} = -1; + end +end +end + diff --git a/src/design/TrimGdel/example1.m b/src/design/TrimGdel/example1.m new file mode 100644 index 0000000000..85b54e5d52 --- /dev/null +++ b/src/design/TrimGdel/example1.m @@ -0,0 +1,19 @@ +function [] = example1() +% example1 calculates the gene deletion strategy for growth coupling +% for succinate in e_coli_core. +% +% USAGE: +% +% function [] = example1() +% +% .. Author: - Takeyuki Tamura, Mar 05, 2025 +% + + +load('e_coli_core.mat'); +model = e_coli_core; + +[gvalue, GR, PR, size1, size2, size3, success] = TrimGdel(model, 'succ_e', 10, 0.1, 0.1) + +end + diff --git a/src/design/TrimGdel/example2.m b/src/design/TrimGdel/example2.m new file mode 100644 index 0000000000..c25df172b6 --- /dev/null +++ b/src/design/TrimGdel/example2.m @@ -0,0 +1,17 @@ +function [outputArg1, outputArg2] = example2() +% example2 calculates the gene deletion strategy for growth coupling +% for biotin in iML1515. +% +% USAGE: +% +% function [] = example2() +% +% .. Author: - Takeyuki Tamura, Mar 06, 2025 +% + +load('iML1515.mat'); +model = iML1515; +[gvalue, GR, PR, size1, size2, size3, success] = TrimGdel(model, 'btn_c', 10, 0.1, 0.1) + +end + diff --git a/src/design/TrimGdel/example3.m b/src/design/TrimGdel/example3.m new file mode 100644 index 0000000000..3fc1788208 --- /dev/null +++ b/src/design/TrimGdel/example3.m @@ -0,0 +1,18 @@ +function [outputArg1, outputArg2] = example3() +% example3 calculates the gene deletion strategy for growth coupling +% for riboflavin in iML1515. +% +% USAGE: +% +% function [] = example3() +% +% .. Author: - Takeyuki Tamura, Mar 06, 2025 +% + +load('iML1515.mat'); +model = iML1515; + +[gvalue, GR, PR, size1, size2, size3, success] = TrimGdel(model, 'ribflv_c', 10, 0.1, 0.1) + +end + diff --git a/src/design/TrimGdel/example4.m b/src/design/TrimGdel/example4.m new file mode 100644 index 0000000000..3b2a6718a9 --- /dev/null +++ b/src/design/TrimGdel/example4.m @@ -0,0 +1,18 @@ +function [outputArg1, outputArg2] = exampl4() +% example4 calculates the gene deletion strategy for growth coupling +% for pantothenate in iML1515. +% +% USAGE: +% +% function [] = example4() +% +% .. Author: - Takeyuki Tamura, Mar 06, 2025 +% + +load('iML1515.mat'); +model = iML1515; + +[gvalue, GR, PR, size1, size2, size3, success] = TrimGdel(model, 'pnto__R_c', 10, 0.1, 0.1) + +end + diff --git a/src/design/TrimGdel/example5.m b/src/design/TrimGdel/example5.m new file mode 100644 index 0000000000..a094ca591f --- /dev/null +++ b/src/design/TrimGdel/example5.m @@ -0,0 +1,18 @@ +function [outputArg1, outputArg2] = example5() +% example5 calculates the gene deletion strategy for growth coupling +% for succinate in iMM904. +% +% USAGE: +% +% function [] = example5() +% +% .. Author: - Takeyuki Tamura, Mar 06, 2025 +% + +load('iMM904.mat'); +model = iMM904; + +[gvalue, GR, PR, size1, size2, size3, success] = TrimGdel(model, 'succ_e', 10, 0.1, 0.1) + +end + diff --git a/src/design/TrimGdel/gDel_minRN.m b/src/design/TrimGdel/gDel_minRN.m new file mode 100644 index 0000000000..55278a0d09 --- /dev/null +++ b/src/design/TrimGdel/gDel_minRN.m @@ -0,0 +1,254 @@ +function [gvalue gr pr it success] = gDel_minRN(model, targetMet, maxLoop, PRLB, GRLB) +% gDel-minRN (Step 1 of TrimGdel) determines gene deletion strategies +% by mixed integer linear programming to achieve growth coupling +% for the target metabolite by repressing the maximum number of reactions +% via gene-protein-reaction relations. +% +% USAGE: +% +% function [vg gr pr it success] +% = gDel_minRN(model, targetMet, maxLoop, PRLB, GRLB) +% +% INPUTS: +% model: COBRA model structure containing the following required fields to perform gDel_minRN. +% +% *.rxns: Rxns in the model +% *.mets: Metabolites in the model +% *.genes: Genes in the model +% *.grRules: Gene-protein-reaction relations in the model +% *.S: Stoichiometric matrix (sparse) +% *.b: RHS of Sv = b (usually zeros) +% *.c: Objective coefficients +% *.lb: Lower bounds for fluxes +% *.ub: Upper bounds for fluxes +% *.rev: Reversibility of fluxes +% +% targetMet: target metabolites (e.g., 'btn_c') +% maxLoop: the maximum number of iterations in gDel_minRN +% PRLB: the minimum required production rates of the target metabolites +% when gDel-minRN searches the gene deletion strategy candidates. +% (But it is not ensured to achieve this minimum required value +% when GR is maximized withoug PRLB.) +% GRLB: the minimum required growth rate when gDel-minRN searches +% the gene deletion strategy candidates. +% +% OUTPUTS: +% gvalue: The first column is the list of genes in the original model. +% The second column contains a 0/1 vector indicating which genes should be deleted. +% 0: indicates genes to be deleted. +% 1: indecates genes to be remained. +% gr: the growth rate obained when the gene deletion strategy is +% applied and the growth rate is maximized. +% pr: the target metabolite production rate obained +% when the gene deletion strategy is applied and the growth rate is maximized. +% it: indicates how many iterations were necessary to obtain the solution. +% success: indicates whether gDel_minRN obained an appropriate gene +% deletion strategy. (1:success, 0:failure) +% +% .. Author: - Takeyuki Tamura, Mar 06, 2025 +% + +tic; +ori_model = model; +n = size(model.rxns, 1); +for i=1:n + if model.ub(i) > 9999 + model.ub(i) = 1000; + end + if model.lb(i) < -9999 + model.lb(i) = -1000; + end +end + +gvalue = []; +gr = -1; pr = -1; it = 0; success = 0; + +params.IntFeasTol = 1e-09; +changeTolerance = 1000; + +sss=sprintf('gDel-minRN.mat'); + +[ori_model, targetRID, extype] = modelSetting(ori_model, targetMet); +[model, targetRID, extype] = modelSetting(model, targetMet); + +m = size(model.mets, 1); +n = size(model.rxns, 1); +g = size(model.genes, 1); +vg = zeros(g, 1); +gid = find(model.c); +pid = targetRID; +model2 = model; +model2.c(gid) = 0; +model2.c(targetRID) = 1; + +gm2.obj = -model2.c; +gm2.A = sparse(model2.S); +gm2.modelsense = 'Min'; +gm2.sense = repmat('=', 1, m); +gm2.lb = model2.lb; +gm2.ub = model2.ub; +optPre = gurobi(gm2); + +if strcmp(optPre.status, 'OPTIMAL') ~= 1 + display('no solution 1') + return; +elseif -optPre.objval < PRLB + display('TMPR < PRLB') + return; +end + +model2 = model; +model.lb(pid) = PRLB; +model.lb(gid) = GRLB; + +gm.obj = -model.c; +gm.A = sparse(model.S); +gm.modelsense = 'Min'; +gm.sense = repmat('=', 1, m); +gm.lb = model.lb; +gm.ub = model.ub; +opt0 = gurobi(gm); + +TMGR = -opt0.objval; +big = TMGR; +[term, ng, nt, nr, nko, reactionKO, reactionKO2term] = readGeneRules(model); + [f, A, b, Aeq, beq, lb, ub, xname] = geneReactionMILP(model, term, ng, nt, nr, nko); + + lp.Aeq = [model.S zeros(m, ng+nt+nko+nko); + zeros(size(Aeq, 1), n) Aeq zeros(size(Aeq, 1),nko); + zeros(nko, nr+ng+nt) changeTolerance*eye(nko) -1*eye(nko)]; +lp.beq = [zeros(m, 1); zeros(size(Aeq, 1),1); zeros(nko, 1)]; +j = 1; +for i=1:size(model.grRules, 1) + if isempty(model.grRules{i, :}) == 0 + ind(1,j) = i; + j = j+1; + end +end +z1 = -diag(model.ub); +z2 = diag(model.lb); +z3 = eye(n); + +lp.A = [zeros(size(A,1), n) A zeros(size(A,1), nko); + z3(ind,:) zeros(size(ind,2), ng+nt) z1(ind, ind) zeros(size(ind, 2), nko); + -z3(ind, :) zeros(size(ind, 2),ng+nt) z2(ind, ind) zeros(size(ind, 2),nko)]; +lp.b = [b; zeros(2*nko, 1)]; +lp.lb = [model.lb; lb; zeros(nko, 1)]; +lp.ub = [model.ub; ub; changeTolerance*ones(nko, 1)]; +lp.f = [-model.c; zeros(ng+nt, 1); big*ones(nko, 1); zeros(nko, 1)]; +for i = 1:n + s1 = repelem('C', n); + s2 = repelem('B', ng+nt+nko); + s3 = repelem('I', nko); + lp.ctype = sprintf('%s%s%s', s1, s2, s3); +end +A2 = lp.A; +b2 = lp.b; + +it = 1; +while it <= maxLoop + it + gr = -1; pr = -1; + + gm.obj = lp.f; + gm.A = sparse([lp.A; lp.Aeq]); + gm.rhs =[lp.b; lp.beq]; + gm.modelsense = 'Min'; + gm.sense = horzcat(repmat('<', 1, size(lp.A, 1)), repmat('=', 1, size(lp.Aeq, 1))); + gm.lb = lp.lb; + gm.ub = lp.ub; + gm.vtype = lp.ctype; + % find a gene deletion strategy candidate + opt = gurobi(gm,params); + + + if strcmp(opt.status, 'OPTIMAL') + for i = 1:n + vx(i, it) = opt.x(i); + result{i, it+1} = opt.x(i); + result{i, 1} = model.rxns{i}; + end + for i=1:ng + vg(i, it) = opt.x(n+i); + result{n+i, 1} = xname{i}; + result{n+i, it+1} = vg(i,it); + gvalue{i, 1} = xname{i}; + gvalue{i, 2} = vg(i, it) > 0.1; + end + for i = 1:nt + vt(i, 1) = opt.x(n+ng+i); + result{n+ng+i, 1} = xname{ng+i}; + result{n+ng+i, it+1} = opt.x(n+ng+i); + end + for i=1:nko + vko(i, it) = opt.x(n+ng+nt+i); + vnewadd(i,it) = opt.x(n+ng+nt+nko+i); + result{n+ng+nt+i, 1} = xname{ng+nt+i}; + result{n+ng+nt+i, it+1} = opt.x(n+ng+nt+i); + end + + else + system('rm -f clone*.log'); + if it == 1 + display('no solution 2') + return; + end + display('no more candidates') + system('rm -f clone*.log'); + return; + end + [grRules] = calculateGR(ori_model, gvalue); + + lb2 = ori_model.lb; + ub2 = ori_model.ub; + for i = 1:nr + if grRules{i, 4} == 0 + lb2(i) = 0; + ub2(i) = 0; + end + end + + gm2.obj = -ori_model.c; + gm2.modelsense = 'Min'; + gm2.A = sparse(model.S); + gm2.sense = repmat('=', 1, m); + gm2.lb = lb2; + gm2.ub = ub2; + gm2.rhs = zeros(m,1); + % validate the gene deletion candidate + opt2 = gurobi(gm2,params); + + gm3 = gm2; + gm3.obj(gid) = 0; + gm3.obj(pid) = 1; + gm3.lb(gid) = opt2.x(gid); + gm3.ub(gid) = opt2.x(gid); + % evaluate the minimum PR under the GR maximazation. + opt3 = gurobi(gm3); + + grprList(it, :) = [opt2.x(gid) opt3.x(pid)]; + gr = opt2.x(gid); pr = opt3.x(pid); + result2(:, it) = opt2.x; + + if (opt2.x(gid) >= GRLB) && (opt3.x(pid) >= PRLB) + [opt2.x(gid) opt3.x(pid)]; + vg(:,it); + success = 1; + time = toc; + system('rm -f clone*.log'); + return; + end + + zeroList(:,it) = vko(:, it) < 0.01; + dA = [zeros(1, nr+ng+nt) -(zeroList(:, it))' zeros(1, nko)]; + db = -1; + lp.A = [lp.A; dA]; + lp.b = [lp.b; db]; + + it = it + 1; + system('rm -f clone*.log'); +end +vg = vg > 0.1; + +end + diff --git a/src/design/TrimGdel/geneReactionMILP.m b/src/design/TrimGdel/geneReactionMILP.m new file mode 100644 index 0000000000..c6d38012ee --- /dev/null +++ b/src/design/TrimGdel/geneReactionMILP.m @@ -0,0 +1,110 @@ +function [f, A, b, Aeq, beq, lb, ub, xname] = geneReactionMILP(model, term, ng, nt, nr, nko) +% geneReactionMILP is a submodule to convert GPR relations to MILP. +% +% USAGE: +% +% function [f, A, b, Aeq, beq, lb, ub, xname] +% = geneReactionMILP(model, term, ng, nt, nr, nko) +% +% INPUTS: +% model: COBRA model structure containing the following required fields to perform gDel_minRN. +% +% *.rxns: Rxns in the model +% *.mets: Metabolites in the model +% *.genes: Genes in the model +% *.grRules: Gene-protein-reaction relations in the model +% *.S: Stoichiometric matrix (sparse) +% *.b: RHS of Sv = b (usually zeros) +% *.c: Objective coefficients +% *.lb: Lower bounds for fluxes +% *.ub: Upper bounds for fluxes +% *.rev: Reversibility of fluxes +% +% term: the list of Boolean functions extracted from the gene-protein-reaction relations +% ng: the number of genes +% nt: the number of internal terms +% nr: the number of reactions +% nko: the number of repressible reactions +% +% OUTPUTS: +% f: the weighted vector for the objective function of the resulting MILP. +% A, b, Aeq, beq, lb, ub: correspond to each matrix included in the +% following part of the MILP constraints. +% +% A*x <= b +% Aeq*x=beq +% lb <= x <= ub +% +% xname: variable names in the resulting MILP. +% +% .. Author: - Takeyuki Tamura, Mar 06, 2025 +% + +n_and = 0; n_or = 0; n_equal = 0; +for i=1:size(term, 2) + switch char(term(i).function) + case 'or' + n_or = n_or+1; + case 'and' + n_and = n_and+1; + case 'equal' + n_equal = n_equal+1; + end +end + +for i=1:ng + xname{i,1} = model.genes{i}; +end +for i=1:nt+nko + xname{ng+i, 1} = term(i).output; +end + +n_column = ng+nt+nko; +n_row = 2*(n_and + n_or); +A = zeros(n_row, n_column); +b = zeros(n_row, 1); +Aeq = zeros(n_equal, n_column); +beq = zeros(n_equal, 1); +lb = zeros(n_column, 1); +ub = ones(n_column, 1); +intcon = [1:n_column]; +f = zeros(n_column, 1); +for i=1:ng+nt + f(i,1) = 1; +end + +jj = 1; kk = 1; +for i=1:size(term, 2) + + k = size(term(i).input, 1); + x = find(strcmp(xname(:,1), term(i).output)); + switch char(term(i).function) + case 'or' + A(jj, x) = -k; + A(jj+1, x) = 1; + for j=1:k + x = find(strcmp(xname(:,1), term(i).input{j})); + A(jj, x) = 1; + A(jj+1, x) = -1; + end + jj = jj + 2; + + case 'and' + A(jj, x) = k; + A(jj+1, x) = -1; + for j=1:k + x = find(strcmp(xname(:,1), term(i).input{j})); + A(jj, x) = -1; + A(jj+1, x) = 1; + b(jj+1, 1) = k-1; + end + jj = jj + 2; + case 'equal' + Aeq(kk, x) = 1; + x = find(strcmp(xname(:, 1), term(i).input{1})); + Aeq(kk, x) = -1; + kk = kk+1; + end +end +end + diff --git a/src/design/TrimGdel/modelSetting.m b/src/design/TrimGdel/modelSetting.m new file mode 100644 index 0000000000..351aeaa768 --- /dev/null +++ b/src/design/TrimGdel/modelSetting.m @@ -0,0 +1,59 @@ +function [model2, targetRID, extype] = modelSetting(model, targetMet) +% modelSetting is a function of gDel_minRN that adds an auxiliary +% exchange reaction for the target metabolite when there is no +% original corresponding exchange reaction. +% +% USAGE: +% +% function [model2,targetRID,extype] = modelSetting(model,targetMet) +% +% INPUTS: +% model: COBRA model structure containing the following required fields to perform gDel_minRN. +% +% *.rxns: Rxns in the model +% *.mets: Metabolites in the model +% *.genes: Genes in the model +% *.grRules: Gene-protein-reaction relations in the model +% *.S: Stoichiometric matrix (sparse) +% *.b: RHS of Sv = b (usually zeros) +% *.c: Objective coefficients +% *.lb: Lower bounds for fluxes +% *.ub: Upper bounds for fluxes +% *.rev: Reversibility of fluxes +% +% targetMet: target metabolites (e.g., 'btn_c') +% +% OUTPUTS: +% model2: the model that has the exchange reaction of the target metabolite +% targetRID: ID of the exchange reaction of the target metabolite +% extype: indicates that an auxiliary exchange reaction was added. +% 1,2: there was the corresponding exchange reaction. +% 3: An auxiliary exchange reaction was added. +% +% .. Author: - Takeyuki Tamura, Mar 06, 2025 +% + +target = findMetIDs(model, targetMet); +m = size(model.mets, 1); +n = size(model.rxns, 1); +if isempty(find(strcmp(model.rxns, strcat('EX_',targetMet))))==0 + targetRID = find(strcmp(model.rxns, strcat('EX_',targetMet))); + model2 = model; + extype = 1; +elseif isempty(find(strcmp(model.rxns, strcat('DM_',targetMet))))==0 + targetRID = find(strcmp(model.rxns, strcat('DM_',targetMet))); + model2 = model; + extype = 2; +else + [model2, rxnIDexists] = addReaction(model, 'Transport', {targetMet}, [-1]); + m = size(model2.mets, 1); + n = size(model2.rxns, 1); + model2.S(target, n) = -1; + model2.ub(n) = 999999; + model2.lb(n) = 0; + model2.rev(n) = 0; + targetRID = n; + extype = 3; +end +end + diff --git a/src/design/TrimGdel/readGeneRules.m b/src/design/TrimGdel/readGeneRules.m new file mode 100644 index 0000000000..0491f5a903 --- /dev/null +++ b/src/design/TrimGdel/readGeneRules.m @@ -0,0 +1,115 @@ +function [term, ng, nt, nr, nko, reactionKO, reactionKO2term] = readGeneRules(model) +% readGeneRules is a function of gDel_minRN that reads +% gene-protein-reaction relations and outputs the necessary +% information for the MILP formalization +% +% USAGE: +% function [term, ng, nt, nr, nko, reactionKO, reactionKO2term] = readGeneRules(model) +% +% INPUTS: +% model: COBRA model structure containing the following required fields to perform gDel_minRN. +% +% *.rxns: Rxns in the model +% *.mets: Metabolites in the model +% *.genes: Genes in the model +% *.grRules: Gene-protein-reaction relations in the model +% *.S: Stoichiometric matrix (sparse) +% *.b: RHS of Sv = b (usually zeros) +% *.c: Objective coefficients +% *.lb: Lower bounds for fluxes +% *.ub: Upper bounds for fluxes +% *.rev: Reversibility of fluxes +% +% OUTPUTS: +% term: the list of Boolean functions extracted from the gene-protein-reaction relations +% ng: the number of genes +% nt: the number of internal terms +% nr: the number of reactions +% nko: the number of repressible reactions +% +% .. Author: - Takeyuki Tamura, Mar 06, 2025 +% + +x = 1; y = 1; reactionKO = 0; qq = 0; ww = 0; +for i=1:size(model.grRules) + GRrelation = model.grRules{i}; + empty = 0; + if isempty(GRrelation) == 0 + if (isempty(strfind(GRrelation, 'or'))==0) || (isempty(strfind(GRrelation, 'and')) == 0) + GRrelation=strcat('(', GRrelation, ')'); + else + GRrelation = strtrim(GRrelation); + term(x).output = sprintf('reactionKO_%d', i); + term(x).function = 'equal'; + term(x).input{1, 1} = GRrelation; + reactionKO2term(i, 1) = x; + x = x+1; + reactionKO = reactionKO + 1; + qq = qq+1; + end + else + empty = 1; + end + flag = 1; + while isempty(strfind(GRrelation, ')')) == 0 + size(GRrelation, 2); + p_tail_list = strfind(GRrelation, ')'); + p_tail = p_tail_list(1); + j = p_tail; + flag = 1; + while flag == 1 + j = j-1; + if GRrelation(j) == '(' + p_head = j; + flag = 0; + end + end + extract = GRrelation(p_head+1:p_tail-1); + extract2 = strtrim(extract); + term(x).output = sprintf('term_%d', y); + term(x).output; + y = y+1; + unit = strsplit(extract2); + for k=1:floor(size(unit, 2)/2) + if strcmp(unit(2), unit(2*k)) == 0 + save('readGeneRules.mat'); + error('Boolean functions are inappropriate.') + end + end + + term(x).function = unit(2); + for k=1:ceil(size(unit,2)/2) + term(x).input{k,1} = unit(2*k-1); + end + str = sprintf('(%s)', extract); + GRrelation = strrep(GRrelation, str, term(x).output); + x = x+1; + end + if (empty == 0 && strcmp(term(x-1).function, 'equal') == 0) + term(x-1).output = sprintf('reactionKO_%d', i); + reactionKO2term(i, 1) = x-1; + reactionKO = reactionKO + 1; + ww = ww + 1; + end +end +ng = size(model.genes, 1); +nr = size(model.rxns, 1); +nko = ww + qq; +nt = y - ww - 1; +j = 1; +for i=1:nt+nko + if contains(term(i).output, 'term') == 1 + term2(j) = term(i); + j = j + 1; + end +end +for i=1:nt+nko + if contains(term(i).output, 'reaction') == 1 + term2(j) = term(i); + j = j + 1; + end +end +term = term2; + +end + diff --git a/src/design/TrimGdel/step2and3.m b/src/design/TrimGdel/step2and3.m new file mode 100644 index 0000000000..8ca2f6b656 --- /dev/null +++ b/src/design/TrimGdel/step2and3.m @@ -0,0 +1,250 @@ +function [gvalue, GR, PR, size1, size2, size3] = step2and3(model, targetMet, givenGvalue) +% The function step2and3 implements Step 2 and Step 3 of TrimGdel. +% Step 2 minimizes the number of deleted genes while maintaining +% which reactions are repressed. +% Step 3 trims unnecessary deleted genes while maintaining GR and PR +% at the maximization of GR. +% +% USAGE: +% +% function [gvalue, finalGRPR, size1, size2, size3] = step2and3(model, targetMet, givenGvalue) +% +% INPUTS: +% +% model: COBRA model structure containing the following required fields to perform gDel_minRN. +% *.rxns: Rxns in the model +% *.mets: Metabolites in the model +% *.genes: Genes in the model +% *.grRules: Gene-protein-reaction relations in the model +% *.S: Stoichiometric matrix (sparse) +% *.b: RHS of Sv = b (usually zeros) +% *.c: Objective coefficients +% *.lb: Lower bounds for fluxes +% *.ub: Upper bounds for fluxes +% *.rev: Reversibility of fluxes +% +% targetMet: target metabolites (e.g., 'btn_c') +% givenGvalue: a large gene deletion strategy (obtained by step1). +% The first column is the list of genes. +% The second column is a 0/1 vector indicating which genes should be deleted. +% 0: indicates genes to be deleted. +% 1: indecates genes to be remained. +% +% OUTPUTS: +% gvalue: a small gene deletion strategy (obtained by TrimGdel). +% The first column is the list of genes. +% The second column is a 0/1 vector indicating which genes should be deleted. +% 0: indicates genes to be deleted. +% 1: indecates genes to be remained. +% GR: the maximum growth rate when the obtained gene deletion +% strategy represented by gvalue is applied. +% PR: the minimum production rate of the target metabolite under +% the maximization of the growth rate when the obtained gene deletion +% strategy represented by gvalue is applied. +% size1: the number of gene deletions after Step1. +% size2: the number of gene deletions after Step2. +% size3: the number of gene deletions after Step3. +% +% .. Author: - Takeyuki Tamura, Mar 06, 2025 +% + +sss = sprintf('step2and3.mat'); +[model, targetRID, extype] = modelSetting(model, targetMet); + +m = size(model.mets, 1); +n = size(model.rxns, 1); +g = size(model.genes, 1); +gid = find(model.c); +pid = targetRID; +model2 = model; + +[grRules0] = calculateGR(model, givenGvalue); +lb2 = model.lb; +ub2 = model.ub; + +for i=1:n + if grRules0{i, 4} == 0 + lb2(i) = 0; + ub2(i) = 0; + end +end + +gm0.obj = -model.c; +gm0.A = sparse(model.S); +gm0.rhs = zeros(m, 1); +gm0.modelsense = 'Min'; +gm0.sense = repmat('=', 1, m); +gm0.lb = lb2; +gm0.ub = ub2; +% derives the initial maximum GR. +opt0 = gurobi(gm0); + +GR0 = -opt0.objval; +lb2(gid) = GR0; +ub2(gid) = GR0; +model2.c(gid) = 0; +model2.c(pid) = 1; + +gm1.obj = -model2.c; +gm1.A = sparse(model.S); +gm1.rhs = zeros(m, 1); +gm1.modelsense = 'Min'; +gm1.sense = repmat('=', 1, m); +gm1.lb = lb2; +gm1.ub = ub2; +% derives the initial minimum PR under the GR maximization. +opt1 = gurobi(gm1); + +GRLB = opt1.x(gid); +PRLB = opt1.x(pid); +[term, ng, nt, nr, nko, reactionKO, reactionKO2term] = readGeneRules(model); +[f, A, b, Aeq, beq, lb, ub, xname] = geneReactionMILP(model, term, ng, nt, nr, nko); + +lp.Aeq = Aeq; +lp.beq = [zeros(size(lp.Aeq, 1), 1)]; +j = 1; +for i=1:size(model.grRules, 1) + if isempty(model.grRules{i, :}) == 0 + ind(1,j) = i; + j = j+1; + end +end +z1 = -diag(model.ub); +z2 = diag(model.lb); +z3 = eye(n); +lp.A = A; +lp.b = b; +lp.lb = lb; +lp.ub = ub; + +for i=1:ng + if givenGvalue{i, 2} == 1 + lp.lb(i, 1) = 1; + lp.ub(i, 1) = 1; + end +end + +[grRules0] = calculateGR(model, givenGvalue); + +j = 1; +for i=1:n + if isempty(model.grRules{i, 1}) == 0 + lp.lb(ng+nt+j, 1) = grRules0{i, 4}; + lp.ub(ng+nt+j, 1) = grRules0{i, 4}; + j = j + 1; + end +end + +lp.f = [-ones(ng, 1); zeros(nt, 1); zeros(nko, 1)]; +for i=1:n + s2 = repelem('B', ng+nt+nko); + lp.ctype = sprintf('%s%s', s2); +end + +gm.obj = lp.f; +gm.A = sparse([lp.A; lp.Aeq]); +gm.rhs = [lp.b; lp.beq]; +gm.modelsense = 'Min'; +gm.sense = horzcat(repmat('<', 1, size(lp.A,1)), repmat('=', 1, size(lp.Aeq, 1))); +gm.lb = lp.lb; +gm.ub = lp.ub; +gm.vtype = lp.ctype; +% MILP for Step 2 +opt = gurobi(gm); + +gvalue = givenGvalue; +if strcmp(opt.status, 'OPTIMAL') + for i=1:ng + vg(i, 1) = opt.x(i); + gvalue{i, 2} = opt.x(i); + end + for i=1:nt + vt(i, 1) = opt.x(ng+i); + end + for i=1:nko + vko(i, 1) = opt.x(ng+nt+i); + end +end +gvalue0 = gvalue; + +trimmed = 1; +model2 = model; +grprlist(1, 1) = opt1.x(gid); +grprlist(1, 2) = opt1.x(pid); +grprlist(1,3) = opt1.x(gid); +grprlist(1,4) = opt1.x(pid); +k = 2; + +while trimmed == 1; + trimmed = 0; + for i=1:ng + i + if gvalue{i, 2} == 0 + gvalue{i, 2} = 1; + [grRules2] = calculateGR(model, gvalue); + lb2 = model.lb; + ub2 = model.ub; + + for j=1:n + if grRules2{j, 4} == 0 + lb2(j) = 0; + ub2(j) = 0; + end + end + + gm2.obj = -model.c; + gm2.A = sparse(model.S); + gm2.rhs = zeros(m, 1); + gm2.modelsense = 'Min'; + gm2.sense = repmat('=', 1, m); + gm2.lb = lb2; + gm2.ub = ub2; + % evaluate the maximum GR when a gene deletion is trimmed. + opt2 = gurobi(gm2); + + + grprlist(k, 1) = opt2.x(gid); + grprlist(k, 2) = opt2.x(pid); + GR2 = -opt2.objval; + lb2(gid) = GR2; + ub2(gid) = GR2; + model2.c(gid) = 0; + model2.c(pid) = 1; + + gm3.obj = model2.c; + gm3.A = sparse(model.S); + gm3.rhs = zeros(m, 1); + gm3.modelsense = 'Min'; + gm3.sense = repmat('=', 1, m); + gm3.lb = lb2; + gm3.ub = ub2; + % evaluate the minimum PR under the GR maximization when a gene + % is trimmed. + opt3 = gurobi(gm3); + + + grprlist(k, 3) = opt3.x(gid); + grprlist(k, 4) = opt3.x(pid); + if ((opt3.x(gid) < 0.999 * GRLB) || (opt3.x(pid) < 0.999*PRLB)) + gvalue{i, 2} = 0; + grprlist(k, :) = grprlist(k-1, :); + + else + trimmed = 1; + + end + GR = grprlist(k, 3); + PR = grprlist(k, 4); + k = k+1; + end + end +end + +gvalueList = horzcat(givenGvalue(:, 2), gvalue0(:, 2), gvalue(:, 2)); +size1 = size(find(cell2mat(givenGvalue(:, 2)) == 0), 1); +size2 = size(find(cell2mat(gvalue0(:, 2)) == 0), 1); +size3 = size(find(cell2mat(gvalue(:, 2)) == 0), 1); + +return; +end + diff --git a/src/design/TrimGdel/testTrimGdel.m b/src/design/TrimGdel/testTrimGdel.m new file mode 100644 index 0000000000..9f686222b7 --- /dev/null +++ b/src/design/TrimGdel/testTrimGdel.m @@ -0,0 +1,73 @@ +function [] = testTrimGdel() +% +% The COBRAToolbox:testTrimGdel.m +% +% This code is to check TrimGdel appropriately works for all target metabolites in e_coli_core. +% +% USAGE: +% +% [] = testTrimGdel() +% +% .. Author: - Takeyuki Tamura, Mar 06, 2025 +% + +initCobraToolbox; +solvers = prepareTest('needsLP', true, 'needsMILP', true, 'useSolversIFAvailable', {'gurobi'}); + +try +assert(strcmp(solvers.LP, 'gurobi')) +catch + display('The MILP solver mulst be gurobi in this version.') + display('The CPLEX version is available on https://github.com/MetNetComp/TrimGdel') + return; +end +try +assert(strcmp(solvers.MILP, 'gurobi')) +catch + display('The MILP solver mulst be gurobi in this version.') + display('The CPLEX version is available on https://github.com/MetNetComp/TrimGdel') + return; +end + + + +load('e_coli_core.mat'); +model = e_coli_core; + + +m = size(model.mets, 1); +for i=1:m + [gvalue, GR, PR, size1, size2, size3, success] = TrimGdel(model, model.mets{i} ,3, 0.1, 0.1); + if success == 1 + try + assert(GR >= 0.001); + catch + display('The GR value is not appropriate.') + return; + end + try + assert(PR >= 0.001); + catch + display('The PR value is not appropriate.') + return; + end + table(i, 1) = success; + table(i, 2) = GR; + table(i, 3) = PR; + table(i, 4) = size1; + table(i, 5) = size2; + table(i, 6) = size3; + + end +end +try + assert(sum(table(:,1)) > 30) +catch + display('It seems that TrimGdel cannot perform optimally in this environment.') + return; +end + +display('The test was successful.') + +end + diff --git a/test/verifiedTests/design/testTrimGdel.m b/test/verifiedTests/design/testTrimGdel.m new file mode 100644 index 0000000000..9f686222b7 --- /dev/null +++ b/test/verifiedTests/design/testTrimGdel.m @@ -0,0 +1,73 @@ +function [] = testTrimGdel() +% +% The COBRAToolbox:testTrimGdel.m +% +% This code is to check TrimGdel appropriately works for all target metabolites in e_coli_core. +% +% USAGE: +% +% [] = testTrimGdel() +% +% .. Author: - Takeyuki Tamura, Mar 06, 2025 +% + +initCobraToolbox; +solvers = prepareTest('needsLP', true, 'needsMILP', true, 'useSolversIFAvailable', {'gurobi'}); + +try +assert(strcmp(solvers.LP, 'gurobi')) +catch + display('The MILP solver mulst be gurobi in this version.') + display('The CPLEX version is available on https://github.com/MetNetComp/TrimGdel') + return; +end +try +assert(strcmp(solvers.MILP, 'gurobi')) +catch + display('The MILP solver mulst be gurobi in this version.') + display('The CPLEX version is available on https://github.com/MetNetComp/TrimGdel') + return; +end + + + +load('e_coli_core.mat'); +model = e_coli_core; + + +m = size(model.mets, 1); +for i=1:m + [gvalue, GR, PR, size1, size2, size3, success] = TrimGdel(model, model.mets{i} ,3, 0.1, 0.1); + if success == 1 + try + assert(GR >= 0.001); + catch + display('The GR value is not appropriate.') + return; + end + try + assert(PR >= 0.001); + catch + display('The PR value is not appropriate.') + return; + end + table(i, 1) = success; + table(i, 2) = GR; + table(i, 3) = PR; + table(i, 4) = size1; + table(i, 5) = size2; + table(i, 6) = size3; + + end +end +try + assert(sum(table(:,1)) > 30) +catch + display('It seems that TrimGdel cannot perform optimally in this environment.') + return; +end + +display('The test was successful.') + +end +