Working towards compatibility with cadnano structures; fixed some issues

CanonicalNucleotideAtoms.py

@@ -5,6 +5,7 @@ from  coords import rotationAboutAxis
 from arbdmodel import Group, PointParticle, ParticleType
 
 seqComplement = dict(A='T',G='C')
+resnames = dict(A="ADE", G="GUA", C="CYY", T="THY")
 for x,y in list(seqComplement.items()): seqComplement[y] = x
 
 
@@ -22,9 +23,18 @@ def stringToIntTuples(string, tupleLen, offset):
 class CanonicalNucleotide(Group):
     # DefaultOrientation = rotationAboutAxis([0,0,1], 68)
     # DefaultOrientation = rotationAboutAxis([0,0,1], 68+180)
-    DefaultOrientation = rotationAboutAxis([0,0,1], 68+180-20)
+    # DefaultOrientation = rotationAboutAxis([0,0,1], 68+180-20)
+    DefaultOrientation = rotationAboutAxis([0,0,1], 90)
+    # DefaultOrientation = rotationAboutAxis([0,0,1], -90)
+    # DefaultOrientation = rotationAboutAxis([0,0,1], 60)
+    # DefaultOrientation = rotationAboutAxis([0,0,1], -90)
+    DefaultOrientation = rotationAboutAxis([0,0,1], -30)
+    DefaultOrientation = rotationAboutAxis([0,0,1], -120)
+    DefaultOrientation = rotationAboutAxis([0,0,1], -150)
+    DefaultOrientation = rotationAboutAxis([0,0,1], 150)
     def __init__(self, prefix, seq):
-        self.sequence = seq
+        self.sequence = seq     # TODO: used?
+        self.resname = resnames[seq]
         Group.__init__(self, orientation = CanonicalNucleotide.DefaultOrientation)
                 
         ## idx, type,   
@@ -49,7 +59,7 @@ class CanonicalNucleotide(Group):
             p = PointParticle( type_ = atomic_types[tn],
                                position = [xs[i],ys[i],zs[i]],
                                name = names[i].decode('UTF-8') )
-            self.children.append( p )
+            self.add( p )
 
         
         atoms = self.children

arbdmodel.py

@@ -317,7 +317,7 @@ class ParticleType():
     
 
 class PointParticle(Transformable, Child):
-    def __init__(self, type_, position, name="A", segname="A", **kwargs):
+    def __init__(self, type_, position, name="A", **kwargs):
         parent = None
         if 'parent' in kwargs:
             parent = kwargs['parent']
@@ -326,7 +326,6 @@ class PointParticle(Transformable, Child):
 
         self.type_    = type_                
         self.idx     = None
-        self.segname = segname
         self.name = name
         self.counter = 0
 
@@ -416,19 +415,29 @@ class PdbModel(Transformable, Parent):
             for p in self.particles:
                 ## http://www.wwpdb.org/documentation/file-format-content/format33/sect9.html#ATOM
                 idx = p.idx+1
-                name = p.type_.name
-                resname = name[:3]
+                try:
+                    resname = p.resname
+                except:
+                    resname = p.name[:3]
+                try: 
+                    resid = p.resid
+                except:
+                    resid = idx
+
+                # name = str(p.name)[:4]
+                name = str(p.name)[:4]
+                resname = resname
                 chain = "A"
                 charge = 0
                 occ = 0
                 beta = 0
                 x,y,z = [x for x in p.collapsedPosition()]
                 
-                assert(idx < 1e5)
-                resid = "{:<4d}".format(idx)
+                assert(resid < 1e5)
+                resid = "{:<4d}".format(resid)
 
                 fh.write( formatString.format(
-                    idx, name[:1], "", resname, chain, resid, x, y, z, occ, beta, "", charge ))
+                    idx, name, "", resname, chain, resid, x, y, z, occ, beta, "", charge ))
         return
         
     def writePsf(self, filename):
@@ -448,13 +457,26 @@ class PdbModel(Transformable, Parent):
                            " {name:7s} {type:7s} {charge:f} {mass:f}\n"
             for p in self.particles:
                 idx = p.idx + 1
+                try:
+                    segname = p.segname
+                except:
+                    segname = "A"
+                try:
+                    resname = p.resname
+                except:
+                    resname = str(p.name)[:3]
+                try: 
+                    resid = p.resid
+                except:
+                    resid = idx
+
                 data = dict(
                     idx     = idx,
-                    segname = "A",
-                    resid   = "%d%c%c" % (idx," "," "), # TODO: work with large indices
-                    name    = p.type_.name[:1],
-                    resname = p.type_.name[:3],
-                    type    = p.type_.name[:1],
+                    segname = segname,
+                    resid   = "%d%c%c" % (resid," "," "), # TODO: work with large indices
+                    name    = str(p.name)[:4],
+                    resname = resname,
+                    type    = p.type_.name[:4],
                     charge  = p.charge,
                     mass    = p.mass
                 )

segmentmodel.py

@@ -18,17 +18,20 @@ import pdb
 TODO:
  + fix handling of crossovers for atomic representation
  + map to atomic representation
-    - add nicks
-    - handle cylcic dna
+    + add nicks
+    - handle circular dna
+    - transform ssDNA nucleotides 
     - shrink ssDNA
     - shrink dsDNA backbone
     - make orientation continuous
     - sequence
+ + ensure crossover bead potentials aren't applied twice 
  - remove performance bottlenecks
  - test for large systems
  - assign sequence
- - document
+ - ENM
  - rework location class Location 
+ - document
 """
 
 class Location():
@@ -45,6 +48,8 @@ class Location():
 
         self.prev_in_strand = None
         self.next_in_strand = None
+        
+        self.combine = None     # some locations might be combined in bead model 
 
     def get_connected_location(self):
         if self.connection is None:
@@ -135,6 +140,7 @@ class ConnectableElement():
 
 class SegmentParticle(PointParticle):
     def __init__(self, type_, position, name="A", segname="A", **kwargs):
+        self.name = name
         self.contour_position = None
         PointParticle.__init__(self, type_, position, name=name, segname=segname, **kwargs)
         self.intrahelical_neighbors = []
@@ -200,6 +206,7 @@ class Segment(ConnectableElement, Group):
         Group.__init__(self, name, children=[])
         ConnectableElement.__init__(self, connection_locations=[], connections=[])
 
+        self.resname = name
         self.start_orientation = None
         self.twist_per_nt = 0
 
@@ -263,6 +270,7 @@ class Segment(ConnectableElement, Group):
     
     def randomize_unset_sequence(self):
         bases = list(seqComplement.keys())
+        bases = ['T']        ## FOR DEBUG
         if self.sequence is None:
             self.sequence = [random.choice(bases) for i in range(self.num_nts)]
         else:
@@ -332,7 +340,6 @@ class Segment(ConnectableElement, Group):
         assert(c*(self.num_nts-1) == round(c*(self.num_nts-1)))
         # TODO? loc = self.Location( address=c, type_=type_, on_fwd_strand=is_fwd )
         loc = Location( self, address=c, type_=type_, on_fwd_strand=on_fwd_strand )
-        print("Appending location ",self.name, loc.type_,c, self.num_nts)
         self.locations.append(loc)
 
     ## TODO? Replace with abstract strand-based model?
@@ -343,8 +350,6 @@ class Segment(ConnectableElement, Group):
                      l.on_fwd_strand == on_fwd_strand] 
         if len(locations) == 0:
             self.add_location(nt,"5prime",on_fwd_strand)
-        else:
-            pdb.set_trace()
 
     def add_3prime(self, nt, on_fwd_strand=True):
         address = nt/(self.num_nts-1) # TODO: put into a fn call
@@ -448,16 +453,14 @@ class Segment(ConnectableElement, Group):
         # assert(b.parent is None)
         if b.parent is not None:
             b.parent.children.remove(b)
-        b.parent = self
-        self.children.append(b)
+        self.add(b)
         self.beads.append(b) # don't add orientation bead
         if "orientation_bead" in b.__dict__: # TODO: think of a cleaner approach
             o = b.orientation_bead
             o.contour_position = b.contour_position
             if o.parent is not None:
                 o.parent.children.remove(o)
-            o.parent = self
-            self.children.append(o)
+            self.add(o)
             self.add_bond(b,o, Segment.orientation_bond, exclude=True)
 
     def _rebuild_children(self, new_children):
@@ -490,7 +493,7 @@ class Segment(ConnectableElement, Group):
 
 
     def _generate_beads(self, bead_model, max_basepairs_per_bead, max_nucleotides_per_bead):
-        
+
         """ Generate beads (positions, types, etcl) and bonds, angles, dihedrals, exclusions """
         ## TODO: decide whether to remove bead_model argument
         ##       (currently unused)
@@ -689,13 +692,13 @@ class DoubleStrandedSegment(Segment):
             opos = pos + orientation.dot( np.array((Segment.orientation_bond.r0,0,0)) )
             o = SegmentParticle( Segment.orientation_particle, opos, nts,
                                  num_nts=nts, parent=self )
-            bead = SegmentParticle( Segment.dsDNA_particle, pos, nts,
+            bead = SegmentParticle( Segment.dsDNA_particle, pos, name="DNA",
                                     num_nts=nts, parent=self, 
                                     orientation_bead=o,
                                     contour_position=contour_position )
 
         else:
-            bead = SegmentParticle( Segment.dsDNA_particle, pos, nts,
+            bead = SegmentParticle( Segment.dsDNA_particle, pos, name="DNA",
                                     num_nts=nts, parent=self,
                                     contour_position=contour_position )
         self._add_bead(bead)
@@ -755,7 +758,8 @@ class SingleStrandedSegment(Segment):
 
     def _generate_one_bead(self, contour_position, nts):
         pos = self.contour_to_position(contour_position)
-        b = SegmentParticle( Segment.ssDNA_particle, pos, nts,
+        b = SegmentParticle( Segment.ssDNA_particle, pos, 
+                             name="NAS",
                              num_nts=nts, parent=self,
                              contour_position=contour_position )
         self._add_bead(b)
@@ -796,20 +800,28 @@ class StrandInSegment(Group):
 
 class Strand(Group):
     """ Class that holds atomic model, maps to segments """
-    
-    def __init__(self):
+    def __init__(self, segname = None):
         Group.__init__(self)
         self.num_nts = 0
         self.children = self.strand_segments = []
-        
+        self.segname = segname
+
     ## TODO disambiguate names of functions
     def add_dna(self, segment, start, end, is_fwd):
         """ start/end should be provided expressed as contour_length, is_fwd tuples """
         s = StrandInSegment( segment, start, end, is_fwd )
-        self.strand_segments.append( s )
+        self.add( s )
         self.num_nts += s.num_nts
 
-    def set_sequence(self):
+    def set_sequence(self,sequence):
+        ## validate input
+        assert( np.all( [i in ('A','T','C','G') for i in sequence] ) )
+
+        ## set sequence on each segment
+        for s in self.children:
+            seg = s.segment
+            ...
+
         ...
 
     # def get_sequence(self):
@@ -836,7 +848,7 @@ class Strand(Group):
                 else:
                     nt = seg._generate_atomic_nucleotide( c, s.is_fwd, "A" )
 
-                s.children.append(nt)
+                s.add(nt)
                 ## Join last basepairs
                 if last is not None:
                     o3,c3,c4,c2,h3 = [last.atoms_by_name[n] 
@@ -1058,7 +1070,34 @@ class SegmentModel(ArbdModel):
         segments = self.segments
         for s in segments:
             s.local_twist = local_twist
-            
+
+        """ Simplify connections """
+        d_nt = dict()           # 
+        for s in segments:
+            d_nt[s] = 1.5/(s.num_nts-1)
+        for s in segments:
+            ## replace consecutive crossovers with
+            cl = sorted( s.get_connections_and_locations("crossover"), key=lambda x: x[1].address )
+            last = None
+            for entry in cl:
+                c,A,B = entry
+                if last is not None and \
+                   (A.address - last[1].address) < d_nt[s]:
+                    same_type = c.type_ == last[0].type_
+                    same_dest_seg = B.container == last[2].container
+                    if same_type and same_dest_seg:
+                        if np.abs(B.address - last[2].address) < d_nt[B.container]:
+                            ## combine
+                            A.combine = last[1]
+                            B.combine = last[2]
+
+                            ...
+                # if last is not None:
+                #     s.bead_locations.append(last)
+                ...
+                last = entry
+        del d_nt
+
         """ Generate beads at intrahelical junctions """
         if self.DEBUG: print( "Adding intrahelical beads at junctions" )
 
@@ -1100,7 +1139,7 @@ class SegmentModel(ArbdModel):
             a1,a2 = [l.address   for l in (A,B)]
 
             b = s1.get_nearest_bead(a1)
-            if b is not None and np.abs(b.contour_position-a1)*s1.num_nts < 1:
+            if b is not None and np.abs(b.contour_position-a1)*s1.num_nts < 1.9:
                 ## combine beads
                 b.contour_position = 0.5*(b.contour_position + a1) # avg position
                 A.particle = b
@@ -1108,7 +1147,7 @@ class SegmentModel(ArbdModel):
                 A.particle = s1._generate_one_bead(a1,0)
 
             b = s2.get_nearest_bead(a2)
-            if b is not None and np.abs(b.contour_position-a2)*s2.num_nts < 1:
+            if b is not None and np.abs(b.contour_position-a2)*s2.num_nts < 1.9:
                 ## combine beads
                 b.contour_position = 0.5*(b.contour_position + a2) # avg position
                 B.particle = b
@@ -1294,9 +1333,15 @@ class SegmentModel(ArbdModel):
             pot = self.get_bond_potential(4,18.5)
             self.add_bond(b1,b2, pot)
 
+        crossover_bead_pots = set()
         for c,A,B in self.get_connections("crossover"):
-            ## TODO: avoid double-counting for double-crossovers
             b1,b2 = [loc.particle for loc in (c.A,c.B)]
+
+            ## Avoid double-counting
+            if (b1,b2) in crossover_bead_pots: continue
+            crossover_bead_pots.add((b1,b2))
+            crossover_bead_pots.add((b2,b1))
+                
             pot = self.get_bond_potential(4,18.5)
             self.add_bond(b1,b2, pot)
 
@@ -1307,13 +1352,13 @@ class SegmentModel(ArbdModel):
                 k = (1.0/2) * 1.5 * kT * (1.0 / (1-np.exp(-float(1)/147))) * 0.00030461742; # kcal_mol/degree^2
                 if 'orientation_bead' in b1.__dict__:
                     t0 = 90 + 60
-                    if A.type_: t0 -= 120
+                    if A.on_fwd_strand: t0 -= 120 # TODO handle antiparallel segments
                     o = b1.orientation_bead
                     pot = self.get_angle_potential(k,t0)
                     self.add_angle( o,b1,b2, pot )
                 else:
                     t0 = 90 + 60
-                    if B.type_: t0 -= 120
+                    if B.on_fwd_strand: t0 -= 120
                     o = b2.orientation_bead
                     pot = self.get_angle_potential(k,t0)
                     self.add_angle( b1,b2,o, pot )
@@ -1326,6 +1371,14 @@ class SegmentModel(ArbdModel):
 
                 k_fn = lambda sep: (1.0/2) * 1.5 * kT * (1.0 / (1-np.exp(-float(sep)/147))) * 0.00030461742; # kcal_mol/degree^2
                 t0 = 90
+                if A.on_fwd_strand: t0 -= 120
+                if B.on_fwd_strand: t0 += 120
+
+                # t0 = t0 % 360
+                t0 = ((t0+180) % 360)-180
+                assert(t0 < 180)
+                assert(t0 >= -180)
+
                 if 'orientation_bead' in b1.__dict__:
                     o1 = b1.orientation_bead
                     if u2 is not None:
@@ -1420,7 +1473,7 @@ class SegmentModel(ArbdModel):
 
         for seg in self.segments:
             locs = seg.get_5prime_locations()
-            print(locs)
+            # print(locs)
             if locs is None: continue
             # for pos, is_fwd in locs:
             for l in locs:
@@ -1428,11 +1481,19 @@ class SegmentModel(ArbdModel):
                 is_fwd = l.on_fwd_strand
                 s = Strand()
                 _recursively_build_strand(s, seg, pos, is_fwd)
-                print("{} {}".format(seg.name,s.num_nts))
+                # print("{} {}".format(seg.name,s.num_nts))
                 strands.append(s)
         
         self.children = self.strands = sorted(strands, key=lambda s:s.num_nts)[::-1] # or something
 
+        ## relabel segname
+        counter = 0
+        for s in self.strands:
+            if s.segname is None:
+                s.segname = "D%03d" % counter
+            counter += 1
+
+
     def _generate_atomic_model(self):
         for s in self.strands:
             s.generate_atomic_model()