Updated method of adding atomistic DNA nucleotides

CanonicalNucleotideAtoms.py

@@ -17,8 +17,50 @@ def stringToIntTuples(string, tupleLen, offset):
 
     return ret
 
+class Nucleotide(Group):
+    def __init__(self, children, orientation, sequence, resname, factory):
+        Group.__init__(self, children=children, orientation = orientation)
+        self.atoms_by_name = {a.name:a for a in self.children}
+        self._factory = factory
+        self.sequence = sequence
+        self.resname = resname
+        self._first_atomic_index = 0
+
+    def _get_atomic_index(self, name):
+        return self._factory._name_to_idx[name] + self._first_atomic_index
+
+    def get_intrahelical_above(self):
+        """ Returns next nucleotide along 5'-to-3' in same strand-side of helix """
+        return self.get_intrahelical_nucleotide(1)
+
+    def get_intrahelical_below(self):
+        """ Returns last nucleotide along 5'-to-3' in same strand-side of helix """
+        return self.get_intrahelical_nucleotide(-1)
+
+    def get_intrahelical_nucleotide(self,offset=1):
+        """ Returns nucleotide in same strand-side of helix with some offset along 5'-to-3'
+
+        TODO: Does not yet cross intrahelical connections
+        """
+
+        strand_piece = self.parent
+        is_fwd = strand_piece.is_fwd
+        lo,hi = sorted([strand_piece.start,strand_piece.end])
+
+        if is_fwd:
+            nt_idx = strand_piece.children.index(self)+lo + offset
+        else:
+            nt_idx = hi-strand_piece.children.index(self) - offset
+
+        seg = strand_piece.segment
+        if nt_idx < 0 or nt_idx >= seg.num_nt:
+            return None
+        else:
+            return seg._get_atomic_nucleotide(nt_idx, is_fwd)
+
 class CanonicalNucleotideFactory(Group):
-    DefaultOrientation = rotationAboutAxis([0,0,1], 180)
+    # DefaultOrientation = rotationAboutAxis([0,0,1], 180)
+    DefaultOrientation = rotationAboutAxis([1,0,0], 180)
     def __init__(self, prefix, seq):
         self.sequence = seq     # TODO: used?
         self.resname = resnames[seq]
@@ -50,20 +92,23 @@ class CanonicalNucleotideFactory(Group):
         
         atoms = self.children
 
+
         # minIdx = np.min(self.indices) - 1
         self.bonds = np.loadtxt("%s.bonds.txt" % prefix, dtype='i4') # - minIdx
         self.angles = np.loadtxt("%s.angles.txt" % prefix, dtype='i4') # - minIdx
         self.dihedrals = np.loadtxt("%s.dihedrals.txt" % prefix, dtype='i4')# - minIdx
         self.impropers = np.loadtxt("%s.impropers.txt" % prefix, dtype='i4')# - minIdx
 
+        self._name_to_idx = { atoms[i].name:i for i in range(len(atoms)) }
+
     def generate(self):
-        nt = Group(orientation = CanonicalNucleotideFactory.DefaultOrientation)
-        nt.__dict__['sequence'] = self.sequence
-        nt.__dict__['resname'] = self.resname              
-        atoms = nt.children = [copy.copy(c) for c in self.children]
-        for a in atoms:
-            a.parent = nt
-        nt.atoms_by_name = {a.name:a for a in nt.children}
+        atoms = [copy.copy(c) for c in self.children]
+        nt = Nucleotide(
+            children = atoms,
+            orientation = CanonicalNucleotideFactory.DefaultOrientation,
+            sequence = self.sequence,
+            resname = self.resname,
+            factory = self)
         
         b = self.bonds
         a = self.angles

segmentmodel.py

@@ -482,10 +482,8 @@ class Segment(ConnectableElement, Group):
         # if self.ntAt5prime is None and self.ntAt3prime is None: key = key+"singlet"
 
         key = seq
-        if not is_fwd:
-            nt_dict = canonicalNtFwd
-        else:
-            nt_dict = canonicalNtRev
+        nt_dict = canonicalNtFwd if is_fwd else canonicalNtRev
+
         atoms = nt_dict[ key ].generate() # TODO: clone?
                         
         atoms.orientation = orientation.dot(atoms.orientation)