diff --git a/mrdna/readers/__init__.py b/mrdna/readers/__init__.py index 31a6cedd5f7f92641d81ac6691280103a97c4e51..4ac260c8886a8bbdc7ca68754fb8b5fe6a8803f2 100644 --- a/mrdna/readers/__init__.py +++ b/mrdna/readers/__init__.py @@ -13,7 +13,6 @@ def read_vhelix(maya_file, **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 def read_list(infile,**model_parameters): diff --git a/mrdna/readers/polygon_mesh.py b/mrdna/readers/polygon_mesh.py index 2eb0015f740581ac76548f6998334a68c5a20dc4..7f96df1509d0cdedcb852f425a775f951cdd384e 100644 --- a/mrdna/readers/polygon_mesh.py +++ b/mrdna/readers/polygon_mesh.py @@ -13,6 +13,8 @@ maya_obj_registry_short = dict() class MayaObj(): """ Class representing a node in a Maya ascii file """ + seq_mapping = {str(i):c for c,i in zip('ATCG',range(4))} + def __init__(self, maya_lines): self.parent = None self.children = dict() @@ -59,6 +61,11 @@ class MayaObj(): 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) @@ -138,6 +145,7 @@ def ParseMayaConnection(line, base_dict): class MayaBase(MayaObj): def __init__(self, maya_lines): + self.sequence = '?' MayaObj.__init__(self, maya_lines) self._parse_maya_lines(maya_lines) self.basepair = None @@ -371,6 +379,10 @@ def convert_maya_bases_to_segment_model(maya_bases, **model_parameters): 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)] + 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() @@ -380,6 +392,15 @@ def convert_maya_bases_to_segment_model(maya_bases, **model_parameters): 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)