Added sequence to vHelix parser

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):

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)