From 613428d9366798973c16a9dcdf2d878a1b62ef36 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Tue, 6 Oct 2026 17:19:18 +0200 Subject: [PATCH] feat: complexPolicy, keepGeneFree and preserveNotes in getModelFromHomology, and no repeated genes complexPolicy decides what happens to an and-linked subunit without an ortholog (flag, keep, drop), keepGeneFree also transfers template reactions without a grRule, and preserveNotes keeps the template's rxnNotes, rxnReferences and rxnConfidenceScores. A gene that several template genes map to appears once in a grRule, and the draft no longer carries the templates' version, date and description. Matches get_model_from_homology in raven-toolbox. --- .../homology/getModelFromHomology.m | 201 +++++++++++++++++- testing/function_tests/tReconstruction.m | 51 +++++ 2 files changed, 241 insertions(+), 11 deletions(-) diff --git a/reconstruction/homology/getModelFromHomology.m b/reconstruction/homology/getModelFromHomology.m index e6cbca09..ca26e70a 100755 --- a/reconstruction/homology/getModelFromHomology.m +++ b/reconstruction/homology/getModelFromHomology.m @@ -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 % ------- @@ -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. @@ -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=[]; @@ -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 @@ -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 @@ -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 '; @@ -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); @@ -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 +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 +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 + kids{end+1}=parts{m}; %#ok + 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 diff --git a/testing/function_tests/tReconstruction.m b/testing/function_tests/tReconstruction.m index dfe9624e..96e3b277 100644 --- a/testing/function_tests/tReconstruction.m +++ b/testing/function_tests/tReconstruction.m @@ -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');