Fixed two issues with ssDNA cadnano reader.

First, ssDNA segments are now broken at crossovers
Second, ssDNA segments aren't allowed to stack intrahelically

mrdna/readers/cadnano_segments.py

@@ -226,7 +226,10 @@ class cadnano_part(SegmentModel):
 
         vh_list = []
         strand_list = []
-        self.xover_list = xover_list = []
+        xover_list = []
+        self.xovers_from = dict()
+        self.xovers_to = dict()
+
         try:
             numHID = part.getMaxIdNum() + 1
         except:
@@ -247,9 +250,16 @@ class cadnano_part(SegmentModel):
                 fwd_ss, rev_ss = part.getStrandSets(id_num)
                 fwd_idxs, fwd_colors  = fwd_ss.dump(xover_list)
                 rev_idxs, rev_colors  = rev_ss.dump(xover_list)
-            
                 strand_list.append((fwd_idxs, rev_idxs))
 
+            self.xovers_from[id_num] = []
+            self.xovers_to[id_num] = []
+
+        for xo in xover_list:
+            h1,f1,z1,h2,f2,z2 = xo
+            self.xovers_from[h1].append(xo)
+            self.xovers_to[h2].append(xo)
+
         ## Get lists of 5/3prime ends
         strands5 = [o.strand5p() for o in part.oligos()]
         strands3 = [o.strand3p() for o in part.oligos()]
@@ -274,7 +284,8 @@ class cadnano_part(SegmentModel):
             ## Build list of tuples containing (idx,length) of insertions/skips
             insertions = sorted( [(i[0],i[1].length()) for i in allInsertions[hid].items()],
                                  key=lambda x: x[0] )
-            
+
+            ## TODO: make the following code (until "regions = ...") more readable
             ## Build list of strand ends and list of mandatory node locations
             ends1,ends2 = self._helixStrandsToEnds(helixStrands)
 
@@ -290,9 +301,28 @@ class cadnano_part(SegmentModel):
             ends1,ends2 = [ [(e[i],e[i+1]) for i in range(0,len(e),2)] for e in (ends1,ends2) ]
 
             regions = combineRegionLists(ends1,ends2)
+
+            ## Split regions in event of ssDNA crossover
+            split_regions = []
+            for zid1,zid2 in regions:
+                zMid = int(0.5*(zid1+zid2))
+                if zMid in strandOccupancies[0] and zMid in strandOccupancies[1]:
+                    split_regions.append( (zid1,zid2) )
+                else:
+                    is_fwd = zMid in strandOccupancies[0]
+                    ends = [z for h,f,z in self._get_crossover_locations( hid, range(zid1+1,zid2), is_fwd )]
+                    z1 = zid1
+                    for z in sorted(ends):
+                        z2 = z
+                        if z2 > z1:
+                            split_regions.append( (z1,z2) )
+                            z1 = z2+1
+                    z2 = zid2
+                    if z2 > z1: split_regions.append( (z1,z2) )
+
             # if hid == 43:
             #     import pdb
-            for zid1,zid2 in regions:
+            for zid1,zid2 in split_regions:
                 zMid = int(0.5*(zid1+zid2))
                 assert( zMid in strandOccupancies[0] or zMid in strandOccupancies[1] )
 
@@ -482,6 +512,8 @@ class cadnano_part(SegmentModel):
             occ = self.strand_occupancies[hid]
             for i in range(len(segs)-1):
                 seg1,seg2 = [segs[j] for j in (i,i+1)]
+                if isinstance(seg1,SingleStrandedSegment) and isinstance(seg2,SingleStrandedSegment):
+                    continue
                 r1,r2 = [self.helix_ranges[hid][j] for j in (i,i+1)]
                 if r1[1]+1 == r2[0]:
                     ## TODO: handle nicks that are at intrahelical connections(?)
@@ -502,16 +534,32 @@ class cadnano_part(SegmentModel):
                         else:
                             seg2.connect_end3(end)
                         
+
+    def _get_crossover_locations(self, helix_idx, nt_idx_range, fwd_strand=None):
+        xos = []
+        def append_if_in_range(h,f,z):
+            if fwd_strand in (None,f) and z in nt_idx_range:
+                xos.append((h,f,z))
+
+        for xo in self.xovers_from[helix_idx]:
+            ## h1,f1,z1,h2,f2,z2 = xo[3:]
+
+            append_if_in_range(*xo[:3])
+        for xo in self.xovers_to[helix_idx]:
+            append_if_in_range(*xo[3:])
+        return xos
+
     def _add_crossovers(self):
-        for h1,f1,z1,h2,f2,z2 in self.xover_list:
-            seg1, nt1 = self._get_segment_nucleotide(h1,z1,not f1)
-            seg2, nt2 = self._get_segment_nucleotide(h2,z2,f2)
-            ## TODO: use different types of crossovers
-            ## fwd?
-            ## 5'-to-3' direction
-            if isinstance(seg1, SingleStrandedSegment): f1 = True
-            if isinstance(seg2, SingleStrandedSegment): f2 = True
-            seg1.add_crossover(nt1,seg2,nt2,[f1,f2])
+        for hid,xos in self.xovers_from.items():
+            for h1,f1,z1,h2,f2,z2 in xos:
+                seg1, nt1 = self._get_segment_nucleotide(h1,z1,not f1)
+                seg2, nt2 = self._get_segment_nucleotide(h2,z2,f2)
+                ## TODO: use different types of crossovers
+                ## fwd?
+                ## 5'-to-3' direction
+                if isinstance(seg1, SingleStrandedSegment): f1 = True
+                if isinstance(seg2, SingleStrandedSegment): f2 = True
+                seg1.add_crossover(nt1,seg2,nt2,[f1,f2])
 
     def _add_prime_ends(self):
         for h,fwd,z in self._5prime_list: