find_stacks routine moved to segmentmodel_from_lists using only coordinates...

find_stacks routine moved to segmentmodel_from_lists using only coordinates and transforms; generate stacks automatically if None

kept old find_stacks routine in segmentmodel_from_pdb as find_stacks_mdanalysis

mrdna/readers/segmentmodel_from_lists.py

@@ -3,11 +3,14 @@
 import pdb
 import numpy as np
 import os,sys
+import scipy
 
 from ..segmentmodel import SegmentModel, SingleStrandedSegment, DoubleStrandedSegment
 from ..coords import quaternion_from_matrix
 from .. import get_resource_path
 
+ref_stack_position = np.array((-2.41851735, -0.259761333, 3.39999978))
+
 def _three_prime_list_to_five_prime(three_prime):
     five_prime = -np.ones(three_prime.shape, dtype=np.int)
     has_three_prime = np.where(three_prime >= 0)[0]
@@ -46,6 +49,31 @@ def _primes_list_to_strands(three_prime, five_prime):
 
     return strands, strand_is_circular
 
+def find_stacks(centers, transforms):
+
+    ## Find orientation and center of each nucleotide
+    expected_stack_positions = []
+    for R,c in zip(transforms,centers):
+        expected_stack_positions.append( c + ref_stack_position.dot(R) )
+
+    expected_stack_positions = np.array(expected_stack_positions, dtype=np.float32)
+
+    dists = scipy.spatial.distance_matrix(expected_stack_positions, centers)
+    dists = dists + 5*np.eye(len(dists))
+    idx1, idx2 = np.where(dists < 3.5)
+
+    ## Convert distances to stacks
+    stacks_above = -np.ones(len(centers), dtype=np.int)
+    _z = np.array((0,0,1))
+    for i in np.unique(idx1):
+        js = idx2[ idx1 == i ]
+        angles = [np.arccos( transforms[j].T.dot( transforms[i].dot(_z) ).dot( _z ) ) for j in js]
+        angles = np.array( angles )
+        tmp = np.argmin(dists[i][js] + 1.0*angles)
+        j = js[tmp]
+        stacks_above[i] = j
+
+    return stacks_above
 
 def basepairs_and_stacks_to_helixmap(basepairs,stacks_above):
 
@@ -224,16 +252,19 @@ def model_from_basepair_stack_3prime(coordinate, basepair, stack, three_prime,
         if orientation.shape != (num_nt,3,3):
             raise ValueError("The 'orientation' array is not properly formatted")
 
-    try:
-        stack = np.array(stack,dtype=np.int)
-    except:
-        raise TypeError("The 'stack' array could not be converted into a numpy integer array")
+    if stack is None:
+        stack = find_stacks(coordinate, orientation)
+    else:
+        try:
+            stack = np.array(stack,dtype=np.int)
+        except:
+            raise TypeError("The 'stack' array could not be converted into a numpy integer array")
 
-    if len(stack.shape) != 1:
-        raise ValueError("The 'stack' array has the wrong dimensionality")
+        if len(stack.shape) != 1:
+            raise ValueError("The 'stack' array has the wrong dimensionality")
 
-    if len(stack) != num_nt:
-        raise ValueError("The length of the 'stack' array does not match other inputs")
+        if len(stack) != num_nt:
+            raise ValueError("The length of the 'stack' array does not match other inputs")
 
     bps = basepair              # alias
 

mrdna/readers/segmentmodel_from_pdb.py

@@ -17,7 +17,7 @@ from ..segmentmodel import SegmentModel, SingleStrandedSegment, DoubleStrandedSe
 from .. import get_resource_path
 from .segmentmodel_from_lists import model_from_basepair_stack_3prime
 from .segmentmodel_from_lists import _three_prime_list_to_five_prime
-
+from .segmentmodel_from_lists import find_stacks
 
 resnames = dict(A='ADE DA DA5 DA3', T='THY DT DT5 DT3',
                 C='CYT DC DC5 DC3', G='GUA DG DG5 DG3')
@@ -178,7 +178,7 @@ def find_basepairs( u, centers, transforms, selection_text='all' ):
 
     return basepairs
 
-def find_stacks( u, centers, transforms, selection_text='all' ):
+def find_stacks_mdanalysis( u, centers, transforms, selection_text='all' ):
 
     ## Find orientation and center of each nucleotide
     expected_stack_positions = []
@@ -274,7 +274,8 @@ def SegmentModelFromPdb(*args,**kwargs):
     ## Find basepairs and base stacks
     centers,transforms = find_base_position_orientation(u)
     bps = find_basepairs(u, centers, transforms)
-    stacks = find_stacks(u, centers, transforms)
+    # stacks = find_stacks_mdanalysis(u, centers, transforms)
+    stacks = find_stacks(centers, transforms)
 
     ## Find three-prime ends
     three_prime = -np.ones(bps.shape, dtype=np.int)