Order of double-stranded segment creation and ds-segment direction should...

Order of double-stranded segment creation and ds-segment direction should match original cadnano reader.
Position and orientation of cadnano bases more consistent with original cadnano reader.

mrdna/readers/segmentmodel_from_cadnano.py

@@ -106,7 +106,7 @@ def get_helix_angle(part, helix_id, indices, fwd=True):
                                                                      '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
+    angle = eulerZ - twist_per_base*indices - 0.5*mgroove
     R = rotationAboutAxis(np.array((0,0,1)),angle)
     if fwd: R = R.dot(rotationAboutAxis((1,0,0),180))
     return R
@@ -142,7 +142,10 @@ def gen_id_series(strand, part, insertion_dict):
                 zids.insert(z_ind,str(insert_base)+"."+str(k))
                 z.insert(z_pos_ind,insert_base+k/(z_val+1))
     df["zid"] = zids
-    df["z"]   = -np.array(z)*3.4
+    df["z"]   = -3.4 * np.linspace( id_lo+0.25, id_hi-0.25, len(df) ) # Note: this spreads
+                                                                      # insertions uniformly through
+                                                                      # the strand, not through
+                                                                      # helical region like original did
 
     L=[(df["vh"][i],df["zid"][i],df["fwd"][i]) for i in df.index]
     if strand.isForward()==True:
@@ -226,8 +229,10 @@ def gen_prop_table(part):
     nt_prop["threeprime"]=tprime
 
     devlogger.info(f'  Finding orientations')
+    # nt_prop["orientation"]=[get_helix_angle(part, helix_id, int(float(indices)), _fwd)
+    #  for helix_id,indices,_fwd in zip(nt_prop["vh"],nt_prop["zid"],nt_prop['fwd'])]
     nt_prop["orientation"]=[get_helix_angle(part, helix_id, int(float(indices)), _fwd)
-                            for helix_id,indices,_fwd in zip(nt_prop["vh"],nt_prop["zid"],nt_prop['fwd'])]
+                            for helix_id,indices,_fwd in zip(nt_prop["vh"],nt_prop["z"]/3.4,nt_prop['fwd'])]
     nt_prop=nt_prop.fillna(-1)
     counter=-1
     devlogger.info(f'  Building bp and bp_map orientations')
@@ -262,9 +267,13 @@ def mrdna_model_from_cadnano(json_file,seq=None,**model_parameters):
     three_prime=np.array((list(nt_prop["threeprime"])))
     stack=np.array((list(nt_prop["stack"])))
     orientation=np.array(list(nt_prop["orientation"]))
+    rank = -np.array(list(nt_prop["z"])) # arrange segments in usual way
+    rank = rank + (np.max(rank)+1)*np.array(list(nt_prop["vh"]))
     devlogger.info(f'Building model from lists')
 
-    model = model_from_basepair_stack_3prime( r, bp, stack, three_prime, seq, orientation, explicitly_unstacked = np.argwhere(stack == -1), **model_parameters )
+    model = model_from_basepair_stack_3prime( r, bp, stack, three_prime, seq, orientation,
+                                              rank = rank, explicitly_unstacked = np.argwhere(stack == -1),
+                                              **model_parameters )
 
     devlogger.info(f'Fixing segment names')
     hmap = model._reader_list_hmap
@@ -285,7 +294,9 @@ def mrdna_model_from_cadnano(json_file,seq=None,**model_parameters):
         devlogger.debug(f'Renaming segs in {vh}')
         for i,(seg,s) in enumerate(sorted(ss, key=lambda x: x[-1])):
             devlogger.debug(f'  seg {i}: {seg} {s} {seg.end3.connection}')
-            seg.name = f'{vh:d}-{i:d}'
+            _newname = f'{vh:d}-{i:d}'
+            devlogger.debug(f'  {seg} -> {_newname}')
+            seg.name = _newname
             seg.segname = seg.name
 
     model._dataframe=nt_prop

mrdna/readers/segmentmodel_from_lists.py

@@ -78,11 +78,15 @@ def find_stacks(centers, transforms):
 
     return stacks_above
 
-def basepairs_and_stacks_to_helixmap(basepairs,stacks_above):
+def basepairs_and_stacks_to_helixmap(basepairs,stacks_above, nucleotide_rank=None):
+    """ Decide how helices should be created from basepair and stack data. Optionally use nucleotide_rank to determine direction of each helix and order of helices """
 
     helixmap = -np.ones(basepairs.shape, dtype=int)
     helixrank = -np.ones(basepairs.shape)
     is_fwd = np.ones(basepairs.shape, dtype=int)
+
+    if nucleotide_rank is None:
+        nucleotide_rank = np.arange(basepairs.shape, dtype=int)
     
     ## Remove stacks with nts lacking a basepairs
     nobp = np.where(basepairs < 0)[0]
@@ -93,7 +97,7 @@ def basepairs_and_stacks_to_helixmap(basepairs,stacks_above):
     end_ids = np.where( (stacks_above < 0)*(basepairs >= 0) )[0]
 
     hid = 0
-    for end in end_ids:
+    for end in sorted(end_ids, key=lambda x: nucleotide_rank[x]):
         if helixmap[end] >= 0:
             continue
         rank = 0
@@ -231,7 +235,7 @@ def set_splines(seg, coordinate, hid, hmap, hrank, fwd, basepair, orientation=No
 
 def model_from_basepair_stack_3prime(coordinate, basepair, stack, three_prime,
                                      sequence=None, orientation=None,
-                                     explicitly_unstacked = tuple(),
+                                     rank=None, explicitly_unstacked = tuple(),
                                      max_basepairs_per_bead = 5,
                                      max_nucleotides_per_bead = 5,
                                      local_twist = False,
@@ -325,7 +329,7 @@ def model_from_basepair_stack_3prime(coordinate, basepair, stack, three_prime,
     five_prime = _three_prime_list_to_five_prime(three_prime)
 
     """ Build map of dsDNA helices and strands """
-    hmap,hrank,fwd,intrahelical = basepairs_and_stacks_to_helixmap(bps,stack)
+    hmap,hrank,fwd,intrahelical = basepairs_and_stacks_to_helixmap(bps,stack,rank)
     double_stranded_helices = np.unique(hmap[hmap >= 0])    
     strands, strand_is_circular = _primes_list_to_strands(three_prime, five_prime)