diff --git a/mrdna/readers/cadnano_segments.py b/mrdna/readers/cadnano_segments.py index f5986074e8c633860c1dd6f041c26c41747c9246..a0e3c6715bb15dc858efe848aa42999b5f578bee 100644 --- a/mrdna/readers/cadnano_segments.py +++ b/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 ): ...