fix irmsd and lrmsd pdb2ql routine

Author: NicoRenaud

pdb2sql/StructureSimilarity.py

@@ -254,8 +254,6 @@ def compute_irmsd_fast(
         # read the izone file
         if izone is None:
             resData = self.compute_izone(cutoff, save_file=False)
-        # elif not os.path.isfile(izone):
-        #     resData = self.compute_izone(cutoff,save_file=True,filename=izone)
         else:
             resData = self.read_zone(izone)
 
@@ -501,37 +499,40 @@ def compute_lrmsd_pdb2sql(self, exportpath=None, method='svd'):
         backbone = ['CA', 'C', 'N', 'O']
 
         # create the sql
-        sql_decoy = pdb2sql(
-            self.decoy, sqlfile='decoy.db', fix_chainID=True)
-        sql_ref = pdb2sql(
-            self.ref, sqlfile='ref.db', fix_chainID=True)
+        sql_decoy = pdb2sql(self.decoy, sqlfile='decoy.db')
+        sql_ref = pdb2sql(self.ref, sqlfile='ref.db')
+
+        # get the chains
+        chains_decoy = sql_decoy.get_chains()
+        chains_ref = sql_ref.get_chains()
+
+        if chains_decoy != chains_ref:
+            raise ValueError(
+                'Chains are different in decoy and reference structure')
+
+        chain1 = chains_decoy[0]
+        chain2 = chains_decoy[1]
 
         # extract the pos of chains A
         xyz_decoy_A = np.array(
-            sql_decoy.get(
-                'x,y,z',
-                chainID='A',
-                name=backbone))
+            sql_decoy.get('x,y,z', chainID=chain1, name=backbone))
         xyz_ref_A = np.array(sql_ref.get(
-            'x,y,z', chainID='A', name=backbone))
+            'x,y,z', chainID=chain1, name=backbone))
 
         # extract the pos of chains B
         xyz_decoy_B = np.array(
-            sql_decoy.get(
-                'x,y,z',
-                chainID='B',
-                name=backbone))
+            sql_decoy.get('x,y,z', chainID=chain2, name=backbone))
         xyz_ref_B = np.array(sql_ref.get(
-            'x,y,z', chainID='B', name=backbone))
+            'x,y,z', chainID=chain2, name=backbone))
 
         # check the lengthes
         if len(xyz_decoy_A) != len(xyz_ref_A):
             xyz_decoy_A, xyz_ref_A = self.get_identical_atoms(
-                sql_decoy, sql_ref, 'A')
+                sql_decoy, sql_ref, chain1)
 
         if len(xyz_decoy_B) != len(xyz_ref_B):
             xyz_decoy_B, xyz_ref_B = self.get_identical_atoms(
-                sql_decoy, sql_ref, 'B')
+                sql_decoy, sql_ref, chain2)
 
         # detect which chain is the longest
         nA, nB = len(xyz_decoy_A), len(xyz_decoy_B)
@@ -688,14 +689,26 @@ def compute_irmsd_pdb2sql(
         """
 
         # create thes sql
-        sql_decoy = interface(self.decoy, fix_chainID=True)
-        sql_ref = interface(self.ref, fix_chainID=True)
+        sql_decoy = interface(self.decoy)
+        sql_ref = interface(self.ref)
+
+        # get the chains
+        chains_decoy = sql_decoy.get_chains()
+        chains_ref = sql_ref.get_chains()
+
+        if chains_decoy != chains_ref:
+            raise ValueError(
+                'Chains are different in decoy and reference structure')
 
         # get the contact atoms
         if izone is None:
+
             contact_ref = sql_ref.get_contact_atoms(
                 cutoff=cutoff,
-                extend_to_residue=True)
+                extend_to_residue=True,
+                chain1=chains_ref[0],
+                chain2=chains_decoy[1])
+
             index_contact_ref = []
             for v in contact_ref.values():
                 index_contact_ref += v

pdb2sql/interface.py

@@ -3,6 +3,7 @@
 import warnings
 from .pdb2sqlcore import pdb2sql
 
+
 class interface(pdb2sql):
 
     def __init__(self, pdb, **kwargs):
@@ -73,6 +74,12 @@ def get_contact_atoms(
         else:
             chainIDs = [chain1, chain2]
 
+        chains = self.get_chains()
+        for c in chainIDs:
+            if c not in chains:
+                raise ValueError(
+                    'chain %s not found in the structure' % c)
+
         xyz = dict()
         index = dict()
         resName = dict()
@@ -265,7 +272,8 @@ def get_contact_residues(
                     residue_contact_pairs[data1] = set()
 
                 # get the res info of the atom in the other chain
-                data2 = self.get('chainID,resSeq,resName', rowID=atoms2)
+                data2 = self.get(
+                    'chainID,resSeq,resName', rowID=atoms2)
 
                 # store that in the dict without double
                 for resData in data2:
@@ -300,4 +308,4 @@ def get_contact_residues(
                 residue_contact[chain] = sorted(
                     set([tuple(resData) for resData in data[chain]]))
 
-            return residue_contact
+            return residue_contact