Skip to content

feat: complexPolicy, keepGeneFree and preserveNotes in getModelFromHomology, and no repeated genes #757

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
edkerk merged 1 commit into from
Oct 7, 2026
Merged

feat: complexPolicy, keepGeneFree and preserveNotes in getModelFromHomology, and no repeated genes #757

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
201 changes: 190 additions & 11 deletions reconstruction/homology/getModelFromHomology.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -56,6 +56,23 @@
% determines how to match genes if not looking at only 1-1 orthologs.
% Either map the new genes to the old or old genes to new. The default
% is to map the new genes (default true).
% complexPolicy : char
% what happens to an "and"-linked subunit without an ortholog:
% "flag" keeps the reaction with the subunit renamed
% OLD_MODELID_geneName, for a curator to review; "keep" drops just that
% subunit; "drop" drops the whole reaction (default "flag").
% keepGeneFree : logical
% also transfer the template reactions that have no grRule (spontaneous,
% transport, exchange, artificial). By default they are dropped, because
% nothing supports them in the new organism. Set it when the template is
% a curated reference whose gene-free reactions belong in every derived
% model. With several templates, a gene-free reaction comes in once per
% template (default false).
% preserveNotes : logical
% keep each transferred reaction's rxnNotes, rxnReferences and
% rxnConfidenceScores from the template. By default the notes say the
% reaction was included by getModelFromHomology and the confidence score
% is 2 (default false).
%
% Returns
% -------
Expand All @@ -77,9 +94,12 @@
% grRules (use expandModel).
%
% The resulting draft model contains only reactions associated with
% orthologous genes. The old (original) genes involved in "and" relations
% orthologous genes, plus the gene-free ones with keepGeneFree. With
% complexPolicy "flag", the old (original) genes involved in "and" relations
% in grRules without any orthologs are still included in the draft model as
% OLD_MODELID_geneName.
% OLD_MODELID_geneName. A gene that several template genes map to appears
% once in a grRule. The draft does not keep the version, date or
% descriptive annotation of the templates.
%
% "to" and "from" means relative to the new organism.

Expand All @@ -90,7 +110,14 @@

p=parseRAVENargs(varargin, {'preferredOrder',[]; 'strictness',1; ...
'onlyGenesInModels',false; 'maxE',10^-30; 'minLen',100; 'minIde',40; ...
'mapNewGenesToOld',true});
'mapNewGenesToOld',true; 'complexPolicy','flag'; 'keepGeneFree',false; ...
'preserveNotes',false});
complexPolicy=lower(char(p.complexPolicy));
if ~ismember(complexPolicy,{'flag','keep','drop'})
error('complexPolicy must be "flag", "keep" or "drop", not "%s"',complexPolicy);
end
keepGeneFree=p.keepGeneFree;
preserveNotes=p.preserveNotes;
preferredOrder=p.preferredOrder;
if isempty(preferredOrder)
preferredOrder=[];
Expand Down Expand Up @@ -194,14 +221,20 @@
blastStructure(i).ppos(~indexes)=[];
end

%Remove all reactions from the models that have no genes encoding for them.
%Also remove all genes that encode for no reactions. There should not be any
%but there might be mistakes
%Remove all reactions from the models that have no genes encoding for them,
%unless keepGeneFree. Also remove all genes that encode for no reactions.
%There should not be any but there might be mistakes
geneFreeRxns=cell(numel(models),1);
geneFreeRxns(:)={{}};
for i=1:numel(models)
[hasGenes, ~]=find(models{i}.rxnGeneMat);
hasNoGenes=1:numel(models{i}.rxns);
hasNoGenes(hasGenes)=[];
models{i}=removeReactions(models{i},hasNoGenes,true,true);
if keepGeneFree
geneFreeRxns{i}=models{i}.rxns(hasNoGenes);
else
models{i}=removeReactions(models{i},hasNoGenes,true,true);
end
end

%Create a structure that contains all genes used in the blasts in any
Expand Down Expand Up @@ -409,6 +442,7 @@
[rxnsToKeep,~] = find(models{useOrderIndexes(i)}.rxnGeneMat(:,a));
rxnsToRemove = repmat(1,numel(models{useOrderIndexes(i)}.rxns),1);
rxnsToRemove(rxnsToKeep) = 0;
rxnsToRemove(ismember(models{useOrderIndexes(i)}.rxns,geneFreeRxns{useOrderIndexes(i)})) = 0;
rxnsToRemove = find(rxnsToRemove);
models{useOrderIndexes(i)}=removeReactions(models{useOrderIndexes(i)},rxnsToRemove,true,true,true);
end
Expand Down Expand Up @@ -571,6 +605,27 @@
draftModel.grRules{oldI}=gprBuildRule(tkns);
end

%Apply the complex policy to the OLD_ genes left in "and" relations, and
%write each grRule without repeated genes or redundant brackets
toRemove=false(numel(draftModel.rxns),1);
for j=1:numel(draftModel.grRules)
rule=draftModel.grRules{j};
if isempty(rule)
continue
end
if contains(rule,'OLD_') && strcmp(complexPolicy,'drop')
toRemove(j)=true;
elseif contains(rule,'OLD_') && strcmp(complexPolicy,'keep')
draftModel.grRules{j}=gprSimplify(rule,'OLD_');
toRemove(j)=isempty(draftModel.grRules{j});
else
draftModel.grRules{j}=gprSimplify(rule,'');
end
end
if any(toRemove)
draftModel=removeReactions(draftModel,draftModel.rxns(toRemove),true,true,true);
end

%Change name of the resulting model
draftModel.id=getModelFor;
name='Generated by getModelFromHomology using ';
Expand All @@ -582,10 +637,21 @@
end
end
draftModel.name=name;
draftModel.rxnNotes=cell(length(draftModel.rxns),1);
draftModel.rxnNotes(:)={'Included by getModelFromHomology'};
draftModel.rxnConfidenceScores=NaN(length(draftModel.rxns),1);
draftModel.rxnConfidenceScores(:)=2;
%The draft is a new model: the templates' version, date and description are
%not its own
draftModel=rmfield(draftModel,intersect({'version','date','description'},fieldnames(draftModel)));
if isfield(draftModel,'annotation')
descriptive={'taxonomy','sourceUrl','note','givenName','familyName','authors','email','organization'};
draftModel.annotation=rmfield(draftModel.annotation,intersect(descriptive,fieldnames(draftModel.annotation)));
end
if ~preserveNotes || ~isfield(draftModel,'rxnNotes')
draftModel.rxnNotes=cell(length(draftModel.rxns),1);
draftModel.rxnNotes(:)={'Included by getModelFromHomology'};
end
if ~preserveNotes || ~isfield(draftModel,'rxnConfidenceScores')
draftModel.rxnConfidenceScores=NaN(length(draftModel.rxns),1);
draftModel.rxnConfidenceScores(:)=2;
end
draftModel=deleteUnusedGenes(draftModel,0);
%Standardize grRules and notify if problematic grRules are found
[draftModel.grRules,draftModel.rxnGeneMat]=standardizeGrRules(draftModel,false);
Expand Down Expand Up @@ -646,3 +712,116 @@
end
end
end


function rule=gprSimplify(rule,removePrefix)
% Rewrite a grRule as a tree: genes starting with removePrefix (if any) are
% dropped from their group, groups of the same operator are flattened,
% repeated genes or groups are removed, and brackets appear only around a
% group of the other operator. A rule with no gene left becomes ''.
tokens=gprTokenize(rule);
if isempty(tokens)
rule='';
return
end
[node,~]=gprParseOr(tokens,1);
node=gprPrune(node,removePrefix);
if isempty(node)
rule='';
else
rule=gprRender(node,'');
end
end


function [node,pos]=gprParseOr(tokens,pos)
[node,pos]=gprParseAnd(tokens,pos);
kids={node};
while pos<=numel(tokens) && strcmpi(tokens{pos},'or')
[node,pos]=gprParseAnd(tokens,pos+1);
kids{end+1}=node; %#ok<AGROW>
end
if numel(kids)>1
node=struct('op','or','gene','','kids',{kids});
end
end


function [node,pos]=gprParseAnd(tokens,pos)
[node,pos]=gprParseAtom(tokens,pos);
kids={node};
while pos<=numel(tokens) && strcmpi(tokens{pos},'and')
[node,pos]=gprParseAtom(tokens,pos+1);
kids{end+1}=node; %#ok<AGROW>
end
if numel(kids)>1
node=struct('op','and','gene','','kids',{kids});
end
end


function [node,pos]=gprParseAtom(tokens,pos)
if pos>numel(tokens)
error('Malformed grRule: it ends with an operator');
end
if strcmp(tokens{pos},'(')
[node,pos]=gprParseOr(tokens,pos+1);
if pos>numel(tokens) || ~strcmp(tokens{pos},')')
error('Malformed grRule: unbalanced brackets');
end
pos=pos+1;
else
node=struct('op','gene','gene',tokens{pos},'kids',{{}});
pos=pos+1;
end
end


function node=gprPrune(node,removePrefix)
if strcmp(node.op,'gene')
if ~isempty(removePrefix) && startsWith(node.gene,removePrefix)
node=[];
end
return
end
kids={};
seen={};
for k=1:numel(node.kids)
kid=gprPrune(node.kids{k},removePrefix);
if isempty(kid)
continue
end
if strcmp(kid.op,node.op)
parts=kid.kids;
else
parts={kid};
end
for m=1:numel(parts)
key=gprRender(parts{m},node.op);
if ~ismember(key,seen)
seen{end+1}=key; %#ok<AGROW>
kids{end+1}=parts{m}; %#ok<AGROW>
end
end
end
if isempty(kids)
node=[];
elseif numel(kids)==1
node=kids{1};
else
node.kids=kids;
end
end


function rule=gprRender(node,parentOp)
if strcmp(node.op,'gene')
rule=node.gene;
return
end
parts=cellfun(@(k) gprRender(k,node.op),node.kids,'UniformOutput',false);
rule=strjoin(parts,[' ' node.op ' ']);
if ~isempty(parentOp) && ~strcmp(parentOp,node.op)
rule=['(' rule ')'];
end
end
51 changes: 51 additions & 0 deletions testing/function_tests/tReconstruction.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,57 @@ function getModelFromHomologyBuildsDraft(testCase)
testCase.verifyClass(draft, 'struct');
end

function getModelFromHomologyComplexPolicyAndOptions(testCase)
% A small template: a complex, an isozyme pair, a gene-free
% reaction, and two template genes that map to one new gene.
t = struct();
t.id = 'tmpl'; t.name = 'tmpl';
t.mets = {'a_c';'b_c';'d_c'}; t.metNames = {'A';'B';'D'};
t.comps = {'c'}; t.compNames = {'cytosol'}; t.metComps = [1;1;1];
t.rxns = {'Rcplx';'Riso';'Rfree';'Rdup'};
t.rxnNames = t.rxns;
t.S = sparse([-1 0 0 -1; 0 -1 0 1; 1 1 -1 0]);
t.lb = [0;0;0;0]; t.ub = [1000;1000;1000;1000]; t.c = [0;0;0;0]; t.b = [0;0;0];
t.rev = [0;0;0;0];
t.grRules = {'tg1 and tg2'; 'tg3 or tg4'; ''; 'tg3 or tg4'};
t.genes = {'tg1';'tg2';'tg3';'tg4'};
t.rxnGeneMat = sparse([1 1 0 0; 0 0 1 1; 0 0 0 0; 0 0 1 1]);
t.rxnNotes = {'n1';'n2';'n3';'n4'};
t.rxnConfidenceScores = [4;4;4;4];
t.version = '9.9.9';
% makeFakeBlastStructure wants at least 10 pairs; the extra ones
% are genes the template does not have
pairs = [{'tg1','ng1'; 'tg3','ng3'; 'tg4','ng3'}; ...
[strcat('x', string(1:7))', strcat('nx', string(1:7))']];
pairs = cellstr(pairs);
bs = makeFakeBlastStructure(pairs, 'tmpl', 'new');
evalc('flag = getModelFromHomology({t}, bs, ''new'', ''minLen'', 0, ''minIde'', 0);');
evalc('keep = getModelFromHomology({t}, bs, ''new'', ''minLen'', 0, ''minIde'', 0, ''complexPolicy'', ''keep'', ''keepGeneFree'', true, ''preserveNotes'', true);');
evalc('drop = getModelFromHomology({t}, bs, ''new'', ''minLen'', 0, ''minIde'', 0, ''complexPolicy'', ''drop'');');
% flag: unmapped subunit kept as OLD_; gene-free reaction dropped
testCase.verifyTrue(contains(flag.grRules{strcmp(flag.rxns,'Rcplx')}, 'OLD_tmpl_tg2'));
testCase.verifyFalse(ismember('Rfree', flag.rxns));
% two template genes mapping to one new gene appear once
testCase.verifyEqual(flag.grRules{strcmp(flag.rxns,'Riso')}, 'ng3');
% keep: subunit dropped, gene-free kept, template notes kept
testCase.verifyEqual(keep.grRules{strcmp(keep.rxns,'Rcplx')}, 'ng1');
testCase.verifyTrue(ismember('Rfree', keep.rxns));
testCase.verifyEqual(keep.rxnNotes{strcmp(keep.rxns,'Rcplx')}, 'n1');
testCase.verifyEqual(keep.rxnConfidenceScores(strcmp(keep.rxns,'Rcplx')), 4);
% drop: the incomplete complex is removed
testCase.verifyFalse(ismember('Rcplx', drop.rxns));
testCase.verifyTrue(ismember('Riso', drop.rxns));
% the template version is not inherited
testCase.verifyFalse(isfield(flag, 'version') && strcmp(flag.version, '9.9.9'));
testCase.verifyEqual(flag.rxnNotes{1}, 'Included by getModelFromHomology');
end

function getModelFromHomologyRejectsUnknownComplexPolicy(testCase)
src = testCase.model; src.id = 'srcModel';
bs = makeFakeBlastStructure([src.genes(1:10), strcat('t_', src.genes(1:10))], 'srcModel', 'tgtOrg');
testCase.verifyError(@() getModelFromHomology({src}, bs, 'tgtOrg', 'complexPolicy', 'maybe'), ?MException);
end

function getBlastRunsWhenAvailable(testCase)
fa1 = fullfile(testCase.ravenRoot,'testing','function_tests','test_data','human_galactosidases.fa');
fa2 = fullfile(testCase.ravenRoot,'testing','function_tests','test_data','yeast_galactosidases.fa');
Expand Down
Loading