From 648b43f254f3c161dec88f9d2c879f638ae6029a Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Fri, 17 Jul 2026 22:43:51 +0200 Subject: [PATCH 1/6] refactor: rearchitect assignCompartments as place -> repair -> certify The previous version coupled a flux model into the placement MILP (Big-M flux gating, a growth-floor row) and set exitFlag=1 from the MILP status alone, materialising outModel without ever solving it. At genome scale that can certify a placement the built model cannot actually grow. Placement is now a flux-free MILP over only the reaction/compartment and gene/compartment binaries, maximising the localization score, so it has nothing for a tolerance-rounded binary to leak into. A solver-free confinement fixpoint keeps the placement connected: a non-transportable metabolite split across compartments forces its reactions to co-locate, while transportable splits get passive transports. The placement is then materialised and that exact model is solved; exitFlag reflects whether it reached the growth floor, and a new fifth output reports certified/status and the real per-medium growths, so a failed certification is honest rather than hidden. Adds a growthConditions argument (certify on extra media). transportCost is now accepted but unused: placement follows the score, and transports are a structural consequence of repair. The three existing tests still pass; new tests cover honest certification reporting and the confinement co-location. --- localization/assignCompartments.m | 511 ++++++++++++------- testing/function_tests/tAssignCompartments.m | 70 +++ 2 files changed, 391 insertions(+), 190 deletions(-) diff --git a/localization/assignCompartments.m b/localization/assignCompartments.m index 37cb08de..115f3cfb 100644 --- a/localization/assignCompartments.m +++ b/localization/assignCompartments.m @@ -1,70 +1,99 @@ -function [outModel, placement, addedTransports, exitFlag] = assignCompartments(model, GSS, reactionsToRelocate, varargin) -% assignCompartments Assign reactions to compartments by a functionality-constrained MILP. +function [outModel, placement, addedTransports, exitFlag, report] = assignCompartments(model, GSS, reactionsToRelocate, varargin) +% assignCompartments Assign reactions to compartments, certified by growth. % -% Deterministic alternative to predictLocalization: a single MILP places the requested reactions -% into compartments to agree with localization scores (GSS) WHILE keeping the model's objective -% (biomass) producible. Because functionality is a hard constraint, a reaction is placed against -% its own top score whenever the network needs it there (pathway coherence emerges from requiring -% flux, not from a heuristic). Mono-localization: each reaction is placed in exactly one -% compartment; a gene still spans compartments when it catalyses reactions placed in different ones. +% Deterministic alternative to predictLocalization. Places the requested +% reactions into the compartments named by localization scores (GSS) and then +% verifies, by an actual FBA on the built model, that the objective (biomass) +% is still producible. The work is done in three phases: % -% The model is merged to a single compartment first (like predictLocalization); the MILP then -% re-distributes the relocatable reactions across the compartments named in GSS, adding passive -% transports (defaultCompartment <-> c) where needed. +% 1. Place - a flux-free MILP that maximises the localization score. +% Its only variables are placement binaries (reaction -> +% compartment, gene -> compartment); it has no flux variables +% and no growth constraint, so a placement can never harvest a +% compartment's score through leaked flux. +% 2. Repair - a solver-free fixpoint that keeps the placement connected: a +% non-transportable metabolite split across compartments forces +% its reactions to co-locate; transportable splits get passive +% transports through the default compartment. +% 3. Certify - the placement is materialised into a compartmentalised model +% and that exact model is solved. exitFlag reflects whether it +% reached the growth floor, so the certificate is the model +% that is returned, never a placement the model cannot grow. +% +% Mono-localization: each reaction is placed in exactly one compartment; a +% gene still spans compartments when it catalyses reactions placed in +% different ones. % % Parameters % ---------- % model : struct -% a RAVEN model with an objective (model.c) and grRules. Multiple compartments are merged. +% a RAVEN model with an objective (model.c) and grRules. Multiple +% compartments are merged before placement. % GSS : struct -% gene scoring structure (genes, compartments, scores) as from parseScores. GSS.compartments -% are the target compartment labels and must include defaultCompartment. +% gene scoring structure (genes, compartments, scores) as from parseScores. +% GSS.compartments are the target compartment labels and must include +% defaultCompartment. % reactionsToRelocate : cell -% reaction ids to (re)place. Boundary reactions and the objective reaction are always pinned. +% reaction ids to (re)place. Boundary reactions and the objective reaction +% are always pinned. % % Name-Value Arguments % -------------------- % defaultCompartment : char -% compartment that transports route through (usually cytosol); must be in GSS.compartments. -% transportCost : double (default 0.5) -% cost per added inter-compartment transport. +% compartment that transports route through (usually cytosol); must be in +% GSS.compartments. % multiCompartmentPenalty : double (default 0.5) -% cost per extra compartment a gene ends up in. +% score cost per extra compartment a gene ends up in. % minGrowth : double (default [] = 10%% of the unconstrained optimum) -% required objective flux. +% required objective flux for certification. % transportable : cell (default [] = all) -% metabolite ids (merged-model names) that may receive transports. Restricting it forces -% functionality-driven placement of reactions whose substrates are confined. -% bigM : double (default 1000) -% Big-M for flux gating. +% metabolite ids (merged-model names) that may receive transports. +% Restricting it forces functionality-driven placement of reactions whose +% substrates are confined. +% growthConditions : struct array (default []) +% extra media the placement must also grow on, each with fields: name, +% medium (struct of exchangeRxnId -> max uptake), and minGrowth. +% maxRounds : double (default 8) +% maximum place/repair rounds before returning the best placement found. +% transportCost : double (default 0.5) +% accepted for backward compatibility but not used: placement is driven by +% the score alone, and transports are a structural consequence of repair. % verbose : logical (default true) % % Returns % ------- % outModel : struct -% the compartmentalised model (objective still producible). +% the compartmentalised model. % placement : struct -% .rxns and .compartment (the assigned compartment id per relocated reaction). +% .rxns and .compartment (the assigned compartment id per relocated +% reaction). % addedTransports : struct -% .mets and .compartment for the transports the MILP added. +% .mets and .compartment for the transports repair added. % exitFlag : double -% 1 = optimal, -1 = infeasible/failed. +% 1 = certified (the built model reaches the growth floor on every +% medium), -1 = could not place or the built model did not certify. +% report : struct +% .certified (logical), .status (char), and .growths (struct of +% medium -> objective flux), so a failed certification reports the real +% growth rather than hiding it. % % See also % -------- % predictLocalization, parseScores, gapFillMILP p = parseRAVENargs(varargin, {'defaultCompartment',[]; 'transportCost',0.5; ... - 'multiCompartmentPenalty',0.5; 'minGrowth',[]; 'transportable',[]; 'bigM',1000; 'verbose',true}); + 'multiCompartmentPenalty',0.5; 'minGrowth',[]; 'transportable',[]; ... + 'growthConditions',[]; 'maxRounds',8; 'verbose',true}); defaultCompartment = char(p.defaultCompartment); -transportCost = p.transportCost; multiPen = p.multiCompartmentPenalty; minGrowth = p.minGrowth; -M = p.bigM; +growthConditions = p.growthConditions; +maxRounds = p.maxRounds; verbose = p.verbose; outModel = model; placement = struct('rxns',{{}},'compartment',{{}}); addedTransports = struct('mets',{{}},'compartment',{{}}); exitFlag = -1; +report = struct('certified',false,'status','not_solved','growths',struct()); if all(model.c == 0) error('RAVEN:badInput','model has no objective (set model.c).'); @@ -74,21 +103,19 @@ error('RAVEN:badInput','defaultCompartment ''%s'' not in GSS.compartments.', defaultCompartment); end -% Merge to a single compartment so every metabolite has one identity; relocations and transports -% are then expressed relative to defaultCompartment. +% Merge to a single compartment so every metabolite has one identity; +% placement and transports are then expressed relative to defaultCompartment. if numel(model.comps) > 1 model = mergeCompartments(model, true, true); end -nMet = numel(model.mets); biomassIdx = find(model.c ~= 0, 1); -% Determine min growth from the merged (single-compartment) model if isempty(minGrowth) sol = solveLP(model); if isempty(sol.f) || sol.f <= 0 error('RAVEN:badInput','merged model does not grow; pass minGrowth.'); end - minGrowth = 0.1 * sol.f; + minGrowth = 0.1 * abs(sol.f); end % ---- scope: movable vs pinned ---- @@ -98,172 +125,142 @@ movableMask(biomassIdx) = false; movIdx = find(movableMask); pinIdx = find(~movableMask); -nMov = numel(movIdx); nPin = numel(pinIdx); nC = numel(comps); +nMov = numel(movIdx); nC = numel(comps); defC = find(strcmp(comps, defaultCompartment)); % Genes in scope: those on movable reactions that have a row in GSS [gInGSS, gssRow] = ismember(model.genes, GSS.genes); geneOnMov = any(model.rxnGeneMat(movIdx, :) ~= 0, 1)'; -scopeMask = gInGSS & geneOnMov; -geneIdx = find(scopeMask); +geneIdx = find(gInGSS & geneOnMov); nGene = numel(geneIdx); -score = zeros(nGene, nC); % gene x compartment score, columns aligned to comps +score = zeros(nGene, nC); for gi = 1:nGene score(gi, :) = GSS.scores(gssRow(geneIdx(gi)), :); end -% Transportable metabolites (touched by movable reactions) +% Transportable metabolites (a base metabolite may be moved between +% compartments by a passive transport during repair). movMetMask = any(model.S(:, movIdx) ~= 0, 2); if isnumeric(p.transportable) && isempty(p.transportable) - transpMet = find(movMetMask); % default ([]): all movable metabolites transportable + transpMet = find(movMetMask); % default: all else - transpMet = find(movMetMask & ismember(model.mets, p.transportable)); % a (possibly empty) set -end -% transport variables only for non-default compartments -trPairs = []; % [metRow, compCol] -for ci = 1:nC - if ci == defC; continue; end - trPairs = [trPairs; [transpMet, repmat(ci, numel(transpMet), 1)]]; %#ok + transpMet = find(movMetMask & ismember(model.mets, p.transportable)); end -nTr = size(trPairs, 1); +isTransp = false(numel(model.mets),1); isTransp(transpMet) = true; if verbose - fprintf('assignCompartments: %d movable, %d pinned reactions, %d genes, %d compartments, %d transports.\n', ... - nMov, nPin, nGene, nC, nTr); -end - -% ---- variable layout (columns) ---- -% fmove[mov,c] : nMov*nC continuous -% fpin[pin] : nPin continuous -% ftr[tr] : nTr continuous -% x[mov,c] : nMov*nC binary -% y[gene,c] : nGene*nC binary -% t[tr] : nTr binary -oFmove = 0; nFmove = nMov*nC; -oFpin = oFmove+nFmove; nFpin = nPin; -oFtr = oFpin+nFpin; % nTr -oX = oFtr+nTr; nX = nMov*nC; -oY = oX+nX; nY = nGene*nC; -oT = oY+nY; % nTr -nVar = oT+nTr; -fmoveCol = @(mi,ci) oFmove + (mi-1)*nC + ci; -xCol = @(mi,ci) oX + (mi-1)*nC + ci; -yCol = @(gi,ci) oY + (gi-1)*nC + ci; - -% ---- bounds ---- -lb = zeros(nVar,1); ub = zeros(nVar,1); -for k = 1:nMov - r = movIdx(k); - for ci = 1:nC - lb(fmoveCol(k,ci)) = min(model.lb(r),0); ub(fmoveCol(k,ci)) = max(model.ub(r),0); - end + fprintf('assignCompartments: %d movable, %d pinned reactions, %d genes, %d compartments.\n', ... + nMov, numel(pinIdx), nGene, nC); end -lb(oFpin+(1:nPin)) = model.lb(pinIdx); ub(oFpin+(1:nPin)) = model.ub(pinIdx); -lb(oFtr+(1:nTr)) = -M; ub(oFtr+(1:nTr)) = M; -lb(oX+(1:nX)) = 0; ub(oX+(1:nX)) = 1; -lb(oY+(1:nY)) = 0; ub(oY+(1:nY)) = 1; -lb(oT+(1:nTr)) = 0; ub(oT+(1:nTr)) = 1; - -% ---- node balance: S over (met, compartment) ---- -% node row index = (metRow-1)*nC + ci -nNode = nMet*nC; -nodeRow = @(mr,ci) (mr-1)*nC + ci; -ri = []; ci_ = []; vv = []; -% movable: metabolites placed in c -for k = 1:nMov - r = movIdx(k); - mrows = find(model.S(:,r) ~= 0); - for ci = 1:nC - for mm = mrows' - ri(end+1)=nodeRow(mm,ci); ci_(end+1)=fmoveCol(k,ci); vv(end+1)=model.S(mm,r); %#ok + +% ---- place / repair rounds ---- +% forced(mi) = compartment index a movable reaction is pinned to (0 = free); +% groups is a cell array of movable-index vectors that must share a +% compartment. Both grow monotonically until the confinement fixpoint holds. +forced = zeros(nMov,1); +groups = {}; +seen = {}; +placeIdx = []; +for round = 1:maxRounds + placeIdx = i_placementMaster(model, movIdx, geneIdx, score, multiPen, ... + nC, forced, groups, verbose); + if isempty(placeIdx) + report.status = 'infeasible'; + if verbose; fprintf('assignCompartments: placement MILP infeasible.\n'); end + return; + end + [newForced, newGroups] = i_diagnoseConfinement(model, movIdx, pinIdx, ... + placeIdx, defC, isTransp); + % Apply only genuinely new tightenings; stop at the fixpoint. + changed = false; + for k = 1:nMov + if newForced(k) ~= 0 && forced(k) == 0 + forced(k) = newForced(k); changed = true; end end -end -% pinned: metabolites in defaultCompartment -for k = 1:nPin - r = pinIdx(k); - mrows = find(model.S(:,r) ~= 0); - for mm = mrows' - ri(end+1)=nodeRow(mm,defC); ci_(end+1)=oFpin+k; vv(end+1)=model.S(mm,r); %#ok + for gi = 1:numel(newGroups) + if ~i_hasGroup(groups, newGroups{gi}) + groups{end+1} = newGroups{gi}; changed = true; %#ok + end end + sig = i_signature(forced, groups); + if ~changed || any(strcmp(seen, sig)) + break; + end + seen{end+1} = sig; %#ok end -% transports: -1 at (met,default), +1 at (met,c) -for k = 1:nTr - mm = trPairs(k,1); cc = trPairs(k,2); - ri(end+1)=nodeRow(mm,defC); ci_(end+1)=oFtr+k; vv(end+1)=-1; %#ok - ri(end+1)=nodeRow(mm,cc); ci_(end+1)=oFtr+k; vv(end+1)=+1; %#ok -end -A_node = sparse(ri, ci_, vv, nNode, nVar); -keepNode = any(A_node ~= 0, 2); % drop empty nodes -A_node = A_node(keepNode, :); -nNodeKept = size(A_node,1); -% ---- flux gating: lb*x <= fmove <= ub*x ---- -giR=[]; giC=[]; giV=[]; row=0; gUb=[]; gLb=[]; -for k=1:nMov - r=movIdx(k); - for ci=1:nC - row=row+1; % ub: fmove - ub*x <= 0 - giR(end+1)=row; giC(end+1)=fmoveCol(k,ci); giV(end+1)=1; %#ok - giR(end+1)=row; giC(end+1)=xCol(k,ci); giV(end+1)=-max(model.ub(r),0); %#ok - gUb(end+1)=0; %#ok +% ---- placement + transports ---- +placeRxns = model.rxns(movIdx); +placeComp = comps(placeIdx); +placement.rxns = placeRxns(:); placement.compartment = placeComp(:); + +[trMet, trComp] = i_splitTransports(model, movIdx, pinIdx, placeIdx, defC, isTransp, nC); +addedTransports.mets = model.mets(trMet); +addedTransports.compartment = comps(trComp); + +% ---- materialise + certify ---- +outModel = i_applyAssignment(model, placement, addedTransports, comps, defaultCompartment); +[certified, growths] = i_certify(outModel, biomassIdx, model.rxns{biomassIdx}, minGrowth, growthConditions); +report.certified = certified; +report.growths = growths; +if certified + exitFlag = 1; report.status = 'certified'; +else + exitFlag = -1; report.status = 'uncertified'; + if verbose + fprintf('assignCompartments: placement did not certify (growth %.4g < %.4g).\n', ... + growths.primary, minGrowth); end end -A_gub = sparse(giR,giC,giV,row,nVar); b_gub=zeros(row,1); -giR=[]; giC=[]; giV=[]; row=0; -for k=1:nMov - r=movIdx(k); - for ci=1:nC - row=row+1; % lb: fmove - lb*x >= 0 - giR(end+1)=row; giC(end+1)=fmoveCol(k,ci); giV(end+1)=1; %#ok - giR(end+1)=row; giC(end+1)=xCol(k,ci); giV(end+1)=-min(model.lb(r),0); %#ok - end end -A_glb = sparse(giR,giC,giV,row,nVar); b_glb=zeros(row,1); -% ---- transport gating: -M*t <= ftr <= M*t ---- -tR=[];tC=[];tV=[]; -for k=1:nTr - tR(end+1)=k; tC(end+1)=oFtr+k; tV(end+1)=1; tR(end+1)=k; tC(end+1)=oT+k; tV(end+1)=-M; %#ok -end -A_tub = sparse(tR,tC,tV,nTr,nVar); % ftr - M*t <= 0 -tR=[];tC=[];tV=[]; -for k=1:nTr - tR(end+1)=k; tC(end+1)=oFtr+k; tV(end+1)=1; tR(end+1)=k; tC(end+1)=oT+k; tV(end+1)=M; %#ok -end -A_tlb = sparse(tR,tC,tV,nTr,nVar); % ftr + M*t >= 0 +% ------------------------------------------------------------- placement MILP +function placeIdx = i_placementMaster(model, movIdx, geneIdx, score, multiPen, nC, forced, groups, verbose) +% Flux-free score MILP. Variables: x[mov,c] and y[gene,c], both binary. +% Returns the chosen compartment index per movable reaction, or [] if +% infeasible. +nMov = numel(movIdx); nGene = numel(geneIdx); +oX = 0; nX = nMov*nC; +oY = oX+nX; nY = nGene*nC; +nVar = oY+nY; +xCol = @(mi,ci) oX + (mi-1)*nC + ci; +yCol = @(gi,ci) oY + (gi-1)*nC + ci; -% ---- placement: sum_c x[mov,c] = 1 ---- +% placement: sum_c x[r,c] = 1 pR=[];pC=[];pV=[]; for k=1:nMov for ci=1:nC; pR(end+1)=k; pC(end+1)=xCol(k,ci); pV(end+1)=1; end %#ok end A_place=sparse(pR,pC,pV,nMov,nVar); b_place=ones(nMov,1); -% ---- gene coupling x[r,c] <= y[g,c]; gene-has y<=sum x; gene1 sum_c y>=1 ---- +% coupling: x[r,c] - y[g,c] <= 0 for each gene g on reaction r cR=[];cC=[];cV=[];row=0; for k=1:nMov gOn=find(model.rxnGeneMat(movIdx(k),:)~=0); for g=gOn gi=find(geneIdx==g,1); if isempty(gi); continue; end for ci=1:nC - row=row+1; cR(end+1)=row; cC(end+1)=xCol(k,ci); cV(end+1)=1; + row=row+1; + cR(end+1)=row; cC(end+1)=xCol(k,ci); cV(end+1)=1; %#ok cR(end+1)=row; cC(end+1)=yCol(gi,ci); cV(end+1)=-1; %#ok end end end A_couple=sparse(cR,cC,cV,row,nVar); b_couple=zeros(row,1); -% gene-has: y[g,c] - sum_{r of g} x[r,c] <= 0 + +% has: y[g,c] - sum_{r of g} x[r,c] <= 0 hR=[];hC=[];hV=[];row=0; for gi=1:nGene g=geneIdx(gi); - rOfG=find(model.rxnGeneMat(movIdx,g)~=0)'; % indices into movIdx + rOfG=find(model.rxnGeneMat(movIdx,g)~=0)'; for ci=1:nC - row=row+1; hR(end+1)=row; hC(end+1)=yCol(gi,ci); hV(end+1)=1; + row=row+1; hR(end+1)=row; hC(end+1)=yCol(gi,ci); hV(end+1)=1; %#ok for k=rOfG; hR(end+1)=row; hC(end+1)=xCol(k,ci); hV(end+1)=-1; end %#ok end end A_has=sparse(hR,hC,hV,row,nVar); b_has=zeros(row,1); + % gene1: sum_c y[g,c] >= 1 g1R=[];g1C=[];g1V=[]; for gi=1:nGene @@ -271,78 +268,212 @@ end A_gene1=sparse(g1R,g1C,g1V,nGene,nVar); b_gene1=ones(nGene,1); -% ---- growth: fpin[biomass] >= minGrowth ---- -pinBio = find(pinIdx==biomassIdx,1); -A_grow = sparse(1, oFpin+pinBio, 1, 1, nVar); b_grow = minGrowth; +% forced pins: x[r,c_force] = 1 +fR=[];fC=[];nForce=0; +for k=1:nMov + if forced(k)~=0 + nForce=nForce+1; fR(end+1)=nForce; fC(end+1)=xCol(k,forced(k)); %#ok + end +end +A_force=sparse(fR,fC,ones(1,nForce),nForce,nVar); b_force=ones(nForce,1); + +% co-location groups: x[a,c] - x[b,c] = 0 for consecutive members +coR=[];coC=[];coV=[];row=0; +for gi=1:numel(groups) + mem=groups{gi}; + for j=1:numel(mem)-1 + a=mem(j); b=mem(j+1); + for ci=1:nC + row=row+1; + coR(end+1)=row; coC(end+1)=xCol(a,ci); coV(end+1)=1; %#ok + coR(end+1)=row; coC(end+1)=xCol(b,ci); coV(end+1)=-1; %#ok + end + end +end +A_colo=sparse(coR,coC,coV,row,nVar); b_colo=zeros(row,1); -% ---- objective: max sum score*y - multiPen*sum y - transportCost*sum t ---- +% objective: max sum score*y - multiPen*sum y c = zeros(nVar,1); for gi=1:nGene - for ci=1:nC; c(yCol(gi,ci)) = c(yCol(gi,ci)) + score(gi,ci) - multiPen; end + for ci=1:nC; c(yCol(gi,ci)) = score(gi,ci) - multiPen; end end -if isscalar(transportCost); tcost = repmat(transportCost,nTr,1); else; tcost = transportCost(trPairs(:,1)); end -c(oT+(1:nTr)) = c(oT+(1:nTr)) - tcost(:); -% ---- assemble ---- -prob.A = [A_node; A_gub; A_glb; A_tub; A_tlb; A_place; A_couple; A_has; A_gene1; A_grow]; +prob.A = [A_place; A_couple; A_has; A_gene1; A_force; A_colo]; prob.a = prob.A; -prob.b = [zeros(nNodeKept,1); b_gub; b_glb; zeros(nTr,1); zeros(nTr,1); b_place; b_couple; b_has; b_gene1; b_grow]; -prob.csense = [repmat('E',1,nNodeKept), repmat('L',1,numel(b_gub)), repmat('G',1,numel(b_glb)), ... - repmat('L',1,nTr), repmat('G',1,nTr), repmat('E',1,nMov), ... - repmat('L',1,numel(b_couple)), repmat('L',1,numel(b_has)), repmat('G',1,nGene), 'G']; -prob.c = -c; % optimizeProb minimises; we maximise c +prob.b = [b_place; b_couple; b_has; b_gene1; b_force; b_colo]; +prob.csense = [repmat('E',1,nMov), repmat('L',1,numel(b_couple)), ... + repmat('L',1,numel(b_has)), repmat('G',1,nGene), ... + repmat('E',1,nForce), repmat('E',1,numel(b_colo))]; +prob.c = -c; % optimizeProb minimises; we maximise c prob.osense = 1; -prob.lb = lb; prob.ub = ub; -prob.vartype = repmat('C',1,nVar); -prob.vartype(oX+(1:nX)) = 'B'; -prob.vartype(oY+(1:nY)) = 'B'; -prob.vartype(oT+(1:nTr)) = 'B'; +prob.lb = zeros(nVar,1); prob.ub = ones(nVar,1); +prob.vartype = repmat('B',1,nVar); params.intTol = 1e-9; params.TimeLimit = 1000; sol = optimizeProb(prob, params, verbose); if ~checkSolution(sol) - if verbose; fprintf('assignCompartments: MILP infeasible or failed.\n'); end - return; + placeIdx = []; return; end - -% ---- extract ---- xval = sol.full(oX+(1:nX)); -placeRxns = {}; placeComp = {}; +placeIdx = zeros(nMov,1); for k=1:nMov [~,ci] = max(xval((k-1)*nC + (1:nC))); - placeRxns{end+1} = model.rxns{movIdx(k)}; placeComp{end+1} = comps{ci}; %#ok + placeIdx(k) = ci; +end end -placement.rxns = placeRxns(:); placement.compartment = placeComp(:); -tval = sol.full(oT+(1:nTr)); -sel = find(tval > 0.5); -addedTransports.mets = model.mets(trPairs(sel,1)); -addedTransports.compartment = comps(trPairs(sel,2)); -exitFlag = 1; -outModel = i_applyAssignment(model, placement, addedTransports, comps, defaultCompartment); +% ------------------------------------------------------ confinement diagnosis +function [forced, groups] = i_diagnoseConfinement(model, movIdx, pinIdx, placeIdx, defC, isTransp) +% Solver-free: find non-transportable metabolites split across compartments +% and decide, per split, whether to force the movers to the pinned +% compartment or to co-locate them. Pinned reactions all live in the default +% compartment (the model was merged first). +nMov = numel(movIdx); +forced = zeros(nMov,1); +groups = {}; + +% used{metRow} maps compartment index -> list of movable local indices. +% pinnedHere(metRow) is true if a pinned reaction touches the metabolite. +nMet = numel(model.mets); +usedComp = cell(nMet,1); +usedMov = cell(nMet,1); +pinnedHere = false(nMet,1); +for k=1:nMov + mrows = find(model.S(:,movIdx(k))~=0)'; + ci = placeIdx(k); + for mm=mrows + if isTransp(mm); continue; end % transportable: handled by a transport + usedComp{mm}(end+1) = ci; + usedMov{mm}(end+1) = k; + end +end +for k=1:numel(pinIdx) + mrows = find(model.S(:,pinIdx(k))~=0)'; + for mm=mrows + if isTransp(mm); continue; end + usedComp{mm}(end+1) = defC; + pinnedHere(mm) = true; + end +end + +for mm=1:nMet + cs = unique(usedComp{mm}); + if numel(cs) <= 1; continue; end % lives in one compartment: fine + movers = unique(usedMov{mm}); + if pinnedHere(mm) + % A pinned reaction anchors this metabolite in the default + % compartment; every movable toucher must join it there. + for k=movers; forced(k) = defC; end + else + % Only movable reactions touch it: co-locate them and let the score + % objective choose the shared compartment. + groups{end+1} = movers(:)'; %#ok + end +end +end + +% ----------------------------------------------------------- split transports +function [trMet, trComp] = i_splitTransports(model, movIdx, pinIdx, placeIdx, defC, isTransp, nC) +% For each transportable base metabolite placed in more than one compartment, +% add a transport between the default compartment and each non-default one. +nMet = numel(model.mets); +compsUsed = cell(nMet,1); +for k=1:numel(movIdx) + mrows = find(model.S(:,movIdx(k))~=0)'; + for mm=mrows; compsUsed{mm}(end+1) = placeIdx(k); end +end +for k=1:numel(pinIdx) + mrows = find(model.S(:,pinIdx(k))~=0)'; + for mm=mrows; compsUsed{mm}(end+1) = defC; end +end +trMet=[]; trComp=[]; +for mm=1:nMet + if ~isTransp(mm); continue; end + cs = unique(compsUsed{mm}); + if numel(cs) <= 1; continue; end + for ci=cs + if ci==defC; continue; end + trMet(end+1)=mm; trComp(end+1)=ci; %#ok + end +end +trMet=trMet(:); trComp=trComp(:); +end + +% ------------------------------------------------------------- certification +function [certified, growths] = i_certify(outModel, biomassIdx, biomassId, minGrowth, growthConditions) +% Solve the materialised model and confirm it reaches minGrowth on the primary +% medium and on every extra growth condition. +tol = 1e-9; +growths = struct(); +bIdx = find(strcmp(outModel.rxns, biomassId), 1); +if isempty(bIdx); bIdx = biomassIdx; end +sol = solveLP(outModel); +primary = 0; +if ~isempty(sol.f); primary = abs(sol.f); end +growths.primary = primary; +certified = primary >= minGrowth - tol; + +if ~isempty(growthConditions) && isstruct(growthConditions) + for i=1:numel(growthConditions) + gc = growthConditions(i); + m2 = i_applyMedium(outModel, gc.medium); + s2 = solveLP(m2); + g = 0; if ~isempty(s2.f); g = abs(s2.f); end + growths.(matlab.lang.makeValidName(gc.name)) = g; + certified = certified && (g >= gc.minGrowth - tol); + end +end +end + +function model = i_applyMedium(model, medium) +% Close all uptake, then open only the listed exchanges. medium is a struct +% of exchangeRxnId -> max uptake. +[~, exchIdx] = getExchangeRxns(model); +model.lb(exchIdx(model.lb(exchIdx) < 0)) = 0; +if ~isempty(medium) && isstruct(medium) + f = fieldnames(medium); + for i=1:numel(f) + r = find(strcmp(model.rxns, f{i}), 1); + if ~isempty(r); model.lb(r) = -abs(medium.(f{i})); end + end +end +end + +% ------------------------------------------------------------------ helpers +function tf = i_hasGroup(groups, g) +tf = false; +g = sort(g(:)'); +for i=1:numel(groups) + if isequal(sort(groups{i}(:)'), g); tf = true; return; end +end +end + +function s = i_signature(forced, groups) +gp = ''; +for i=1:numel(groups); gp = [gp '|' num2str(sort(groups{i}(:)'))]; end %#ok +s = [num2str(forced(:)') '#' gp]; end % ------------------------------------------------------------------- apply function outModel = i_applyAssignment(model, placement, addedTransports, comps, defaultCompartment) -% Build the compartmentalised model: relabel each placed reaction's metabolites into its -% compartment (creating per-compartment metabolite copies) and add the requested transports. +% Build the compartmentalised model: relabel each placed reaction's metabolites +% into its compartment (creating per-compartment metabolite copies) and add the +% requested transports. outModel = model; -% ensure all target compartments exist for ci = 1:numel(comps) if ~ismember(comps{ci}, outModel.comps) outModel.comps{end+1,1} = comps{ci}; if isfield(outModel,'compNames'); outModel.compNames{end+1,1} = comps{ci}; end end end -% metabolite (baseName, compartment) -> index, seeded with existing -metKey = strcat(model.metNames, '###', model.comps(model.metComps)); for k = 1:numel(placement.rxns) r = find(strcmp(outModel.rxns, placement.rxns{k}), 1); comp = placement.compartment{k}; - mrows = find(model.S(:, find(strcmp(model.rxns,placement.rxns{k}),1)) ~= 0); %#ok<*FNDSB> + srcR = find(strcmp(model.rxns,placement.rxns{k}),1); + mrows = find(model.S(:, srcR) ~= 0); %#ok<*FNDSB> for mm = mrows' [outModel, newMet] = i_metInComp(outModel, model, mm, comp); - outModel.S(newMet, r) = model.S(mm, find(strcmp(model.rxns,placement.rxns{k}),1)); + outModel.S(newMet, r) = model.S(mm, srcR); if ~isequal(newMet, mm); outModel.S(mm, r) = 0; end end end @@ -355,7 +486,7 @@ ci = find(strcmp(outModel.comps, comp), 1); name = model.metNames{srcMetRow}; cand = find(strcmp(outModel.metNames, name) & outModel.metComps == ci, 1); -if ~isempty(cand); idx = cand; return; end % reuse existing per-compartment metabolite +if ~isempty(cand); idx = cand; return; end newId = [model.mets{srcMetRow} '_' comp]; outModel.mets{end+1,1} = newId; outModel.metNames{end+1,1} = name; diff --git a/testing/function_tests/tAssignCompartments.m b/testing/function_tests/tAssignCompartments.m index 756352c2..eb59a6a5 100644 --- a/testing/function_tests/tAssignCompartments.m +++ b/testing/function_tests/tAssignCompartments.m @@ -45,6 +45,54 @@ function infeasibleGrowthFloorReported(testCase) testCase.verifyEqual(eFlag, -1); % growth floor unreachable end + function certificationReportsRealGrowth(testCase) + % A certified placement reports certified=true and the achieved + % growth; an unreachable floor reports the real (short) growth + % rather than hiding it behind an optimal-MILP status. + testCase.assumeMILPSolver(); + model = tAssignCompartments.toy(); + GSS = tAssignCompartments.gss(); + evalc(['[~, ~, ~, ok, rep] = assignCompartments(model, GSS, {''r1''}, ' ... + '''defaultCompartment'', ''c'', ''transportable'', {}, ''verbose'', false);']); + testCase.verifyEqual(ok, 1); + testCase.verifyTrue(rep.certified); + testCase.verifyEqual(rep.status, 'certified'); + testCase.verifyGreaterThan(rep.growths.primary, 1e-6); + + evalc(['[~, ~, ~, bad, badRep] = assignCompartments(model, GSS, {''r1''}, ' ... + '''defaultCompartment'', ''c'', ''minGrowth'', 1e6, ''verbose'', false);']); + testCase.verifyEqual(bad, -1); + testCase.verifyFalse(badRep.certified); + testCase.verifyEqual(badRep.status, 'uncertified'); + % the real growth is small and reported, not concealed + testCase.verifyGreaterThan(badRep.growths.primary, 1e-6); + testCase.verifyLessThan(badRep.growths.primary, 1e6); + end + + function confinementColocatesMovableReactions(testCase) + % Two movable reactions share a non-transportable intermediate X. + % Their scores pull them to different compartments, but X cannot + % be transported, so the confinement repair must co-locate them. + testCase.assumeMILPSolver(); + model = tAssignCompartments.chainToy(); + GSS = tAssignCompartments.chainGss(); + % X_c non-transportable, S_c/P_c transportable (so the biomass + % path is not what forces co-location). + evalc(['[~, place, ~, ok] = assignCompartments(model, GSS, {''r1'';''r2''}, ' ... + '''defaultCompartment'', ''c'', ''transportable'', {''S_c'';''P_c''}, ''verbose'', false);']); + testCase.verifyEqual(ok, 1); + c1 = place.compartment{strcmp(place.rxns,'r1')}; + c2 = place.compartment{strcmp(place.rxns,'r2')}; + testCase.verifyEqual(c1, c2); % co-located + + % With X_c transportable too, nothing forces co-location and the + % opposing scores split the two reactions across compartments. + evalc(['[~, sp] = assignCompartments(model, GSS, {''r1'';''r2''}, ' ... + '''defaultCompartment'', ''c'', ''transportable'', {''S_c'';''X_c'';''P_c''}, ''verbose'', false);']); + testCase.verifyNotEqual(sp.compartment{strcmp(sp.rxns,'r1')}, ... + sp.compartment{strcmp(sp.rxns,'r2')}); + end + end methods (Static) @@ -60,5 +108,27 @@ function infeasibleGrowthFloorReported(testCase) function GSS = gss() GSS.genes={'g1'}; GSS.compartments={'c';'m'}; GSS.scores=[0.4 0.9]; end + + function model = chainToy() + % EX_S -> S -> [r1] -> X -> [r2] -> P -> bio. X is the shared, + % non-transportable intermediate of the two movable reactions. + model.id='chain'; model.comps={'c'}; model.compNames={'cytoplasm'}; + model.mets={'S_c';'X_c';'P_c'}; model.metNames={'S';'X';'P'}; model.metComps=[1;1;1]; + % EX_S r1 r2 bio + model.S=sparse([ 1 -1 0 0; % S + 0 1 -1 0; % X + 0 0 1 -1]); % P + model.rxns={'EX_S';'r1';'r2';'bio'}; model.rxnNames=model.rxns; + model.lb=[-10;0;0;0]; model.ub=[1000;1000;1000;1000]; model.rev=[1;0;0;0]; + model.c=[0;0;0;1]; model.b=zeros(3,1); + model.genes={'g1';'g2'}; model.grRules={'';'g1';'g2';''}; + model.rxnGeneMat=sparse([0 0; 1 0; 0 1; 0 0]); + end + function GSS = chainGss() + % g1 prefers 'c' strongly, g2 prefers 'm'; co-located in 'c' wins + % on combined score (0.9+0.6 > 0.1+0.9). + GSS.genes={'g1';'g2'}; GSS.compartments={'c';'m'}; + GSS.scores=[0.9 0.1; 0.6 0.9]; + end end end From 62ae554e0cd73151a4a4b23e1b4f1a5731e7e5da Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 18 Jul 2026 01:04:34 +0200 Subject: [PATCH 2/6] refactor: handle merged multi-compartment input and follow the score at genome scale Benchmarking against raven-toolbox on yeast-GEM (identical scores, relocate set and curated ground truth) surfaced two issues in the first version: - On a real multi-compartment model, the internal merge produced a 'System' compartment while pinned reactions were treated as being in the default, disconnecting the merged metabolites. The merged compartment is now relabelled to the default so placement and transports are consistent. - The score objective rewarded only the gene placement (y), leaving the reaction placement (x) underdetermined: any placement consistent with the gene assignment was optimal, so the solver scattered reactions arbitrarily among each gene's compartments. Reaction agreement with curated yeast-GEM was 51.9% (following the score argmax only 57.5% of the time). A small score-consistent tie-break on x now pulls each reaction into the compartment its own genes score highest, raising agreement to 70.8% (Python 72.0%) and direct placement overlap with Python to 81.4%, without changing the gene assignment. The model still certifies at the same growth. --- localization/assignCompartments.m | 23 ++++++++++++++++++++++- 1 file changed, 22 insertions(+), 1 deletion(-) diff --git a/localization/assignCompartments.m b/localization/assignCompartments.m index 115f3cfb..be606a33 100644 --- a/localization/assignCompartments.m +++ b/localization/assignCompartments.m @@ -105,8 +105,15 @@ % Merge to a single compartment so every metabolite has one identity; % placement and transports are then expressed relative to defaultCompartment. +% After merging, treat that single compartment AS the default: pinned +% reactions are considered to sit there, and placed reactions and transports +% are expressed relative to it. if numel(model.comps) > 1 model = mergeCompartments(model, true, true); + model.comps = {defaultCompartment}; + if isfield(model,'compNames'); model.compNames = {defaultCompartment}; end + if isfield(model,'compOutside'); model.compOutside = {''}; end + model.metComps = ones(numel(model.mets),1); end biomassIdx = find(model.c ~= 0, 1); @@ -292,11 +299,25 @@ end A_colo=sparse(coR,coC,coV,row,nVar); b_colo=zeros(row,1); -% objective: max sum score*y - multiPen*sum y +% objective: max sum score*y - multiPen*sum y, with a small score-consistent +% tie-break on the reaction placement. The gene reward on y decides which +% compartment each gene occupies, but leaves the reaction binaries x +% underdetermined (any placement consistent with the gene assignment is +% optimal); without the tie-break the solver picks one arbitrarily. The tiny +% reward on x pulls each reaction into the compartment its own genes score +% highest, without being large enough to change the gene assignment. +xTie = 1e-3; c = zeros(nVar,1); for gi=1:nGene for ci=1:nC; c(yCol(gi,ci)) = score(gi,ci) - multiPen; end end +for k=1:nMov + gOn = find(model.rxnGeneMat(movIdx(k),:)~=0); + for g=gOn + gi=find(geneIdx==g,1); if isempty(gi); continue; end + for ci=1:nC; c(xCol(k,ci)) = c(xCol(k,ci)) + xTie*score(gi,ci); end + end +end prob.A = [A_place; A_couple; A_has; A_gene1; A_force; A_colo]; prob.a = prob.A; From fe381e7e227eea29a8f7de5c4229d2d019105112 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 18 Jul 2026 01:33:10 +0200 Subject: [PATCH 3/6] refactor: make assignCompartments a faithful port of the raven-toolbox method Reworked assignCompartments to mirror raven_toolbox.localization .assign_compartments function-for-function, removing the earlier MATLAB-only divergences: - no internal merge: each reaction keeps its own compartment (pinned stay put), matching Python; the caller merges the draft first if desired. - placement objective is score*y - penalty*y only, dropping the MATLAB-only score tie-break on the reaction binaries. - transport pruning (pruneTransports, on by default): transports blocked in every certification medium are dropped, via haveFlux. - growth-gap feedback: a stranded biomass precursor's producer is pinned to biomass's compartment and the placement re-solved, with the protected/ gap-pinned mechanism so confinement and the growth-gap diagnoser cannot oscillate. - confinement is the full force / co-locate / relax logic, and relaxed metabolites flow into the split transports. - materialisation matches Python's id conventions and reconciles all model fields; the compartment set is the union of model and score compartments; transportable is keyed by metabolite name. Benchmarked against Python on yeast-GEM (identical scores/relocate/ground truth): both certify at growth 0.1426; transports 1117 vs 1001. Reaction agreement is 51.9% vs Python's 72.0% -- the objective rewards the gene placement and leaves the reaction placement underdetermined, so the residual is solver tie-breaking between the two Gurobi stacks, not a method difference. --- localization/assignCompartments.m | 646 +++++++++++-------- testing/function_tests/tAssignCompartments.m | 21 +- 2 files changed, 381 insertions(+), 286 deletions(-) diff --git a/localization/assignCompartments.m b/localization/assignCompartments.m index be606a33..b24472eb 100644 --- a/localization/assignCompartments.m +++ b/localization/assignCompartments.m @@ -1,63 +1,58 @@ function [outModel, placement, addedTransports, exitFlag, report] = assignCompartments(model, GSS, reactionsToRelocate, varargin) % assignCompartments Assign reactions to compartments, certified by growth. % -% Deterministic alternative to predictLocalization. Places the requested -% reactions into the compartments named by localization scores (GSS) and then -% verifies, by an actual FBA on the built model, that the objective (biomass) -% is still producible. The work is done in three phases: +% Faithful MATLAB port of raven_toolbox.localization.assign_compartments. Places +% reactions into subcellular compartments from soft gene x compartment +% localization scores while keeping the result functional, WITHOUT ever putting +% a flux model inside the placement optimisation: % -% 1. Place - a flux-free MILP that maximises the localization score. -% Its only variables are placement binaries (reaction -> -% compartment, gene -> compartment); it has no flux variables -% and no growth constraint, so a placement can never harvest a -% compartment's score through leaked flux. -% 2. Repair - a solver-free fixpoint that keeps the placement connected: a -% non-transportable metabolite split across compartments forces -% its reactions to co-locate; transportable splits get passive -% transports through the default compartment. -% 3. Certify - the placement is materialised into a compartmentalised model -% and that exact model is solved. exitFlag reflects whether it -% reached the growth floor, so the certificate is the model -% that is returned, never a placement the model cannot grow. +% 1. a flux-free placement master (maximise localization score, +% mono-localisation): with no flux variable and no growth constraint, a +% tolerance-rounded binary has nothing to leak into; +% 2. a structural confinement repair: reactions sharing a non-transportable +% metabolite are co-located (or pinned to a fixed reaction's compartment), +% and transportable pools a placement splits get star-topology transports; +% 3. materialised-FBA certification: the placement is confirmed functional by a +% real solveLP on the model i_applyAssignment builds, over the primary +% medium and every growth condition; +% 4. feedback on real failure: for a genuine placement gap, tighten the +% placement and re-solve until certified or the round budget is hit. % -% Mono-localization: each reaction is placed in exactly one compartment; a -% gene still spans compartments when it catalyses reactions placed in -% different ones. +% The model is NOT merged: each reaction keeps its current compartment, movable +% reactions are re-placed, and pinned reactions stay where they are (merge the +% model first, as a caller, to place onto a single-compartment draft). % % Parameters % ---------- % model : struct -% a RAVEN model with an objective (model.c) and grRules. Multiple -% compartments are merged before placement. +% a RAVEN model with an objective (model.c) and grRules. % GSS : struct % gene scoring structure (genes, compartments, scores) as from parseScores. -% GSS.compartments are the target compartment labels and must include -% defaultCompartment. % reactionsToRelocate : cell -% reaction ids to (re)place. Boundary reactions and the objective reaction -% are always pinned. +% reaction ids to (re)place. Boundary, multi-compartment and objective +% reactions are always pinned. % % Name-Value Arguments % -------------------- % defaultCompartment : char -% compartment that transports route through (usually cytosol); must be in -% GSS.compartments. +% compartment transports route through (usually cytosol); must be in the +% union of model and GSS compartments. % multiCompartmentPenalty : double (default 0.5) -% score cost per extra compartment a gene ends up in. +% score cost per compartment a gene ends up in. % minGrowth : double (default [] = 10%% of the unconstrained optimum) % required objective flux for certification. -% transportable : cell (default [] = all) -% metabolite ids (merged-model names) that may receive transports. -% Restricting it forces functionality-driven placement of reactions whose -% substrates are confined. +% transportable : cell (default [] = all movable base metabolites) +% metabolite NAMES that may receive transports. % growthConditions : struct array (default []) -% extra media the placement must also grow on, each with fields: name, -% medium (struct of exchangeRxnId -> max uptake), and minGrowth. +% extra media to certify on, each with fields name, medium +% (struct exchangeRxnId -> max uptake) and minGrowth. % maxRounds : double (default 8) -% maximum place/repair rounds before returning the best placement found. +% budget on placement-tightening rounds. +% pruneTransports : logical (default true) +% drop transports blocked in every certification medium. % transportCost : double (default 0.5) -% accepted for backward compatibility but not used: placement is driven by -% the score alone, and transports are a structural consequence of repair. +% accepted for signature compatibility but NOT used (placement is +% flux-free and score-only). % verbose : logical (default true) % % Returns @@ -65,17 +60,13 @@ % outModel : struct % the compartmentalised model. % placement : struct -% .rxns and .compartment (the assigned compartment id per relocated -% reaction). +% .rxns and .compartment (cell of compartment ids per placed reaction). % addedTransports : struct -% .mets and .compartment for the transports repair added. +% .mets (base metabolite names) and .compartment for the added transports. % exitFlag : double -% 1 = certified (the built model reaches the growth floor on every -% medium), -1 = could not place or the built model did not certify. +% 1 = certified on every medium, -1 = could not place or did not certify. % report : struct -% .certified (logical), .status (char), and .growths (struct of -% medium -> objective flux), so a failed certification reports the real -% growth rather than hiding it. +% .certified, .status, .growths (struct medium -> flux), .unplaced. % % See also % -------- @@ -83,105 +74,69 @@ p = parseRAVENargs(varargin, {'defaultCompartment',[]; 'transportCost',0.5; ... 'multiCompartmentPenalty',0.5; 'minGrowth',[]; 'transportable',[]; ... - 'growthConditions',[]; 'maxRounds',8; 'verbose',true}); + 'growthConditions',[]; 'maxRounds',8; 'pruneTransports',true; 'verbose',true}); defaultCompartment = char(p.defaultCompartment); multiPen = p.multiCompartmentPenalty; minGrowth = p.minGrowth; growthConditions = p.growthConditions; maxRounds = p.maxRounds; +pruneTransports = p.pruneTransports; verbose = p.verbose; -outModel = model; placement = struct('rxns',{{}},'compartment',{{}}); -addedTransports = struct('mets',{{}},'compartment',{{}}); exitFlag = -1; -report = struct('certified',false,'status','not_solved','growths',struct()); +outModel = model; +placement = struct('rxns',{{}},'compartment',{{}}); +addedTransports = struct('mets',{{}},'compartment',{{}}); +exitFlag = -1; +report = struct('certified',false,'status','not_solved','growths',struct(),'unplaced',{{}}); if all(model.c == 0) error('RAVEN:badInput','model has no objective (set model.c).'); end -comps = GSS.compartments(:); +% compartments = union of model and score compartments +comps = unique([model.comps(:); GSS.compartments(:)], 'stable'); +comps = sort(comps); if ~ismember(defaultCompartment, comps) - error('RAVEN:badInput','defaultCompartment ''%s'' not in GSS.compartments.', defaultCompartment); -end - -% Merge to a single compartment so every metabolite has one identity; -% placement and transports are then expressed relative to defaultCompartment. -% After merging, treat that single compartment AS the default: pinned -% reactions are considered to sit there, and placed reactions and transports -% are expressed relative to it. -if numel(model.comps) > 1 - model = mergeCompartments(model, true, true); - model.comps = {defaultCompartment}; - if isfield(model,'compNames'); model.compNames = {defaultCompartment}; end - if isfield(model,'compOutside'); model.compOutside = {''}; end - model.metComps = ones(numel(model.mets),1); -end -biomassIdx = find(model.c ~= 0, 1); - -if isempty(minGrowth) - sol = solveLP(model); - if isempty(sol.f) || sol.f <= 0 - error('RAVEN:badInput','merged model does not grow; pass minGrowth.'); - end - minGrowth = 0.1 * abs(sol.f); + error('RAVEN:badInput','defaultCompartment ''%s'' not in the model/score compartments.', defaultCompartment); end -% ---- scope: movable vs pinned ---- -isBoundary = (sum(model.S ~= 0, 1)' == 1); % one-metabolite reactions -relocate = ismember(model.rxns, reactionsToRelocate); -movableMask = relocate & ~isBoundary; -movableMask(biomassIdx) = false; -movIdx = find(movableMask); -pinIdx = find(~movableMask); -nMov = numel(movIdx); nC = numel(comps); -defC = find(strcmp(comps, defaultCompartment)); - -% Genes in scope: those on movable reactions that have a row in GSS -[gInGSS, gssRow] = ismember(model.genes, GSS.genes); -geneOnMov = any(model.rxnGeneMat(movIdx, :) ~= 0, 1)'; -geneIdx = find(gInGSS & geneOnMov); -nGene = numel(geneIdx); -score = zeros(nGene, nC); -for gi = 1:nGene - score(gi, :) = GSS.scores(gssRow(geneIdx(gi)), :); -end - -% Transportable metabolites (a base metabolite may be moved between -% compartments by a passive transport during repair). -movMetMask = any(model.S(:, movIdx) ~= 0, 2); -if isnumeric(p.transportable) && isempty(p.transportable) - transpMet = find(movMetMask); % default: all -else - transpMet = find(movMetMask & ismember(model.mets, p.transportable)); -end -isTransp = false(numel(model.mets),1); isTransp(transpMet) = true; +% ---- scope: biomass, growth floor, movable/pinned, genes, transportable ---- +sc = i_prepareScope(model, GSS, reactionsToRelocate, comps, defaultCompartment, minGrowth, p.transportable, verbose); +minGrowth = sc.minGrowth; +report.unplaced = sc.unplaced; if verbose fprintf('assignCompartments: %d movable, %d pinned reactions, %d genes, %d compartments.\n', ... - nMov, numel(pinIdx), nGene, nC); + numel(sc.movIdx), numel(sc.pinIdx), numel(sc.geneIdx), numel(comps)); end -% ---- place / repair rounds ---- -% forced(mi) = compartment index a movable reaction is pinned to (0 = free); -% groups is a cell array of movable-index vectors that must share a -% compartment. Both grow monotonically until the confinement fixpoint holds. -forced = zeros(nMov,1); +% ---- place -> repair -> certify (-> tighten) loop ---- +forced = containers.Map('KeyType','double','ValueType','double'); % movable local idx -> comp idx groups = {}; +gapPinned = containers.Map('KeyType','double','ValueType','logical'); seen = {}; -placeIdx = []; +best = struct('placeIdx',[],'trMet',[],'trComp',[],'growths',struct('primary',-1), ... + 'certified',false,'status','uncertified'); + for round = 1:maxRounds - placeIdx = i_placementMaster(model, movIdx, geneIdx, score, multiPen, ... - nC, forced, groups, verbose); + sig = i_signature(forced, groups); + if any(strcmp(seen, sig)) + break; % configuration already tried: diagnosers are cycling + end + seen{end+1} = sig; %#ok + + placeIdx = i_placementMaster(model, sc, comps, multiPen, forced, groups, verbose); if isempty(placeIdx) report.status = 'infeasible'; - if verbose; fprintf('assignCompartments: placement MILP infeasible.\n'); end return; end - [newForced, newGroups] = i_diagnoseConfinement(model, movIdx, pinIdx, ... - placeIdx, defC, isTransp); - % Apply only genuinely new tightenings; stop at the fixpoint. + + % confinement fixpoint + [newForced, newGroups, relaxed] = i_diagnoseConfinement(model, sc, comps, placeIdx, gapPinned); changed = false; - for k = 1:nMov - if newForced(k) ~= 0 && forced(k) == 0 + fk = keys(newForced); + for i = 1:numel(fk) + k = fk{i}; + if (~isKey(forced,k) || forced(k) ~= newForced(k)) && ~isKey(gapPinned,k) forced(k) = newForced(k); changed = true; end end @@ -190,58 +145,128 @@ groups{end+1} = newGroups{gi}; changed = true; %#ok end end - sig = i_signature(forced, groups); - if ~changed || any(strcmp(seen, sig)) - break; + if changed + continue; end - seen{end+1} = sig; %#ok + + % transports + materialise + certify + [trMet, trComp] = i_splitTransports(model, sc, comps, defaultCompartment, placeIdx, relaxed); + outModel = i_applyAssignment(model, sc, comps, defaultCompartment, placeIdx, trMet, trComp); + [ok, growths] = i_certify(outModel, sc.biomassId, minGrowth, growthConditions); + + if ok + if pruneTransports && ~isempty(trMet) + [trMet, trComp] = i_usableTransports(outModel, trMet, trComp, comps, sc.biomassId, growthConditions); + outModel = i_applyAssignment(model, sc, comps, defaultCompartment, placeIdx, trMet, trComp); + end + exitFlag = 1; report.status = 'certified'; report.certified = true; report.growths = growths; + placement.rxns = model.rxns(sc.movIdx); placement.compartment = comps(placeIdx); + addedTransports.mets = sc.metNames(trMet); addedTransports.compartment = comps(trComp); + return; + end + + % keep best partial (largest primary growth) + if growths.primary > best.growths.primary + best = struct('placeIdx',placeIdx,'trMet',trMet,'trComp',trComp,'growths',growths, ... + 'certified',false,'status','uncertified'); + end + + % growth-gap feedback: pin the sole producer of a stranded biomass precursor + fb = i_diagnoseGrowthGap(outModel, model, sc, comps, placeIdx); + if ~isempty(fb) && (~isKey(forced,fb(1)) || forced(fb(1)) ~= fb(2)) + forced(fb(1)) = fb(2); gapPinned(fb(1)) = true; + continue; + end + break; % no tightening available end -% ---- placement + transports ---- -placeRxns = model.rxns(movIdx); -placeComp = comps(placeIdx); -placement.rxns = placeRxns(:); placement.compartment = placeComp(:); +% honest uncertified result: return the best partial found +if ~isempty(best.placeIdx) + outModel = i_applyAssignment(model, sc, comps, defaultCompartment, best.placeIdx, best.trMet, best.trComp); + placement.rxns = model.rxns(sc.movIdx); placement.compartment = comps(best.placeIdx); + addedTransports.mets = sc.metNames(best.trMet); addedTransports.compartment = comps(best.trComp); + report.growths = best.growths; +end +report.status = 'uncertified'; +if verbose + fprintf('assignCompartments: did not certify (growth %.4g < %.4g).\n', best.growths.primary, minGrowth); +end +end -[trMet, trComp] = i_splitTransports(model, movIdx, pinIdx, placeIdx, defC, isTransp, nC); -addedTransports.mets = model.mets(trMet); -addedTransports.compartment = comps(trComp); +% ============================================================ scope +function sc = i_prepareScope(model, GSS, relocate, comps, defaultCompartment, minGrowth, transportable, verbose) %#ok +biomassIdx = find(model.c ~= 0, 1); +sc.biomassId = model.rxns{biomassIdx}; -% ---- materialise + certify ---- -outModel = i_applyAssignment(model, placement, addedTransports, comps, defaultCompartment); -[certified, growths] = i_certify(outModel, biomassIdx, model.rxns{biomassIdx}, minGrowth, growthConditions); -report.certified = certified; -report.growths = growths; -if certified - exitFlag = 1; report.status = 'certified'; -else - exitFlag = -1; report.status = 'uncertified'; - if verbose - fprintf('assignCompartments: placement did not certify (growth %.4g < %.4g).\n', ... - growths.primary, minGrowth); +if isempty(minGrowth) + sol = solveLP(model); + if isempty(sol.f) || abs(sol.f) <= 0 + error('RAVEN:badInput','the draft model does not grow; pass minGrowth.'); end + minGrowth = 0.1 * abs(sol.f); end +sc.minGrowth = minGrowth; + +% each reaction's single compartment ('' if boundary/multi-compartment) +sc.rxnComp = i_reactionCompartments(model); + +toRelocate = ismember(model.rxns, relocate); +isBoundary = (sum(model.S ~= 0, 1)' == 1); +movable = toRelocate & ~isBoundary & ~cellfun(@isempty, sc.rxnComp); +movable(biomassIdx) = false; +sc.movIdx = find(movable); +sc.pinIdx = find(~movable); + +% genes on movable reactions that are scored +[gInGSS, gssRow] = ismember(model.genes, GSS.genes); +geneOnMov = any(model.rxnGeneMat(sc.movIdx, :) ~= 0, 1)'; +scope = gInGSS & geneOnMov; +sc.geneIdx = find(scope); +sc.score = zeros(numel(sc.geneIdx), numel(comps)); +% align score columns (GSS.compartments) to comps +[~, colOf] = ismember(GSS.compartments, comps); +for gi = 1:numel(sc.geneIdx) + row = GSS.scores(gssRow(sc.geneIdx(gi)), :); + sc.score(gi, colOf(colOf>0)) = row(colOf>0); +end + +% unplaced: movable reactions with genes but none scored +sc.unplaced = {}; +for k = 1:numel(sc.movIdx) + g = find(model.rxnGeneMat(sc.movIdx(k),:) ~= 0); + if ~isempty(g) && ~any(ismember(g, sc.geneIdx)) + sc.unplaced{end+1,1} = model.rxns{sc.movIdx(k)}; %#ok + end end -% ------------------------------------------------------------- placement MILP -function placeIdx = i_placementMaster(model, movIdx, geneIdx, score, multiPen, nC, forced, groups, verbose) -% Flux-free score MILP. Variables: x[mov,c] and y[gene,c], both binary. -% Returns the chosen compartment index per movable reaction, or [] if -% infeasible. -nMov = numel(movIdx); nGene = numel(geneIdx); +% base metabolite = metabolite name; transportable base names +sc.metNames = model.metNames; +movMetMask = any(model.S(:, sc.movIdx) ~= 0, 2); +if isnumeric(transportable) && isempty(transportable) + sc.isTransp = movMetMask; % default: all movable bases +else + sc.isTransp = movMetMask & ismember(model.metNames, transportable); +end +end + +% ============================================================ placement MILP +function placeIdx = i_placementMaster(model, sc, comps, multiPen, forced, groups, verbose) +% Flux-free score MILP. Variables x[mov,c] and y[gene,c], both binary. +% Objective: max sum score*y - multiPen*sum y (no reward on x). +movIdx = sc.movIdx; geneIdx = sc.geneIdx; score = sc.score; +nMov = numel(movIdx); nGene = numel(geneIdx); nC = numel(comps); oX = 0; nX = nMov*nC; oY = oX+nX; nY = nGene*nC; nVar = oY+nY; xCol = @(mi,ci) oX + (mi-1)*nC + ci; yCol = @(gi,ci) oY + (gi-1)*nC + ci; -% placement: sum_c x[r,c] = 1 pR=[];pC=[];pV=[]; for k=1:nMov for ci=1:nC; pR(end+1)=k; pC(end+1)=xCol(k,ci); pV(end+1)=1; end %#ok end A_place=sparse(pR,pC,pV,nMov,nVar); b_place=ones(nMov,1); -% coupling: x[r,c] - y[g,c] <= 0 for each gene g on reaction r cR=[];cC=[];cV=[];row=0; for k=1:nMov gOn=find(model.rxnGeneMat(movIdx(k),:)~=0); @@ -256,7 +281,6 @@ end A_couple=sparse(cR,cC,cV,row,nVar); b_couple=zeros(row,1); -% has: y[g,c] - sum_{r of g} x[r,c] <= 0 hR=[];hC=[];hV=[];row=0; for gi=1:nGene g=geneIdx(gi); @@ -268,23 +292,19 @@ end A_has=sparse(hR,hC,hV,row,nVar); b_has=zeros(row,1); -% gene1: sum_c y[g,c] >= 1 g1R=[];g1C=[];g1V=[]; for gi=1:nGene for ci=1:nC; g1R(end+1)=gi; g1C(end+1)=yCol(gi,ci); g1V(end+1)=1; end %#ok end A_gene1=sparse(g1R,g1C,g1V,nGene,nVar); b_gene1=ones(nGene,1); -% forced pins: x[r,c_force] = 1 -fR=[];fC=[];nForce=0; -for k=1:nMov - if forced(k)~=0 - nForce=nForce+1; fR(end+1)=nForce; fC(end+1)=xCol(k,forced(k)); %#ok - end +fR=[];fC=[];nForce=0; fkeys=keys(forced); +for i=1:numel(fkeys) + k=fkeys{i}; + nForce=nForce+1; fR(end+1)=nForce; fC(end+1)=xCol(k,forced(k)); %#ok end A_force=sparse(fR,fC,ones(1,nForce),nForce,nVar); b_force=ones(nForce,1); -% co-location groups: x[a,c] - x[b,c] = 0 for consecutive members coR=[];coC=[];coV=[];row=0; for gi=1:numel(groups) mem=groups{gi}; @@ -299,25 +319,10 @@ end A_colo=sparse(coR,coC,coV,row,nVar); b_colo=zeros(row,1); -% objective: max sum score*y - multiPen*sum y, with a small score-consistent -% tie-break on the reaction placement. The gene reward on y decides which -% compartment each gene occupies, but leaves the reaction binaries x -% underdetermined (any placement consistent with the gene assignment is -% optimal); without the tie-break the solver picks one arbitrarily. The tiny -% reward on x pulls each reaction into the compartment its own genes score -% highest, without being large enough to change the gene assignment. -xTie = 1e-3; c = zeros(nVar,1); for gi=1:nGene for ci=1:nC; c(yCol(gi,ci)) = score(gi,ci) - multiPen; end end -for k=1:nMov - gOn = find(model.rxnGeneMat(movIdx(k),:)~=0); - for g=gOn - gi=find(geneIdx==g,1); if isempty(gi); continue; end - for ci=1:nC; c(xCol(k,ci)) = c(xCol(k,ci)) + xTie*score(gi,ci); end - end -end prob.A = [A_place; A_couple; A_has; A_gene1; A_force; A_colo]; prob.a = prob.A; @@ -325,8 +330,7 @@ prob.csense = [repmat('E',1,nMov), repmat('L',1,numel(b_couple)), ... repmat('L',1,numel(b_has)), repmat('G',1,nGene), ... repmat('E',1,nForce), repmat('E',1,numel(b_colo))]; -prob.c = -c; % optimizeProb minimises; we maximise c -prob.osense = 1; +prob.c = -c; prob.osense = 1; prob.lb = zeros(nVar,1); prob.ub = ones(nVar,1); prob.vartype = repmat('B',1,nVar); @@ -343,112 +347,103 @@ end end -% ------------------------------------------------------ confinement diagnosis -function [forced, groups] = i_diagnoseConfinement(model, movIdx, pinIdx, placeIdx, defC, isTransp) -% Solver-free: find non-transportable metabolites split across compartments -% and decide, per split, whether to force the movers to the pinned -% compartment or to co-locate them. Pinned reactions all live in the default -% compartment (the model was merged first). -nMov = numel(movIdx); -forced = zeros(nMov,1); +% ============================================================ confinement +function [forced, groups, relaxed] = i_diagnoseConfinement(model, sc, comps, placeIdx, gapPinned) +% Find non-transportable metabolites a placement splits across compartments; +% force to a pinned compartment, co-locate all-movable sets, or relax +% (transport anyway) when pinned into >=2 comps or shared with a protected pin. +forced = containers.Map('KeyType','double','ValueType','double'); groups = {}; +relaxed = containers.Map('KeyType','double','ValueType','logical'); + +placedComp = i_placedCompartments(model, sc, comps, placeIdx); % per reaction, comp idx or 0 +movLocal = containers.Map('KeyType','double','ValueType','double'); % rxn idx -> movable local idx +for k=1:numel(sc.movIdx); movLocal(sc.movIdx(k)) = k; end -% used{metRow} maps compartment index -> list of movable local indices. -% pinnedHere(metRow) is true if a pinned reaction touches the metabolite. nMet = numel(model.mets); -usedComp = cell(nMet,1); -usedMov = cell(nMet,1); -pinnedHere = false(nMet,1); -for k=1:nMov - mrows = find(model.S(:,movIdx(k))~=0)'; - ci = placeIdx(k); - for mm=mrows - if isTransp(mm); continue; end % transportable: handled by a transport - usedComp{mm}(end+1) = ci; - usedMov{mm}(end+1) = k; - end -end -for k=1:numel(pinIdx) - mrows = find(model.S(:,pinIdx(k))~=0)'; - for mm=mrows - if isTransp(mm); continue; end - usedComp{mm}(end+1) = defC; - pinnedHere(mm) = true; +usedComp = cell(nMet,1); usedRxn = cell(nMet,1); +allRxn = [sc.movIdx; sc.pinIdx]; +for r = allRxn' + comp = placedComp(r); + if comp == 0; continue; end % multi-compartment reaction: bridges pools + for mm = find(model.S(:,r)~=0)' + if sc.isTransp(mm); continue; end + usedComp{mm}(end+1) = comp; usedRxn{mm}(end+1) = r; end end for mm=1:nMet cs = unique(usedComp{mm}); - if numel(cs) <= 1; continue; end % lives in one compartment: fine - movers = unique(usedMov{mm}); - if pinnedHere(mm) - % A pinned reaction anchors this metabolite in the default - % compartment; every movable toucher must join it there. - for k=movers; forced(k) = defC; end - else - % Only movable reactions touch it: co-locate them and let the score - % objective choose the shared compartment. - groups{end+1} = movers(:)'; %#ok + if numel(cs) <= 1; continue; end + touching = usedRxn{mm}; + isMov = arrayfun(@(r) isKey(movLocal,r), touching); + movers = touching(isMov); + % pinned compartments = comps where a pinned (non-movable) reaction touches mm + pinnedComps = unique(arrayfun(@(i) placedComp(touching(i)), find(~isMov))); + anyProtected = any(arrayfun(@(r) isKey(movLocal,r) && isKey(gapPinned,movLocal(r)), touching)); + if numel(pinnedComps) > 1 || anyProtected + relaxed(mm) = true; + elseif ~isempty(pinnedComps) + for r = movers; forced(movLocal(r)) = pinnedComps(1); end + elseif ~isempty(movers) + g = sort(unique(arrayfun(@(r) movLocal(r), movers))); + groups{end+1} = g(:)'; %#ok end end end -% ----------------------------------------------------------- split transports -function [trMet, trComp] = i_splitTransports(model, movIdx, pinIdx, placeIdx, defC, isTransp, nC) -% For each transportable base metabolite placed in more than one compartment, -% add a transport between the default compartment and each non-default one. +% ============================================================ split transports +function [trMet, trComp] = i_splitTransports(model, sc, comps, defaultCompartment, placeIdx, relaxed) +defC = find(strcmp(comps, defaultCompartment), 1); +placedComp = i_placedCompartments(model, sc, comps, placeIdx); +allRxn = [sc.movIdx; sc.pinIdx]; nMet = numel(model.mets); -compsUsed = cell(nMet,1); -for k=1:numel(movIdx) - mrows = find(model.S(:,movIdx(k))~=0)'; - for mm=mrows; compsUsed{mm}(end+1) = placeIdx(k); end -end -for k=1:numel(pinIdx) - mrows = find(model.S(:,pinIdx(k))~=0)'; - for mm=mrows; compsUsed{mm}(end+1) = defC; end +used = cell(nMet,1); +for r = allRxn' + comp = placedComp(r); + if comp == 0; continue; end + for mm = find(model.S(:,r)~=0)' + if sc.isTransp(mm) || (isKey(relaxed,mm)) + used{mm}(end+1) = comp; + end + end end trMet=[]; trComp=[]; for mm=1:nMet - if ~isTransp(mm); continue; end - cs = unique(compsUsed{mm}); + cs = unique(used{mm}); if numel(cs) <= 1; continue; end - for ci=cs - if ci==defC; continue; end - trMet(end+1)=mm; trComp(end+1)=ci; %#ok + for ci = sort(cs) + if ci ~= defC; trMet(end+1)=mm; trComp(end+1)=ci; end %#ok end end trMet=trMet(:); trComp=trComp(:); end -% ------------------------------------------------------------- certification -function [certified, growths] = i_certify(outModel, biomassIdx, biomassId, minGrowth, growthConditions) -% Solve the materialised model and confirm it reaches minGrowth on the primary -% medium and on every extra growth condition. +% ============================================================ certification +function [ok, growths] = i_certify(outModel, biomassId, minGrowth, growthConditions) tol = 1e-9; growths = struct(); -bIdx = find(strcmp(outModel.rxns, biomassId), 1); -if isempty(bIdx); bIdx = biomassIdx; end -sol = solveLP(outModel); -primary = 0; -if ~isempty(sol.f); primary = abs(sol.f); end -growths.primary = primary; -certified = primary >= minGrowth - tol; - +growths.primary = i_growOn(outModel, biomassId, []); +ok = growths.primary >= minGrowth - tol; if ~isempty(growthConditions) && isstruct(growthConditions) for i=1:numel(growthConditions) gc = growthConditions(i); - m2 = i_applyMedium(outModel, gc.medium); - s2 = solveLP(m2); - g = 0; if ~isempty(s2.f); g = abs(s2.f); end + g = i_growOn(outModel, biomassId, gc.medium); growths.(matlab.lang.makeValidName(gc.name)) = g; - certified = certified && (g >= gc.minGrowth - tol); + ok = ok && (g >= gc.minGrowth - tol); end end end +function g = i_growOn(model, biomassId, medium) +if ~isempty(medium); model = i_applyMedium(model, medium); end +bIdx = find(strcmp(model.rxns, biomassId), 1); +model.c = zeros(numel(model.rxns),1); model.c(bIdx) = 1; +sol = solveLP(model); +if isempty(sol.f); g = 0; else; g = abs(sol.f); end +end + function model = i_applyMedium(model, medium) -% Close all uptake, then open only the listed exchanges. medium is a struct -% of exchangeRxnId -> max uptake. [~, exchIdx] = getExchangeRxns(model); model.lb(exchIdx(model.lb(exchIdx) < 0)) = 0; if ~isempty(medium) && isstruct(medium) @@ -460,26 +455,64 @@ end end -% ------------------------------------------------------------------ helpers -function tf = i_hasGroup(groups, g) -tf = false; -g = sort(g(:)'); -for i=1:numel(groups) - if isequal(sort(groups{i}(:)'), g); tf = true; return; end +% ============================================================ transport pruning +function [trMet, trComp] = i_usableTransports(outModel, trMet, trComp, comps, biomassId, growthConditions) +% Keep only transports that can carry flux in at least one certification +% medium (sound: an unusable reaction's removal cannot change any of those +% FBAs). +trId = cell(numel(trMet),1); +for i=1:numel(trMet) + trId{i} = ['tr_' num2str(i-1) '_' comps{trComp(i)}]; end +media = {[]}; +if ~isempty(growthConditions) && isstruct(growthConditions) + for i=1:numel(growthConditions); media{end+1} = growthConditions(i).medium; end %#ok +end +keep = false(numel(trMet),1); +for mi=1:numel(media) + m = outModel; + if ~isempty(media{mi}); m = i_applyMedium(m, media{mi}); end + bIdx = find(strcmp(m.rxns, biomassId), 1); + m.c = zeros(numel(m.rxns),1); m.c(bIdx) = 1; + idx = zeros(numel(trId),1); + for i=1:numel(trId); j=find(strcmp(m.rxns,trId{i}),1); if ~isempty(j); idx(i)=j; end; end + valid = idx>0; + fl = false(numel(trId),1); + fl(valid) = haveFlux(m, 'rxns', idx(valid)); + keep = keep | fl; +end +trMet = trMet(keep); trComp = trComp(keep); +end + +% ============================================================ growth-gap feedback +function fb = i_diagnoseGrowthGap(outModel, model, sc, comps, placeIdx) +% If a biomass precursor is blocked, pin its sole movable producer to biomass's +% compartment. Returns [movableLocalIdx, compIdx] or []. +fb = []; +bIdx = find(strcmp(outModel.rxns, sc.biomassId), 1); +if isempty(bIdx); return; end +bioComp = i_reactionCompartments(outModel); bioComp = bioComp{bIdx}; +if isempty(bioComp); return; end +bcIdx = find(strcmp(comps, bioComp),1); +% precursors: base names biomass consumes +precRows = find(outModel.S(:,bIdx) < 0); +precNames = unique(outModel.metNames(precRows)); +% blocked reactions in the materialised model +canFlux = haveFlux(outModel); +placedComp = i_placedCompartments(model, sc, comps, placeIdx); +for k=1:numel(sc.movIdx) + r = sc.movIdx(k); + ri = find(strcmp(outModel.rxns, model.rxns{r}),1); + if isempty(ri) || canFlux(ri); continue; end + touchNames = model.metNames(model.S(:,r)~=0); + if any(ismember(touchNames, precNames)) && placedComp(r) ~= bcIdx + fb = [k, bcIdx]; return; + end end - -function s = i_signature(forced, groups) -gp = ''; -for i=1:numel(groups); gp = [gp '|' num2str(sort(groups{i}(:)'))]; end %#ok -s = [num2str(forced(:)') '#' gp]; end -% ------------------------------------------------------------------- apply -function outModel = i_applyAssignment(model, placement, addedTransports, comps, defaultCompartment) -% Build the compartmentalised model: relabel each placed reaction's metabolites -% into its compartment (creating per-compartment metabolite copies) and add the -% requested transports. +% ============================================================ materialisation +function outModel = i_applyAssignment(model, sc, comps, defaultCompartment, placeIdx, trMet, trComp) outModel = model; for ci = 1:numel(comps) if ~ismember(comps{ci}, outModel.comps) @@ -487,19 +520,40 @@ if isfield(outModel,'compNames'); outModel.compNames{end+1,1} = comps{ci}; end end end -for k = 1:numel(placement.rxns) - r = find(strcmp(outModel.rxns, placement.rxns{k}), 1); - comp = placement.compartment{k}; - srcR = find(strcmp(model.rxns,placement.rxns{k}),1); - mrows = find(model.S(:, srcR) ~= 0); %#ok<*FNDSB> +for k = 1:numel(sc.movIdx) + r = sc.movIdx(k); comp = comps{placeIdx(k)}; + mrows = find(model.S(:, r) ~= 0); for mm = mrows' [outModel, newMet] = i_metInComp(outModel, model, mm, comp); - outModel.S(newMet, r) = model.S(mm, srcR); + outModel.S(newMet, r) = model.S(mm, r); if ~isequal(newMet, mm); outModel.S(mm, r) = 0; end end end -for k = 1:numel(addedTransports.mets) - outModel = i_addTransport(outModel, model, addedTransports.mets{k}, addedTransports.compartment{k}, defaultCompartment); +for k = 1:numel(trMet) + outModel = i_addTransport(outModel, model, trMet(k), comps{trComp(k)}, defaultCompartment, k-1); +end +outModel = i_reconcileFields(outModel); +end + +function model = i_reconcileFields(model) +% Pad every registered rxn/met/gene/comp-indexed field to the current entity +% count, so the materialised model stays consistent for downstream functions +% (haveFlux, removeReactions, ...). Only the core fields are grown as +% reactions/metabolites are added; the rest are reconciled here. +reg = ravenModelFields(); +counts = struct('rxn',numel(model.rxns),'met',numel(model.mets), ... + 'gene',numel(model.genes),'comp',numel(model.comps)); +for i=1:numel(reg) + f = reg(i).name; + if ~isfield(model,f) || strcmp(f,'rxns') || strcmp(f,'mets') || ... + strcmp(f,'genes') || strcmp(f,'comps'); continue; end + n = counts.(reg(i).type); + cur = size(model.(f),1); + if cur >= n; continue; end + def = reg(i).default; + for j = cur+1:n + if iscell(model.(f)); model.(f){j,1} = def; else; model.(f)(j,1) = def; end + end end end @@ -508,7 +562,7 @@ name = model.metNames{srcMetRow}; cand = find(strcmp(outModel.metNames, name) & outModel.metComps == ci, 1); if ~isempty(cand); idx = cand; return; end -newId = [model.mets{srcMetRow} '_' comp]; +newId = [model.mets{srcMetRow} '__' comp]; outModel.mets{end+1,1} = newId; outModel.metNames{end+1,1} = name; outModel.metComps(end+1,1) = ci; @@ -517,12 +571,10 @@ idx = numel(outModel.mets); end -function outModel = i_addTransport(outModel, model, metId, comp, defaultCompartment) -srcRow = find(strcmp(model.mets, metId), 1); -if isempty(srcRow); return; end +function outModel = i_addTransport(outModel, model, srcRow, comp, defaultCompartment, n) [outModel, dRow] = i_metInComp(outModel, model, srcRow, defaultCompartment); [outModel, cRow] = i_metInComp(outModel, model, srcRow, comp); -trId = ['tr_' model.metNames{srcRow} '_' comp]; +trId = ['tr_' num2str(n) '_' comp]; if any(strcmp(outModel.rxns, trId)); return; end outModel.rxns{end+1,1} = trId; outModel.S(:, end+1) = 0; outModel.S(dRow, end) = -1; outModel.S(cRow, end) = 1; @@ -533,3 +585,43 @@ if isfield(outModel,'grRules'); outModel.grRules{end+1,1} = ''; end if isfield(outModel,'rxnGeneMat'); outModel.rxnGeneMat(end+1,:) = 0; end end + +% ============================================================ helpers +function comp = i_reactionCompartments(model) +% Per reaction: the single compartment id of its metabolites, '' if boundary +% or multi-compartment. +comp = cell(numel(model.rxns),1); +for r=1:numel(model.rxns) + cc = unique(model.metComps(model.S(:,r)~=0)); + if numel(cc)==1; comp{r} = model.comps{cc}; else; comp{r} = ''; end +end +end + +function pc = i_placedCompartments(model, sc, comps, placeIdx) +% Per reaction index: compartment index. Movable -> placeIdx; pinned -> its own +% compartment; 0 if boundary/multi-compartment. +pc = zeros(numel(model.rxns),1); +for k=1:numel(sc.movIdx); pc(sc.movIdx(k)) = placeIdx(k); end +for r = sc.pinIdx' + if ~isempty(sc.rxnComp{r}) + ci = find(strcmp(comps, sc.rxnComp{r}),1); + if ~isempty(ci); pc(r) = ci; end + end +end +end + +function tf = i_hasGroup(groups, g) +tf = false; g = sort(g(:)'); +for i=1:numel(groups) + if isequal(sort(groups{i}(:)'), g); tf = true; return; end +end +end + +function s = i_signature(forced, groups) +fk = sort(cell2mat(keys(forced))); +fs = ''; +for i=1:numel(fk); fs = [fs sprintf('%d:%d,', fk(i), forced(fk(i)))]; end %#ok +gs = ''; +for i=1:numel(groups); gs = [gs '|' num2str(sort(groups{i}(:)'))]; end %#ok +s = [fs '#' gs]; +end diff --git a/testing/function_tests/tAssignCompartments.m b/testing/function_tests/tAssignCompartments.m index eb59a6a5..c0236fd5 100644 --- a/testing/function_tests/tAssignCompartments.m +++ b/testing/function_tests/tAssignCompartments.m @@ -76,19 +76,19 @@ function confinementColocatesMovableReactions(testCase) testCase.assumeMILPSolver(); model = tAssignCompartments.chainToy(); GSS = tAssignCompartments.chainGss(); - % X_c non-transportable, S_c/P_c transportable (so the biomass - % path is not what forces co-location). + % X non-transportable (S/P transportable, keyed by metabolite + % name), so the biomass path is not what forces co-location. evalc(['[~, place, ~, ok] = assignCompartments(model, GSS, {''r1'';''r2''}, ' ... - '''defaultCompartment'', ''c'', ''transportable'', {''S_c'';''P_c''}, ''verbose'', false);']); + '''defaultCompartment'', ''c'', ''transportable'', {''S'';''P''}, ''verbose'', false);']); testCase.verifyEqual(ok, 1); c1 = place.compartment{strcmp(place.rxns,'r1')}; c2 = place.compartment{strcmp(place.rxns,'r2')}; testCase.verifyEqual(c1, c2); % co-located - % With X_c transportable too, nothing forces co-location and the - % opposing scores split the two reactions across compartments. + % With X transportable too, nothing forces co-location and the + % dominant scores split the two reactions across compartments. evalc(['[~, sp] = assignCompartments(model, GSS, {''r1'';''r2''}, ' ... - '''defaultCompartment'', ''c'', ''transportable'', {''S_c'';''X_c'';''P_c''}, ''verbose'', false);']); + '''defaultCompartment'', ''c'', ''transportable'', {''S'';''X'';''P''}, ''verbose'', false);']); testCase.verifyNotEqual(sp.compartment{strcmp(sp.rxns,'r1')}, ... sp.compartment{strcmp(sp.rxns,'r2')}); end @@ -125,10 +125,13 @@ function confinementColocatesMovableReactions(testCase) model.rxnGeneMat=sparse([0 0; 1 0; 0 1; 0 0]); end function GSS = chainGss() - % g1 prefers 'c' strongly, g2 prefers 'm'; co-located in 'c' wins - % on combined score (0.9+0.6 > 0.1+0.9). + % Each gene has a single dominant compartment (the other score is + % below the 0.5 multi-compartment penalty), so placement is + % determined without a tie-break: g1 -> 'c', g2 -> 'm'. When X is + % transportable they split; when X is non-transportable the shared + % pool co-locates them, in 'c' (combined 0.9+0.3 > 0.1+0.9). GSS.genes={'g1';'g2'}; GSS.compartments={'c';'m'}; - GSS.scores=[0.9 0.1; 0.6 0.9]; + GSS.scores=[0.9 0.1; 0.3 0.9]; end end end From a7139a016a77c6dee5873e72c9c94c5499a2eb39 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 18 Jul 2026 08:32:28 +0200 Subject: [PATCH 4/6] refactor: close faithfulness gaps found in the line-by-line review Reconciled the divergences from raven_toolbox.localization.assign_compartments surfaced by a function-by-function comparison: - movable reactions are now processed in sorted-reaction-id order, matching the Python master's variable order (a faithfulness fix; it does not change the genome-scale result, so the placement difference vs Python is solver tie-breaking between the two Gurobi stacks, not variable ordering). - biomass reaction is the largest objective-coefficient reaction (was the first non-zero), and a biomassReaction argument can name it explicitly. - a baseMetabolite argument selects the compartment-agnostic key (default 'name', the RAVEN convention; the Python default strips the id's compartment suffix, which is cobra's convention). - confinement and split transports are keyed by base metabolite, not by per-compartment row, so a species split across compartments is detected on a non-merged model as it is in Python (identical on a merged draft). - minimizeTransports (parsimonious-FBA transport pruning), off by default, is ported. Not ported: universal-model gap-fill, and multi-localization (which needs a loopless FVA RAVEN does not provide). Genome-scale benchmark unchanged (certified, growth 0.1426, 1117 transports). --- localization/assignCompartments.m | 165 +++++++++++++++++++++--------- 1 file changed, 118 insertions(+), 47 deletions(-) diff --git a/localization/assignCompartments.m b/localization/assignCompartments.m index b24472eb..8ab57ed2 100644 --- a/localization/assignCompartments.m +++ b/localization/assignCompartments.m @@ -50,11 +50,24 @@ % budget on placement-tightening rounds. % pruneTransports : logical (default true) % drop transports blocked in every certification medium. +% minimizeTransports : logical (default false) +% after certifying, prune transports to those carrying flux in a +% parsimonious-FBA solution, re-certifying (falls back if a needed +% transport is pruned). +% biomassReaction : char (default [] = the largest objective-coefficient +% reaction) +% the reaction whose flux is the growth objective. +% baseMetabolite : char (default 'name') +% the compartment-agnostic metabolite key: 'name' (metNames, the RAVEN +% convention) or 'id' (model.mets). % transportCost : double (default 0.5) % accepted for signature compatibility but NOT used (placement is % flux-free and score-only). % verbose : logical (default true) % +% Not ported from the raven-toolbox reference: universal-model gap-fill, and +% multi-localization (which requires a loopless FVA RAVEN does not provide). +% % Returns % ------- % outModel : struct @@ -74,13 +87,15 @@ p = parseRAVENargs(varargin, {'defaultCompartment',[]; 'transportCost',0.5; ... 'multiCompartmentPenalty',0.5; 'minGrowth',[]; 'transportable',[]; ... - 'growthConditions',[]; 'maxRounds',8; 'pruneTransports',true; 'verbose',true}); + 'growthConditions',[]; 'maxRounds',8; 'pruneTransports',true; ... + 'minimizeTransports',false; 'biomassReaction',[]; 'baseMetabolite',[]; 'verbose',true}); defaultCompartment = char(p.defaultCompartment); multiPen = p.multiCompartmentPenalty; minGrowth = p.minGrowth; growthConditions = p.growthConditions; maxRounds = p.maxRounds; pruneTransports = p.pruneTransports; +minimizeTransports = p.minimizeTransports; verbose = p.verbose; outModel = model; @@ -100,7 +115,7 @@ end % ---- scope: biomass, growth floor, movable/pinned, genes, transportable ---- -sc = i_prepareScope(model, GSS, reactionsToRelocate, comps, defaultCompartment, minGrowth, p.transportable, verbose); +sc = i_prepareScope(model, GSS, reactionsToRelocate, comps, minGrowth, p.transportable, p.biomassReaction, p.baseMetabolite); minGrowth = sc.minGrowth; report.unplaced = sc.unplaced; @@ -114,7 +129,7 @@ groups = {}; gapPinned = containers.Map('KeyType','double','ValueType','logical'); seen = {}; -best = struct('placeIdx',[],'trMet',[],'trComp',[],'growths',struct('primary',-1), ... +best = struct('placeIdx',[],'trBase',[],'trComp',[],'growths',struct('primary',-1), ... 'certified',false,'status','uncertified'); for round = 1:maxRounds @@ -150,24 +165,29 @@ end % transports + materialise + certify - [trMet, trComp] = i_splitTransports(model, sc, comps, defaultCompartment, placeIdx, relaxed); - outModel = i_applyAssignment(model, sc, comps, defaultCompartment, placeIdx, trMet, trComp); + [trBase, trComp] = i_splitTransports(model, sc, comps, defaultCompartment, placeIdx, relaxed); + outModel = i_applyAssignment(model, sc, comps, defaultCompartment, placeIdx, trBase, trComp); [ok, growths] = i_certify(outModel, sc.biomassId, minGrowth, growthConditions); if ok - if pruneTransports && ~isempty(trMet) - [trMet, trComp] = i_usableTransports(outModel, trMet, trComp, comps, sc.biomassId, growthConditions); - outModel = i_applyAssignment(model, sc, comps, defaultCompartment, placeIdx, trMet, trComp); + if pruneTransports && ~isempty(trBase) + keep = i_usableTransports(outModel, trComp, comps, sc.biomassId, growthConditions); + trBase = trBase(keep); trComp = trComp(keep); + outModel = i_applyAssignment(model, sc, comps, defaultCompartment, placeIdx, trBase, trComp); + end + if minimizeTransports && ~isempty(trBase) + [trBase, trComp, outModel] = i_minimizeTransports(model, sc, comps, ... + defaultCompartment, placeIdx, trBase, trComp, minGrowth, growthConditions); end exitFlag = 1; report.status = 'certified'; report.certified = true; report.growths = growths; placement.rxns = model.rxns(sc.movIdx); placement.compartment = comps(placeIdx); - addedTransports.mets = sc.metNames(trMet); addedTransports.compartment = comps(trComp); + addedTransports.mets = sc.baseNames(trBase); addedTransports.compartment = comps(trComp); return; end % keep best partial (largest primary growth) if growths.primary > best.growths.primary - best = struct('placeIdx',placeIdx,'trMet',trMet,'trComp',trComp,'growths',growths, ... + best = struct('placeIdx',placeIdx,'trBase',trBase,'trComp',trComp,'growths',growths, ... 'certified',false,'status','uncertified'); end @@ -182,9 +202,9 @@ % honest uncertified result: return the best partial found if ~isempty(best.placeIdx) - outModel = i_applyAssignment(model, sc, comps, defaultCompartment, best.placeIdx, best.trMet, best.trComp); + outModel = i_applyAssignment(model, sc, comps, defaultCompartment, best.placeIdx, best.trBase, best.trComp); placement.rxns = model.rxns(sc.movIdx); placement.compartment = comps(best.placeIdx); - addedTransports.mets = sc.metNames(best.trMet); addedTransports.compartment = comps(best.trComp); + addedTransports.mets = sc.baseNames(best.trBase); addedTransports.compartment = comps(best.trComp); report.growths = best.growths; end report.status = 'uncertified'; @@ -194,8 +214,16 @@ end % ============================================================ scope -function sc = i_prepareScope(model, GSS, relocate, comps, defaultCompartment, minGrowth, transportable, verbose) %#ok -biomassIdx = find(model.c ~= 0, 1); +function sc = i_prepareScope(model, GSS, relocate, comps, minGrowth, transportable, biomassReaction, baseMetabolite) +% biomass reaction: the largest-objective-coefficient reaction (or a named one) +if isempty(biomassReaction) + [~, biomassIdx] = max(model.c); +else + biomassIdx = find(strcmp(model.rxns, biomassReaction), 1); + if isempty(biomassIdx) + error('RAVEN:badInput','biomassReaction ''%s'' not in model.', biomassReaction); + end +end sc.biomassId = model.rxns{biomassIdx}; if isempty(minGrowth) @@ -207,6 +235,13 @@ end sc.minGrowth = minGrowth; +% base metabolite key (identifies the same species across compartments). +% Default 'name': RAVEN gives a metabolite a different id per compartment but a +% shared metName, so the name is the compartment-agnostic key (the equivalent +% of the Python default's compartment-suffix strip on cobra ids). +if isempty(baseMetabolite); baseMetabolite = 'name'; end +if strcmp(baseMetabolite,'name'); sc.base = model.metNames; else; sc.base = model.mets; end + % each reaction's single compartment ('' if boundary/multi-compartment) sc.rxnComp = i_reactionCompartments(model); @@ -214,7 +249,11 @@ isBoundary = (sum(model.S ~= 0, 1)' == 1); movable = toRelocate & ~isBoundary & ~cellfun(@isempty, sc.rxnComp); movable(biomassIdx) = false; +% movable in sorted-reaction-id order (matching the Python master's variable +% order, so an equal-cost placement is broken the same way). sc.movIdx = find(movable); +[~, ord] = sort(model.rxns(sc.movIdx)); +sc.movIdx = sc.movIdx(ord); sc.pinIdx = find(~movable); % genes on movable reactions that are scored @@ -239,13 +278,20 @@ end end -% base metabolite = metabolite name; transportable base names +% base metabolites: the same species across compartments shares a base key, so +% confinement and transports are keyed by base (not by per-compartment row), +% matching Python. baseId maps each metabolite row to its base index. sc.metNames = model.metNames; -movMetMask = any(model.S(:, sc.movIdx) ~= 0, 2); +[sc.baseNames, baseRep, sc.baseId] = unique(sc.base, 'stable'); +sc.nBase = numel(sc.baseNames); +sc.baseRep = baseRep; % a met row per base +movBase = false(sc.nBase,1); +mrows = find(any(model.S(:, sc.movIdx) ~= 0, 2)); +movBase(sc.baseId(mrows)) = true; % bases touched by movable rxns if isnumeric(transportable) && isempty(transportable) - sc.isTransp = movMetMask; % default: all movable bases + sc.baseTransp = movBase; % default: all movable bases else - sc.isTransp = movMetMask & ismember(model.metNames, transportable); + sc.baseTransp = movBase & ismember(sc.baseNames, transportable); end end @@ -360,19 +406,20 @@ movLocal = containers.Map('KeyType','double','ValueType','double'); % rxn idx -> movable local idx for k=1:numel(sc.movIdx); movLocal(sc.movIdx(k)) = k; end -nMet = numel(model.mets); -usedComp = cell(nMet,1); usedRxn = cell(nMet,1); +% used{base} = compartments (and touching reactions) a non-transportable base +% appears in, keyed by base index. +usedComp = cell(sc.nBase,1); usedRxn = cell(sc.nBase,1); allRxn = [sc.movIdx; sc.pinIdx]; for r = allRxn' comp = placedComp(r); if comp == 0; continue; end % multi-compartment reaction: bridges pools - for mm = find(model.S(:,r)~=0)' - if sc.isTransp(mm); continue; end - usedComp{mm}(end+1) = comp; usedRxn{mm}(end+1) = r; + for b = unique(sc.baseId(model.S(:,r)~=0))' + if sc.baseTransp(b); continue; end + usedComp{b}(end+1) = comp; usedRxn{b}(end+1) = r; end end -for mm=1:nMet +for mm=1:sc.nBase cs = unique(usedComp{mm}); if numel(cs) <= 1; continue; end touching = usedRxn{mm}; @@ -393,30 +440,31 @@ end % ============================================================ split transports -function [trMet, trComp] = i_splitTransports(model, sc, comps, defaultCompartment, placeIdx, relaxed) +function [trBase, trComp] = i_splitTransports(model, sc, comps, defaultCompartment, placeIdx, relaxed) +% Returns (baseIndex, compIndex) transports for transportable/relaxed bases a +% placement splits across compartments (routed through defaultCompartment). defC = find(strcmp(comps, defaultCompartment), 1); placedComp = i_placedCompartments(model, sc, comps, placeIdx); allRxn = [sc.movIdx; sc.pinIdx]; -nMet = numel(model.mets); -used = cell(nMet,1); +used = cell(sc.nBase,1); for r = allRxn' comp = placedComp(r); if comp == 0; continue; end - for mm = find(model.S(:,r)~=0)' - if sc.isTransp(mm) || (isKey(relaxed,mm)) - used{mm}(end+1) = comp; + for b = unique(sc.baseId(model.S(:,r)~=0))' + if sc.baseTransp(b) || isKey(relaxed,b) + used{b}(end+1) = comp; end end end -trMet=[]; trComp=[]; -for mm=1:nMet - cs = unique(used{mm}); +trBase=[]; trComp=[]; +for b=1:sc.nBase + cs = unique(used{b}); if numel(cs) <= 1; continue; end for ci = sort(cs) - if ci ~= defC; trMet(end+1)=mm; trComp(end+1)=ci; end %#ok + if ci ~= defC; trBase(end+1)=b; trComp(end+1)=ci; end %#ok end end -trMet=trMet(:); trComp=trComp(:); +trBase=trBase(:); trComp=trComp(:); end % ============================================================ certification @@ -456,32 +504,55 @@ end % ============================================================ transport pruning -function [trMet, trComp] = i_usableTransports(outModel, trMet, trComp, comps, biomassId, growthConditions) -% Keep only transports that can carry flux in at least one certification +function keep = i_usableTransports(outModel, trComp, comps, biomassId, growthConditions) +% Logical mask of transports that can carry flux in at least one certification % medium (sound: an unusable reaction's removal cannot change any of those % FBAs). -trId = cell(numel(trMet),1); -for i=1:numel(trMet) +nTr = numel(trComp); +trId = cell(nTr,1); +for i=1:nTr trId{i} = ['tr_' num2str(i-1) '_' comps{trComp(i)}]; end media = {[]}; if ~isempty(growthConditions) && isstruct(growthConditions) for i=1:numel(growthConditions); media{end+1} = growthConditions(i).medium; end %#ok end -keep = false(numel(trMet),1); +keep = false(nTr,1); for mi=1:numel(media) m = outModel; if ~isempty(media{mi}); m = i_applyMedium(m, media{mi}); end bIdx = find(strcmp(m.rxns, biomassId), 1); m.c = zeros(numel(m.rxns),1); m.c(bIdx) = 1; - idx = zeros(numel(trId),1); - for i=1:numel(trId); j=find(strcmp(m.rxns,trId{i}),1); if ~isempty(j); idx(i)=j; end; end + idx = zeros(nTr,1); + for i=1:nTr; j=find(strcmp(m.rxns,trId{i}),1); if ~isempty(j); idx(i)=j; end; end valid = idx>0; - fl = false(numel(trId),1); + fl = false(nTr,1); fl(valid) = haveFlux(m, 'rxns', idx(valid)); keep = keep | fl; end -trMet = trMet(keep); trComp = trComp(keep); +end + +% ============================================================ transport minimisation +function [trBase, trComp, outModel] = i_minimizeTransports(model, sc, comps, defaultCompartment, placeIdx, trBase, trComp, minGrowth, growthConditions) +% Prune the transports to those carrying flux in a parsimonious-FBA solution, +% then re-certify. Sound: a zero-flux reaction can be removed without changing +% the solution. Falls back to the input set if pFBA is infeasible or the +% pruned set regresses any medium. +outModel = i_applyAssignment(model, sc, comps, defaultCompartment, placeIdx, trBase, trComp); +bIdx = find(strcmp(outModel.rxns, sc.biomassId), 1); +m = outModel; m.c = zeros(numel(m.rxns),1); m.c(bIdx) = 1; +sol = solveLP(m, 1); +if isempty(sol.x); return; end +carry = false(numel(trBase),1); +for i=1:numel(trBase) + j = find(strcmp(outModel.rxns, ['tr_' num2str(i-1) '_' comps{trComp(i)}]), 1); + if ~isempty(j) && abs(sol.x(j)) > 1e-7; carry(i) = true; end +end +newBase = trBase(carry); newComp = trComp(carry); +minModel = i_applyAssignment(model, sc, comps, defaultCompartment, placeIdx, newBase, newComp); +ok = i_certify(minModel, sc.biomassId, minGrowth, growthConditions); +if ~ok; return; end % a needed transport was pruned: keep the safe set +trBase = newBase; trComp = newComp; outModel = minModel; end % ============================================================ growth-gap feedback @@ -512,7 +583,7 @@ end % ============================================================ materialisation -function outModel = i_applyAssignment(model, sc, comps, defaultCompartment, placeIdx, trMet, trComp) +function outModel = i_applyAssignment(model, sc, comps, defaultCompartment, placeIdx, trBase, trComp) outModel = model; for ci = 1:numel(comps) if ~ismember(comps{ci}, outModel.comps) @@ -529,8 +600,8 @@ if ~isequal(newMet, mm); outModel.S(mm, r) = 0; end end end -for k = 1:numel(trMet) - outModel = i_addTransport(outModel, model, trMet(k), comps{trComp(k)}, defaultCompartment, k-1); +for k = 1:numel(trBase) + outModel = i_addTransport(outModel, model, sc.baseRep(trBase(k)), comps{trComp(k)}, defaultCompartment, k-1); end outModel = i_reconcileFields(outModel); end From 588b4fd08026c03f7517842cc89aad9dff635b12 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 18 Jul 2026 08:56:42 +0200 Subject: [PATCH 5/6] feat: loopless FVA, plus universal gap-fill and multi-localization in assignCompartments Adds the two pieces the faithfulness review flagged as unported, and the loopless FVA RAVEN was missing to make one of them possible: - cycleFreeFlux: remove thermodynamically infeasible internal loops from a flux distribution (Desouki et al. 2015) by fixing the exchanges and minimising internal flux. - looplessFVA: loop-free flux variability, via the loop-law MILP of Schellenberger et al. (2011) with metabolite potentials, so a reaction that reaches a bound only around an internal cycle reports the flux it can actually carry loop-free. (cycleFreeFlux on an FBA solution is insufficient for correct bounds; the loop-law MILP is used instead of the cobra cycleFreeFlux variant the Python reference calls.) - assignCompartments now accepts a universal model and gap-fills from it (via fillGaps) when a placement fails to certify, and a multiLocalize option that proposes a second compartment per reaction (gene score >= threshold), materialises the duplicate, and keeps it only if looplessFVA shows it can carry flux at the growth floor. Both re-certify and fall back on failure. Both off by default; the default path is unchanged. Tests for cycleFreeFlux/looplessFVA (a futile-loop model) and for the multi-localization keep path. --- analysis/cycleFreeFlux.m | 65 ++++++++++ analysis/looplessFVA.m | 106 ++++++++++++++++ localization/assignCompartments.m | 127 ++++++++++++++++++- testing/function_tests/tAnalysis.m | 39 ++++++ testing/function_tests/tAssignCompartments.m | 32 +++++ 5 files changed, 365 insertions(+), 4 deletions(-) create mode 100644 analysis/cycleFreeFlux.m create mode 100644 analysis/looplessFVA.m diff --git a/analysis/cycleFreeFlux.m b/analysis/cycleFreeFlux.m new file mode 100644 index 00000000..76cc5a86 --- /dev/null +++ b/analysis/cycleFreeFlux.m @@ -0,0 +1,65 @@ +function fluxOut = cycleFreeFlux(model, fluxes) +% cycleFreeFlux Remove thermodynamically infeasible internal loops from a flux +% distribution. +% +% Given a steady-state flux distribution, return an equivalent one that has the +% same exchange (boundary) fluxes but the smallest possible total internal flux +% and no net flux around internal loops. Method of Desouki et al. (2015), +% Bioinformatics 31:2159 (the LP-based loop removal cobra uses for its +% loopless='cycleFreeFlux' option): fix the boundary reactions and the sign of +% every internal reaction, then minimise the sum of absolute internal fluxes. +% A flux carried only to close a loop is removed because it adds to that sum +% without changing any exchange. +% +% Parameters +% ---------- +% model : struct +% a RAVEN model structure. +% fluxes : double +% a flux vector (one value per reaction) to make cycle-free. +% +% Returns +% ------- +% fluxOut : double +% the loop-free flux vector, with the same boundary fluxes as the input. +% +% See also +% -------- +% looplessFVA, getExchangeRxns, solveLP + +tol = 1e-9; +nRxns = numel(model.rxns); +fluxes = fluxes(:); + +m = model; +% Fix the boundary (exchange) reactions to their input values so the network's +% exchange with the environment is preserved. +[~, exchIdx] = getExchangeRxns(m); +isExch = false(nRxns,1); isExch(exchIdx) = true; +m.lb(exchIdx) = fluxes(exchIdx); +m.ub(exchIdx) = fluxes(exchIdx); + +% Fix the sign of every internal reaction (so the cleaned solution keeps the +% same directions), and cap the magnitude at the input flux. +internal = find(~isExch); +for i = internal' + v = fluxes(i); + if v > tol + m.lb(i) = max(m.lb(i), 0); m.ub(i) = min(m.ub(i), v); + elseif v < -tol + m.ub(i) = min(m.ub(i), 0); m.lb(i) = max(m.lb(i), v); + else + m.lb(i) = 0; m.ub(i) = 0; + end +end + +% Minimise the total internal flux (no objective; solveLP's minFlux pass +% minimises the sum of absolute fluxes, and the boundary fluxes are fixed). +m.c = zeros(nRxns,1); +sol = solveLP(m, 1); +if isempty(sol.x) + fluxOut = fluxes; % infeasible under the fixings: return input +else + fluxOut = sol.x; +end +end diff --git a/analysis/looplessFVA.m b/analysis/looplessFVA.m new file mode 100644 index 00000000..6eafeb03 --- /dev/null +++ b/analysis/looplessFVA.m @@ -0,0 +1,106 @@ +function [minFlux, maxFlux] = looplessFVA(model, rxns, minGrowth) +% looplessFVA Loop-free flux variability analysis. +% +% Minimum and maximum flux of each listed reaction over the steady-state flux +% space, excluding any flux that could exist only around a thermodynamically +% infeasible internal loop. Uses the loop-law formulation of Schellenberger et +% al. (2011), Biophys J 100:544: a metabolite-potential vector mu assigns each +% internal reaction an energy G = S'*mu that must oppose its flux direction, so +% no internal cycle can carry net flux. Each reaction's min/max is then a MILP. +% +% Parameters +% ---------- +% model : struct +% a RAVEN model structure. +% rxns : cell or double +% reaction ids or indexes to analyse (default all reactions). +% minGrowth : double +% if given, the model objective is held at or above this value while the +% variability is computed (default []: the objective is not constrained). +% +% Returns +% ------- +% minFlux : double +% the loop-free minimum flux of each reaction in rxns. +% maxFlux : double +% the loop-free maximum flux of each reaction in rxns. +% +% See also +% -------- +% cycleFreeFlux, getAllowedBounds, solveLP + +if nargin < 2 || isempty(rxns) + idx = (1:numel(model.rxns))'; +else + idx = getIndexes(model, rxns, 'rxns'); +end +if nargin < 3; minGrowth = []; end + +M = 1000; +nRxns = numel(model.rxns); +nMets = numel(model.mets); +% Internal reactions get the loop law; exchange reactions cannot be in an +% internal loop, so they are exempt. +[~, exchIdx] = getExchangeRxns(model); +isExch = false(nRxns,1); isExch(exchIdx) = true; +intIdx = find(~isExch); +nInt = numel(intIdx); + +% Variables: v (nRxns), mu (nMets), z (nInt binary). +oV = 0; oMu = nRxns; oZ = nRxns + nMets; nVar = nRxns + nMets + nInt; + +% --- constraints --- +% S v = 0 +A_S = [model.S, sparse(nMets, nMets + nInt)]; +% per-internal loop-law rows +rA = []; cA = []; vA = []; bDir = []; senseDir = ''; row = 0; +for k = 1:nInt + gi = intIdx(k); zc = oZ + k; + lbi = model.lb(gi); ubi = model.ub(gi); + % (a) v_gi - ub*z <= 0 + row=row+1; rA(end+1)=row; cA(end+1)=gi; vA(end+1)=1; rA(end+1)=row; cA(end+1)=zc; vA(end+1)=-ubi; bDir(end+1)=0; senseDir(end+1)='L'; %#ok + % (b) v_gi + lb*z >= lb + row=row+1; rA(end+1)=row; cA(end+1)=gi; vA(end+1)=1; rA(end+1)=row; cA(end+1)=zc; vA(end+1)=lbi; bDir(end+1)=lbi; senseDir(end+1)='G'; %#ok + % (c) S(:,gi)'mu + M*z <= M-1 + mrows = find(model.S(:,gi)); + for mm=mrows'; row0=row+1; rA(end+1)=row0; cA(end+1)=oMu+mm; vA(end+1)=model.S(mm,gi); end %#ok + row=row+1; rA(end+1)=row; cA(end+1)=zc; vA(end+1)=M; bDir(end+1)=M-1; senseDir(end+1)='L'; %#ok + % (d) S(:,gi)'mu + M*z >= 1 + for mm=mrows'; row0=row+1; rA(end+1)=row0; cA(end+1)=oMu+mm; vA(end+1)=model.S(mm,gi); end %#ok + row=row+1; rA(end+1)=row; cA(end+1)=zc; vA(end+1)=M; bDir(end+1)=1; senseDir(end+1)='G'; %#ok +end +A_dir = sparse(rA, cA, vA, row, nVar); + +A = [A_S; A_dir]; +b = [zeros(nMets,1); bDir(:)]; +csense = [repmat('E',1,nMets), senseDir]; + +% objective (biomass) floor +if ~isempty(minGrowth) + objIdx = find(model.c ~= 0); + for o=objIdx' + A(end+1,:) = sparse(1, oV+o, 1, 1, nVar); %#ok + b(end+1,1) = minGrowth; csense(end+1) = 'G'; %#ok + end +end + +lb = [model.lb(:); -M*ones(nMets,1); zeros(nInt,1)]; +ub = [model.ub(:); M*ones(nMets,1); ones(nInt,1)]; +vartype = [repmat('C',1,nRxns+nMets), repmat('B',1,nInt)]; + +prob.a = A; prob.A = A; prob.b = b; prob.csense = csense; +prob.lb = lb; prob.ub = ub; prob.vartype = vartype; prob.osense = 1; +params.intTol = 1e-9; params.TimeLimit = 1000; + +minFlux = zeros(numel(idx),1); maxFlux = zeros(numel(idx),1); +for k = 1:numel(idx) + r = idx(k); + c = zeros(nVar,1); c(oV+r) = 1; + prob.c = -c; % maximise v_r + sol = optimizeProb(prob, params, false); + if checkSolution(sol); maxFlux(k) = sol.full(oV+r); end + prob.c = c; % minimise v_r + sol = optimizeProb(prob, params, false); + if checkSolution(sol); minFlux(k) = sol.full(oV+r); end +end +end diff --git a/localization/assignCompartments.m b/localization/assignCompartments.m index 8ab57ed2..3af4f17f 100644 --- a/localization/assignCompartments.m +++ b/localization/assignCompartments.m @@ -60,14 +60,23 @@ % baseMetabolite : char (default 'name') % the compartment-agnostic metabolite key: 'name' (metNames, the RAVEN % convention) or 'id' (model.mets). +% universal : struct (default []) +% a template model; if a placement fails to certify because the primary +% medium falls short, gaps are filled from it (via fillGaps). +% multiLocalize : logical (default false) +% after certifying, add multi-compartment placements: for each placed +% reaction, propose a second compartment its genes score highest for, +% materialise the duplicate, and keep it only if a loopless FVA +% (looplessFVA) shows it can carry flux in a biomass-supporting solution. +% multiLocalizeThreshold : double (default 0.7) +% the gene score a second compartment needs for a duplicate to be proposed. +% multiLocalizeEps : double (default 1e-6) +% the loopless flux a proposed duplicate must reach to be kept. % transportCost : double (default 0.5) % accepted for signature compatibility but NOT used (placement is % flux-free and score-only). % verbose : logical (default true) % -% Not ported from the raven-toolbox reference: universal-model gap-fill, and -% multi-localization (which requires a loopless FVA RAVEN does not provide). -% % Returns % ------- % outModel : struct @@ -88,7 +97,9 @@ p = parseRAVENargs(varargin, {'defaultCompartment',[]; 'transportCost',0.5; ... 'multiCompartmentPenalty',0.5; 'minGrowth',[]; 'transportable',[]; ... 'growthConditions',[]; 'maxRounds',8; 'pruneTransports',true; ... - 'minimizeTransports',false; 'biomassReaction',[]; 'baseMetabolite',[]; 'verbose',true}); + 'minimizeTransports',false; 'biomassReaction',[]; 'baseMetabolite',[]; ... + 'universal',[]; 'multiLocalize',false; 'multiLocalizeThreshold',0.7; ... + 'multiLocalizeEps',1e-6; 'verbose',true}); defaultCompartment = char(p.defaultCompartment); multiPen = p.multiCompartmentPenalty; minGrowth = p.minGrowth; @@ -96,6 +107,10 @@ maxRounds = p.maxRounds; pruneTransports = p.pruneTransports; minimizeTransports = p.minimizeTransports; +universal = p.universal; +multiLocalize = p.multiLocalize; +mlThreshold = p.multiLocalizeThreshold; +mlEps = p.multiLocalizeEps; verbose = p.verbose; outModel = model; @@ -169,6 +184,14 @@ outModel = i_applyAssignment(model, sc, comps, defaultCompartment, placeIdx, trBase, trComp); [ok, growths] = i_certify(outModel, sc.biomassId, minGrowth, growthConditions); + % feedback: gap-fill from the universal model if the primary medium fell short + if ~ok && ~isempty(universal) && growths.primary < minGrowth - 1e-9 + gf = i_gapfill(outModel, universal, sc.biomassId, minGrowth); + if ~isempty(gf.rxns) + outModel = gf; [ok, growths] = i_certify(outModel, sc.biomassId, minGrowth, growthConditions); + end + end + if ok if pruneTransports && ~isempty(trBase) keep = i_usableTransports(outModel, trComp, comps, sc.biomassId, growthConditions); @@ -179,7 +202,13 @@ [trBase, trComp, outModel] = i_minimizeTransports(model, sc, comps, ... defaultCompartment, placeIdx, trBase, trComp, minGrowth, growthConditions); end + multiLoc = {}; + if multiLocalize + [outModel, multiLoc] = i_enrichMultiloc(model, sc, comps, defaultCompartment, ... + placeIdx, outModel, minGrowth, mlThreshold, mlEps, growthConditions); + end exitFlag = 1; report.status = 'certified'; report.certified = true; report.growths = growths; + report.multiLocalized = multiLoc; placement.rxns = model.rxns(sc.movIdx); placement.compartment = comps(placeIdx); addedTransports.mets = sc.baseNames(trBase); addedTransports.compartment = comps(trComp); return; @@ -532,6 +561,96 @@ end end +% ============================================================ universal gap-fill +function outModel = i_gapfill(outModel, universal, biomassId, minGrowth) +% Fill gaps from the universal model to restore the growth floor, using +% RAVEN's fillGaps. The added reactions are validated by the caller's +% certification FBA. On any failure the model is returned unchanged. +outModel.rxns = outModel.rxns; % ensure struct is returned even on failure +try + bIdx = find(strcmp(outModel.rxns, biomassId), 1); + m = outModel; m.lb(bIdx) = minGrowth; % require growth while filling + [~, ~, addedRxns, newModel] = evalc('fillGaps(m, {universal}, ''useModelConstraints'', true, ''minGrowth'', minGrowth, ''verbose'', false)'); + if ~isempty(addedRxns) + newModel.lb(strcmp(newModel.rxns, biomassId)) = outModel.lb(bIdx); % restore biomass bound + outModel = newModel; + end +catch + % gap-fill infeasible or unavailable: leave the model unchanged +end +end + +% ============================================================ multi-localisation +function [outModel, multiLoc] = i_enrichMultiloc(model, sc, comps, defaultCompartment, placeIdx, outModel, minGrowth, threshold, eps, growthConditions) +% Add sound multi-compartment placements to a certified mono model: propose, +% per placed reaction, a second compartment its genes score >= threshold for; +% materialise every candidate as a duplicate; keep only those a loopless FVA +% shows can carry flux >= eps at the growth floor. Re-certifies; on failure or +% no survivor, returns the mono model. +multiLoc = {}; +cand = []; % [movableLocalIdx, compIdx] +for k=1:numel(sc.movIdx) + genes = find(model.rxnGeneMat(sc.movIdx(k),:) ~= 0); + gl = arrayfun(@(g) find(sc.geneIdx==g,1), genes, 'UniformOutput', false); + gl = [gl{:}]; + if isempty(gl); continue; end + for ci=1:numel(comps) + if ci == placeIdx(k); continue; end + if max(sc.score(gl, ci)) >= threshold + cand = [cand; k, ci]; %#ok + end + end +end +if isempty(cand); return; end + +% materialise all candidate duplicates onto the certified model +trial = outModel; dupId = cell(size(cand,1),1); +for i=1:size(cand,1) + r = sc.movIdx(cand(i,1)); comp = comps{cand(i,2)}; + [trial, dupId{i}] = i_duplicateReaction(trial, model, r, comp); +end +trial = i_reconcileFields(trial); + +% loopless FVA over the duplicates at the growth floor; keep flux-carriers +present = find(~cellfun(@isempty, dupId)); +if isempty(present); return; end +ids = dupId(present); +[lo, hi] = looplessFVA(trial, ids, minGrowth); +keep = present(max(abs(lo), abs(hi)) >= eps); +if isempty(keep); return; end + +% rebuild the model keeping only surviving duplicates, re-certify +enriched = outModel; +for i=keep' + r = sc.movIdx(cand(i,1)); comp = comps{cand(i,2)}; + enriched = i_duplicateReaction(enriched, model, r, comp); + multiLoc{end+1,1} = {model.rxns{r}, comp}; %#ok +end +enriched = i_reconcileFields(enriched); +if i_certify(enriched, sc.biomassId, minGrowth, growthConditions) + outModel = enriched; +else + multiLoc = {}; % safety net: duplicates only add capability, keep mono +end +end + +function [outModel, dupId] = i_duplicateReaction(outModel, model, r, comp) +dupId = [model.rxns{r} '_' comp]; +if any(strcmp(outModel.rxns, dupId)); dupId = ''; return; end +outModel.rxns{end+1,1} = dupId; +outModel.S(:, end+1) = 0; +for mm = find(model.S(:,r) ~= 0)' + [outModel, newMet] = i_metInComp(outModel, model, mm, comp); + outModel.S(newMet, end) = model.S(mm, r); +end +outModel.lb(end+1,1) = model.lb(r); outModel.ub(end+1,1) = model.ub(r); +if isfield(outModel,'rev'); outModel.rev(end+1,1) = model.rev(r); end +if isfield(outModel,'c'); outModel.c(end+1,1) = 0; end +if isfield(outModel,'rxnNames'); outModel.rxnNames{end+1,1} = dupId; end +if isfield(outModel,'grRules'); outModel.grRules{end+1,1} = model.grRules{r}; end +if isfield(outModel,'rxnGeneMat'); outModel.rxnGeneMat(end+1,:) = 0; end +end + % ============================================================ transport minimisation function [trBase, trComp, outModel] = i_minimizeTransports(model, sc, comps, defaultCompartment, placeIdx, trBase, trComp, minGrowth, growthConditions) % Prune the transports to those carrying flux in a parsimonious-FBA solution, diff --git a/testing/function_tests/tAnalysis.m b/testing/function_tests/tAnalysis.m index 6d27df31..baedd67c 100644 --- a/testing/function_tests/tAnalysis.m +++ b/testing/function_tests/tAnalysis.m @@ -21,6 +21,29 @@ function haveFluxReturnsLogical(testCase) testCase.verifyNumElements(I, numel(testCase.model.rxns)); end + function cycleFreeFluxRemovesLoop(testCase) + % R1/R2 form a futile a<->b loop; the input carries 1000 around it + % with a net through-flux of 10. cycleFreeFlux must strip the loop + % (R1 -> 10, R2 -> 0) while preserving the exchange fluxes. + m = tAnalysis.loopModel(); + cf = cycleFreeFlux(m, [10; 1000; 990; 10]); + testCase.verifyEqual(cf(1), 10, 'AbsTol', 1e-6); % R_in preserved + testCase.verifyEqual(cf(4), 10, 'AbsTol', 1e-6); % R_out preserved + testCase.verifyEqual(cf(2), 10, 'AbsTol', 1e-6); % R1 net only + testCase.verifyLessThan(abs(cf(3)), 1e-6); % R2 loop flux gone + end + + function looplessFVAExcludesLoopFlux(testCase) + % Standard FVA lets R1 reach 1000 around the loop; loopless FVA + % must cap it at the net 10, and R2 (only ever loop flux) at 0. + testCase.assumeMILPSolver(); + m = tAnalysis.loopModel(); + [lo, hi] = looplessFVA(m, {'R1';'R2'}); + testCase.verifyEqual(hi(1), 10, 'AbsTol', 1e-4); + testCase.verifyLessThan(abs(hi(2)), 1e-4); + testCase.verifyGreaterThanOrEqual(lo(1), -1e-4); + end + function getMinNrFluxesReturnsFlux(testCase) testCase.assumeMILPSolver(); evalc('[x, I, exitFlag] = getMinNrFluxes(testCase.model, testCase.model.rxns);'); @@ -229,4 +252,20 @@ function runSimpleOptKnockRuns(testCase) end end + + methods (Static, Access = private) + function m = loopModel() + % R_in -> a; R1: a->b; R2: b->a (futile loop); R_out: b ->. + m = struct(); + m.id='loop'; m.comps={'c'}; m.compNames={'cyt'}; + m.mets={'a';'b'}; m.metNames={'a';'b'}; m.metComps=[1;1]; + m.S = sparse([ 1 -1 1 0; + 0 1 -1 -1]); + m.rxns={'R_in';'R1';'R2';'R_out'}; m.rxnNames=m.rxns; + m.lb=[0;0;0;0]; m.ub=[10;1000;1000;1000]; m.rev=[0;0;0;0]; + m.c=[0;0;0;1]; m.b=zeros(2,1); + m.genes={}; m.grRules={'';'';'';''}; m.rxnGeneMat=sparse(4,0); + m.metFormulas={'C';'C'}; + end + end end diff --git a/testing/function_tests/tAssignCompartments.m b/testing/function_tests/tAssignCompartments.m index c0236fd5..d4ef7deb 100644 --- a/testing/function_tests/tAssignCompartments.m +++ b/testing/function_tests/tAssignCompartments.m @@ -93,6 +93,20 @@ function confinementColocatesMovableReactions(testCase) sp.compartment{strcmp(sp.rxns,'r2')}); end + function multiLocalizeKeepsFluxCarryingDuplicate(testCase) + % A reaction whose gene scores high in a second compartment, and + % which can carry flux there (independent uptake/drain), is + % duplicated into it; the model stays certified. + testCase.assumeMILPSolver(); + model = tAssignCompartments.dualToy(); + GSS = struct('genes',{{'g1'}},'compartments',{{'c';'m'}},'scores',[0.9 0.9]); + evalc(['[oM, ~, ~, ok, rep] = assignCompartments(model, GSS, {''R''}, ' ... + '''defaultCompartment'', ''c'', ''multiLocalize'', true, ''verbose'', false);']); + testCase.verifyEqual(ok, 1); + testCase.verifyNumElements(rep.multiLocalized, 1); + testCase.verifyTrue(any(strcmp(oM.rxns, 'R_m'))); + end + end methods (Static) @@ -124,6 +138,24 @@ function confinementColocatesMovableReactions(testCase) model.genes={'g1';'g2'}; model.grRules={'';'g1';'g2';''}; model.rxnGeneMat=sparse([0 0; 1 0; 0 1; 0 0]); end + function model = dualToy() + % A/B exist in c and m; R (A->B) can run in either compartment, + % each fed and drained independently, so a duplicate carries flux. + model = struct(); model.id='dual'; + model.comps={'c';'m'}; model.compNames={'c';'m'}; + model.mets={'A_c';'B_c';'A_m';'B_m'}; model.metNames={'A';'B';'A';'B'}; + model.metComps=[1;1;2;2]; + % EXAc EXAm R bio EXBm + model.S = sparse([ 1 0 -1 0 0; % A_c + 0 0 1 -1 0; % B_c + 0 1 0 0 0; % A_m + 0 0 0 0 -1]); % B_m + model.rxns={'EX_A_c';'EX_A_m';'R';'bio';'EX_B_m'}; model.rxnNames=model.rxns; + model.lb=[0;0;0;0;0]; model.ub=[10;10;1000;1000;1000]; model.rev=[0;0;0;0;0]; + model.c=[0;0;0;1;0]; model.b=zeros(4,1); + model.genes={'g1'}; model.grRules={'';'';'g1';'';''}; + model.rxnGeneMat=sparse([0;0;1;0;0]); model.metFormulas={'C';'C';'C';'C'}; + end function GSS = chainGss() % Each gene has a single dominant compartment (the other score is % below the 0.5 multi-compartment penalty), so placement is From 735e8505707a7d5f262d0c9f922e9bd08cd68771 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 18 Jul 2026 22:12:51 +0200 Subject: [PATCH 6/6] refactor: make the placement master deterministic and matched to the port Iterate genes in sorted-id order and emit the placement MILP constraint rows in the same interleaved order as the raven-toolbox master (per movable: its place row then its couple rows; per gene: its gene1 row then its has rows). Pin the solve to a single thread, a fixed seed and a zero MIP gap so a degenerate score objective resolves to the same co-optimal placement on every run. Column order, row order and parameters now match the Python master, which return byte-identical placements on yeast-GEM. --- localization/assignCompartments.m | 89 ++++++++++++++++--------------- 1 file changed, 45 insertions(+), 44 deletions(-) diff --git a/localization/assignCompartments.m b/localization/assignCompartments.m index 3af4f17f..6e30346e 100644 --- a/localization/assignCompartments.m +++ b/localization/assignCompartments.m @@ -290,6 +290,10 @@ geneOnMov = any(model.rxnGeneMat(sc.movIdx, :) ~= 0, 1)'; scope = gInGSS & geneOnMov; sc.geneIdx = find(scope); +% genes in sorted-id order, so the placement MILP is built in a canonical +% order (identical to the Python master) and the solver sees the same problem. +[~, go] = sort(model.genes(sc.geneIdx)); +sc.geneIdx = sc.geneIdx(go); sc.score = zeros(numel(sc.geneIdx), numel(comps)); % align score columns (GSS.compartments) to comps [~, colOf] = ismember(GSS.compartments, comps); @@ -336,80 +340,77 @@ xCol = @(mi,ci) oX + (mi-1)*nC + ci; yCol = @(gi,ci) oY + (gi-1)*nC + ci; -pR=[];pC=[];pV=[]; +% Constraint rows are emitted in the exact order the Python master builds them, +% because the row order (like the column order) selects which of the degenerate +% co-optimal placements the solver returns -- a permutation of the rows changes +% the vertex. Order: per movable, its place row then its couple rows (in-scope +% genes on the reaction, in sorted-id order, x compartments); then per gene, its +% gene1 row then its has rows (x compartments); then the forced and colocation +% rows. Assembled as one triplet list with a running row index. +aR=[];aC=[];aV=[]; bVec=[]; cs=''; row=0; for k=1:nMov - for ci=1:nC; pR(end+1)=k; pC(end+1)=xCol(k,ci); pV(end+1)=1; end %#ok -end -A_place=sparse(pR,pC,pV,nMov,nVar); b_place=ones(nMov,1); - -cR=[];cC=[];cV=[];row=0; -for k=1:nMov - gOn=find(model.rxnGeneMat(movIdx(k),:)~=0); - for g=gOn - gi=find(geneIdx==g,1); if isempty(gi); continue; end - for ci=1:nC + row=row+1; % place: sum_c x[k,c] = 1 + for ci=1:nC; aR(end+1)=row; aC(end+1)=xCol(k,ci); aV(end+1)=1; end %#ok + bVec(end+1)=1; cs(end+1)='E'; %#ok + gOn=find(model.rxnGeneMat(movIdx(k),geneIdx)~=0); % in-scope genes on k, sorted-id order + for gi=gOn + for ci=1:nC % couple: x[k,c] - y[gi,c] <= 0 row=row+1; - cR(end+1)=row; cC(end+1)=xCol(k,ci); cV(end+1)=1; %#ok - cR(end+1)=row; cC(end+1)=yCol(gi,ci); cV(end+1)=-1; %#ok + aR(end+1)=row; aC(end+1)=xCol(k,ci); aV(end+1)=1; %#ok + aR(end+1)=row; aC(end+1)=yCol(gi,ci); aV(end+1)=-1; %#ok + bVec(end+1)=0; cs(end+1)='L'; %#ok end end end -A_couple=sparse(cR,cC,cV,row,nVar); b_couple=zeros(row,1); - -hR=[];hC=[];hV=[];row=0; for gi=1:nGene - g=geneIdx(gi); - rOfG=find(model.rxnGeneMat(movIdx,g)~=0)'; - for ci=1:nC - row=row+1; hR(end+1)=row; hC(end+1)=yCol(gi,ci); hV(end+1)=1; %#ok - for k=rOfG; hR(end+1)=row; hC(end+1)=xCol(k,ci); hV(end+1)=-1; end %#ok + row=row+1; % gene1: sum_c y[gi,c] >= 1 + for ci=1:nC; aR(end+1)=row; aC(end+1)=yCol(gi,ci); aV(end+1)=1; end %#ok + bVec(end+1)=1; cs(end+1)='G'; %#ok + rOfG=find(model.rxnGeneMat(movIdx,geneIdx(gi))~=0)'; % movable reactions on gi, movable order + for ci=1:nC % has: y[gi,c] - sum_r x[r,c] <= 0 + row=row+1; aR(end+1)=row; aC(end+1)=yCol(gi,ci); aV(end+1)=1; %#ok + for k=rOfG; aR(end+1)=row; aC(end+1)=xCol(k,ci); aV(end+1)=-1; end %#ok + bVec(end+1)=0; cs(end+1)='L'; %#ok end end -A_has=sparse(hR,hC,hV,row,nVar); b_has=zeros(row,1); - -g1R=[];g1C=[];g1V=[]; -for gi=1:nGene - for ci=1:nC; g1R(end+1)=gi; g1C(end+1)=yCol(gi,ci); g1V(end+1)=1; end %#ok -end -A_gene1=sparse(g1R,g1C,g1V,nGene,nVar); b_gene1=ones(nGene,1); - -fR=[];fC=[];nForce=0; fkeys=keys(forced); +fkeys=keys(forced); for i=1:numel(fkeys) - k=fkeys{i}; - nForce=nForce+1; fR(end+1)=nForce; fC(end+1)=xCol(k,forced(k)); %#ok + k=fkeys{i}; % force: x[k,forced(k)] = 1 + row=row+1; aR(end+1)=row; aC(end+1)=xCol(k,forced(k)); aV(end+1)=1; %#ok + bVec(end+1)=1; cs(end+1)='E'; %#ok end -A_force=sparse(fR,fC,ones(1,nForce),nForce,nVar); b_force=ones(nForce,1); - -coR=[];coC=[];coV=[];row=0; for gi=1:numel(groups) mem=groups{gi}; for j=1:numel(mem)-1 a=mem(j); b=mem(j+1); - for ci=1:nC + for ci=1:nC % colo: x[a,c] - x[b,c] = 0 row=row+1; - coR(end+1)=row; coC(end+1)=xCol(a,ci); coV(end+1)=1; %#ok - coR(end+1)=row; coC(end+1)=xCol(b,ci); coV(end+1)=-1; %#ok + aR(end+1)=row; aC(end+1)=xCol(a,ci); aV(end+1)=1; %#ok + aR(end+1)=row; aC(end+1)=xCol(b,ci); aV(end+1)=-1; %#ok + bVec(end+1)=0; cs(end+1)='E'; %#ok end end end -A_colo=sparse(coR,coC,coV,row,nVar); b_colo=zeros(row,1); c = zeros(nVar,1); for gi=1:nGene for ci=1:nC; c(yCol(gi,ci)) = score(gi,ci) - multiPen; end end -prob.A = [A_place; A_couple; A_has; A_gene1; A_force; A_colo]; +prob.A = sparse(aR,aC,aV,row,nVar); prob.a = prob.A; -prob.b = [b_place; b_couple; b_has; b_gene1; b_force; b_colo]; -prob.csense = [repmat('E',1,nMov), repmat('L',1,numel(b_couple)), ... - repmat('L',1,numel(b_has)), repmat('G',1,nGene), ... - repmat('E',1,nForce), repmat('E',1,numel(b_colo))]; +prob.b = bVec(:); +prob.csense = cs; prob.c = -c; prob.osense = 1; prob.lb = zeros(nVar,1); prob.ub = ones(nVar,1); prob.vartype = repmat('B',1,nVar); params.intTol = 1e-9; params.TimeLimit = 1000; +% deterministic, solver-independent solve: single thread and fixed seed remove +% thread-count/seed nondeterminism, MIPGap 0 forces the exact optimum rather +% than an early heuristic stop. With a canonically ordered problem this makes +% the placement reproducible and identical to the Python master. +params.Threads = 1; params.Seed = 0; params.MIPGap = 0; sol = optimizeProb(prob, params, verbose); if ~checkSolution(sol) placeIdx = []; return;