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 37cb08de..6e30346e 100644 --- a/localization/assignCompartments.m +++ b/localization/assignCompartments.m @@ -1,353 +1,750 @@ -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. +% 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: % -% 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. 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. +% +% 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. +% 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. +% gene scoring structure (genes, compartments, scores) as from parseScores. % 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. -% transportCost : double (default 0.5) -% cost per added inter-compartment transport. +% compartment transports route through (usually cytosol); must be in the +% union of model and GSS compartments. % multiCompartmentPenalty : double (default 0.5) -% 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. -% 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. +% required objective flux for certification. +% transportable : cell (default [] = all movable base metabolites) +% metabolite NAMES that may receive transports. +% growthConditions : struct array (default []) +% extra media to certify on, each with fields name, medium +% (struct exchangeRxnId -> max uptake) and minGrowth. +% maxRounds : double (default 8) +% 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). +% 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) % % 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 (cell of compartment ids per placed reaction). % addedTransports : struct -% .mets and .compartment for the transports the MILP added. +% .mets (base metabolite names) and .compartment for the added transports. % exitFlag : double -% 1 = optimal, -1 = infeasible/failed. +% 1 = certified on every medium, -1 = could not place or did not certify. +% report : struct +% .certified, .status, .growths (struct medium -> flux), .unplaced. % % 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; 'pruneTransports',true; ... + 'minimizeTransports',false; 'biomassReaction',[]; 'baseMetabolite',[]; ... + 'universal',[]; 'multiLocalize',false; 'multiLocalizeThreshold',0.7; ... + 'multiLocalizeEps',1e-6; 'verbose',true}); defaultCompartment = char(p.defaultCompartment); -transportCost = p.transportCost; multiPen = p.multiCompartmentPenalty; minGrowth = p.minGrowth; -M = p.bigM; +growthConditions = p.growthConditions; +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; placement = struct('rxns',{{}},'compartment',{{}}); -addedTransports = struct('mets',{{}},'compartment',{{}}); exitFlag = -1; +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); + error('RAVEN:badInput','defaultCompartment ''%s'' not in the model/score compartments.', defaultCompartment); end -% Merge to a single compartment so every metabolite has one identity; relocations and transports -% are then expressed relative to defaultCompartment. -if numel(model.comps) > 1 - model = mergeCompartments(model, true, true); +% ---- scope: biomass, growth floor, movable/pinned, genes, transportable ---- +sc = i_prepareScope(model, GSS, reactionsToRelocate, comps, minGrowth, p.transportable, p.biomassReaction, p.baseMetabolite); +minGrowth = sc.minGrowth; +report.unplaced = sc.unplaced; + +if verbose + fprintf('assignCompartments: %d movable, %d pinned reactions, %d genes, %d compartments.\n', ... + numel(sc.movIdx), numel(sc.pinIdx), numel(sc.geneIdx), numel(comps)); 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.'); +% ---- 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 = {}; +best = struct('placeIdx',[],'trBase',[],'trComp',[],'growths',struct('primary',-1), ... + 'certified',false,'status','uncertified'); + +for round = 1:maxRounds + sig = i_signature(forced, groups); + if any(strcmp(seen, sig)) + break; % configuration already tried: diagnosers are cycling end - minGrowth = 0.1 * sol.f; -end + seen{end+1} = sig; %#ok -% ---- 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); nPin = numel(pinIdx); nC = numel(comps); -defC = find(strcmp(comps, defaultCompartment)); + placeIdx = i_placementMaster(model, sc, comps, multiPen, forced, groups, verbose); + if isempty(placeIdx) + report.status = 'infeasible'; + return; + end -% 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); -nGene = numel(geneIdx); -score = zeros(nGene, nC); % gene x compartment score, columns aligned to comps -for gi = 1:nGene - score(gi, :) = GSS.scores(gssRow(geneIdx(gi)), :); -end - -% Transportable metabolites (touched by movable reactions) -movMetMask = any(model.S(:, movIdx) ~= 0, 2); -if isnumeric(p.transportable) && isempty(p.transportable) - transpMet = find(movMetMask); % default ([]): all movable metabolites transportable -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 -end -nTr = size(trPairs, 1); + % confinement fixpoint + [newForced, newGroups, relaxed] = i_diagnoseConfinement(model, sc, comps, placeIdx, gapPinned); + changed = false; + 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 + for gi = 1:numel(newGroups) + if ~i_hasGroup(groups, newGroups{gi}) + groups{end+1} = newGroups{gi}; changed = true; %#ok + end + end + if changed + continue; + end -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); + % transports + materialise + certify + [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); + + % 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 -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 + + if ok + 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 + 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; 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 + + % keep best partial (largest primary growth) + if growths.primary > best.growths.primary + best = struct('placeIdx',placeIdx,'trBase',trBase,'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 + +% honest uncertified result: return the best partial found +if ~isempty(best.placeIdx) + 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.baseNames(best.trBase); 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 -% 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 +% ============================================================ scope +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 -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 +sc.biomassId = model.rxns{biomassIdx}; + +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 -A_glb = sparse(giR,giC,giV,row,nVar); b_glb=zeros(row,1); +sc.minGrowth = minGrowth; -% ---- 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 +% 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); + +toRelocate = ismember(model.rxns, relocate); +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 +[gInGSS, gssRow] = ismember(model.genes, GSS.genes); +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); +for gi = 1:numel(sc.geneIdx) + row = GSS.scores(gssRow(sc.geneIdx(gi)), :); + sc.score(gi, colOf(colOf>0)) = row(colOf>0); 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 + +% 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 -A_tlb = sparse(tR,tC,tV,nTr,nVar); % ftr + M*t >= 0 -% ---- placement: sum_c x[mov,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 +% 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; +[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.baseTransp = movBase; % default: all movable bases +else + sc.baseTransp = movBase & ismember(sc.baseNames, transportable); end -A_place=sparse(pR,pC,pV,nMov,nVar); b_place=ones(nMov,1); +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; -% ---- gene coupling x[r,c] <= y[g,c]; gene-has y<=sum x; gene1 sum_c y>=1 ---- -cR=[];cC=[];cV=[];row=0; +% 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 - 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; - cR(end+1)=row; cC(end+1)=yCol(gi,ci); cV(end+1)=-1; %#ok + 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; + 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); -% gene-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 - for ci=1:nC - row=row+1; hR(end+1)=row; hC(end+1)=yCol(gi,ci); hV(end+1)=1; - 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); -% 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 +fkeys=keys(forced); +for i=1:numel(fkeys) + 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 +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 % colo: x[a,c] - x[b,c] = 0 + row=row+1; + 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_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; -% ---- objective: max sum score*y - multiPen*sum y - transportCost*sum t ---- 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 = sparse(aR,aC,aV,row,nVar); 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.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.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) - 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 + +% ============================================================ 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{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 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:sc.nBase + cs = unique(usedComp{mm}); + 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 [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]; +used = cell(sc.nBase,1); +for r = allRxn' + comp = placedComp(r); + if comp == 0; continue; end + for b = unique(sc.baseId(model.S(:,r)~=0))' + if sc.baseTransp(b) || isKey(relaxed,b) + used{b}(end+1) = comp; + end + end +end +trBase=[]; trComp=[]; +for b=1:sc.nBase + cs = unique(used{b}); + if numel(cs) <= 1; continue; end + for ci = sort(cs) + if ci ~= defC; trBase(end+1)=b; trComp(end+1)=ci; end %#ok + end +end +trBase=trBase(:); trComp=trComp(:); +end + +% ============================================================ certification +function [ok, growths] = i_certify(outModel, biomassId, minGrowth, growthConditions) +tol = 1e-9; +growths = struct(); +growths.primary = i_growOn(outModel, biomassId, []); +ok = growths.primary >= minGrowth - tol; +if ~isempty(growthConditions) && isstruct(growthConditions) + for i=1:numel(growthConditions) + gc = growthConditions(i); + g = i_growOn(outModel, biomassId, gc.medium); + growths.(matlab.lang.makeValidName(gc.name)) = g; + 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) +[~, 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 + +% ============================================================ transport pruning +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). +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(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(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(nTr,1); + fl(valid) = haveFlux(m, 'rxns', idx(valid)); + keep = keep | fl; +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 -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; +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 -outModel = i_applyAssignment(model, placement, addedTransports, comps, defaultCompartment); +% 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, +% 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 +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 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, trBase, trComp) 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> +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, find(strcmp(model.rxns,placement.rxns{k}),1)); + 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(trBase) + outModel = i_addTransport(outModel, model, sc.baseRep(trBase(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 @@ -355,8 +752,8 @@ 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 -newId = [model.mets{srcMetRow} '_' comp]; +if ~isempty(cand); idx = cand; return; end +newId = [model.mets{srcMetRow} '__' comp]; outModel.mets{end+1,1} = newId; outModel.metNames{end+1,1} = name; outModel.metComps(end+1,1) = ci; @@ -365,12 +762,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; @@ -381,3 +776,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/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 756352c2..d4ef7deb 100644 --- a/testing/function_tests/tAssignCompartments.m +++ b/testing/function_tests/tAssignCompartments.m @@ -45,6 +45,68 @@ 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 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'';''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 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'';''X'';''P''}, ''verbose'', false);']); + testCase.verifyNotEqual(sp.compartment{strcmp(sp.rxns,'r1')}, ... + 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) @@ -60,5 +122,48 @@ 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 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 + % 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.3 0.9]; + end end end