Repository: haddocking/haddock3
Issue: correct lrmsd and ilrmsd calculations in caprieval

There is a bug in how the caprieval module calculates the `lrmsd` and `ilrmsd` quantities. `irmsd` values should be correct.

For lrmsd:

- (minor) the fitting on the receptor is including two atoms of the ligand (indexing error)
- (minor) the rmsd calculation on the ligand is missing those two atoms
- (major) translations are removed by using the centroid of the full complex instead of the centroid of the receptor alone. This gives a wrong rotation matrix from the kabsch alignment, thus resulting in a wrong value of lrmsd

For ilrmsd I just re-wrote the function after realising the output value did not agree with that of an external software 

Changed files:
src/haddock/modules/analysis/caprieval/capri.py, tests/test_module_caprieval.py

diff --git a/src/haddock/modules/analysis/caprieval/capri.py b/src/haddock/modules/analysis/caprieval/capri.py
index f41de3db..978be76d 100644
--- a/src/haddock/modules/analysis/caprieval/capri.py
+++ b/src/haddock/modules/analysis/caprieval/capri.py
@@ -146,9 +146,8 @@ class CAPRI:
             if chain not in chain_ranges:
                 chain_ranges[chain] = []
             chain_ranges[chain].append(i)
-
-        chain_ranges = make_range(chain_ranges)
         
+        chain_ranges = make_range(chain_ranges)
         obs_chains = list(chain_ranges.keys())  # observed chains
         if len(obs_chains) < 2:
             log.warning("Not enough chains for calculating lrmsd")
@@ -167,30 +166,35 @@ class CAPRI:
             Q = np.asarray(Q)
             P = np.asarray(P)
 
-            # write_coord_dic("ref.pdb", ref_coord_dic)
-            # write_coord_dic("model.pdb", mod_coord_dic)
+            # write_coords("ref_first.pdb", Q)
+            # write_coords("model_first.pdb", P)
 
-            # write_coords("ref.pdb", Q)
-            # write_coords("model.pdb", P)
+            # get receptor and ligand coordinates
+            Q_r_first = Q[r_start: r_end + 1]
+            P_r_first = P[r_start: r_end + 1]
+            # write_coords("ref_r_first.pdb", Q_r_first)
+            # write_coords("model_r_first.pdb", P_r_first)
+            # Q_l_first = Q[l_start: l_end + 1]
+            # P_l_first = P[l_start: l_end + 1]
+            # write_coords("ref_l_first.pdb", Q_l_first)
+            # write_coords("model_l_first.pdb", P_l_first)
 
-            # move to the origin
-            Q = Q - centroid(Q)
-            P = P - centroid(P)
+            # move to the origin of the receptor
+            
+            Q = Q - centroid(Q_r_first)
+            P = P - centroid(P_r_first)
 
             # get receptor coordinates
-            Q_r = Q[r_start: r_end - 1]
-            P_r = P[r_start: r_end - 1]
-
+            Q_r = Q[r_start: r_end + 1]
+            P_r = P[r_start: r_end + 1]
             # Center receptors and get rotation matrix
             # Q_r = Q_r - centroid(Q_r)
             # P_r = P_r - centroid(P_r)
+            # write_coords("ref_r_centr.pdb", Q_r)
+            # write_coords("model_r_centr.pdb", P_r)
 
             U_r = kabsch(P_r, Q_r)
 
-            # Center complexes at receptor centroids
-            Q = Q - centroid(Q_r)
-            P = P - centroid(P_r)
-
             # Apply rotation to complex
             #  - complex are now aligned by the receptor
             P = np.dot(P, U_r)
@@ -199,8 +203,8 @@ class CAPRI:
             # write_coords("model.pdb", P)
 
             # Identify the ligand coordinates
-            Q_l = Q[l_start: l_end - 1]
-            P_l = P[l_start: l_end - 1]
+            Q_l = Q[l_start: l_end + 1]
+            P_l = P[l_start: l_end + 1]
 
             # write_coords("ref_l.pdb", Q_l)
             # write_coords("model_l.pdb", P_l)
@@ -208,9 +212,6 @@ class CAPRI:
             # Calculate the RMSD of the ligands
             self.lrmsd = calc_rmsd(P_l, Q_l)
 
-            # write_coords("ref.pdb", Q)
-            # write_coords("model.pdb", P)
-
     def calc_ilrmsd(self, cutoff=10.0):
         """
         Calculate the Interface Ligand RMSD.
@@ -222,20 +223,12 @@ class CAPRI:
         """
         # Identify interface
         ref_interface_resdic = self.identify_interface(self.reference, cutoff)
-
         # Load interface coordinates
-        ref_coord_dic, _ = load_coords(self.reference, self.atoms)
 
         ref_int_coord_dic, _ = load_coords(
             self.reference, self.atoms, ref_interface_resdic
             )
 
-        mod_coord_dic, _ = load_coords(
-            self.model,
-            self.atoms,
-            numbering_dic=self.model2ref_numbering
-            )
-
         mod_int_coord_dic, _ = load_coords(
             self.model,
             self.atoms,
@@ -262,20 +255,15 @@ class CAPRI:
 
         # write_coords("ref.pdb", Q_int)
         # write_coords("model.pdb", P_int)
-
-        # find atoms present in both molecules
-        Q = []
-        P = []
-        intersection = sorted(ref_coord_dic.keys() & mod_coord_dic.keys())
+        
         chain_ranges = {}
-        for i, segment in enumerate(intersection):
+        for i, segment in enumerate(sorted(common_keys)):
             chain, _, _ = segment
             if chain not in chain_ranges:
                 chain_ranges[chain] = []
             chain_ranges[chain].append(i)
 
         chain_ranges = make_range(chain_ranges)
-
         obs_chains = list(chain_ranges.keys())  # observed chains
         if len(obs_chains) < 2:
             log.warning("Not enough chains for calculating ilrmsd")
@@ -283,64 +271,39 @@ class CAPRI:
             r_chain, l_chain = self.check_chains(obs_chains)
 
             l_start, l_end = chain_ranges[l_chain]
-
-            for k in sorted(ref_coord_dic.keys() & mod_coord_dic.keys()):
-                ref_xyz = ref_coord_dic[k]
-                mod_xyz = mod_coord_dic[k]
-
-                Q.append(ref_xyz)
-                P.append(mod_xyz)
-
-            Q = np.asarray(Q)
-            P = np.asarray(P)
+            r_start, r_end = chain_ranges[r_chain]
 
             # write_coords("ref.pdb", Q)
             # write_coords("model.pdb", P)
 
-            # put system at origin
-            Q_int = Q_int - centroid(Q_int)
-            P_int = P_int - centroid(P_int)
+            # put system at origin of the receptor interface
+            Q_r_int = Q_int[r_start: r_end + 1]
+            P_r_int = P_int[r_start: r_end + 1]
 
+            Q_int = Q_int - centroid(Q_r_int)
+            P_int = P_int - centroid(P_r_int)
             # put interfaces at the origin
-            # Q_int = Q_int - centroid(Q_int)
-            # P_int = P_int - centroid(P_int)
 
-            # find the rotation that minimizes the interface rmsd
-            U_int = kabsch(P_int, Q_int)
+            # find the rotation that minimizes the receptor interface rmsd
+            Q_r_int = Q_int[r_start: r_end + 1]
+            P_r_int = P_int[r_start: r_end + 1]
+            
+            U_int = kabsch(P_r_int, Q_r_int)
             P_int = np.dot(P_int, U_int)
-
-            # write_coords("ref.pdb", Q_int)
-            # write_coords("model.pdb", P_int)
-
-            # # move the system to the centroid of the interfaces
-            Q = Q - centroid(Q)
-            P = P - centroid(P)
-
-            # write_coords("ref_1.pdb", Q)
-            # write_coords("model_1.pdb", P)
-
-            Q = Q - centroid(Q_int)
-            P = P - centroid(P_int)
-
-            # write_coords("ref_2.pdb", Q)
-            # write_coords("model_2.pdb", P)
-
-            # apply this rotation to the model
-            #  - complexes are now aligned by the interfaces
-            P = np.dot(P, U_int)
-
-            # write_coords("ref_i.pdb", Q)
-            # write_coords("model_i.pdb", P)
-
-            # Calculate the rmsd of the ligand
-            Q_l = Q[l_start: l_end - 1]
-            P_l = P[l_start: l_end - 1]
-
-            # write_coords("ref_l.pdb", Q_l)
-            # write_coords("model_l.pdb", P_l)
-
-            # this will be the interface-ligand-rmsd
-            self.ilrmsd = calc_rmsd(P_l, Q_l)
+            # just for checks.
+            # the interface rmsd for the rec interface should be almost zero
+            # Q_r_int = Q_int[r_start: r_end + 1]
+            # P_r_int = P_int[r_start: r_end + 1]
+            # r_rmsd = calc_rmsd(Q_r_int, P_int[r_start: r_end + 1])
+            # print(r_rmsd)
+
+            Q_l_int = Q_int[l_start: l_end + 1]
+            P_l_int = P_int[l_start: l_end + 1]
+            # write_coords("ref_l_int_fin.pdb", Q_l_int)
+            # write_coords("mod_l_int_fin.pdb", P_l_int)
+
+            # # this will be the interface-ligand-rmsd
+            self.ilrmsd = calc_rmsd(P_l_int, Q_l_int)
             
     def calc_fnat(self, cutoff=5.0):
         """
@@ -447,7 +410,7 @@ class CAPRI:
         if self.params["irmsd"]:
             log.debug(f"id {self.identificator}, calculating I-RMSD")
             irmsd_cutoff = self.params["irmsd_cutoff"]
-            log.debug(f" cutoff: {irmsd_cutoff}A")
+            log.info(f" cutoff: {irmsd_cutoff}A")
             self.calc_irmsd(cutoff=irmsd_cutoff)
 
         if self.params["lrmsd"]:
diff --git a/tests/test_module_caprieval.py b/tests/test_module_caprieval.py
index cda7dc2b..a326ab96 100644
--- a/tests/test_module_caprieval.py
+++ b/tests/test_module_caprieval.py
@@ -158,6 +158,10 @@ def protprot_caprimodule_parallel(protprot_input_list):
 
 def test_protprot_irmsd(protprot_caprimodule):
     """Test protein-protein i-rmsd calculation."""
+    # using standard cutoff
+    protprot_caprimodule.calc_irmsd(cutoff=10.0)
+    assert round_two_dec(protprot_caprimodule.irmsd) == 8.33
+    # using default cutoff = 5.0
     protprot_caprimodule.calc_irmsd()
     assert round_two_dec(protprot_caprimodule.irmsd) == 7.38
 
@@ -165,13 +169,13 @@ def test_protprot_irmsd(protprot_caprimodule):
 def test_protprot_lrmsd(protprot_caprimodule):
     """Test protein-protein l-rmsd calculation."""
     protprot_caprimodule.calc_lrmsd()
-    assert round_two_dec(protprot_caprimodule.lrmsd) == 15.9
+    assert round_two_dec(protprot_caprimodule.lrmsd) == 20.94
 
 
 def test_protprot_ilrmsd(protprot_caprimodule):
     """Test protein-protein i-l-rmsd calculation."""
     protprot_caprimodule.calc_ilrmsd()
-    assert round_two_dec(protprot_caprimodule.ilrmsd) == 9.67
+    assert round_two_dec(protprot_caprimodule.ilrmsd) == 18.25
 
 
 def test_protprot_fnat(protprot_caprimodule):
@@ -198,13 +202,13 @@ def test_protlig_irmsd(protlig_caprimodule):
 def test_protlig_lrmsd(protlig_caprimodule):
     """Test protein-ligand l-rmsd calculation."""
     protlig_caprimodule.calc_lrmsd()
-    assert round_two_dec(protlig_caprimodule.lrmsd) == 0.51
+    assert round_two_dec(protlig_caprimodule.lrmsd) == 0.49
 
 
 def test_protlig_ilrmsd(protlig_caprimodule):
     """Test protein-ligand i-l-rmsd calculation."""
     protlig_caprimodule.calc_ilrmsd()
-    assert round_two_dec(protlig_caprimodule.ilrmsd) == 0.5
+    assert round_two_dec(protlig_caprimodule.ilrmsd) == 0.49
 
 
 def test_protlig_fnat(protlig_caprimodule):
@@ -222,13 +226,13 @@ def test_protdna_irmsd(protdna_caprimodule):
 def test_protdna_lrmsd(protdna_caprimodule):
     """Test protein-dna l-rmsd calculation."""
     protdna_caprimodule.calc_lrmsd()
-    assert round_two_dec(protdna_caprimodule.lrmsd) == 4.19
+    assert round_two_dec(protdna_caprimodule.lrmsd) == 6.13
 
 
 def test_protdna_ilrmsd(protdna_caprimodule):
     """Test protein-dna i-l-rmsd calculation."""
     protdna_caprimodule.calc_ilrmsd()
-    assert round_two_dec(protdna_caprimodule.ilrmsd) == 1.89
+    assert round_two_dec(protdna_caprimodule.ilrmsd) == 5.97
 
 
 def test_protdna_fnat(protdna_caprimodule, protdna_input_list):