diff --git a/CHANGELOG.md b/CHANGELOG.md
index d5200eb8..6892c2d4 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -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
diff --git a/src/reconstruction/homology/getModelFromHomology.m b/src/reconstruction/homology/getModelFromHomology.m
index 4dba5084..2768b5fb 100755
--- a/src/reconstruction/homology/getModelFromHomology.m
+++ b/src/reconstruction/homology/getModelFromHomology.m
@@ -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):
 %
@@ -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
 
@@ -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
@@ -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
@@ -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;
diff --git a/tests/function_tests/tReconstruction.m b/tests/function_tests/tReconstruction.m
index 494b5aaf..607a3007 100644
--- a/tests/function_tests/tReconstruction.m
+++ b/tests/function_tests/tReconstruction.m
@@ -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.