Added functionality for create rotlibs and RotamerEnsembles for macrmolecules with arbitrary backbones

README.md

@@ -1,3 +1,16 @@
+<p align="center">
+   <a href="https://stolllab.github.io/chiLife/">
+      <img src="https://img.shields.io/badge/docs-latest-brightgreen.svg" height="60">
+   </a>
+</p>
+
+
+|                       :exclamation:   News   :exclamation:                                |                       
+|-------------------------------------------------------------------------------------------|
+| chiLife now supports arbitrary backbone attachments including DNA and RNA labels and more! Check out our new nucleic acid lables at https://github.com/mtessmer/chilife_rotlibs|
+
+
+
 # chiLife
 chiLife is a Python package for modeling non-canonical amino acid side chain ensembles and using those ensembles to
 predict experimental observables. Currently, it is focused primarily on site-directed spin labels (SDSLs). The goal of

scripts/update_rotlib.py

@@ -102,7 +102,11 @@ def get_data(filename, description = None, comment = None, reference = None):
         if (key not in data) or (val is not None):
             data[key] = val
 
-    data['format_version'] = 1.3
+    if 'backbone_atoms' not in data:
+        data['backbone_atoms'] = ["H", "N", "CA", "HA", "C", "O"]
+        data['aln_atoms'] = ['N', 'CA', 'C']
+
+    data['format_version'] = 1.4
 
     return data
 

src/chilife/MolSysIC.py

@@ -529,6 +529,7 @@ def get_dihedral(self, resi: int, atom_list: ArrayLike, chain: Union[int, str] =
                                    f'considers dihedrals as directional so you may want to try the reverse dihedral.')
 
             dihedral_idxs.append(list(self.topology.dihedrals_by_resnum[tag]))
+
         dihedral_idxs = np.array(dihedral_idxs).T
         dihedral_values = self.coords[dihedral_idxs]
         dihedrals = get_dihedrals(*dihedral_values)

src/chilife/RotamerEnsemble.py

@@ -22,6 +22,7 @@
 from .numba_utils import get_sasa as nu_getsasa
 from .alignment_methods import alignment_methods, parse_backbone, local_mx, global_mx
 from .protein_utils import FreeAtom, guess_mobile_dihedrals, get_dihedral, get_angle
+from .pdb_utils import get_bb_candidates, get_backbone_atoms
 from .MolSys import MolSys, MolecularSystemBase
 from .MolSysIC import MolSysIC
 
@@ -103,8 +104,6 @@ class RotamerEnsemble:
        """
 
 
-    backbone_atoms = ["H", "N", "CA", "HA", "C", "O"]
-
     def __init__(self, res, site=None, protein=None, chain=None, rotlib=None, **kwargs):
 
 
@@ -165,9 +164,9 @@ def __init__(self, res, site=None, protein=None, chain=None, rotlib=None, **kwar
             raise RuntimeError('The kwarg `forcefield` must be a string or ForceField object.')
 
         # Parse important indices
-        self.backbone_idx = np.argwhere(np.isin(self.atom_names, ["N", "CA", "C"]))
+        self.backbone_idx = np.squeeze(np.argwhere(np.isin(self.atom_names, self.aln_atoms)))
         self.side_chain_idx = np.argwhere(
-            np.isin(self.atom_names, RotamerEnsemble.backbone_atoms, invert=True)
+            np.isin(self.atom_names, self.backbone_atoms, invert=True)
         ).flatten()
 
         self._graph = ig.Graph(edges=self.bonds)
@@ -352,6 +351,23 @@ def from_trajectory(cls, traj, site, chain=None, energy=None, burn_in=0, **kwarg
         dihedrals = np.array([ic.get_dihedral(1, ddefs) for ic in ICs])
         sigmas = kwargs.get('sigmas', np.array([]))
 
+        bb_candidates = get_bb_candidates(res.names, resname)
+        aln_atoms = kwargs.get('ln_atoms', [])
+        if not bb_candidates and not aln_atoms:
+            raise RuntimeError('No suitable backbone candidates were found. Please use from_trajectory on residues with'
+                               'standard backbone atom names or specify alignment backbone atoms with the aln_atoms '
+                               'keyword argument.')
+        elif aln_atoms:
+            graph = res.topology.graph
+            root_idx = np.argwhere(res.names == aln_atoms[1]).flat[0]
+            aname_lst = ICs.atom_names.tolist()
+            neighbor_idx = [aname_lst.index(a) for a in aln_atoms[::2]]
+            backbone_atoms = get_backbone_atoms(graph, root_idx, neighbor_idx)
+
+        else:
+            backbone_atoms = [b for b in bb_candidates if b in res.names]
+            aln_atoms = ['N', 'CA', 'C'] if np.all(np.isin(['N', 'CA', 'C'], backbone_atoms)) else ["O4'", "C1'", "C2'"]
+
         lib = {'rotlib': f'{resname}_from_traj',
                'resname': resname,
                'coords': np.atleast_3d(coords),
@@ -361,6 +377,8 @@ def from_trajectory(cls, traj, site, chain=None, energy=None, burn_in=0, **kwarg
                'atom_names': res.names.copy(),
                'dihedral_atoms': ddefs,
                'dihedrals': np.rad2deg(dihedrals),
+               'aln_atoms': aln_atoms,
+               'backbone_atoms': backbone_atoms,
                '_dihedrals': dihedrals.copy(),
                '_rdihedrals': dihedrals,
                'sigmas': sigmas,
@@ -496,12 +514,14 @@ def to_rotlib(self,
                'atom_names': self.atom_names.copy(),
                'dihedral_atoms': self.dihedral_atoms,
                'dihedrals': self.dihedrals,
+               'aln_atoms': self.aln_atoms,
+               'backbone_atoms': self.backbone_atoms,
                'sigmas': self.sigmas,
                'type': 'chilife rotamer library',
                'description': description,
                'comment': comment,
                'reference': reference,
-               'format_version': 1.3}
+               'format_version': 1.4}
 
         if hasattr(self, 'spin_atoms'):
             lib['spin_atoms'] = self.spin_atoms
@@ -581,15 +601,20 @@ def to_site(self, site_pos: ArrayLike = None) -> None:
             3x3 array of ordered backbone atom coordinates of new site (N CA C) (Default value = None)
         """
         if site_pos is None:
-            N, CA, C = parse_backbone(self, kind="global")
+            bbs = parse_backbone(self, kind="global")
         else:
-            N, CA, C = site_pos
+            bbs = site_pos
+
+        if len(bbs) < 3 :
+            raise RuntimeError(f'The residue/rotlib you have selected {self.res}, does not share a compatible backbone '
+                               f'with the residue you are trying to label. Check the site and rotlib and try again.\n'
+                               f'The label backbone atoms: {self.backbone_atoms}')
 
-        mx, ori = global_mx(N, CA, C, method=self.alignment_method)
+        mx, ori = global_mx(*bbs, method=self.alignment_method)
 
         # if self.alignment_method not in {'fit', 'rosetta'}:
-        N, CA, C = parse_backbone(self, kind="local")
-        old_ori, ori_mx = local_mx(N, CA, C, method=self.alignment_method)
+        bbs = parse_backbone(self, kind="local")
+        old_ori, ori_mx = local_mx(*bbs, method=self.alignment_method)
         self._coords -= old_ori
         cmx = ori_mx @ mx
 
@@ -659,7 +684,10 @@ def ICs_to_site(self, cori, cmx):
     def backbone_to_site(self):
         """Modify backbone atoms to match the site that the RotamerEnsemble is attached to """
         # Keep protein backbone dihedrals for oxygen and hydrogen atoms
-        for atom in ["H", "O"]:
+
+        for atom in self.backbone_atoms:
+
+            if atom in self.aln_atoms: continue  # TODO update test to change this
             mask = self.atom_names == atom
             if any(mask) and self.protein is not None:
 

src/chilife/Topology.py

@@ -9,6 +9,7 @@
 import igraph as ig
 
 from .globals import bond_hmax_dict
+from .math_utils import simple_cycle_vertices
 
 class Topology:
     """
@@ -63,17 +64,7 @@ def __init__(self, mol, bonds, **kwargs):
 
     @property
     def ring_idxs(self):
-        found_cycle_edges = self.graph.fundamental_cycles()
-        found_cycle_nodes = []
-
-        for cycle in found_cycle_edges:
-            cyverts = set()
-            for edge in self.graph.es(cycle):
-                cyverts.update(edge.tuple)
-
-            found_cycle_nodes.append(sorted(cyverts))
-
-        return found_cycle_nodes
+        return simple_cycle_vertices(self.graph)
 
     @property
     def has_rings(self):
@@ -283,7 +274,7 @@ def neighbors(edges, node):
     return nbs
 
 
-def bfs_edges(edges, root):
+def modified_bfs_edges(edges, root, bb_idxs):
     """
     Breadth first search of nodes given a set of edges
     Parameters
@@ -307,18 +298,22 @@ def bfs_edges(edges, root):
 
     n = len(nodes)
     depth = 0
-    next_parents_children = [(root, neighbors(edges, root))]
-
-    while next_parents_children and depth < depth_limit:
-        this_parents_children = next_parents_children
-        next_parents_children = []
-        for parent, children in this_parents_children:
-            for child in children:
-                if child not in seen:
-                    seen.add(child)
-                    next_parents_children.append((child, neighbors(edges, child)))
-                    yield parent, child
-            if len(seen) == n:
-                return
-        depth += 1
-
+    neigh = neighbors(edges, root)
+    # Prioritize side chains
+    neigh1 = [n for n in neigh if n not in bb_idxs]
+    neigh2 = [n for n in neigh if n in bb_idxs]
+
+    for neigh in neigh1, neigh2:
+        next_parents_children = [(root, neigh)]
+        while next_parents_children and depth < depth_limit:
+            this_parents_children = next_parents_children
+            next_parents_children = []
+            for parent, children in this_parents_children:
+                for child in children:
+                    if child not in seen:
+                        seen.add(child)
+                        next_parents_children.append((child, neighbors(edges, child)))
+                        yield parent, child
+                if len(seen) == n:
+                    return
+            depth += 1

src/chilife/__init__.py

@@ -26,4 +26,4 @@
 # SpinLabel = SpinLabel.SpinLabel
 # dSpinLabel = dSpinLabel.dSpinLabel
 
-__version__ = '0.3.0'
+__version__ = '0.4.0'

src/chilife/alignment_methods.py

@@ -230,13 +230,13 @@ def parse_backbone(rotamer_ensemble, kind):
         with the first coordinate as the rotamer ensemble backbone and the second as the protein site backbone.
     """
     method = rotamer_ensemble.alignment_method
-
+    aln_atoms = " ".join(rotamer_ensemble.aln_atoms)
     if method.__name__ == "fit_alignment":
         N1, CA1, C1 = rotamer_ensemble.backbone
         N2, CA2, C2 = rotamer_ensemble.protein.select_atoms(
             f"segid {rotamer_ensemble.chain} and "
             f"resnum {rotamer_ensemble.site} "
-            f"and name N CA C and not altloc B"
+            f"and name {aln_atoms} and not altloc B"
         ).positions
         return np.array([[N1, N2], [CA1, CA2], [C1, C2]])
 
@@ -247,7 +247,7 @@ def parse_backbone(rotamer_ensemble, kind):
         return rotamer_ensemble.protein.select_atoms(
             f"segid {rotamer_ensemble.chain} and "
             f"resnum {rotamer_ensemble.site} "
-            f"and name N CA C and not altloc B"
+            f"and name {aln_atoms} and not altloc B"
         ).positions
 
 
@@ -322,4 +322,5 @@ def global_mx(*p: ArrayLike, method: Union[str, callable] = "bisect") -> Tuple[A
         p = [pi[::-1] for pi in p]
 
     rotation_matrix, origin = method(*p)
+
     return rotation_matrix, origin