Skip to content
Snippets Groups Projects
Commit 74d1ddf7 authored by cmaffeo2's avatar cmaffeo2
Browse files

Fixed polygon_mesh reader so that it works with user-generated .ma files, not... 

Fixed polygon_mesh reader so that it works with user-generated .ma files, not just autogenerated ones
parent 7b417c83
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
mrdna/readers/__init__.py
+ 3
− 7
View file @ 74d1ddf7
...@@ -8,15 +8,11 @@ def read_cadnano(json_file, **model_parameters): ...@@ -8,15 +8,11 @@ def read_cadnano(json_file, **model_parameters):
return model_from_cadnano_json(data, **model_parameters) return model_from_cadnano_json(data, **model_parameters)
def read_vhelix(maya_file, **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 ..model.dna_sequence import m13 as m13seq
from ..segmentmodel import SegmentModel
data = parse_maya_file(maya_file) maya_bases = parse_maya_file(maya_file)
segments, dsSegmentDict = convert_maya_to_segments( data ) model = convert_maya_bases_to_segment_model( maya_bases, **model_parameters )
if 'dimensions' not in model_parameters:
model_parameters['dimensions']=(5000,5000,5000)
model = SegmentModel( segments,**model_parameters )
model.set_sequence(m13seq*10) model.set_sequence(m13seq*10)
return model return model
......
mrdna/readers/polygon_mesh.py
+ 236
− 418
View file @ 74d1ddf7
import numpy as np import numpy as np
import sys import sys
import re import re
from ..coords import rotationAboutAxis from ..coords import rotationAboutAxis, minimizeRmsd
from ..version import maintainer
from ..segmentmodel import SegmentModel, SingleStrandedSegment, DoubleStrandedSegment from ..segmentmodel import SegmentModel, SingleStrandedSegment, DoubleStrandedSegment
DEBUG=False
maya_obj_registry = dict()
maya_obj_registry_short = dict()
class MayaObj(): class MayaObj():
""" Class representing a node in a Maya ascii file """
def __init__(self, maya_lines): 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.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): def _parse_first_maya_line(l):
## createNode vHelix -n "helix_1"; """ Extracts name and parent from maya file """
name = parent = None name = parent = None
vals = l.split() vals = l.split()
for i in range(len(vals)): for i in range(len(vals)):
...@@ -19,10 +45,59 @@ class MayaObj(): ...@@ -19,10 +45,59 @@ class MayaObj():
if vals[i] == "-p": if vals[i] == "-p":
parent = MayaObj._strip_maya_string( vals[i+1] ) parent = MayaObj._strip_maya_string( vals[i+1] )
return name,parent 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): 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(): class MayaConnection():
def __init__(self, helix1, base1, suff1, helix2, base2, suff2): def __init__(self, helix1, base1, suff1, helix2, base2, suff2):
...@@ -33,8 +108,6 @@ class MayaConnection(): ...@@ -33,8 +108,6 @@ class MayaConnection():
self.base2 = base2 self.base2 = base2
self.suff2 = suff2 self.suff2 = suff2
# MayaConnections = dict()
def ParseMayaConnection(line, base_dict): def ParseMayaConnection(line, base_dict):
words = line.split() words = line.split()
if len(words) != 3: if len(words) != 3:
...@@ -63,79 +136,14 @@ def ParseMayaConnection(line, base_dict): ...@@ -63,79 +136,14 @@ def ParseMayaConnection(line, base_dict):
return MayaConnection(helix1,base1,suff1,helix2,base2,suff2) return MayaConnection(helix1,base1,suff1,helix2,base2,suff2)
return None 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): class MayaBase(MayaObj):
def __init__(self, maya_lines): def __init__(self, maya_lines):
MayaObj.__init__(self, maya_lines) MayaObj.__init__(self, maya_lines)
self._parse_maya_lines(maya_lines) self._parse_maya_lines(maya_lines)
self.basepair = None self.basepair = None
self.parent = None
self.end3 = None self.end3 = None
self.end5 = 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): def add_basepair(self, base):
assert( self.parent == base.parent ) assert( self.parent == base.parent )
assert( self.basepair is None or self.basepair == base ) assert( self.basepair is None or self.basepair == base )
...@@ -151,27 +159,22 @@ setAttr ".t" -type "double3" 0.36767788771440857 -0.78848777472188547 -6.014 ; ...@@ -151,27 +159,22 @@ setAttr ".t" -type "double3" 0.36767788771440857 -0.78848777472188547 -6.014 ;
def add_end5(self, base): def add_end5(self, base):
base.add_end3(self) 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): 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() helices = dict()
bases = [] bases = []
aimConstraints = [] aim_constraints = []
connections = [] connections = []
## Parse through .ma file """ Parse .ma file """
with open(maya_file) as fh: with open(maya_file) as fh:
linesBuffer = [] linesBuffer = []
for line in fh: for line in fh:
...@@ -182,12 +185,15 @@ def parse_maya_file(maya_file): ...@@ -182,12 +185,15 @@ def parse_maya_file(maya_file):
if len(linesBuffer) > 0: if len(linesBuffer) > 0:
objType = linesBuffer[0].split()[1] objType = linesBuffer[0].split()[1]
if objType == "vHelix": if objType == "vHelix":
h = MayaHelix( linesBuffer ) h = MayaObj( linesBuffer )
helices[h.name] = h helices[h.get_full_name()] = h
elif objType == "HelixBase": elif objType == "HelixBase":
bases.append( MayaBase( linesBuffer ) ) bases.append( MayaBase( linesBuffer ) )
# elif objType == "aimConstraint": elif objType == "transform":
# aimConstraints.append( MayaAimConstraint( linesBuffer ) ) o = MayaObj( linesBuffer )
objects.append(o)
elif objType == "aimConstraint":
aim_constraints.append( MayaObj( linesBuffer ) )
## Clear lines buffer ## Clear lines buffer
linesBuffer = [] linesBuffer = []
elif words[0] == "connectAttr": elif words[0] == "connectAttr":
...@@ -195,13 +201,8 @@ def parse_maya_file(maya_file): ...@@ -195,13 +201,8 @@ def parse_maya_file(maya_file):
## Extend lines buffer ## Extend lines buffer
linesBuffer.append(line) 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 = [] new_connections = []
base_dict = { b.name:b for b in bases } base_dict = { b.name:b for b in bases }
for line in connections: for line in connections:
...@@ -210,7 +211,7 @@ def parse_maya_file(maya_file): ...@@ -210,7 +211,7 @@ def parse_maya_file(maya_file):
new_connections.append(conn) new_connections.append(conn)
connections = new_connections connections = new_connections
## Assign base connectivity """ Assign base connectivity """
for c in connections: for c in connections:
h1 = c.helix1 h1 = c.helix1
bn1 = c.base1 bn1 = c.base1
...@@ -219,349 +220,166 @@ def parse_maya_file(maya_file): ...@@ -219,349 +220,166 @@ def parse_maya_file(maya_file):
bn2 = c.base2 bn2 = c.base2
s2 = c.suff2 s2 = c.suff2
b1 = helices[h1].bases[bn1] b1 = helices[h1].children[bn1]
b2 = helices[h2].bases[bn2] b2 = helices[h2].children[bn2]
if s1 == "lb" and s2 == "lb": if s1 == "lb" and s2 == "lb":
b1.add_basepair(b2) 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": elif s1 == "fw" and s2 == "bw":
b1.add_end3(b2)
elif s1 == "bw" and s2 == "fw":
b1.add_end5(b2) b1.add_end5(b2)
elif s1 == "bw" and s2 == "fw":
b1.add_end3(b2)
else: else:
raise Exception("Unkown connection %s %s" % (s1,s2)) raise Exception("Unrecognized connection %s %s" % (s1,s2))
# # base_positions = np.zeros( [ len(bases), 3 ] ) """ Find local basis of each base so that
# base_positions = np.array( [b.position for b in bases] ) model_from_basepair_stack_3prime can determine stacks """
# for aim in aimConstraints:
# base = base_dict[ aim.parent_name ]
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): shift = np.array((1.19, -3.52, 0.08))
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
def get_end5_ssDNA(nt): if DEBUG:
ret = [] write_pdb_psf(bases, 'before')
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
""" Build dsDNA helices """ def update_geometry(b,beads):
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)
## Find helical axis """ Update the position and orientation of bead b using the 5
# if len(paired_bases) < 6: neighboring bases in "beads" as a reference """
# 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)]
if top not in connectionToSingleStrand and top.end3 is not None: def _get_beads_index(i):
b2 = top.end3 j=0
h2 = b2.parent.name while j <= i:
s2 = segments[h2] if beads[j] is None:
if b2 == segmentsBot[h2] or b2.basepair == segmentsBot[h2]: i+=1
seg.connect_end3(s2.start5,type_="terminal_crossover") j+=1
elif b2 == segmentsTop[h2] or b2.basepair == segmentsTop[h2]: return j-1
seg.connect_end3(s2.end5,type_="terminal_crossover")
else: """ Get position and filter out beads == None """
raise Exception("BUG") 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: """ Assign position and orientation to base """
b2 = top_bp.end5 if b.parent is None:
h2 = b2.parent.name o_parent = np.eye(3)
s2 = segments[h2] else:
if b2 == segmentsBot[h2] or b2.basepair == segmentsBot[h2]: o_parent = b.parent.get_orientation()
seg.connect_end5(s2.start3,type_="terminal_crossover")
elif b2 == segmentsTop[h2] or b2.basepair == segmentsTop[h2]: b.orientation = o_parent.T.dot( R )
seg.connect_end5(s2.end3, type_="terminal_crossover") b.position = b.position + b.orientation.dot( shift )
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")
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: def convert_maya_bases_to_segment_model(maya_bases, **model_parameters):
fh.write("%d\n\n" % nbeads)
for hn,h in helices.items():
h.write_xyz(fh)
# h.find_double_stranded_regions
""" Converts bases parse from .ma file into mrdna an model """
if __name__ == "__main__": from .segmentmodel_from_lists import model_from_basepair_stack_3prime
maya_helices = parse_maya_file("/home/cmaffeo2/Downloads/bunny.ma")
segments,dsSegmentDict = convert_maya_to_segments(maya_helices) base_to_index = {b:i for i,b in zip(range(len(maya_bases)),maya_bases)}
coord = np.zeros((len(maya_bases),3))
model = SegmentModel( segments, orientation = np.zeros((len(maya_bases),3,3))
local_twist = True, bp,end3 = [-np.ones((len(maya_bases)),dtype=np.int) for i in range(2)]
max_basepairs_per_bead = 1,
max_nucleotides_per_bead = 1,
dimensions=(5000,5000,5000), DEBUG=True )
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)
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment