Made cadnano reader set azimuthal angle more accurately; some problems could...

Made cadnano reader set azimuthal angle more accurately; some problems could arise with cadnano2 files that are saved in cadnano2.5 due to incorrect assignment of eulerZ in cadnano2.5

mrdna/readers/cadnano_segments.py

@@ -183,10 +183,7 @@ class cadnano_part(SegmentModel):
                  **kwargs
     ):
         self.part = part
-        try:
-            self.lattice_type = part.getGridType().name.lower()
-        except:
-            print("WARNING: unable to set cadnano lattice type")
+        self.lattice_type = _get_lattice(part)
 
         self._cadnano_part_to_segments(part)
         # SegmentModel.__init__(self,...)
@@ -198,6 +195,31 @@ class cadnano_part(SegmentModel):
         SegmentModel.__init__(self, self.segments,
                               **kwargs)
 
+    def _get_helix_angle(self, helix_id, indices):
+        """ Get "start_orientation" for helix """
+        # import ipdb
+        # ipdb.set_trace()
+
+        """ FROM CADNANO2.5
+        + angle is CCW
+        - angle is CW
+        Right handed DNA rotates clockwise from 5' to 3'
+        we use the convention the 5' end starts at 0 degrees
+        and it's pair is minor_groove_angle degrees away
+        direction, hence the minus signs.  eulerZ
+        """
+
+        hp, bpr, tpr, eulerZ, mgroove = self.part.vh_properties.loc[helix_id,
+                                                                    ['helical_pitch',
+                                                                     'bases_per_repeat',
+                                                                     'turns_per_repeat',
+                                                                     'eulerZ',
+                                                                     'minor_groove_angle']]
+        twist_per_base = tpr*360./bpr
+        # angle = eulerZ - twist_per_base*indices + 0.5*mgroove + 180
+        angle = eulerZ + twist_per_base*indices - 0.5*mgroove
+        return angle
+
     def _cadnano_part_to_segments(self,part):
         try:
             from cadnano.cnenum import PointType
@@ -351,7 +373,8 @@ class cadnano_part(SegmentModel):
                 ## TODO get sequence from cadnano api
                 if zMid in strandOccupancies[0] and zMid in strandOccupancies[1]:
                     kwargs['num_bp'] = numBps
-                    start_orientation = rotationAboutAxis(np.array((0,0,1)), (zid1-3.25)*34)
+                    _angle = self._get_helix_angle(hid, zid1)
+                    start_orientation = rotationAboutAxis(np.array((0,0,1)), _angle)
                     seg = DoubleStrandedSegment(**kwargs,**posargs1, start_orientation = start_orientation)
                 elif zMid in strandOccupancies[0]:
                     kwargs['num_nt'] = numBps
@@ -612,16 +635,41 @@ def decode_cadnano_part(json_data):
     try:
         doc = Document()
         cadnano.fileio.v3decode.decode(doc, json_data)
+        decoder = 3
     except:
         doc = Document()
         cadnano.fileio.v2decode.decode(doc, json_data)
+        decoder = 2
 
     parts = [p for p in doc.getParts()]
     if len(parts) != 1:
         raise Exception("Only documents containing a single cadnano part are implemented at this time.")
     part = parts[0]
+
+    if decoder == 2:
+        """ It seems cadnano2.5 (as of ce6ff019) does not set the EulerZ for square lattice structures correctly, doing so here """
+        l = _get_lattice(part)
+        if l == 'square':
+            for id_num in part.getIdNums():
+                if part.vh_properties.loc[id_num,'eulerZ'] == 0:
+                    part.vh_properties.loc[id_num,'eulerZ'] = 360*(6/10.5)
+
     return part
 
+def _get_lattice(part):
+    lattice_type = None
+    _gt = part.getGridType()
+    try:
+        lattice_type = _gt.name.lower()
+    except:
+        if _gt == 1:
+            lattice_type = 'square'
+        elif _gt == 2:
+            lattice_type = 'honeycomb'
+        else:
+            print("WARNING: unable to determine cadnano part lattice type")
+    return lattice_type
+
 def package_archive( name, directory ):
     ...