htsworkflow/pipelines/bustard.py

   1 """
   2 Extract configuration from Illumina Bustard Directory.
   3
   4 This includes the version number, run date, bustard executable parameters, and
   5 phasing estimates.
   6 """
   7 from copy import copy
   8 from datetime import date
   9 from glob import glob
  10 import logging
  11 import os
  12 import re
  13 import sys
  14 import time
  15
  16 from htsworkflow.pipelines.runfolder import \
  17    ElementTree, \
  18    VERSION_RE, \
  19    EUROPEAN_STRPTIME
  20
  21 # make epydoc happy
  22 __docformat__ = "restructuredtext en"
  23
  24 LANE_LIST = range(1,9)
  25
  26 class Phasing(object):
  27     PHASING = 'Phasing'
  28     PREPHASING = 'Prephasing'
  29
  30     def __init__(self, fromfile=None, xml=None):
  31         self.lane = None
  32         self.phasing = None
  33         self.prephasing = None
  34
  35         if fromfile is not None:
  36             self._initialize_from_file(fromfile)
  37         elif xml is not None:
  38             self.set_elements(xml)
  39
  40     def _initialize_from_file(self, pathname):
  41         path, name = os.path.split(pathname)
  42         basename, ext = os.path.splitext(name)
  43         # the last character of the param base filename should be the
  44         # lane number
  45         tree = ElementTree.parse(pathname).getroot()
  46         self.set_elements(tree)
  47         self.lane = int(basename[-1])
  48
  49     def get_elements(self):
  50         root = ElementTree.Element(Phasing.PHASING, {'lane': str(self.lane)})
  51         root.tail = os.linesep
  52         phasing = ElementTree.SubElement(root, Phasing.PHASING)
  53         phasing.text = str(self.phasing)
  54         phasing.tail = os.linesep
  55         prephasing = ElementTree.SubElement(root, Phasing.PREPHASING)
  56         prephasing.text = str(self.prephasing)
  57         prephasing.tail = os.linesep
  58         return root
  59
  60     def set_elements(self, tree):
  61         if tree.tag not in ('Phasing', 'Parameters'):
  62             raise ValueError('exptected Phasing or Parameters')
  63         lane = tree.attrib.get('lane', None)
  64         if lane is not None:
  65             self.lane = int(lane)
  66         for element in list(tree):
  67             if element.tag == Phasing.PHASING:
  68                 self.phasing = float(element.text)
  69             elif element.tag == Phasing.PREPHASING:
  70                 self.prephasing = float(element.text)
  71
  72 class CrosstalkMatrix(object):
  73     CROSSTALK = "MatrixElements"
  74     BASE = 'Base'
  75     ELEMENT = 'Element'
  76
  77     def __init__(self, fromfile=None, xml=None):
  78         self.base = {}
  79
  80         if fromfile is not None:
  81             self._initialize_from_file(fromfile)
  82         elif xml is not None:
  83             self.set_elements(xml)
  84
  85     def _initialize_from_file(self, pathname):
  86         data = open(pathname).readlines()
  87         auto_header = '# Auto-generated frequency response matrix'
  88         if data[0].strip() != auto_header or len(data) != 9:
  89             raise RuntimeError("matrix file %s is unusual" % (pathname,))
  90         # skip over lines 1,2,3,4 which contain the 4 bases
  91         self.base['A'] = [ float(v) for v in data[5].split() ]
  92         self.base['C'] = [ float(v) for v in data[6].split() ]
  93         self.base['G'] = [ float(v) for v in data[7].split() ]
  94         self.base['T'] = [ float(v) for v in data[8].split() ]
  95
  96     def get_elements(self):
  97         root = ElementTree.Element(CrosstalkMatrix.CROSSTALK)
  98         root.tail = os.linesep
  99         base_order = ['A','C','G','T']
 100         for b in base_order:
 101             base_element = ElementTree.SubElement(root, CrosstalkMatrix.BASE)
 102             base_element.text = b
 103             base_element.tail = os.linesep
 104         for b in base_order:
 105             for value in self.base[b]:
 106                 crosstalk_value = ElementTree.SubElement(root, CrosstalkMatrix.ELEMENT)
 107                 crosstalk_value.text = unicode(value)
 108                 crosstalk_value.tail = os.linesep
 109
 110         return root
 111
 112     def set_elements(self, tree):
 113         if tree.tag != CrosstalkMatrix.CROSSTALK:
 114             raise ValueError('Invalid run-xml exptected '+CrosstalkMatrix.CROSSTALK)
 115         base_order = []
 116         current_base = None
 117         current_index = 0
 118         for element in tree.getchildren():
 119             # read in the order of the bases
 120             if element.tag == 'Base':
 121                 base_order.append(element.text)
 122             elif element.tag == 'Element':
 123                 # we're done reading bases, now its just the 4x4 matrix
 124                 # written out as a list of elements
 125                 # if this is the first element, make a copy of the list
 126                 # to play with and initialize an empty list for the current base
 127                 if current_base is None:
 128                     current_base = copy(base_order)
 129                     self.base[current_base[0]] = []
 130                 # we found (probably) 4 bases go to the next base
 131                 if current_index == len(base_order):
 132                     current_base.pop(0)
 133                     current_index = 0
 134                     self.base[current_base[0]] = []
 135                 value = float(element.text)
 136                 self.base[current_base[0]].append(value)
 137
 138                 current_index += 1
 139             else:
 140                 raise RuntimeError("Unrecognized tag in run xml: %s" %(element.tag,))
 141
 142 def crosstalk_matrix_from_bustard_config(bustard_path, bustard_config_tree):
 143     """
 144     Analyze the bustard config file and try to find the crosstalk matrix.
 145     """
 146     bustard_run = bustard_config_tree[0]
 147     if bustard_run.tag != 'Run':
 148         raise RuntimeError('Expected Run tag, got %s' % (bustard_run.tag,))
 149
 150     call_parameters = bustard_run.find('BaseCallParameters')
 151     if call_parameters is None:
 152         raise RuntimeError('Missing BaseCallParameters section')
 153
 154     matrix = call_parameters.find('Matrix')
 155     if matrix is None:
 156         raise RuntimeError('Expected to find Matrix in Bustard BaseCallParameters')
 157
 158     matrix_auto_flag = int(matrix.find('AutoFlag').text)
 159     matrix_auto_lane = int(matrix.find('AutoLane').text)
 160
 161     crosstalk = None
 162     if matrix_auto_flag:
 163         # we estimated the matrix from something in this run.
 164         # though we don't really care which lane it was
 165         matrix_path = os.path.join(bustard_path, 'Matrix', 's_02_matrix.txt')
 166         crosstalk = CrosstalkMatrix(matrix_path)
 167     else:
 168         matrix_elements = call_parameters.find('MatrixElements')
 169         # the matrix was provided
 170         if matrix_elements is not None:
 171             crosstalk = CrosstalkMatrix(xml=matrix_elements)
 172         else:
 173             # we have no crosstalk matrix?
 174             pass
 175
 176     return crosstalk
 177
 178 class Bustard(object):
 179     XML_VERSION = 2
 180
 181     # Xml Tags
 182     BUSTARD = 'Bustard'
 183     SOFTWARE_VERSION = 'version'
 184     DATE = 'run_time'
 185     USER = 'user'
 186     PARAMETERS = 'Parameters'
 187     BUSTARD_CONFIG = 'BaseCallAnalysis'
 188
 189     def __init__(self, xml=None):
 190         self.version = None
 191         self.date = None
 192         self.user = None
 193         self.phasing = {}
 194         self.crosstalk = None
 195         self.pathname = None
 196         self.bustard_config = None
 197
 198         if xml is not None:
 199             self.set_elements(xml)
 200
 201     def _get_time(self):
 202         if self.date is None:
 203             return None
 204         return time.mktime(self.date.timetuple())
 205     time = property(_get_time, doc='return run time as seconds since epoch')
 206
 207     def dump(self):
 208         #print ElementTree.tostring(self.get_elements())
 209         ElementTree.dump(self.get_elements())
 210
 211     def get_elements(self):
 212         root = ElementTree.Element('Bustard',
 213                                    {'version': str(Bustard.XML_VERSION)})
 214         version = ElementTree.SubElement(root, Bustard.SOFTWARE_VERSION)
 215         version.text = self.version
 216         if self.time is not None:
 217             run_date = ElementTree.SubElement(root, Bustard.DATE)
 218             run_date.text = str(self.time)
 219         if self.user is not None:
 220             user = ElementTree.SubElement(root, Bustard.USER)
 221             user.text = self.user
 222         params = ElementTree.SubElement(root, Bustard.PARAMETERS)
 223
 224         # add phasing parameters
 225         for lane in LANE_LIST:
 226             if self.phasing.has_key(lane):
 227                 params.append(self.phasing[lane].get_elements())
 228
 229         # add crosstalk matrix if it exists
 230         if self.crosstalk is not None:
 231             root.append(self.crosstalk.get_elements())
 232
 233         # add bustard config if it exists
 234         if self.bustard_config is not None:
 235             root.append(self.bustard_config)
 236         return root
 237
 238     def set_elements(self, tree):
 239         if tree.tag != Bustard.BUSTARD:
 240             raise ValueError('Expected "Bustard" SubElements')
 241         xml_version = int(tree.attrib.get('version', 0))
 242         if xml_version > Bustard.XML_VERSION:
 243             logging.warn('Bustard XML tree is a higher version than this class')
 244         for element in list(tree):
 245             if element.tag == Bustard.SOFTWARE_VERSION:
 246                 self.version = element.text
 247             elif element.tag == Bustard.DATE:
 248                 self.date = date.fromtimestamp(float(element.text))
 249             elif element.tag == Bustard.USER:
 250                 self.user = element.text
 251             elif element.tag == Bustard.PARAMETERS:
 252                 for param in element:
 253                     p = Phasing(xml=param)
 254                     self.phasing[p.lane] = p
 255             elif element.tag == CrosstalkMatrix.CROSSTALK:
 256                 self.crosstalk = CrosstalkMatrix(xml=element)
 257             elif element.tag == Bustard.BUSTARD_CONFIG:
 258                 self.bustard_config = element
 259             else:
 260                 raise ValueError("Unrecognized tag: %s" % (element.tag,))
 261
 262 def bustard(pathname):
 263     """
 264     Construct a Bustard object by analyzing an Illumina Bustard directory.
 265
 266     :Parameters:
 267       - `pathname`: A bustard directory
 268
 269     :Return:
 270       Fully initialized Bustard object.
 271     """
 272     b = Bustard()
 273     pathname = os.path.abspath(pathname)
 274     path, name = os.path.split(pathname)
 275     groups = name.split("_")
 276     if groups[0].lower().startswith('bustard'):
 277         version = re.search(VERSION_RE, groups[0])
 278         b.version = version.group(1)
 279         t = time.strptime(groups[1], EUROPEAN_STRPTIME)
 280         b.date = date(*t[0:3])
 281         b.user = groups[2]
 282     elif groups[0] == 'BaseCalls':
 283         # stub values
 284         b.version = None
 285         b.date = None
 286         b.user = None
 287
 288     b.pathname = pathname
 289     bustard_config_filename = os.path.join(pathname, 'config.xml')
 290     paramfiles = glob(os.path.join(pathname, "params?.xml"))
 291     for paramfile in paramfiles:
 292         phasing = Phasing(paramfile)
 293         assert (phasing.lane >= 1 and phasing.lane <= 8)
 294         b.phasing[phasing.lane] = phasing
 295     # I only found these in Bustard1.9.5/1.9.6 directories
 296     if b.version in ('1.9.5', '1.9.6'):
 297         # at least for our runfolders for 1.9.5 and 1.9.6 matrix[1-8].txt are always the same
 298         crosstalk_file = os.path.join(pathname, "matrix1.txt")
 299         b.crosstalk = CrosstalkMatrix(crosstalk_file)
 300     # for version 1.3.2 of the pipeline the bustard version number went down
 301     # to match the rest of the pipeline. However there's now a nifty
 302     # new (useful) bustard config file.
 303     elif os.path.exists(bustard_config_filename):
 304         bustard_config_root = ElementTree.parse(bustard_config_filename)
 305         b.bustard_config = bustard_config_root.getroot()
 306         b.crosstalk = crosstalk_matrix_from_bustard_config(b.pathname, b.bustard_config)
 307         software = bustard_config_root.find('*/Software')
 308         b.version = software.attrib['Version']
 309         #b.version = software.attrib['Name'] + "-" + software.attrib['Version']
 310
 311     return b
 312
 313 def fromxml(tree):
 314     """
 315     Reconstruct a htsworkflow.pipelines.Bustard object from an xml block
 316     """
 317     b = Bustard()
 318     b.set_elements(tree)
 319     return b
 320
 321 def make_cmdline_parser():
 322     from optparse import OptionParser
 323     parser = OptionParser('%prog: bustard_directory')
 324     return parser
 325
 326 def main(cmdline):
 327     parser = make_cmdline_parser()
 328     opts, args = parser.parse_args(cmdline)
 329
 330     for bustard_dir in args:
 331         print u'analyzing bustard directory: ' + unicode(bustard_dir)
 332         bustard_object = bustard(bustard_dir)
 333         bustard_object.dump()
 334
 335         bustard_object2 = Bustard(xml=bustard_object.get_elements())
 336         print ('-------------------------------------')
 337         bustard_object2.dump()
 338         print ('=====================================')
 339         b1_tree = bustard_object.get_elements()
 340         b1 = ElementTree.tostring(b1_tree).split(os.linesep)
 341         b2_tree = bustard_object2.get_elements()
 342         b2 = ElementTree.tostring(b2_tree).split(os.linesep)
 343         for line1, line2 in zip(b1, b2):
 344             if b1 != b2:
 345                 print "b1: ", b1
 346                 print "b2: ", b2
 347
 348 if __name__ == "__main__":
 349     main(sys.argv[1:])