Skip to content

Treat a gene taken by an earlier model as a subunit without an ortholog #26

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

Open
edkerk wants to merge 3 commits into
base: develop
Choose a base branch
from
Open

Treat a gene taken by an earlier model as a subunit without an ortholog #26

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
a037baa
Drop the complexes without an ortholog under the drop policy
edkerk Nov 5, 2026
ded3620
Treat a gene taken by an earlier model as a subunit without an ortholog
edkerk Nov 7, 2026
d21e997
Merge develop into the removed-genes branch
edkerk Nov 7, 2026
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
1 change: 1 addition & 0 deletions CHANGELOG.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@

### Breaking changes

* `getModelFromHomology` with `preferredOrder` treats a gene of a later model that maps to a new gene an earlier model already maps as a subunit without an ortholog, and `complexPolicy` decides what happens to it (was: the gene and every part of the rule that needs it were removed). For the rule `(a and b) or x` with `a` taken by an earlier model, `flag` gives `(OLD_MODELID_a and nb) or nx`, `keep` gives `nb or nx` and `drop` gives `nx`. A reaction that has genes, none of which keeps a mapping, is not included.
* `getModelFromHomology` with `complexPolicy` `drop` removes each complex that has a subunit without an ortholog, and removes the reaction only when no complex remains (was: the reaction was removed as soon as one subunit had no ortholog). For the rule `(A and B) or (C and D)` without an ortholog of `B`, the draft has the orthologs of `C` and `D`.

### Fixes
Expand Down
42 changes: 30 additions & 12 deletions src/reconstruction/homology/getModelFromHomology.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,12 @@
% the order in which reactions are added from the models, as a cell array
% of model IDs. If not supplied, reactions are included from all models,
% otherwise one gene results in reactions from one model only
% (default {}).
% (default {}). A gene of a later model that maps to a new gene which an
% earlier model already maps loses its mappings and counts as a subunit
% without an ortholog, so complexPolicy decides what happens to it. For
% the rule (a and b) or x with a taken by an earlier model, "flag" gives
% (OLD_MODELID_a and nb) or nx, "keep" gives nb or nx and "drop" gives nx.
% A reaction that has genes, none of which keeps a mapping, is not included.
% strictness : double
% integer that specifies which reactions are included (default 1):
%
Expand Down Expand Up @@ -503,16 +508,24 @@
[~, genesToDelete]=find(finalMappings{i}(usedGenes(a),:)); %IMPORTANT! IS it really correct to remove all genes that map?
genesToDelete=unique(genesToDelete); %Maybe not needed, but for clarity if nothing else

%Remove all the genes that were already found and add the other
%ones to allUsedGenes
if ~isempty(genesToDelete)
models{useOrderIndexes(i)}=removeGenes(models{useOrderIndexes(i)},allGenes{i+1}(genesToDelete),true,true,false);
end
%The genes that were already found lose their mappings and are
%treated as genes without an ortholog. Add the other ones to
%allUsedGenes
allUsedGenes(usedGenes)=true;

%Remove the deleted genes from finalMappings and allGenes.
finalMappings{i}(:,genesToDelete)=[];
allGenes{i+1}(genesToDelete)=[];

%Remove the reactions that have no gene with a mapping left
if ~isempty(genesToDelete)
a=ismember(models{useOrderIndexes(i)}.genes,allGenes{i+1});
[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;
models{useOrderIndexes(i)}=removeReactions(models{useOrderIndexes(i)},find(rxnsToRemove),true,true,true);
end
end
end

Expand Down Expand Up @@ -606,7 +619,7 @@
else
models{useOrderIndexes(i)}.genes=allGenes{1}(unique(newGenes));
end
if isfield(models{useOrderIndexes(i)},'geneComps')
if isfield(models{useOrderIndexes(i)},'geneComps') && ~isempty(models{useOrderIndexes(i)}.geneComps)
geneComps=models{useOrderIndexes(i)}.geneComps(1);
models{useOrderIndexes(i)}.geneComps=zeros(numel(models{useOrderIndexes(i)}.genes),1);
%Assume that all genes are in the same compartment, and this
Expand All @@ -616,7 +629,12 @@
end

%Now merge the models. All information should be correct except for "or"
%complexes
%complexes. A model that has lost all its reactions adds nothing to the draft
modelIds=cellfun(@(m) m.id,models,'UniformOutput',false);
hasRxns=cellfun(@(m) ~isempty(m.rxns),models);
if any(hasRxns)
models=models(hasRxns);
end
draftModel=mergeModels(models,'metNames');

%Remove OLD_ genes that are a member of an isozyme group (an "or" group) with
Expand Down Expand Up @@ -659,11 +677,11 @@
%Change name of the resulting model
draftModel.id=getModelFor;
name='Generated by getModelFromHomology using ';
for i=1:numel(models)
if i<numel(models)
name=[name models{i}.id ', '];
for i=1:numel(modelIds)
if i<numel(modelIds)
name=[name modelIds{i} ', '];
else
name=[name models{i}.id];
name=[name modelIds{i}];
end
end
draftModel.name=name;
Expand Down
53 changes: 53 additions & 0 deletions tests/function_tests/tReconstruction.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -152,6 +152,59 @@ function getModelFromHomologyPreferredOrderDropsSharedGenes(testCase)
testCase.verifyEqual(sort(draft.genes), {'N1'; 'N2'});
end

function getModelFromHomologyPreferredOrderRetainsTheOrphanSubunit(testCase)
% A gene taken by an earlier template is a subunit without an ortholog:
% complexPolicy decides what happens to the rest of its complex.
t1 = tReconstruction.homologyTemplate('T1', {'R_T1'}, {'Y'});
t2 = tReconstruction.homologyTemplate('T2', {'R_T2', 'R_T3'}, {'X or Z', 'X and Z'});
pairs = {{'Y', 'N1'}, {'X', 'N1'; 'Z', 'N3'}};
draft = tReconstruction.homologyDraft({t1, t2}, pairs, 'flag', {'T1', 'T2'});
testCase.verifyEqual(draft.grRules(strcmp(draft.rxns, 'R_T3')), {'OLD_T2_X and N3'});
testCase.verifyEqual(draft.grRules(strcmp(draft.rxns, 'R_T2')), {'N3'});
draft = tReconstruction.homologyDraft({t1, t2}, pairs, 'keep', {'T1', 'T2'});
testCase.verifyEqual(draft.grRules(strcmp(draft.rxns, 'R_T3')), {'N3'});
draft = tReconstruction.homologyDraft({t1, t2}, pairs, 'drop', {'T1', 'T2'});
testCase.verifyFalse(ismember('R_T3', draft.rxns));
testCase.verifyTrue(ismember('R_T2', draft.rxns));
% a taken gene in a complex beside an isozyme
t3 = tReconstruction.homologyTemplate('T2', {'R_T2'}, {'(a and b) or x'});
pairs = {{'Y', 'N1'}, {'a', 'N1'; 'b', 'Nb'; 'x', 'Nx'}};
rules = cell(1, 3);
policies = {'flag', 'keep', 'drop'};
for k = 1:3
draft = tReconstruction.homologyDraft({t1, t3}, pairs, policies{k}, {'T1', 'T2'});
rules{k} = draft.grRules{strcmp(draft.rxns, 'R_T2')};
end
testCase.verifyEqual(rules, {'(OLD_T2_a and Nb) or Nx', 'Nb or Nx', 'Nx'});
% a reaction whose genes were all taken is not included
t4 = tReconstruction.homologyTemplate('T2', {'R_T2'}, {'a and b'});
draft = tReconstruction.homologyDraft({t1, t4}, {{'Y', 'N1'}, {'a', 'N1'}}, 'flag', {'T1', 'T2'});
testCase.verifyEqual(draft.rxns, {'R_T1'});
end

function getModelFromHomologyPreferredOrderEmptiesATemplateWithGeneComps(testCase)
% A template with a geneComps field that loses all its genes in the preferred-order step
t1 = tReconstruction.homologyTemplate('T1', {'R_T1'}, {'Y'});
t2 = tReconstruction.homologyTemplate('T2', {'R_T2'}, {'X'});
t2.geneComps = ones(numel(t2.genes), 1);
draft = tReconstruction.homologyDraft({t1, t2}, {{'Y', 'N1'}, {'X', 'N1'}}, 'flag', {'T1', 'T2'});
testCase.verifyEqual(draft.rxns, {'R_T1'});
end

function getModelFromHomologyPreferredOrderChainsTheConflicts(testCase)
% Three templates: a gene of the third that maps to a gene taken by the first, and a
% gene that maps to a gene taken by the second, follow complexPolicy.
t1 = tReconstruction.homologyTemplate('T1', {'R_T1'}, {'Y'});
t2 = tReconstruction.homologyTemplate('T2', {'R_T2'}, {'X'});
t3 = tReconstruction.homologyTemplate('T3', {'R_T3'}, {'(p and q) or r'});
pairs = {{'Y', 'N1'}, {'X', 'N2'}, {'p', 'N1'; 'q', 'N2'; 'r', 'N3'}};
draft = tReconstruction.homologyDraft({t1, t2, t3}, pairs, 'flag', {'T1', 'T2', 'T3'});
testCase.verifyEqual(draft.grRules(strcmp(draft.rxns, 'R_T3')), {'(OLD_T3_p and OLD_T3_q) or N3'});
draft = tReconstruction.homologyDraft({t1, t2, t3}, pairs, 'keep', {'T1', 'T2', 'T3'});
testCase.verifyEqual(draft.grRules(strcmp(draft.rxns, 'R_T3')), {'N3'});
testCase.verifyEqual(draft.name, 'Generated by getModelFromHomology using T1, T2, T3');
end

function getModelFromHomologyPreferredOrderWithoutSharedGenes(testCase)
% Templates with the same reaction id and no shared gene: the order does not
% change the result, and no empty gene list reaches removeGenes.
Expand Down
Loading