+class CrosstalkMatrix(object):
+ CROSSTALK = "MatrixElements"
+ BASE = 'Base'
+ ELEMENT = 'Element'
+
+ def __init__(self, fromfile=None, xml=None):
+ self.base = {}
+
+ if fromfile is not None:
+ self._initialize_from_file(fromfile)
+ elif xml is not None:
+ self.set_elements(xml)
+
+ def _initialize_from_file(self, pathname):
+ data = open(pathname).readlines()
+ auto_header = '# Auto-generated frequency response matrix'
+ if data[0].strip() != auto_header or len(data) != 9:
+ raise RuntimeError("matrix file %s is unusual" % (pathname,))
+ # skip over lines 1,2,3,4 which contain the 4 bases
+ self.base['A'] = [ float(v) for v in data[5].split() ]
+ self.base['C'] = [ float(v) for v in data[6].split() ]
+ self.base['G'] = [ float(v) for v in data[7].split() ]
+ self.base['T'] = [ float(v) for v in data[8].split() ]
+
+ def get_elements(self):
+ root = ElementTree.Element(CrosstalkMatrix.CROSSTALK)
+ root.tail = os.linesep
+ base_order = ['A','C','G','T']
+ for b in base_order:
+ base_element = ElementTree.SubElement(root, CrosstalkMatrix.BASE)
+ base_element.text = b
+ base_element.tail = os.linesep
+ for b in base_order:
+ for value in self.base[b]:
+ crosstalk_value = ElementTree.SubElement(root, CrosstalkMatrix.ELEMENT)
+ crosstalk_value.text = unicode(value)
+ crosstalk_value.tail = os.linesep
+
+ return root
+
+ def set_elements(self, tree):
+ if tree.tag != CrosstalkMatrix.CROSSTALK:
+ raise ValueError('Invalid run-xml exptected '+CrosstalkMatrix.CROSSTALK)
+ base_order = []
+ current_base = None
+ current_index = 0
+ for element in tree.getchildren():
+ # read in the order of the bases
+ if element.tag == 'Base':
+ base_order.append(element.text)
+ elif element.tag == 'Element':
+ # we're done reading bases, now its just the 4x4 matrix
+ # written out as a list of elements
+ # if this is the first element, make a copy of the list
+ # to play with and initialize an empty list for the current base
+ if current_base is None:
+ current_base = copy(base_order)
+ self.base[current_base[0]] = []
+ # we found (probably) 4 bases go to the next base
+ if current_index == len(base_order):
+ current_base.pop(0)
+ current_index = 0
+ self.base[current_base[0]] = []
+ value = float(element.text)
+ self.base[current_base[0]].append(value)
+
+ current_index += 1
+ else:
+ raise RuntimeError("Unrecognized tag in run xml: %s" %(element.tag,))
+
+def crosstalk_matrix_from_bustard_config(bustard_path, bustard_config_tree):
+ """
+ Analyze the bustard config file and try to find the crosstalk matrix.
+ """
+ bustard_run = bustard_config_tree[0]
+ if bustard_run.tag != 'Run':
+ raise RuntimeError('Expected Run tag, got %s' % (bustard_run.tag,))
+
+ call_parameters = bustard_run.find('BaseCallParameters')
+ if call_parameters is None:
+ raise RuntimeError('Missing BaseCallParameters section')
+
+ matrix = call_parameters.find('Matrix')
+ if matrix is None:
+ raise RuntimeError('Expected to find Matrix in Bustard BaseCallParameters')
+
+ matrix_auto_flag = int(matrix.find('AutoFlag').text)
+ matrix_auto_lane = int(matrix.find('AutoLane').text)
+
+ if matrix_auto_flag:
+ # we estimated the matrix from something in this run.
+ # though we don't really care which lane it was
+ matrix_path = os.path.join(bustard_path, 'Matrix', 's_02_matrix.txt')
+ matrix = CrosstalkMatrix(matrix_path)
+ else:
+ # the matrix was provided
+ matrix_elements = call_parameters.find('MatrixElements')
+ if matrix_elements is None:
+ raise RuntimeError('Expected to find MatrixElements in Bustard BaseCallParameters')
+ matrix = CrosstalkMatrix(xml=matrix_elements)
+
+ print "matrix:", matrix
+ return matrix
+
projects / htsworkflow.git / commitdiff