From ddfa2ab5c2ca56fe82f9e2a26efd7a0070526e16 Mon Sep 17 00:00:00 2001
From: Gabriel Studer <gabriel.studer@unibas.ch>
Date: Mon, 12 Feb 2024 13:23:15 +0100
Subject: [PATCH] chain mapping: replace GlobalAlign with SemiGlobalAlign in
 chain alignment

Terminal residues can be weird in GlobalAlign. E.g.

AAABAAB
AAA---B

instead of

AAABAAB
AAAB---

SemiGlobalAlign produces the latter
---
 modules/mol/alg/pymod/chain_mapping.py | 12 ++++++------
 modules/mol/alg/tests/test_qsscore.py  | 16 ++++++++--------
 2 files changed, 14 insertions(+), 14 deletions(-)

diff --git a/modules/mol/alg/pymod/chain_mapping.py b/modules/mol/alg/pymod/chain_mapping.py
index 5e90c49dc..466679ec2 100644
--- a/modules/mol/alg/pymod/chain_mapping.py
+++ b/modules/mol/alg/pymod/chain_mapping.py
@@ -1687,13 +1687,13 @@ class _Aligner:
         :returns: Alignment with s1 as first and s2 as second sequence 
         """
         if chem_type == mol.ChemType.AMINOACIDS:
-            return seq.alg.GlobalAlign(s1, s2, self.pep_subst_mat,
-                                       gap_open=self.pep_gap_open,
-                                       gap_ext=self.pep_gap_ext)[0]
+            return seq.alg.SemiGlobalAlign(s1, s2, self.pep_subst_mat,
+                                           gap_open=self.pep_gap_open,
+                                           gap_ext=self.pep_gap_ext)[0]
         elif chem_type == mol.ChemType.NUCLEOTIDES:
-            return seq.alg.GlobalAlign(s1, s2, self.nuc_subst_mat,
-                                       gap_open=self.nuc_gap_open,
-                                       gap_ext=self.nuc_gap_ext)[0]
+            return seq.alg.SemiGlobalAlign(s1, s2, self.nuc_subst_mat,
+                                           gap_open=self.nuc_gap_open,
+                                           gap_ext=self.nuc_gap_ext)[0]
         else:
             raise RuntimeError("Invalid ChemType")
         return aln
diff --git a/modules/mol/alg/tests/test_qsscore.py b/modules/mol/alg/tests/test_qsscore.py
index 5aed15ee8..a9b3c6131 100644
--- a/modules/mol/alg/tests/test_qsscore.py
+++ b/modules/mol/alg/tests/test_qsscore.py
@@ -130,8 +130,8 @@ class TestQSScore(unittest.TestCase):
         res = mapper.GetRigidMapping(ent_2, strategy="greedy_iterative_rmsd")
         qs_scorer = QSScorer.FromMappingResult(res)
         score_result = qs_scorer.Score(res.mapping)
-        self.assertAlmostEqual(score_result.QS_global, 0.3131, 2)
-        self.assertAlmostEqual(score_result.QS_best, 0.941, 2)
+        self.assertAlmostEqual(score_result.QS_global, 0.3191, 2)
+        self.assertAlmostEqual(score_result.QS_best, 0.9781, 2)
 
     def test_hetero_case_2_switched_order(self):
         # different stoichiometry
@@ -141,8 +141,8 @@ class TestQSScore(unittest.TestCase):
         res = mapper.GetRigidMapping(ent_2, strategy="greedy_iterative_rmsd")
         qs_scorer = QSScorer.FromMappingResult(res)
         score_result = qs_scorer.Score(res.mapping)
-        self.assertAlmostEqual(score_result.QS_global, 0.3131, 2)
-        self.assertAlmostEqual(score_result.QS_best, 0.941, 2)
+        self.assertAlmostEqual(score_result.QS_global, 0.3191, 2)
+        self.assertAlmostEqual(score_result.QS_best, 0.9781, 2)
 
     def test_hetero_case_3(self):
         # more chains
@@ -196,8 +196,8 @@ class TestQSScore(unittest.TestCase):
         res = mapper.GetRigidMapping(model, strategy="greedy_iterative_rmsd")
         qs_scorer = QSScorer.FromMappingResult(res)
         score_result = qs_scorer.Score(res.mapping)
-        self.assertAlmostEqual(score_result.QS_global, 0.323, 2)
-        self.assertAlmostEqual(score_result.QS_best, 0.921, 2)
+        self.assertAlmostEqual(score_result.QS_global, 0.321, 2)
+        self.assertAlmostEqual(score_result.QS_best, 0.932, 2)
 
     def test_hetero_model_switched_order(self):
         # same as above but with switched order to test for symmetric behaviour
@@ -208,8 +208,8 @@ class TestQSScore(unittest.TestCase):
         res = mapper.GetRigidMapping(model, strategy="greedy_iterative_rmsd")
         qs_scorer = QSScorer.FromMappingResult(res)
         score_result = qs_scorer.Score(res.mapping)
-        self.assertAlmostEqual(score_result.QS_global, 0.323, 2)
-        self.assertAlmostEqual(score_result.QS_best, 0.921, 2)
+        self.assertAlmostEqual(score_result.QS_global, 0.321, 2)
+        self.assertAlmostEqual(score_result.QS_best, 0.932, 2)
 
     def test_homo_1(self):
         # different stoichiometry SOD
-- 
GitLab