Source code for mrdna.readers.polygon_mesh

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)