From 8f6b23368ca66dd29436b87cbf8a0844773e894a Mon Sep 17 00:00:00 2001
From: Eduard Kerkhoven <eduardk@chalmers.se>
Date: Fri, 28 Aug 2026 09:43:00 +0200
Subject: [PATCH] fix: gene field alignment in removeLowScoreGenes

Ports the fix from #670 (merged to develop) to develop3, which carried the
identical bug in its reworked removeLowScoreGenes.

removeLowScoreGenes regenerated model.genes from getGenesFromGrRules, which
returns a sorted list, but trimmed the gene-associated fields
(geneShortNames, proteins, geneMiriams, geneFrom, geneComps) with a mask in
the original gene order. When model.genes was not already sorted, this left
every annotation field shifted relative to model.genes, corrupting the gene
ID to gene symbol mapping in ftINIT-reconstructed models (removeGenes=true).

Retain the remaining genes in their original order and reorder the
rxnGeneMat columns to match, so the trimmed annotation fields stay aligned.

Reproduced on develop3 with an unsorted-genes model: removing one isozyme
gene paired G2->short3 and G3->short2; after the fix G2->short2, G3->short3.
Adds removeLowScoreGenesKeepsFieldsAligned to tINIT covering exactly this.
---
 INIT/removeLowScoreGenes.m     | 15 +++++++++++----
 testing/function_tests/tINIT.m | 32 ++++++++++++++++++++++++++++++++
 2 files changed, 43 insertions(+), 4 deletions(-)

diff --git a/INIT/removeLowScoreGenes.m b/INIT/removeLowScoreGenes.m
index e945d2a3..df2378d8 100644
--- a/INIT/removeLowScoreGenes.m
+++ b/INIT/removeLowScoreGenes.m
@@ -101,13 +101,20 @@
 
 % regenerate "genes" and "rxnGeneMat" model fields
 [genes,rxnGeneMat] = getGenesFromGrRules(newModel.grRules);
-newModel.genes = genes;
-newModel.rxnGeneMat = rxnGeneMat;
 
-% update other gene-related fields
-remInd = ~ismember(model.genes,newModel.genes);
+% determine which of the original genes were removed
+remInd = ~ismember(model.genes,genes);
 remGenes = model.genes(remInd);
 
+% Keep the retained genes in their original order rather than the sorted
+% order returned by getGenesFromGrRules. Gene removal never introduces new
+% genes, so the remaining genes are model.genes(~remInd). Preserving this
+% order ensures the gene-associated fields trimmed below (geneShortNames,
+% proteins, etc.), which are indexed by remInd, stay aligned with genes.
+newModel.genes = model.genes(~remInd);
+[~,reorderInd] = ismember(newModel.genes,genes);
+newModel.rxnGeneMat = rxnGeneMat(:,reorderInd);
+
 if isfield(newModel,'geneShortNames')
     newModel.geneShortNames(remInd) = [];
 end
diff --git a/testing/function_tests/tINIT.m b/testing/function_tests/tINIT.m
index a0eec8e4..476d9761 100644
--- a/testing/function_tests/tINIT.m
+++ b/testing/function_tests/tINIT.m
@@ -29,6 +29,38 @@ function removeLowScoreGenesRuns(testCase)
             testCase.verifyClass(m2, 'struct');
         end
 
+        function removeLowScoreGenesKeepsFieldsAligned(testCase)
+            % getGenesFromGrRules returns a sorted gene list, but the
+            % annotation fields are trimmed with a mask in the original gene
+            % order. If newModel.genes took the sorted order, every annotation
+            % field would shift relative to it. Use deliberately unsorted genes
+            % so sorted != original, and drop one isozyme gene.
+            m = struct();
+            m.id = 'test';
+            m.rxns = {'R1'}; m.rxnNames = {'R1'};
+            m.mets = {'A';'B'}; m.metNames = {'A';'B'}; m.metComps = [1;1];
+            m.comps = {'c'}; m.compNames = {'c'};
+            m.S = sparse([-1;1]); m.lb = 0; m.ub = 1000; m.rev = 0; m.c = 0;
+            m.b = [0;0];
+            m.genes = {'G3';'G1';'G2'};                      % unsorted
+            m.geneShortNames = {'short3';'short1';'short2'}; % aligned to genes
+            m.grRules = {'G3 or G1 or G2'};
+            m.rxnGeneMat = sparse([1 1 1]);
+            scores = [1; -1; 1];   % G1 (negative) is dropped from the isozyme
+            evalc('[mm, removed] = removeLowScoreGenes(m, scores);');
+
+            testCase.verifyEqual(removed, {'G1'});
+            % Each surviving gene keeps its own short name.
+            for k = 1:numel(mm.genes)
+                want = strrep(mm.genes{k}, 'G', 'short');
+                testCase.verifyEqual(mm.geneShortNames{k}, want);
+            end
+            % rxnGeneMat columns must correspond to mm.genes, and R1 still uses
+            % both surviving genes.
+            testCase.verifyEqual(numel(mm.genes), size(mm.rxnGeneMat, 2));
+            testCase.verifyEqual(full(mm.rxnGeneMat), ones(1, numel(mm.genes)));
+        end
+
         function reverseRxnsRuns(testCase)
             % R1 = '=> a[s]'; R3 = 'a[c] <=> b[c] + c[c]'
             testModel = getTstModel();