Improve pdb and list readers

mrdna/RELEASE-VERSION

-1.0a.dev55
+1.0a.dev57

mrdna/readers/segmentmodel_from_lists.py

@@ -6,7 +6,7 @@ import os,sys
 import scipy
 
 from ..segmentmodel import SegmentModel, SingleStrandedSegment, DoubleStrandedSegment
-from ..arbdmodel.coords import quaternion_from_matrix
+from ..arbdmodel.coords import quaternion_from_matrix, rotationAboutAxis
 from .. import get_resource_path
 
 ref_stack_position = np.array((-2.41851735, -0.259761333, 3.39999978))
@@ -162,7 +162,7 @@ def basepairs_and_stacks_to_helixmap(basepairs,stacks_above):
     return helixmap, helixrank, is_fwd, intrahelical
 
 
-def set_splines(seg, coordinate, hid, hmap, hrank, fwd, orientation=None):
+def set_splines(seg, coordinate, hid, hmap, hrank, fwd, basepair, orientation=None):
     maxrank = np.max( hrank[hmap==hid] )
     if maxrank == 0:
         ids = np.where((hmap == hid))[0]
@@ -189,8 +189,20 @@ def set_splines(seg, coordinate, hid, hmap, hrank, fwd, orientation=None):
                 q = quaternion_from_matrix( orientation[ids[0]] )
                 if last_q is not None and last_q[1:].dot(q[1:]) < 0:
                     q = -1*q
+
+                ## Average quaterion with reverse direction
+                bp = basepair[ids[0]]
+                if bp >= 0:
+                    q2 = quaternion_from_matrix( orientation[bp].dot(rotationAboutAxis(np.array((1,0,0)),180)) )
+                    if last_q is None:
+                        if q[1:].dot(q2[1:]) < 0:
+                            q2 = -1*q2
+                    else:
+                        if last_q[1:].dot(q2[1:]) < 0:
+                            q2 = -1*q2
+                    q = (q + q2)*0.5
                 quats.append(q)
-                ## TODO: average quaterion with reverse direction
+                last_q = q
 
     coords = np.array(coords)
     seg.set_splines(contours,coords)
@@ -326,7 +338,7 @@ def model_from_basepair_stack_3prime(coordinate, basepair, stack, three_prime,
     for hid in double_stranded_helices:
         seg = DoubleStrandedSegment(name=str(hid),
                                     num_bp = np.sum(hmap==hid)//2)
-        set_splines(seg, coordinate, hid, hmap, hrank, fwd, orientation)
+        set_splines(seg, coordinate, hid, hmap, hrank, fwd, basepair, orientation)
 
         assert(hid == len(doubleSegments))
         doubleSegments.append(seg)
@@ -336,7 +348,7 @@ def model_from_basepair_stack_3prime(coordinate, basepair, stack, three_prime,
     for hid in single_stranded_helices:
         seg = SingleStrandedSegment(name=str(hid),
                                     num_nt = np.sum(hmap==hid))
-        set_splines(seg, coordinate, hid, hmap, hrank, fwd, orientation)
+        set_splines(seg, coordinate, hid, hmap, hrank, fwd, basepair, orientation)
 
         assert(hid == len(doubleSegments) + len(singleSegments))
         singleSegments.append(seg)
@@ -422,6 +434,16 @@ def model_from_basepair_stack_3prime(coordinate, basepair, stack, three_prime,
                           **model_parameters )
 
 
+    model._reader_list_coordinates = coordinate
+    model._reader_list_basepair = basepair
+    model._reader_list_stack = stack
+    model._reader_list_three_prime = three_prime
+    model._reader_list_sequence = sequence
+    model._reader_list_orientation = orientation
+    model._reader_list_hmap = hmap
+    model._reader_list_fwd = fwd
+    model._reader_list_hrank = hrank
+
     if sequence is None:
         for s in model.segments:
             s.randomize_unset_sequence()

mrdna/readers/segmentmodel_from_pdb.py

@@ -38,8 +38,7 @@ for key,val in refUnis.items():
     ref_res = val.residues[0]
     names = ref_res.atoms.select_atoms( stack_text ).atoms.names
     ref_name_ids_inv[key] = np.argsort(np.argsort(names))
-
-    
+    val.atoms.positions = val.atoms.positions.dot(CanonicalNucleotideFactory.DefaultOrientation) # return to default orientation
 
 ## TODO: possibly improve ref_stack_position
 # avg_stack_position = [ru.residues[0].atoms.select_atoms(stack_text).center_of_mass()
@@ -297,7 +296,10 @@ def SegmentModelFromPdb(*args,**kwargs):
     three_prime = -np.ones(bps.shape, dtype=np.int)
     for s in u.segments:
         r = s.atoms.select_atoms("name C1'").resindices
+        if u.segments[0] == s: print(r)
+        # three_prime[r[:-1]] = r[1:]
         three_prime[r[:-1]] = r[1:]
+        if u.segments[0] == s: print(three_prime[r])
 
         """ ## Remove connections between resids that are very far away
         dist = centers[r[1:]]-centers[r[:-1]]
@@ -321,17 +323,12 @@ def SegmentModelFromPdb(*args,**kwargs):
     seq = [resname_to_key[rn] 
            for rn in u.select_atoms("name C1'").resnames]
 
-    ## Get orientations
-    orientations = []
-    for R in transforms:
-        orientations.append( R.T.dot(CanonicalNucleotideFactory.DefaultOrientation) )
-
     return model_from_basepair_stack_3prime(centers,
                                             bps,
                                             stacks,
                                             three_prime,
                                             seq,
-                                            np.array(orientations)
+                                            np.array([t.T for t in transforms])
     )
 
 if __name__ == "__main__":