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 LOGGER = logging.getLogger(__name__)
  22
  23 # make epydoc happy
  24 __docformat__ = "restructuredtext en"
  25
  26 LANE_LIST = range(1,9)
  27
  28 class Phasing(object):
  29     PHASING = 'Phasing'
  30     PREPHASING = 'Prephasing'
  31
  32     def __init__(self, fromfile=None, xml=None):
  33         self.lane = None
  34         self.phasing = None
  35         self.prephasing = None
  36
  37         if fromfile is not None:
  38             self._initialize_from_file(fromfile)
  39         elif xml is not None:
  40             self.set_elements(xml)
  41
  42     def _initialize_from_file(self, pathname):
  43         path, name = os.path.split(pathname)
  44         basename, ext = os.path.splitext(name)
  45         # the last character of the param base filename should be the
  46         # lane number
  47         tree = ElementTree.parse(pathname).getroot()
  48         self.set_elements(tree)
  49         self.lane = int(basename[-1])
  50
  51     def get_elements(self):
  52         root = ElementTree.Element(Phasing.PHASING, {'lane': str(self.lane)})
  53         root.tail = os.linesep
  54         phasing = ElementTree.SubElement(root, Phasing.PHASING)
  55         phasing.text = str(self.phasing)
  56         phasing.tail = os.linesep
  57         prephasing = ElementTree.SubElement(root, Phasing.PREPHASING)
  58         prephasing.text = str(self.prephasing)
  59         prephasing.tail = os.linesep
  60         return root
  61
  62     def set_elements(self, tree):
  63         if tree.tag not in ('Phasing', 'Parameters'):
  64             raise ValueError('exptected Phasing or Parameters')
  65         lane = tree.attrib.get('lane', None)
  66         if lane is not None:
  67             self.lane = int(lane)
  68         for element in list(tree):
  69             if element.tag == Phasing.PHASING:
  70                 self.phasing = float(element.text)
  71             elif element.tag == Phasing.PREPHASING:
  72                 self.prephasing = float(element.text)
  73
  74 class CrosstalkMatrix(object):
  75     CROSSTALK = "MatrixElements"
  76     BASE = 'Base'
  77     ELEMENT = 'Element'
  78
  79     def __init__(self, fromfile=None, xml=None):
  80         self.base = {}
  81
  82         if fromfile is not None:
  83             self._initialize_from_file(fromfile)
  84         elif xml is not None:
  85             self.set_elements(xml)
  86
  87     def _initialize_from_file(self, pathname):
  88         data = open(pathname).readlines()
  89         auto_header = '# Auto-generated frequency response matrix'
  90         if data[0].strip() == auto_header and len(data) == 9:
  91             # skip over lines 1,2,3,4 which contain the 4 bases
  92             self.base['A'] = [ float(v) for v in data[5].split() ]
  93             self.base['C'] = [ float(v) for v in data[6].split() ]
  94             self.base['G'] = [ float(v) for v in data[7].split() ]
  95             self.base['T'] = [ float(v) for v in data[8].split() ]
  96         elif len(data) == 16:
  97             self.base['A'] = [ float(v) for v in data[:4] ]
  98             self.base['C'] = [ float(v) for v in data[4:8] ]
  99             self.base['G'] = [ float(v) for v in data[8:12] ]
 100             self.base['T'] = [ float(v) for v in data[12:16] ]
 101         else:
 102             raise RuntimeError("matrix file %s is unusual" % (pathname,))
 103     def get_elements(self):
 104         root = ElementTree.Element(CrosstalkMatrix.CROSSTALK)
 105         root.tail = os.linesep
 106         base_order = ['A','C','G','T']
 107         for b in base_order:
 108             base_element = ElementTree.SubElement(root, CrosstalkMatrix.BASE)
 109             base_element.text = b
 110             base_element.tail = os.linesep
 111         for b in base_order:
 112             for value in self.base[b]:
 113                 crosstalk_value = ElementTree.SubElement(root, CrosstalkMatrix.ELEMENT)
 114                 crosstalk_value.text = unicode(value)
 115                 crosstalk_value.tail = os.linesep
 116
 117         return root
 118
 119     def set_elements(self, tree):
 120         if tree.tag != CrosstalkMatrix.CROSSTALK:
 121             raise ValueError('Invalid run-xml exptected '+CrosstalkMatrix.CROSSTALK)
 122         base_order = []
 123         current_base = None
 124         current_index = 0
 125         for element in tree.getchildren():
 126             # read in the order of the bases
 127             if element.tag == 'Base':
 128                 base_order.append(element.text)
 129             elif element.tag == 'Element':
 130                 # we're done reading bases, now its just the 4x4 matrix
 131                 # written out as a list of elements
 132                 # if this is the first element, make a copy of the list
 133                 # to play with and initialize an empty list for the current base
 134                 if current_base is None:
 135                     current_base = copy(base_order)
 136                     self.base[current_base[0]] = []
 137                 # we found (probably) 4 bases go to the next base
 138                 if current_index == len(base_order):
 139                     current_base.pop(0)
 140                     current_index = 0
 141                     self.base[current_base[0]] = []
 142                 value = float(element.text)
 143                 self.base[current_base[0]].append(value)
 144
 145                 current_index += 1
 146             else:
 147                 raise RuntimeError("Unrecognized tag in run xml: %s" %(element.tag,))
 148
 149 def crosstalk_matrix_from_bustard_config(bustard_path, bustard_config_tree):
 150     """
 151     Analyze the bustard config file and try to find the crosstalk matrix.
 152     """
 153     bustard_run = bustard_config_tree[0]
 154     if bustard_run.tag != 'Run':
 155         raise RuntimeError('Expected Run tag, got %s' % (bustard_run.tag,))
 156
 157     call_parameters = bustard_run.find('BaseCallParameters')
 158     if call_parameters is None:
 159         raise RuntimeError('Missing BaseCallParameters section')
 160
 161     matrix = call_parameters.find('Matrix')
 162     if matrix is None:
 163         raise RuntimeError('Expected to find Matrix in Bustard BaseCallParameters')
 164
 165     matrix_auto_flag = int(matrix.find('AutoFlag').text)
 166     matrix_auto_lane = int(matrix.find('AutoLane').text)
 167
 168     crosstalk = None
 169     if matrix_auto_flag:
 170         # we estimated the matrix from something in this run.
 171         # though we don't really care which lane it was
 172         if matrix_auto_lane == 0:
 173             # its defaulting to all of the lanes, so just pick one
 174             auto_lane_fragment = "_1"
 175         else:
 176             auto_lane_fragment = "_%d" % ( matrix_auto_lane,)
 177
 178         for matrix_name in ['s%s_02_matrix.txt' % (auto_lane_fragment,),
 179                             's%s_1_matrix.txt' % (auto_lane_fragment,),
 180                             ]:
 181             matrix_path = os.path.join(bustard_path, 'Matrix', matrix_name)
 182             if os.path.exists(matrix_path):
 183                 break
 184         else:
 185             raise RuntimeError("Couldn't find matrix for lane %d" % \
 186                                (matrix_auto_lane,))
 187
 188         crosstalk = CrosstalkMatrix(matrix_path)
 189     else:
 190         matrix_elements = call_parameters.find('MatrixElements')
 191         # the matrix was provided
 192         if matrix_elements is not None:
 193             crosstalk = CrosstalkMatrix(xml=matrix_elements)
 194         else:
 195             # we have no crosstalk matrix?
 196             pass
 197
 198     return crosstalk
 199
 200 class Bustard(object):
 201     XML_VERSION = 2
 202
 203     # Xml Tags
 204     BUSTARD = 'Bustard'
 205     SOFTWARE_VERSION = 'version'
 206     DATE = 'run_time'
 207     USER = 'user'
 208     PARAMETERS = 'Parameters'
 209     BUSTARD_CONFIG = 'BaseCallAnalysis'
 210
 211     def __init__(self, xml=None):
 212         self._path_version = None # version number from directory name
 213         self.date = None
 214         self.user = None
 215         self.phasing = {}
 216         self.crosstalk = None
 217         self.pathname = None
 218         self.bustard_config = None
 219
 220         if xml is not None:
 221             self.set_elements(xml)
 222
 223     def update_attributes_from_pathname(self):
 224         """Update version, date, user from bustard directory names
 225         Obviously this wont work for BaseCalls or later runfolders
 226         """
 227         if self.pathname is None:
 228             raise ValueError(
 229                 "Set pathname before calling update_attributes_from_pathname")
 230         path, name = os.path.split(self.pathname)
 231
 232         if not re.match('bustard', name, re.IGNORECASE):
 233             return
 234
 235         groups = name.split("_")
 236         version = re.search(VERSION_RE, groups[0])
 237         self._path_version = version.group(1)
 238         t = time.strptime(groups[1], EUROPEAN_STRPTIME)
 239         self.date = date(*t[0:3])
 240         self.user = groups[2]
 241
 242     def _get_sequence_format(self):
 243         """Guess sequence format"""
 244         projects = glob(os.path.join(self.pathname, 'Project_*'))
 245         if len(projects) > 0:
 246             # Hey we look like a demultiplexed run
 247             return 'fastq'
 248         return 'qseq'
 249         #qseqs = glob(os.path.join(self.pathname, '*_qseq.txt'))
 250         #if len(qseqs) > 0:
 251         #    return 'qseq'
 252         #seqs = glob(os.path.join(self.pathname, '*_seq.txt'))
 253         #if len(seqs) > 0:
 254         #    return 'srf'
 255         return None
 256     sequence_format = property(_get_sequence_format)
 257
 258     def _get_software_version(self):
 259         """return software name, version tuple"""
 260         if self.bustard_config is None:
 261             if self._path_version is not None:
 262                 return 'Bustard', self._path_version
 263             else:
 264                 return None
 265         software_nodes = self.bustard_config.xpath('Run/Software')
 266         if len(software_nodes) == 0:
 267             return None
 268         elif len(software_nodes) > 1:
 269             raise RuntimeError("Too many software XML elements for bustard.py")
 270         else:
 271             return (software_nodes[0].attrib['Name'],
 272                     software_nodes[0].attrib['Version'])
 273
 274     def _get_software(self):
 275         """Return software name"""
 276         software_version = self._get_software_version()
 277         return software_version[0] if software_version is not None else None
 278     software = property(_get_software)
 279
 280     def _get_version(self):
 281         """Return software name"""
 282         software_version = self._get_software_version()
 283         return software_version[1] if software_version is not None else None
 284     version = property(_get_version)
 285
 286
 287     def _get_time(self):
 288         if self.date is None:
 289             return None
 290         return time.mktime(self.date.timetuple())
 291     time = property(_get_time, doc='return run time as seconds since epoch')
 292
 293     def dump(self):
 294         #print ElementTree.tostring(self.get_elements())
 295         ElementTree.dump(self.get_elements())
 296
 297     def get_elements(self):
 298         root = ElementTree.Element('Bustard',
 299                                    {'version': str(Bustard.XML_VERSION)})
 300         version = ElementTree.SubElement(root, Bustard.SOFTWARE_VERSION)
 301         version.text = self.version
 302         if self.time is not None:
 303             run_date = ElementTree.SubElement(root, Bustard.DATE)
 304             run_date.text = str(self.time)
 305         if self.user is not None:
 306             user = ElementTree.SubElement(root, Bustard.USER)
 307             user.text = self.user
 308         params = ElementTree.SubElement(root, Bustard.PARAMETERS)
 309
 310         # add phasing parameters
 311         for lane in LANE_LIST:
 312             if self.phasing.has_key(lane):
 313                 params.append(self.phasing[lane].get_elements())
 314
 315         # add crosstalk matrix if it exists
 316         if self.crosstalk is not None:
 317             root.append(self.crosstalk.get_elements())
 318
 319         # add bustard config if it exists
 320         if self.bustard_config is not None:
 321             root.append(self.bustard_config)
 322         return root
 323
 324     def set_elements(self, tree):
 325         if tree.tag != Bustard.BUSTARD:
 326             raise ValueError('Expected "Bustard" SubElements')
 327         xml_version = int(tree.attrib.get('version', 0))
 328         if xml_version > Bustard.XML_VERSION:
 329             LOGGER.warn('Bustard XML tree is a higher version than this class')
 330         for element in list(tree):
 331             if element.tag == Bustard.SOFTWARE_VERSION:
 332                 self._path_version = element.text
 333             elif element.tag == Bustard.DATE:
 334                 self.date = date.fromtimestamp(float(element.text))
 335             elif element.tag == Bustard.USER:
 336                 self.user = element.text
 337             elif element.tag == Bustard.PARAMETERS:
 338                 for param in element:
 339                     p = Phasing(xml=param)
 340                     self.phasing[p.lane] = p
 341             elif element.tag == CrosstalkMatrix.CROSSTALK:
 342                 self.crosstalk = CrosstalkMatrix(xml=element)
 343             elif element.tag == Bustard.BUSTARD_CONFIG:
 344                 self.bustard_config = element
 345             else:
 346                 raise ValueError("Unrecognized tag: %s" % (element.tag,))
 347
 348 def bustard(pathname):
 349     """
 350     Construct a Bustard object by analyzing an Illumina Bustard directory.
 351
 352     :Parameters:
 353       - `pathname`: A bustard directory
 354
 355     :Return:
 356       Fully initialized Bustard object.
 357     """
 358     pathname = os.path.abspath(pathname)
 359     bustard_filename = os.path.join(pathname, 'config.xml')
 360     demultiplexed_filename = os.path.join(pathname,
 361                                           'DemultiplexedBustardConfig.xml')
 362
 363     if os.path.exists(demultiplexed_filename):
 364         b = bustard_from_hiseq(pathname, demultiplexed_filename)
 365     elif os.path.exists(bustard_filename):
 366         b = bustard_from_ga2(pathname, bustard_filename)
 367     else:
 368         b = bustard_from_ga1(pathname)
 369
 370     return b
 371
 372 def bustard_from_ga1(pathname):
 373     """Initialize bustard class from ga1 era runfolders.
 374     """
 375     path, name = os.path.split(pathname)
 376
 377     groups = name.split("_")
 378     if len(groups) < 3:
 379         msg = "Not enough information to create attributes"\
 380               " from directory name: %s"
 381         LOGGER.error(msg % (pathname,))
 382         return None
 383
 384     b = Bustard()
 385     b.pathname = pathname
 386     b.update_attributes_from_pathname()
 387     version = re.search(VERSION_RE, groups[0])
 388     b._path_version = version.group(1)
 389     t = time.strptime(groups[1], EUROPEAN_STRPTIME)
 390     b.date = date(*t[0:3])
 391     b.user = groups[2]
 392
 393     # I only found these in Bustard1.9.5/1.9.6 directories
 394     if b.version in ('1.9.5', '1.9.6'):
 395         # at least for our runfolders for 1.9.5 and 1.9.6 matrix[1-8].txt are always the same
 396         crosstalk_file = os.path.join(pathname, "matrix1.txt")
 397         b.crosstalk = CrosstalkMatrix(crosstalk_file)
 398
 399     add_phasing(b)
 400     return b
 401
 402
 403 def bustard_from_ga2(pathname, config_filename):
 404     """Initialize bustard class from ga2-era runfolder
 405     Actually I don't quite remember if it is exactly the GA2s, but
 406     its after the original runfolder style and before the HiSeq.
 407     """
 408     # for version 1.3.2 of the pipeline the bustard version number went down
 409     # to match the rest of the pipeline. However there's now a nifty
 410     # new (useful) bustard config file.
 411
 412     # stub values
 413     b = Bustard()
 414     b.pathname = pathname
 415     b.update_attributes_from_pathname()
 416     bustard_config_root = ElementTree.parse(config_filename)
 417     b.bustard_config = bustard_config_root.getroot()
 418     b.crosstalk = crosstalk_matrix_from_bustard_config(b.pathname,
 419                                                        b.bustard_config)
 420     add_phasing(b)
 421
 422     return b
 423
 424 def bustard_from_hiseq(pathname, config_filename):
 425     b = Bustard()
 426     b.pathname = pathname
 427     bustard_config_root = ElementTree.parse(config_filename)
 428     b.bustard_config = bustard_config_root.getroot()
 429     add_phasing(b)
 430     return b
 431
 432 def add_phasing(bustard_obj):
 433     paramfiles = glob(os.path.join(bustard_obj.pathname,
 434                                    "params?.xml"))
 435     for paramfile in paramfiles:
 436         phasing = Phasing(paramfile)
 437         assert (phasing.lane >= 1 and phasing.lane <= 8)
 438         bustard_obj.phasing[phasing.lane] = phasing
 439
 440 def fromxml(tree):
 441     """
 442     Reconstruct a htsworkflow.pipelines.Bustard object from an xml block
 443     """
 444     b = Bustard()
 445     b.set_elements(tree)
 446     return b
 447
 448 def make_cmdline_parser():
 449     from optparse import OptionParser
 450     parser = OptionParser('%prog: bustard_directory')
 451     return parser
 452
 453 def main(cmdline):
 454     parser = make_cmdline_parser()
 455     opts, args = parser.parse_args(cmdline)
 456
 457     for bustard_dir in args:
 458         print u'analyzing bustard directory: ' + unicode(bustard_dir)
 459         bustard_object = bustard(bustard_dir)
 460         bustard_object.dump()
 461
 462         bustard_object2 = Bustard(xml=bustard_object.get_elements())
 463         print ('-------------------------------------')
 464         bustard_object2.dump()
 465         print ('=====================================')
 466         b1_tree = bustard_object.get_elements()
 467         b1 = ElementTree.tostring(b1_tree).split(os.linesep)
 468         b2_tree = bustard_object2.get_elements()
 469         b2 = ElementTree.tostring(b2_tree).split(os.linesep)
 470         for line1, line2 in zip(b1, b2):
 471             if b1 != b2:
 472                 print "b1: ", b1
 473                 print "b2: ", b2
 474
 475 if __name__ == "__main__":
 476     main(sys.argv[1:])