import numpy as np
import sys
import re
from ..arbdmodel.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()
[docs]
class MayaObj():
""" Class representing a node in a Maya ascii file """
seq_mapping = {str(i):c for c,i in zip('ATCG',range(4))}
[docs]
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):
""" Extracts name and parent from maya file """
name = parent = None
vals = l.split()
for i in range(len(vals)):
if vals[i] == "-n":
name = MayaObj._strip_maya_string( vals[i+1] )
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))
elif vals[1] == '".lb"':
try:
self.sequence = MayaObj.seq_mapping[vals[2][0]]
except:
pass
def _strip_maya_string(string):
return re.match('"\|?([a-zA-Z_0-9|]+)";?$', string).group(1)
[docs]
def add(self, obj):
self.children[obj.name] = obj
obj.parent = self
[docs]
def get_position(self):
if self.parent is not None:
return self.parent.get_orientation().dot(self.position) + self.parent.get_position()
else:
return self.position
[docs]
def get_orientation(self):
if self.parent is not None:
return self.parent.get_orientation().dot(self.orientation)
else:
return self.orientation
[docs]
def get_full_name(self):
ret = self.name
if self.parent is not None:
ret = '{}|{}'.format(self.parent.get_full_name(), ret)
return ret
[docs]
class MayaConnection():
[docs]
def __init__(self, helix1, base1, suff1, helix2, base2, suff2):
self.helix1 = helix1
self.base1 = base1
self.suff1 = suff1
self.helix2 = helix2
self.base2 = base2
self.suff2 = suff2
[docs]
def ParseMayaConnection(line, base_dict):
words = line.split()
if len(words) != 3:
return None
cmd,obj1,obj2 = words
def parse_obj(obj):
## Typically, base name is not unique, so parent name used too
result = re.match('"(.*(forw)?(backw)?_[0-9]+)\.([a-z]{2})"', obj)
if result is not None:
name, tmp, tmp, suff = result.groups()
if name in base_dict:
return base_dict[name].parent.name, name, suff
else:
result = re.match('\|(.*)\|(.*)', name)
if result is not None:
helix, name = result.groups()
return helix, name, suff
return None, None, None
helix1,base1,suff1 = parse_obj(obj1)
if helix1 is not None:
helix2,base2,suff2 = parse_obj(obj2)
if helix2 is not None:
return MayaConnection(helix1,base1,suff1,helix2,base2,suff2)
return None
[docs]
class MayaBase(MayaObj):
[docs]
def __init__(self, maya_lines):
self.sequence = '?'
MayaObj.__init__(self, maya_lines)
self._parse_maya_lines(maya_lines)
self.basepair = None
self.end3 = None
self.end5 = None
[docs]
def add_basepair(self, base):
assert( self.parent == base.parent )
assert( self.basepair is None or self.basepair == base )
assert( base.basepair is None or base.basepair == self )
self.basepair = base
base.basepair = self
[docs]
def add_end3(self, base):
assert( self.end3 is None or self.end3 == base )
assert( base.end5 is None or base.end5 == self )
self.end3 = base
base.end5 = self
[docs]
def add_end5(self, base):
base.add_end3(self)
def __repr__(self):
return "Base {} {:.2f}".format(self.name, self.position[2])
[docs]
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 = []
aim_constraints = []
connections = []
""" Parse .ma file """
with open(maya_file) as fh:
linesBuffer = []
for line in fh:
words = line.split()
if words[0] == "createNode":
## Create a new object
if len(linesBuffer) > 0:
objType = linesBuffer[0].split()[1]
if objType == "vHelix":
h = MayaObj( linesBuffer )
helices[h.get_full_name()] = h
elif objType == "HelixBase":
bases.append( MayaBase( 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":
connections.append(line)
## Extend lines buffer
linesBuffer.append(line)
""" Parse connections """
new_connections = []
base_dict = { b.name:b for b in bases }
for line in connections:
conn = ParseMayaConnection(line, base_dict)
if conn is not None:
new_connections.append(conn)
connections = new_connections
""" Assign base connectivity """
for c in connections:
h1 = c.helix1
bn1 = c.base1
s1 = c.suff1
h2 = c.helix2
bn2 = c.base2
s2 = c.suff2
b1 = helices[h1].children[bn1]
b2 = helices[h2].children[bn2]
if s1 == "lb" and s2 == "lb":
b1.add_basepair(b2)
elif s1 == "fw" and s2 == "bw":
b1.add_end5(b2)
elif s1 == "bw" and s2 == "fw":
b1.add_end3(b2)
else:
raise Exception("Unrecognized connection %s %s" % (s1,s2))
""" Find local basis of each base so that
model_from_basepair_stack_3prime can determine stacks """
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))
ref_positions = np.array((ref_below_position,
ref_stack_position,
ref_bp_position,
ref_bp_below_position,
ref_bp_stack_position))
shift = np.array((1.19, -3.52, 0.08))
if DEBUG:
write_pdb_psf(bases, 'before')
def update_geometry(b,beads):
""" Update the position and orientation of bead b using the 5
neighboring bases in "beads" as a reference """
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)
""" 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 )
""" 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
[docs]
def write_pdb_psf(bases, prefix):
""" Function for debugging .ma parser """
from ..arbdmodel import ArbdModel,ParticleType, PointParticle, Group
from ..arbdmodel.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')
[docs]
def convert_maya_bases_to_segment_model(maya_bases, **model_parameters):
""" Converts bases parse from .ma file into mrdna an model """
from .segmentmodel_from_lists import model_from_basepair_stack_3prime
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=int) for i in range(2)]
sequence = []
seq_map = {'A':'T', 'C':'G', 'T':'A', 'G':'C'}
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 ]
try:
s = b.sequence
if s == '?':
s = seq_map[s.basepair.sequence]
except:
s = 'T'
sequence.append( s )
if 'sequence' not in model_parameters:
model_parameters['sequence'] = sequence
return model_from_basepair_stack_3prime(coord, bp, None, end3,
orientation=orientation,
**model_parameters)