Merge pull request #191 from UC-Davis-molecular-computing/dev

Dev

.gitignore

@@ -1,3 +1,4 @@
+*.DS_Store
 *.pyc
 __pyache__/*
 /venv/

scadnano/scadnano.py

@@ -54,7 +54,7 @@
 # commented out for now to support Py3.6, which does not support this feature
 # from __future__ import annotations
 
-__version__ = "0.16.2"  # version line; WARNING: do not remove or change this line or comment
+__version__ = "0.16.3"  # version line; WARNING: do not remove or change this line or comment
 
 import dataclasses
 from abc import abstractmethod, ABC, ABCMeta
@@ -82,7 +82,7 @@ def _pairwise(iterable: Iterable) -> Iterable:
     """s -> (s0,s1), (s1,s2), (s2, s3), ..."""
     a, b = itertools.tee(iterable)
     next(b, None)
-    return zip(a, b) 
+    return zip(a, b)
 
 
 # for putting evaluated expressions in docstrings
@@ -3862,7 +3862,7 @@ def from_scadnano_file(filename: str) -> 'Design':  # remove quotes when Py3.6 s
         """
         Loads a :any:`Design` from the file with the given name.
 
-        :param filename: name of the file with the design. Should be a JSON file ending in .dna
+        :param filename: name of the file with the design. Should be a JSON file ending in .sc
         :return: Design described in the file
         """
         with open(filename) as f:
@@ -4627,14 +4627,25 @@ def _cadnano_v2_place_crossover(helix_from_dct: Dict[str, Any], helix_to_dct: Di
         start_from, end_from, forward_from = domain_from.start, domain_from.end, domain_from.forward
 
         helix_to = helix_to_dct['num']
-        start_to, end_to = domain_to.start, domain_to.end
+        start_to, end_to, forward_to = domain_to.start, domain_to.end, domain_to.forward
 
-        if forward_from:
+        # Because of paranemic crossovers it is possible
+        # to crossover to a strand that goes in the same direction
+        # In total there are four cases corresponding to
+        # [forward_from, not forward_from] x [forward_to, not forward_to]
+        if forward_from and not forward_to:
             helix_from_dct[strand_type][end_from - 1][2:] = [helix_to, end_to - 1]
             helix_to_dct[strand_type][end_to - 1][:2] = [helix_from, end_from - 1]
-        else:
+        elif not forward_from and forward_to:
             helix_from_dct[strand_type][start_from][2:] = [helix_to, start_to]
             helix_to_dct[strand_type][start_to][:2] = [helix_from, start_from]
+        elif forward_from and forward_to:
+            helix_from_dct[strand_type][end_from - 1][2:] = [helix_to, start_to]
+            helix_to_dct[strand_type][end_to - 1][:2] = [helix_from, start_from]
+        elif not forward_from and not forward_to:
+            helix_from_dct[strand_type][start_from][2:] = [helix_to, end_to - 1]
+            helix_to_dct[strand_type][start_to][:2] = [helix_from, end_from - 1]
+
 
     @staticmethod
     def _cadnano_v2_color_of_stap(color: Color, domain: Domain) -> List[int]:
@@ -4706,7 +4717,8 @@ def _cadnano_v2_fill_blank(self, dct: dict, num_bases: int, design_grid: Grid) -
         """Creates blank cadnanov2 helices in and initialized all their fields.
         """
         helices_ids_reverse = {}
-        for i, helix in self.helices.items():
+        i = 0
+        for _, helix in self.helices.items():
             helix_dct: Dict[str, Any] = OrderedDict()
             helix_dct['num'] = helix.idx
 
@@ -4732,6 +4744,7 @@ def _cadnano_v2_fill_blank(self, dct: dict, num_bases: int, design_grid: Grid) -
 
             helices_ids_reverse[helix_dct['num']] = i
             dct['vstrands'].append(helix_dct)
+            i += 1
         return helices_ids_reverse
 
     def to_cadnano_v2(self) -> Dict[str, Any]:
@@ -5639,7 +5652,6 @@ def _write_plates_default(self, directory: str, filename: Optional[str], strands
         filename_plate, workbook = self._setup_excel_file(directory, filename)
         worksheet = self._add_new_excel_plate_sheet(f'plate{plate}', workbook)
 
-
         num_strands_per_plate = plate_type.num_wells_per_plate()
         num_plates_needed = len(strands) // num_strands_per_plate
         if len(strands) % num_strands_per_plate != 0:
@@ -5676,8 +5688,8 @@ def _write_plates_default(self, directory: str, filename: Optional[str], strands
             # So if we would have fewer than that many on the last plate, 
             # shift some from the penultimate plate.
             if not on_final_plate and \
-                final_plate_less_than_min_required and \
-                num_strands_remaining == min_strands_per_plate:
+                    final_plate_less_than_min_required and \
+                    num_strands_remaining == min_strands_per_plate:
                 plate_coord.advance_to_next_plate()
             else:
                 plate_coord.increment()
@@ -5696,7 +5708,7 @@ def to_oxdna_format(self) -> Tuple[str, str]:
         """Exports to oxdna format.
 
         :return:
-            two strings that are the contents of the .conf and .top file
+            two strings that are the contents of the .dat and .top file
             suitable for reading by oxdna (https://sulcgroup.github.io/oxdna-viewer/)
         """
         system = _convert_design_to_oxdna_system(self)
@@ -5706,21 +5718,21 @@ def write_oxdna_files(self, directory: str = '.', filename_no_extension: Optiona
         """Write text file representing this :any:`Design`,
         suitable for reading by oxdna (https://sulcgroup.github.io/oxdna-viewer/),
         with the output files having the same name as the running script but with ``.py`` changed to
-        ``.conf`` and ``.top``,
+        ``.dat`` and ``.top``,
         unless `filename_no_extension` is explicitly specified.
         For instance, if the script is named ``my_origami.py``,
-        then the design will be written to ``my_origami.conf`` and ``my_origami.top``.
+        then the design will be written to ``my_origami.dat`` and ``my_origami.top``.
 
         The strings written are those returned by :meth:`Design.to_oxdna_format`.
 
         :param directory:
             directory in which to put file (default: current working directory)
         :param filename_no_extension:
-            filename without extension (default: name of script without ``.py``).
+            filename without extension (default: name of running script without ``.py``).
         """
-        conf, top = self.to_oxdna_format()
+        dat, top = self.to_oxdna_format()
 
-        _write_file_same_name_as_running_python_script(conf, 'conf', directory, filename_no_extension)
+        _write_file_same_name_as_running_python_script(dat, 'dat', directory, filename_no_extension)
         _write_file_same_name_as_running_python_script(top, 'top', directory, filename_no_extension)
 
     # @_docstring_parameter was used to substitute sc in for the filename extension, but it is
@@ -6413,10 +6425,11 @@ def _write_file_same_name_as_running_python_script(contents: str, extension: str
 
 
 def _get_filename_same_name_as_running_python_script(directory: str, extension: str,
-                                                     filename: Optional[str]) -> str:
-    if filename is None:
-        filename = _name_of_this_script() + f'.{extension}'
-    relative_filename = _create_directory_and_set_filename(directory, filename)
+                                                     filename_no_extension: Optional[str]) -> str:
+    if filename_no_extension is None:
+        filename_no_extension = _name_of_this_script()
+    filename_with_extension = f'{filename_no_extension}.{extension}'
+    relative_filename = _create_directory_and_set_filename(directory, filename_with_extension)
     return relative_filename
 
 
@@ -6492,29 +6505,36 @@ def rotate(self, angle: float, axis: "_OxdnaVector") -> "_OxdnaVector":
 _OXDNA_ORIGIN = _OxdnaVector(0, 0, 0)
 
 
-# r, b, and n represent the oxDNA conf file vectors that describe a nucleotide
+# r, b, and n represent the oxDNA configuration (.dat) file vectors that describe a nucleotide
 # r is the position of the base
 # b is the backbone-base vector (in documentation as versor: more info on versors here https://eater.net/quaternions)
 # n is the forward direction of the helix
 @dataclass(frozen=True)
 class _OxdnaNucleotide:
-    center: _OxdnaVector   # center of the slice of the helix to which this nucleotide belongs
-    normal: _OxdnaVector   # unit vector from the backbone to the center of the helix, same as negated b from oxDNA configuration file
-    forward: _OxdnaVector  # unit vector pointing in direction from 3' to 5' ends of the helix, same as oxDNA n vector
+    center: _OxdnaVector
+    # center of the slice of the helix to which this nucleotide belongs
+
+    normal: _OxdnaVector
+    # unit vector from the backbone to the center of the helix, same as negated b from oxDNA dat file
+
+    forward: _OxdnaVector
+    # unit vector pointing in direction from 3' to 5' ends of the helix, same as oxDNA n vector
+
     base: str
 
-    v: _OxdnaVector = field(init=False, default=_OXDNA_ORIGIN)  # velocity for oxDNA conf file
-    L: _OxdnaVector = field(init=False, default=_OXDNA_ORIGIN)  # angular velocity for oxDNA conf file
+    v: _OxdnaVector = field(init=False, default=_OXDNA_ORIGIN)  # velocity for oxDNA dat file
+    L: _OxdnaVector = field(init=False, default=_OXDNA_ORIGIN)  # angular velocity for oxDNA dat file
 
     # https://dna.physics.ox.ac.uk/index.php/Documentation
-    # r, b, and n represent the oxDNA conf file vectors that describe a nucleotide
+    # r, b, and n represent the oxDNA dat file vectors that describe a nucleotide
 
-    # r is the position of the center of mass of the oxDNA binding sites (stacking, hydrogen bonding, backbone-repulsion https://dna.physics.ox.ac.uk/index.php/Documentation)
+    # r is the position of the center of mass of the oxDNA binding sites
+    # (stacking, hydrogen bonding, backbone-repulsion https://dna.physics.ox.ac.uk/index.php/Documentation)
     # offset from the center of the helix by _BASE_DIST
     @property
     def r(self) -> _OxdnaVector:
         return self.center - self.b * _BASE_DIST
-    
+
     # b is the backbone-base vector (in documentation as versor: more info on versors here https://eater.net/quaternions)
     @property
     def b(self) -> _OxdnaVector:
@@ -6525,6 +6545,7 @@ def b(self) -> _OxdnaVector:
     def n(self) -> _OxdnaVector:
         return self.forward
 
+
 @dataclass(frozen=True)
 class _OxdnaStrand:
     nucleotides: List[_OxdnaNucleotide] = field(default_factory=list)
@@ -6556,7 +6577,7 @@ def compute_bounding_box(self, cubic: bool = True) -> _OxdnaVector:
             # 5 is arbitrarily chosen so that the box has a bit of wiggle room
             # 1.5 multiplier is to make all crossovers appear (advice from Oxdna authors)
             box = 1.5 * (max_vec - min_vec + _OxdnaVector(5, 5, 5))
-            if cubic: # oxDNA requires cubic bounding box with default simulation options
+            if cubic:  # oxDNA requires cubic bounding box with default simulation options
                 max_side = max(box.x, box.y, box.z)
                 box = _OxdnaVector(max_side, max_side, max_side)
             return box
@@ -6625,7 +6646,8 @@ def _oxdna_get_helix_vectors(design: Design, helix: Helix) -> Tuple[_OxdnaVector
     forward = forward.rotate(design.yaw_of_helix(helix), normal)
     forward = forward.rotate(-design.pitch_of_helix(helix), _OxdnaVector(1, 0, 0))
     normal = normal.rotate(-design.pitch_of_helix(helix), _OxdnaVector(1, 0, 0))
-    normal = normal.rotate(helix.roll, forward)
+
+    normal = normal.rotate(-helix.roll, forward)
 
     x: float = 0.0
     y: float = 0.0
@@ -6719,7 +6741,7 @@ def _convert_design_to_oxdna_system(design: Design) -> _OxdnaSystem:
                 # we apply the opposite rotation so that we get the expected vector from scadnano in oxDNA
                 groove_gamma_correction = _GROOVE_GAMMA if domain.forward else -_GROOVE_GAMMA
                 normal = normal.rotate(groove_gamma_correction, forward).normalize()
-                    
+
                 # dict / set for insertions / deletions to make lookup cheaper when there are lots of them
                 deletions = set(domain.deletions)
                 insertions = dict(domain.insertions)
@@ -6740,14 +6762,14 @@ def _convert_design_to_oxdna_system(design: Design) -> _OxdnaSystem:
                                 r = origin_ + forward * (
                                         offset + mod - num + i) * geometry.rise_per_base_pair * NM_TO_OX_UNITS
                                 b = normal.rotate(step_rot * (offset + mod - num + i), forward)
-                                n = -forward if domain.forward else forward # note oxDNA n vector points 3' to 5' opposite of scadnano forward vector
+                                n = -forward if domain.forward else forward  # note oxDNA n vector points 3' to 5' opposite of scadnano forward vector
                                 nuc = _OxdnaNucleotide(r, b, n, seq[index])
                                 dom_strand.nucleotides.append(nuc)
                                 index += 1
 
                         r = origin_ + forward * (offset + mod) * geometry.rise_per_base_pair * NM_TO_OX_UNITS
                         b = normal.rotate(step_rot * (offset + mod), forward)
-                        n = -forward if domain.forward else forward # note oxDNA n vector points 3' to 5' opposite of scadnano forward vector
+                        n = -forward if domain.forward else forward  # note oxDNA n vector points 3' to 5' opposite of scadnano forward vector
                         nuc = _OxdnaNucleotide(r, b, n, seq[index])
                         dom_strand.nucleotides.append(nuc)
                         index += 1

tests/scadnano_tests.py

@@ -610,6 +610,13 @@ def test_circular_auto_staple_hex(self) -> None:
         design.write_scadnano_file(directory=self.output_path,
                                    filename=f'{file_name}.{sc.default_scadnano_file_extension}')
 
+    def test_paranemic_crossover(self) -> None:
+        file_name = "test_paranemic_crossover"
+        design = sc.Design.from_cadnano_v2(directory=self.input_path,
+                                           filename=file_name + ".json")
+        design.write_scadnano_file(directory=self.output_path,
+                                   filename=f'{file_name}.{sc.default_scadnano_file_extension}')
+
 
 class TestExportDNASequences(unittest.TestCase):
 
@@ -1034,6 +1041,12 @@ def test_big_circular_staples(self) -> None:
         design.export_cadnano_v2(directory=self.output_path,
                                  filename='test_big_circular_staples.json')
 
+    def test_paranemic_crossover(self) -> None:
+        design = sc.Design.from_scadnano_file(
+            os.path.join(self.input_path, f'test_paranemic_crossover.{self.ext}'))
+        design.export_cadnano_v2(directory=self.output_path,
+                                 filename='test_paranemic_crossover.json')
+
     def test_bad_cases(self) -> None:
         """ We do not handle Loopouts and design where the parity of the helix
         does not correspond to the direction.