diff --git a/mrdna/readers/__init__.py b/mrdna/readers/__init__.py
index 5e9a828d00e815cd50694180856df4f5aa083f73..31a6cedd5f7f92641d81ac6691280103a97c4e51 100644
--- a/mrdna/readers/__init__.py
+++ b/mrdna/readers/__init__.py
@@ -8,15 +8,11 @@ def read_cadnano(json_file, **model_parameters):
return model_from_cadnano_json(data, **model_parameters)
def read_vhelix(maya_file, **model_parameters):
- from .polygon_mesh import parse_maya_file, convert_maya_to_segments
+ from .polygon_mesh import parse_maya_file, convert_maya_bases_to_segment_model
from ..model.dna_sequence import m13 as m13seq
- from ..segmentmodel import SegmentModel
- data = parse_maya_file(maya_file)
- segments, dsSegmentDict = convert_maya_to_segments( data )
- if 'dimensions' not in model_parameters:
- model_parameters['dimensions']=(5000,5000,5000)
- model = SegmentModel( segments,**model_parameters )
+ maya_bases = parse_maya_file(maya_file)
+ model = convert_maya_bases_to_segment_model( maya_bases, **model_parameters )
model.set_sequence(m13seq*10)
return model
diff --git a/mrdna/readers/polygon_mesh.py b/mrdna/readers/polygon_mesh.py
index 700ad09fd34295430323d766254195bd4e6b4519..2eb0015f740581ac76548f6998334a68c5a20dc4 100644
--- a/mrdna/readers/polygon_mesh.py
+++ b/mrdna/readers/polygon_mesh.py
@@ -1,16 +1,42 @@
import numpy as np
import sys
import re
-from ..coords import rotationAboutAxis
+from ..coords import rotationAboutAxis, minimizeRmsd
+from ..version import maintainer
from ..segmentmodel import SegmentModel, SingleStrandedSegment, DoubleStrandedSegment
+DEBUG=False
+
+maya_obj_registry = dict()
+maya_obj_registry_short = dict()
class MayaObj():
+ """ Class representing a node in a Maya ascii file """
+
def __init__(self, maya_lines):
+ self.parent = None
+ self.children = dict()
+ self.position = np.array((0,0,0))
+ self.orientation = np.eye(3)
+
self.name, self.parent_name = MayaObj._parse_first_maya_line( maya_lines[0] )
+ self._parse_maya_lines(maya_lines)
+
+ if self.parent_name is not None:
+ if self.parent_name in maya_obj_registry:
+ p = maya_obj_registry[self.parent_name]
+ p.add(self)
+ elif self.parent_name in maya_obj_registry_short:
+ p = maya_obj_registry_short[self.parent_name]
+ p.add(self)
+ else:
+ raise Exception("Parent for object not yet defined")
+
+ maya_obj_registry[self.get_full_name()] = self
+ maya_obj_registry_short[self.name] = self
def _parse_first_maya_line(l):
- ## createNode vHelix -n "helix_1";
+ """ Extracts name and parent from maya file """
name = parent = None
vals = l.split()
for i in range(len(vals)):
@@ -19,10 +45,59 @@ class MayaObj():
if vals[i] == "-p":
parent = MayaObj._strip_maya_string( vals[i+1] )
return name,parent
-
+
+ def _parse_maya_lines(self, maya_lines):
+ """ Sets position and orientation from attributes in vHelix maya files """
+ for l in maya_lines:
+ vals = l.split()
+ if vals[0] == "setAttr":
+ if vals[1] == '".t"':
+ self.position = 10*np.array([float(a) for a in vals[4:7]])
+ elif vals[1] in ('".r"','".rsrr"'):
+ ax,ay,az = [float(a) for a in vals[4:7]]
+ Rx = rotationAboutAxis( [1,0,0], ax, normalizeAxis=False )
+ Ry = rotationAboutAxis( [0,1,0], ay, normalizeAxis=False )
+ Rz = rotationAboutAxis( [0,0,1], az, normalizeAxis=False )
+ self.orientation = Rz.dot(Ry.dot(Rx))
+
def _strip_maya_string(string):
- return re.match('"([a-zA-Z_[0-9]+)";?$', string).group(1)
+ return re.match('"\|?([a-zA-Z_0-9|]+)";?$', string).group(1)
+ def add(self, obj):
+ self.children[obj.name] = obj
+ obj.parent = self
+
+ def get_position(self):
+ if self.parent is not None:
+ return self.parent._transform_child_position(self.position)
+ else:
+ return self.position
+
+ def get_orientation(self):
+ if self.parent is not None:
+ return self.parent._transform_child_orientation(self.orientation)
+ else:
+ return self.orientation
+
+ def _transform_child_position(self, position):
+ if self.parent is not None:
+ p = self.parent._transform_child_position( position )
+ else:
+ p = position
+ return self.orientation.dot(p)+self.position
+
+ def _transform_child_orientation(self, orientation):
+ if self.parent is not None:
+ o = self.parent._transform_child_orientation( orientation )
+ else:
+ o = orientation
+ return self.orientation.dot(o)
+
+ def get_full_name(self):
+ ret = self.name
+ if self.parent is not None:
+ ret = '{}|{}'.format(self.parent.get_full_name(), ret)
+ return ret
class MayaConnection():
def __init__(self, helix1, base1, suff1, helix2, base2, suff2):
@@ -33,8 +108,6 @@ class MayaConnection():
self.base2 = base2
self.suff2 = suff2
-# MayaConnections = dict()
-
def ParseMayaConnection(line, base_dict):
words = line.split()
if len(words) != 3:
@@ -63,79 +136,14 @@ def ParseMayaConnection(line, base_dict):
return MayaConnection(helix1,base1,suff1,helix2,base2,suff2)
return None
-
-class MayaHelix(MayaObj):
- def __init__(self, maya_lines):
- MayaObj.__init__(self, maya_lines)
- # self.bases = []
- self.bases = dict()
- self.position = None
- self.orientation = None
-
- self._parse_maya_lines(maya_lines)
-
- def _parse_maya_lines(self, maya_lines):
- """
-createNode vHelix -n "helix_30";
-setAttr ".t" -type "double3" 9.24809 21.4638 -12.465 ;
-setAttr ".r" -type "double3" 28.99971393816978 2.2177410772603801 87.637950315471372 ;
-setAttr ".dh" yes;
- """
- for l in maya_lines:
- vals = l.split()
- if vals[0] == "setAttr":
- if vals[1] == '".t"':
- self.position = 10*np.array([float(a) for a in vals[4:7]])
- elif vals[1] == '".r"':
- # self.eulerAngles = np.array([float(a) for a in vals[4:7]]) # defined here
- ax,ay,az = [float(a) for a in vals[4:7]]
- Rx = rotationAboutAxis( [1,0,0], ax, normalizeAxis=False )
- Ry = rotationAboutAxis( [0,1,0], ay, normalizeAxis=False )
- Rz = rotationAboutAxis( [0,0,1], az, normalizeAxis=False )
- self.orientation = Rz.dot(Ry.dot(Rx)) # TODO: check
- # rotXYZ = np.array()
-
- def add_base(self, base):
- # self.bases.append(base)
- self.bases[base.name] = base
- base.parent = self
-
- def find_double_helix_regions(self):
- double_helices = []
- for bn,b1 in self.bases.items():
- b2 = b1.basepair
- if b2 is None: continue
- d = np.linalg.norm( b1.get_position()-b2.get_position() )
- print( d )
-
-
- def write_xyz(self,file_handle = sys.stdout):
- for bn,b in self.bases.items():
- position = b.get_position()
- # position = b.position.dot( self.orientation ) + self.position # TODO: double check
- file_handle.write("%s %f %f %f\n" % (b.name[0], position[0], position[1], position[2]))
-
-
class MayaBase(MayaObj):
def __init__(self, maya_lines):
MayaObj.__init__(self, maya_lines)
self._parse_maya_lines(maya_lines)
self.basepair = None
- self.parent = None
self.end3 = None
self.end5 = None
- def _parse_maya_lines(self, maya_lines):
- """
-createNode HelixBase -n "backw_1" -p "helix_30";
-setAttr ".t" -type "double3" 0.36767788771440857 -0.78848777472188547 -6.014 ;
- """
- for l in maya_lines:
- vals = l.split()
- if vals[0] == "setAttr":
- if vals[1] == '".t"':
- self.position = 10*np.array([float(a) for a in vals[4:7]])
-
def add_basepair(self, base):
assert( self.parent == base.parent )
assert( self.basepair is None or self.basepair == base )
@@ -151,27 +159,22 @@ setAttr ".t" -type "double3" 0.36767788771440857 -0.78848777472188547 -6.014 ;
def add_end5(self, base):
base.add_end3(self)
-
- def get_position(self):
- h = self.parent
- return h.orientation.dot(self.position) + h.position # TODO: double check
-
-
-# class MayaAimConstraint(MayaObj):
-# ...
-
-
-# set ch [open /home/cmaffeo2/Downloads/ball.ma]
+ def __repr__(self):
+ return "Base {} {:.2f}".format(self.name, self.position[2])
def parse_maya_file(maya_file):
+ """ Function to parse vHelix maya ascii file, extracting useful information about base """
+
+ objects = []
+ objects_by_full_name = dict()
helices = dict()
bases = []
- aimConstraints = []
+ aim_constraints = []
connections = []
- ## Parse through .ma file
+ """ Parse .ma file """
with open(maya_file) as fh:
linesBuffer = []
for line in fh:
@@ -182,12 +185,15 @@ def parse_maya_file(maya_file):
if len(linesBuffer) > 0:
objType = linesBuffer[0].split()[1]
if objType == "vHelix":
- h = MayaHelix( linesBuffer )
- helices[h.name] = h
+ h = MayaObj( linesBuffer )
+ helices[h.get_full_name()] = h
elif objType == "HelixBase":
bases.append( MayaBase( linesBuffer ) )
- # elif objType == "aimConstraint":
- # aimConstraints.append( MayaAimConstraint( linesBuffer ) )
+ elif objType == "transform":
+ o = MayaObj( linesBuffer )
+ objects.append(o)
+ elif objType == "aimConstraint":
+ aim_constraints.append( MayaObj( linesBuffer ) )
## Clear lines buffer
linesBuffer = []
elif words[0] == "connectAttr":
@@ -195,13 +201,8 @@ def parse_maya_file(maya_file):
## Extend lines buffer
linesBuffer.append(line)
-
- ## Crawl through bases, attach to parent helix
- for b in bases:
- helix_name = b.parent_name
- helices[helix_name].add_base( b ) # function updates base position?
- ## Parse connections
+ """ Parse connections """
new_connections = []
base_dict = { b.name:b for b in bases }
for line in connections:
@@ -210,7 +211,7 @@ def parse_maya_file(maya_file):
new_connections.append(conn)
connections = new_connections
- ## Assign base connectivity
+ """ Assign base connectivity """
for c in connections:
h1 = c.helix1
bn1 = c.base1
@@ -219,349 +220,166 @@ def parse_maya_file(maya_file):
bn2 = c.base2
s2 = c.suff2
- b1 = helices[h1].bases[bn1]
- b2 = helices[h2].bases[bn2]
+ b1 = helices[h1].children[bn1]
+ b2 = helices[h2].children[bn2]
if s1 == "lb" and s2 == "lb":
b1.add_basepair(b2)
- # elif (bn1[0] == "f" and s1 == "fw" and s2 == "bw") or \
- # (bn1[0] == "b" and s1 == "bw" and s2 == "fw"):
- # b1.add_end3(b2)
- # elif (bn1[0] == "f" and s1 == "bw" and s2 == "fw") or \
- # (bn1[0] == "b" and s1 == "fw" and s2 == "bw"):
- # b1.add_end5(b2)
elif s1 == "fw" and s2 == "bw":
- b1.add_end3(b2)
- elif s1 == "bw" and s2 == "fw":
b1.add_end5(b2)
+ elif s1 == "bw" and s2 == "fw":
+ b1.add_end3(b2)
else:
- raise Exception("Unkown connection %s %s" % (s1,s2))
+ raise Exception("Unrecognized connection %s %s" % (s1,s2))
- # # base_positions = np.zeros( [ len(bases), 3 ] )
- # base_positions = np.array( [b.position for b in bases] )
- # for aim in aimConstraints:
- # base = base_dict[ aim.parent_name ]
+ """ Find local basis of each base so that
+ model_from_basepair_stack_3prime can determine stacks """
- return helices
+ ref_below_position = np.array((5.11, -0.4, -3.34))
+ ref_stack_position = np.array((-1.5, -4.7, 3.34))
+ ref_bp_position = np.array((0.042, -16.98, 0.0))
+ ref_bp_below_position = np.array((4.2, -12.9, 3.34))
+ ref_bp_stack_position = np.array((-4.42, -14.46, -3.34))
-def convert_maya_to_segments(helices):
+ ref_positions = np.array((ref_below_position,
+ ref_stack_position,
+ ref_bp_position,
+ ref_bp_below_position,
+ ref_bp_stack_position))
- def get_end3_ssDNA(nt):
- ret = []
- if nt.end3 is not None:
- # print( nt.parent.name, nt.end3.parent.name, nt.end3.name )
- nt2 = nt.end3
- if nt2.basepair is None:
- ret.append(nt2)
- ret.extend(get_end3_ssDNA(nt2))
- return ret
+ shift = np.array((1.19, -3.52, 0.08))
- def get_end5_ssDNA(nt):
- ret = []
- if nt.end5 is not None:
- nt2 = nt.end5
- if nt2.basepair is None:
- ret.append(nt2)
- ret.extend(get_end5_ssDNA(nt2))
- return ret
+ if DEBUG:
+ write_pdb_psf(bases, 'before')
- """ Build dsDNA helices """
- segments = dict()
- segmentsBot = dict()
- segmentsTop = dict()
- for hn,h in helices.items():
- paired_bases = [b for bn,b in h.bases.items() if b.basepair is not None]
- upaired = [b for bn,b in h.bases.items() if b.basepair is None]
- # for b in upaired:
- # print(hn, b.name, b.end3, b.end5)
+ def update_geometry(b,beads):
- ## Find helical axis
- # if len(paired_bases) < 6:
- # raise Exception("Helices shorter than 3 bp are not supported")
- # coords = np.array([b.position for b in paired_bases])
- # u,s,vh = np.linalg.svd(coords)
- # axis = vh[0] / np.linalg.norm(vh[0])
-
- fwd_bases = list(filter(lambda b:b.name[0] == 'f', paired_bases))
- # rank = [np.dot(b.position,axis) for b in fwd_bases]
- rank = [b.position[2] for b in fwd_bases]
- order = np.argsort(rank)
- ordered_fwd_bases = [fwd_bases[o] for o in order]
- next_fwd_is_3end = [b1.end3==b2 for b1,b2 in
- zip(ordered_fwd_bases[:-1],ordered_fwd_bases[1:])]
-
- ## Check if axis points the wrong way
- assert( np.mean(next_fwd_is_3end) > 0.5 )
- # if np.mean( next_fwd_is_3end ) < 0.5:
- # # axis = axis * -1
- # rank = rank[order[-1]]-rank
- # order = reversed(order)
- # ordered_fwd_bases = reversed(ordered_fwd_bases)
-
- top = ordered_fwd_bases[-1]
- bot = ordered_fwd_bases[0].basepair
-
- assert( bot.name[0] == 'b' )
- assert( top.name[0] == 'f' )
-
- botPos = 0.5 * (bot.get_position() + bot.basepair.get_position())
- topPos = 0.5 * (top.get_position() + top.basepair.get_position())
-
- segname = hn
- if segname[:6] == "helix_":
- segname = 'D'+segname[6:]
- segments[hn] = DoubleStrandedSegment(segname, num_bp = len(fwd_bases),
- start_position = botPos,
- end_position = topPos)
-
- ## Add 3prime/5prime ends
- nt = 0
- def add_primes(b,nt,fwd):
- if b.end3 is None:
- segments[hn].add_3prime(nt,fwd)
- if b.end5 is None:
- segments[hn].add_5prime(nt,fwd)
-
- for b in ordered_fwd_bases:
- add_primes(b,nt,True)
- add_primes(b.basepair,nt,False)
- nt+=1
-
- segmentsBot[hn] = bot
- segmentsTop[hn] = top
- assert( bot.end3 not in paired_bases )
- assert( top.end3 not in paired_bases )
-
- """ Sanity check """
- def end_crosses(end,direction):
- assert(direction in ("end3","end5"))
- other = end.__dict__[direction]
- if other is None or end.parent != other.parent:
- return 1
- else:
- return 0
-
- hs = segments.keys()
- botCrossings1 = [end_crosses(segmentsBot[h],"end3") for h in hs]
- botCrossings2 = [end_crosses(segmentsBot[h].basepair,"end5") for h in hs]
- topCrossings1 = [end_crosses(segmentsTop[h],"end3") for h in hs]
- topCrossings2 = [end_crosses(segmentsTop[h].basepair,"end5") for h in hs]
- # print( [np.mean(a) for a in (botCrossings1,botCrossings2,topCrossings1,topCrossings2)] )
-
- # for a in (botCrossings1,botCrossings2,topCrossings1,topCrossings2:
- # print(np.mean(a))
- # assert( np.mean(a) > 0.5 )
-
- """ Build ssDNA segments """
- ss_segments = []
- connectionToSingleStrand = dict()
- for hn,seg in segments.items():
- bot = segmentsBot[hn]
- top = segmentsTop[hn]
- # if top.end3 is not None:
- # print(top.parent.name,top.name, top.end3.parent.name, top.end3.name)
- # if bot.end3 is not None:
- # print(bot.parent.name,bot.name, bot.end3.parent.name, bot.end3.name)
-
- ## TODO: clean this up
- def get_end_positions(ssSeg, is_fwd):
- start_idx, end_idx = (0,-1)
- if not is_fwd:
- start_idx, end_idx = (-1,0)
-
- r0 = ssSeg[0].get_position()
- if ssSeg[0].end5 is not None:
- r0 = 0.75*r0 + 0.25*ssSeg[0].end5.get_position()
- r1 = ssSeg[-1].get_position()
- if ssSeg[-1].end3 is not None:
- r1 = 0.75*r1 + 0.25*ssSeg[-1].end3.get_position()
- return r0,r1
-
- if top not in connectionToSingleStrand:
- ssSeg = get_end3_ssDNA(top)
- num_nt = len(ssSeg)
- if num_nt > 0:
- r0,r1 = get_end_positions(ssSeg,is_fwd=True)
- seg2 = SingleStrandedSegment("S{:03d}".format(len(ss_segments)),
- start_position = r0,
- end_position = r1,
- num_nt = num_nt)
- seg.connect_end3(seg2)
- connectionToSingleStrand[top] = seg2
- other = ssSeg[-1].end3
- if other is not None:
- hn3 = other.parent.name
- seg3 = segments[hn3]
- bot3 = segmentsBot[hn3]
- top3 = segmentsTop[hn3]
- # if bot3 == other or bot3.basepair == other:
- # seg3.connect_start5(seg2)
- # elif top3 == other or top3.basepair == other:
- # seg3.connect_end5(seg2)
- if bot3.basepair == other:
- seg3.connect_start5(seg2)
- elif top3.basepair == other:
- seg3.connect_end5(seg2)
- else:
- raise Exception("BUG")
- connectionToSingleStrand[other] = seg2
- ss_segments.append(seg2)
-
- if top.basepair not in connectionToSingleStrand:
- ssSeg = get_end5_ssDNA(top.basepair)[::-1]
- if len(ssSeg) > 0:
- r0,r1 = get_end_positions(ssSeg,is_fwd=False)
- seg2 = SingleStrandedSegment("S{:03d}".format(len(ss_segments)),
- start_position = r0,
- end_position = r1,
- num_nt = len(ssSeg))
- seg.connect_end5(seg2)
- connectionToSingleStrand[top.basepair] = seg2
- other = ssSeg[0].end5
- if other is not None:
- hn3 = other.parent.name
- seg3 = segments[hn3]
- bot3 = segmentsBot[hn3]
- top3 = segmentsTop[hn3]
- # if bot3 == other or bot3.basepair == other:
- # seg3.connect_start5(seg2)
- # elif top3 == other or top3.basepair == other:
- # seg3.connect_end5(seg2)
- if bot3 == other:
- seg3.connect_start3(seg2)
- elif top3 == other:
- seg3.connect_end3(seg2)
- else:
- raise Exception("BUG")
- connectionToSingleStrand[other] = seg2
- ss_segments.append(seg2)
-
- if bot not in connectionToSingleStrand:
- ssSeg = get_end3_ssDNA(bot)
- if len(ssSeg) > 0:
- r0,r1 = get_end_positions(ssSeg,is_fwd=True)
- seg2 = SingleStrandedSegment("S{:03d}".format(len(ss_segments)),
- start_position = r0,
- end_position = r1,
- num_nt = len(ssSeg))
- seg.connect_start3(seg2)
- connectionToSingleStrand[bot] = seg2
- other = ssSeg[-1].end3
- if other is not None:
- hn3 = other.parent.name
- seg3 = segments[hn3]
- bot3 = segmentsBot[hn3]
- top3 = segmentsTop[hn3]
- if bot3.basepair == other:
- seg3.connect_start5(seg2)
- elif top3.basepair == other:
- seg3.connect_end5(seg2)
- else:
- raise Exception("BUG")
- connectionToSingleStrand[other] = seg2
- ss_segments.append(seg2)
-
- if bot.basepair not in connectionToSingleStrand:
- ssSeg = get_end5_ssDNA(bot.basepair)[::-1]
- if len(ssSeg) > 0:
- r0,r1 = get_end_positions(ssSeg,is_fwd=False)
- seg2 = SingleStrandedSegment("S{:03d}".format(len(ss_segments)),
- start_position = r0,
- end_position = r1,
- num_nt = len(ssSeg))
- seg.connect_start5(seg2)
- connectionToSingleStrand[bot.basepair] = seg2
- other = ssSeg[0].end5
- if other is not None:
- hn3 = other.parent.name
- seg3 = segments[hn3]
- bot3 = segmentsBot[hn3]
- top3 = segmentsTop[hn3]
- if bot3 == other:
- seg3.connect_start3(seg2)
- elif top3 == other:
- seg3.connect_end3(seg2)
- else:
- raise Exception("BUG")
- connectionToSingleStrand[other] = seg2
- ss_segments.append(seg2)
-
- ## Connect remaining dsDNA ends
- no_connection_count = 0
- for hn,seg in segments.items():
- bot = segmentsBot[hn]
- top = segmentsTop[hn]
-
- top_bp, bot_bp = [b.basepair for b in (top,bot)]
+ """ Update the position and orientation of bead b using the 5
+ neighboring bases in "beads" as a reference """
- if top not in connectionToSingleStrand and top.end3 is not None:
- b2 = top.end3
- h2 = b2.parent.name
- s2 = segments[h2]
- if b2 == segmentsBot[h2] or b2.basepair == segmentsBot[h2]:
- seg.connect_end3(s2.start5,type_="terminal_crossover")
- elif b2 == segmentsTop[h2] or b2.basepair == segmentsTop[h2]:
- seg.connect_end3(s2.end5,type_="terminal_crossover")
- else:
- raise Exception("BUG")
+ def _get_beads_index(i):
+ j=0
+ while j <= i:
+ if beads[j] is None:
+ i+=1
+ j+=1
+ return j-1
+
+ """ Get position and filter out beads == None """
+ r0 = b.get_position()
+ tmp_ref_positions = ref_positions[[b2 is not None for b2 in beads]]
+ positions = np.array([b2.get_position()-r0 for b2 in beads if b2 is not None])
+ dist = np.linalg.norm(positions,axis=-1)
+
+ """ Discard neighboring bases that are very far from current base """
+ if np.any(dist > 20) and len(positions) > 3:
+ i = _get_beads_index(np.argmax(dist))
+ beads[i] = None
+ return update_geometry(b, beads)
+
+ """ Find transform that fits neighborhood around base to a reference """
+ R,comB,comA = minimizeRmsd(positions, tmp_ref_positions)
+
+ dr = np.array([R.T.dot(p-comB)-p0+comA for p,p0 in zip(positions, tmp_ref_positions)])
+ dr2 = (dr**2).sum(axis=-1)
+ msd = dr2.sum()/len(positions)
+
+ """ Discard neighboring bases that are very far from their expected position """
+ if msd > 17 and len(positions) > 3:
+ # print("WARNING:",b,msd)
+ # print(dr2)
+ i = _get_beads_index(np.argmax(dr2))
+ beads[i] = None
+ return update_geometry(b,beads)
- if top_bp not in connectionToSingleStrand and top_bp.end5 is not None:
- b2 = top_bp.end5
- h2 = b2.parent.name
- s2 = segments[h2]
- if b2 == segmentsBot[h2] or b2.basepair == segmentsBot[h2]:
- seg.connect_end5(s2.start3,type_="terminal_crossover")
- elif b2 == segmentsTop[h2] or b2.basepair == segmentsTop[h2]:
- seg.connect_end5(s2.end3, type_="terminal_crossover")
- else:
- raise Exception("BUG")
-
- if bot not in connectionToSingleStrand and bot.end3 is not None:
- b2 = bot.end3
- h2 = b2.parent.name
- s2 = segments[h2]
- if b2 == segmentsBot[h2] or b2.basepair == segmentsBot[h2]:
- seg.connect_start3(s2.start5,type_="terminal_crossover")
- elif b2 == segmentsTop[h2] or b2.basepair == segmentsTop[h2]:
- seg.connect_start3(s2.end5,type_="terminal_crossover")
- else:
- raise Exception("BUG")
-
- if bot_bp not in connectionToSingleStrand and bot_bp.end5 is not None:
- b2 = bot_bp.end5
- h2 = b2.parent.name
- s2 = segments[h2]
- if b2 == segmentsBot[h2] or b2.basepair == segmentsBot[h2]:
- seg.connect_start5(s2.start3,type_="terminal_crossover")
- elif b2 == segmentsTop[h2] or b2.basepair == segmentsTop[h2]:
- seg.connect_start5(s2.end3, type_="terminal_crossover")
- else:
- raise Exception("BUG")
+ """ Assign position and orientation to base """
+ if b.parent is None:
+ o_parent = np.eye(3)
+ else:
+ o_parent = b.parent.get_orientation()
+
+ b.orientation = o_parent.T.dot( R )
+ b.position = b.position + b.orientation.dot( shift )
- return [s for n,s in segments.items()] + ss_segments, segments
+ """ Loop over bases updating their positions and orientations """
+ for b in bases:
+ if b.basepair is not None:
+ beads = [b.end5, b.end3, b.basepair, b.basepair.end5, b.basepair.end3]
+ update_geometry(b,beads)
+
+ if DEBUG:
+ write_pdb_psf(bases, 'after')
+
+ return bases
+
+
+def write_pdb_psf(bases, prefix):
+
+ """ Function for debugging .ma parser """
+
+ from ..model.arbdmodel import ArbdModel,ParticleType, PointParticle, Group
+ from ..model.nonbonded import HarmonicBond
+
+ types = {c:ParticleType(c.upper(), mass=1,radius=1) for c in "bfuxyzs"}
+ bond=HarmonicBond(k=0,r0=1)
+
+ base_to_particle = dict()
+ seg = []
+ count=0
+ for b in bases:
+ r = b.get_position()
+ o = b.get_orientation()
+
+ type_character = 'f' if re.match('.*forw.*',b.name) else 'b' if re.match('.*bac.*',b.name) else 'u'
+ name = type_character + b.name.split("_")[-1]
+ p = PointParticle( name=name, type_=types[type_character], position = r )
+ particles = [p]
+ base_to_particle[b] = p
+ for axis,c in zip(np.eye(3),"xyz"):
+ particles.append( PointParticle( name=c, type_=types[c],
+ position = r + o.T.dot(axis) ))
+
+ ref_stack_position = np.array((-2.41851735, -0.259761333, 3.39999978))
+ particles.append( PointParticle( name=c, type_=types[c],
+ position = r + o.dot(ref_stack_position) ))
+
+ res = Group(name='r'+str(count), children=particles, resid=count)
+ count+=1
+ for p2 in particles[1:]:
+ res.add_bond(p,p2, bond=bond)
+ seg.append(res)
+
+ model = ArbdModel(seg)
+ for b in bases:
+ if b.end3 is not None:
+ model.add_bond( base_to_particle[b], base_to_particle[b.end3], bond=bond )
+ if b.basepair is not None:
+ model.add_bond( base_to_particle[b], base_to_particle[b.basepair], bond=bond )
+ model._countParticleTypes()
+ model._updateParticleOrder()
+ model.writePsf(prefix+'.psf')
+ model.writePdb(prefix+'.pdb')
-def write_xyz(filename, helices):
- nbeads = 0
- for hn,h in helices.items():
- nbeads += len(h.bases)
- with open("test.xyz","w") as fh:
- fh.write("%d\n\n" % nbeads)
- for hn,h in helices.items():
- h.write_xyz(fh)
- # h.find_double_stranded_regions
+def convert_maya_bases_to_segment_model(maya_bases, **model_parameters):
+ """ Converts bases parse from .ma file into mrdna an model """
-if __name__ == "__main__":
- maya_helices = parse_maya_file("/home/cmaffeo2/Downloads/bunny.ma")
+ from .segmentmodel_from_lists import model_from_basepair_stack_3prime
- segments,dsSegmentDict = convert_maya_to_segments(maya_helices)
-
- model = SegmentModel( segments,
- local_twist = True,
- max_basepairs_per_bead = 1,
- max_nucleotides_per_bead = 1,
- dimensions=(5000,5000,5000), DEBUG=True )
+ base_to_index = {b:i for i,b in zip(range(len(maya_bases)),maya_bases)}
+ coord = np.zeros((len(maya_bases),3))
+ orientation = np.zeros((len(maya_bases),3,3))
+ bp,end3 = [-np.ones((len(maya_bases)),dtype=np.int) for i in range(2)]
- model.simulate( outputPrefix = 'bunny2', outputPeriod=1e4, numSteps=1e8, gpu=1 )
+ for i,b in zip(range(len(maya_bases)),maya_bases):
+ coord[i] = b.get_position()
+ orientation[i] = b.get_orientation().T
+ if b.basepair is not None:
+ bp[i] = base_to_index[ b.basepair ]
+ if b.end3 is not None:
+ end3[i] = base_to_index[ b.end3 ]
+ return model_from_basepair_stack_3prime(coord, bp, None, end3,
+ orientation=orientation,
+ **model_parameters)