SCHWED-877: add special treatment for peptide ligands

sidechain/pymod/_reconstruct_sidechains.py

@@ -15,6 +15,46 @@ def _GetRotamerIDs(res_list):
         rotamer_ids[i] = rot_id
     return rotamer_ids
 
+def _GetPhiAngle(r):
+    '''Extract phi angle for residue r.'''
+    # def. fallback = helix
+    phi = -1.0472
+    # try to get phi from torsion angles
+    tor = r.GetPhiTorsion()
+    if tor.IsValid():
+        phi = tor.GetAngle()
+    else:
+        r_prev = r.handle.prev
+        if r_prev.IsValid() and mol.InSequence(r_prev, r.handle):
+            c_prev = r_prev.FindAtom("C")
+            n = r.FindAtom("N")
+            ca = r.FindAtom("CA")
+            c = r.FindAtom("C")
+            if c_prev.IsValid() and n.IsValid() and ca.IsValid() and c.IsValid():
+                phi = geom.DihedralAngle(c_prev.GetPos(),n.GetPos(),
+                                         ca.GetPos(),c.GetPos())
+    return phi
+
+def _GetPsiAngle(r):
+    '''Extract psi angle for residue r.'''
+    # def. fallback = helix 
+    psi = -0.7854
+    # try to get psi from torsion angles
+    tor = r.GetPsiTorsion()
+    if tor.IsValid():
+        psi = tor.GetAngle()
+    else:
+        r_next = r.handle.next
+        if r_next.IsValid() and mol.InSequence(r.handle, r_next):
+            n = r.FindAtom("N")
+            ca = r.FindAtom("CA")
+            c = r.FindAtom("C")
+            n_next = r_next.FindAtom("N")
+            if n.IsValid() and ca.IsValid() and c.IsValid() and n_next.IsValid():
+                psi = geom.DihedralAngle(n.GetPos(), ca.GetPos(),
+                                         c.GetPos(), n_next.GetPos())
+    return psi
+
 def _GetDihedrals(res_list):
     '''Extract dihedral angles for all residues.
     Returns phi and psi angles as 2 lists with length = len(res_list).
@@ -22,42 +62,8 @@ def _GetDihedrals(res_list):
     phi_angles = [0.0] * len(res_list)
     psi_angles = [0.0] * len(res_list)
     for i,r in enumerate(res_list):
-        # def. fallback = helix
-        phi = -1.0472 
-        psi = -0.7854
-
-        # try to get phi from torsion angles
-        tor = r.GetPhiTorsion()
-        if tor.IsValid():
-            phi = tor.GetAngle()
-        elif i > 0:
-            r_prev = res_list[i-1]
-            if mol.InSequence(r_prev.handle,r.handle):
-                c_prev = r_prev.FindAtom("C")
-                n = r.FindAtom("N")
-                ca = r.FindAtom("CA")
-                c = r.FindAtom("C")
-                if c_prev.IsValid() and n.IsValid() and ca.IsValid() and c.IsValid():
-                    phi = geom.DihedralAngle(c_prev.GetPos(),n.GetPos(),
-                                             ca.GetPos(),c.GetPos())
-        
-        # try to get psi from torsion angles
-        tor = r.GetPsiTorsion()
-        if tor.IsValid():
-            psi = tor.GetAngle()
-        elif i < (len(res_list) - 1):
-            r_next = res_list[i+1]
-            if mol.InSequence(r.handle,r_next.handle):
-                n = r.FindAtom("N")
-                ca = r.FindAtom("CA")
-                c = r.FindAtom("C")
-                n_next = r_next.FindAtom("N")
-                if n.IsValid() and ca.IsValid() and c.IsValid() and n_next.IsValid():
-                    psi = geom.DihedralAngle(n.GetPos(), ca.GetPos(),
-                                             c.GetPos(), n_next.GetPos())
-        # store
-        phi_angles[i] = phi
-        psi_angles[i] = psi
+        phi_angles[i] = _GetPhiAngle(r)
+        psi_angles[i] = _GetPsiAngle(r)
     return phi_angles, psi_angles
 
 def _AddBackboneFrameResidues(frame_residues, res_list, rotamer_ids,
@@ -76,14 +82,46 @@ def _AddLigandFrameResidues(frame_residues, ent_lig, rotamer_settings, offset):
     '''Update frame_residues (list) with FrameResidues for res. in ent_lig.
     Set offset >= number of non-ligand residues (used for residue_index).
     '''
-    # TODO? Split peptides and non-peptides?
-    # for now: add as static carbons
-    for i,r in enumerate(ent_lig.residues):
-        try:
-            frame_residue = sidechain.ConstructFrameResidue(r.handle, offset+i)
-            frame_residues.append(frame_residue)
-        except:
-            continue
+    # parse ligand residues
+    for i, res in enumerate(ent_lig.residues):
+        res_idx = offset + i
+        is_done = False
+        # special treatment for peptides
+        if res.IsPeptideLinking():
+            rot_id = sidechain.TLCToRotID(res.GetName())
+            if rot_id != sidechain.XXX:
+                # get more info
+                phi = _GetPhiAngle(res)
+                r_prev = res.handle.prev
+                n_ter = not r_prev.IsValid() \
+                        or not mol.InSequence(r_prev, res.handle)
+                r_next = res.handle.next
+                c_ter = not r_next.IsValid() \
+                        or not mol.InSequence(res.handle, r_next)
+                # try to add frame residues (ignore exceptions)
+                try:
+                    fr1 = sidechain.ConstructBackboneFrameResidue(\
+                              res.handle, rot_id, res_idx, rotamer_settings,
+                              phi, n_ter, c_ter)
+                    if rot_id != sidechain.ALA and rot_id != sidechain.GLY:
+                        fr2 = sidechain.ConstructSidechainFrameResidue(\
+                                  res.handle, rot_id, res_idx, rotamer_settings)
+                        frame_residues.extend([fr1,fr2])
+                    else:
+                        frame_residues.append(fr1)
+                except:
+                    pass   # ignore peptide treatment and treat below
+                else:
+                    is_done = True
+        # if it failed, treat it as an unknown entity
+        if not is_done:
+            # try to add frame residues (skip exceptions)
+            try:
+                # TODO: allow non-carbon treatment
+                fr = sidechain.ConstructFrameResidue(res.handle, res_idx)
+                frame_residues.append(fr)
+            except:
+                continue
 
 def _AddSidechainFrameResidues(frame_residues, incomplete_sidechains,
                                keep_sidechains, res_list, rotamer_ids,