f2687671c2e95278a532a02e21dfd29788b2d7b2
[htsworkflow.git] / gaworkflow / pipeline / gerald.py
1 """
2 Provide access to information stored in the GERALD directory.
3 """
4 from datetime import datetime, date
5 from glob import glob
6 import logging
7 import os
8 import stat
9 import time
10 import types
11
12 from gaworkflow.pipeline.runfolder import \
13    ElementTree, \
14    EUROPEAN_STRPTIME, \
15    LANES_PER_FLOWCELL, \
16    VERSION_RE
17 from gaworkflow.util.ethelp import indent, flatten
18 from gaworkflow.util.opener import autoopen
19
20 class Gerald(object):
21     """
22     Capture meaning out of the GERALD directory
23     """
24     XML_VERSION = 1
25     GERALD='Gerald'
26     RUN_PARAMETERS='RunParameters'
27     SUMMARY='Summary'
28
29     class LaneParameters(object):
30         """
31         Make it easy to access elements of LaneSpecificRunParameters from python
32         """
33         def __init__(self, gerald, key):
34             self._gerald = gerald
35             self._key = key
36         
37         def __get_attribute(self, xml_tag):
38             subtree = self._gerald.tree.find('LaneSpecificRunParameters')
39             container = subtree.find(xml_tag)
40             if container is None:
41                 return None
42             if len(container.getchildren()) > LANES_PER_FLOWCELL:
43                 raise RuntimeError('GERALD config.xml file changed')
44             lanes = [x.tag.split('_')[1] for x in container.getchildren()]
45             index = lanes.index(self._key)
46             element = container[index]
47             return element.text
48         def _get_analysis(self):
49             return self.__get_attribute('ANALYSIS')
50         analysis = property(_get_analysis)
51
52         def _get_eland_genome(self):
53             genome = self.__get_attribute('ELAND_GENOME')
54             # default to the chipwide parameters if there isn't an
55             # entry in the lane specific paramaters
56             if genome is None:
57                 subtree = self._gerald.tree.find('ChipWideRunParameters')
58                 container = subtree.find('ELAND_GENOME')
59                 genome = container.text
60             return genome
61         eland_genome = property(_get_eland_genome)
62
63         def _get_read_length(self):
64             return self.__get_attribute('READ_LENGTH')
65         read_length = property(_get_read_length)
66
67         def _get_use_bases(self):
68             return self.__get_attribute('USE_BASES')
69         use_bases = property(_get_use_bases)
70
71     class LaneSpecificRunParameters(object):
72         """
73         Provide access to LaneSpecificRunParameters
74         """
75         def __init__(self, gerald):
76             self._gerald = gerald
77             self._keys = None
78         def __getitem__(self, key):
79             return Gerald.LaneParameters(self._gerald, key)
80         def keys(self):
81             if self._keys is None:
82                 tree = self._gerald.tree
83                 analysis = tree.find('LaneSpecificRunParameters/ANALYSIS')
84                 # according to the pipeline specs I think their fields 
85                 # are sampleName_laneID, with sampleName defaulting to s
86                 # since laneIDs are constant lets just try using 
87                 # those consistently.
88                 self._keys = [ x.tag.split('_')[1] for x in analysis]
89             return self._keys
90         def values(self):
91             return [ self[x] for x in self.keys() ]
92         def items(self):
93             return zip(self.keys(), self.values())
94         def __len__(self):
95             return len(self.keys())
96
97     def __init__(self, xml=None):
98         self.pathname = None
99         self.tree = None
100
101         # parse lane parameters out of the config.xml file
102         self.lanes = Gerald.LaneSpecificRunParameters(self)
103
104         self.summary = None
105         self.eland_results = None
106
107         if xml is not None:
108             self.set_elements(xml)
109
110     def _get_date(self):
111         if self.tree is None:
112             return datetime.today()
113         timestamp = self.tree.findtext('ChipWideRunParameters/TIME_STAMP')
114         epochstamp = time.mktime(time.strptime(timestamp, '%c'))
115         return datetime.fromtimestamp(epochstamp)
116     date = property(_get_date)
117
118     def _get_time(self):
119         return time.mktime(self.date.timetuple())
120     time = property(_get_time, doc='return run time as seconds since epoch')
121
122     def _get_version(self):
123         if self.tree is None:
124             return None
125         return self.tree.findtext('ChipWideRunParameters/SOFTWARE_VERSION')
126     version = property(_get_version)
127
128     def dump(self):
129         """
130         Debugging function, report current object
131         """
132         print 'Gerald version:', self.version
133         print 'Gerald run date:', self.date
134         print 'Gerald config.xml:', self.tree
135         self.summary.dump()
136
137     def get_elements(self):
138         if self.tree is None or self.summary is None:
139             return None
140
141         gerald = ElementTree.Element(Gerald.GERALD, 
142                                      {'version': unicode(Gerald.XML_VERSION)})
143         gerald.append(self.tree)
144         gerald.append(self.summary.get_elements())
145         if self.eland_results:
146             gerald.append(self.eland_results.get_elements())
147         return gerald
148
149     def set_elements(self, tree):
150         if tree.tag !=  Gerald.GERALD:
151             raise ValueError('exptected GERALD')
152         xml_version = int(tree.attrib.get('version', 0))
153         if xml_version > Gerald.XML_VERSION:
154             logging.warn('XML tree is a higher version than this class')
155         for element in list(tree):
156             tag = element.tag.lower()
157             if tag == Gerald.RUN_PARAMETERS.lower():
158                 self.tree = element
159             elif tag == Gerald.SUMMARY.lower():
160                 self.summary = Summary(xml=element)
161             elif tag == ELAND.ELAND.lower():
162                 self.eland_results = ELAND(xml=element)
163             else:
164                 logging.warn("Unrecognized tag %s" % (element.tag,))
165         
166
167 def gerald(pathname):
168     g = Gerald()
169     g.pathname = pathname
170     path, name = os.path.split(pathname)
171     config_pathname = os.path.join(pathname, 'config.xml')
172     g.tree = ElementTree.parse(config_pathname).getroot()
173
174     # parse Summary.htm file
175     summary_pathname = os.path.join(pathname, 'Summary.htm')
176     g.summary = Summary(summary_pathname)
177     # parse eland files
178     g.eland_results = eland(g.pathname, g)
179     return g
180
181 def tonumber(v):
182     """
183     Convert a value to int if its an int otherwise a float.
184     """
185     try:
186         v = int(v)
187     except ValueError, e:
188         v = float(v)
189     return v
190
191 def parse_mean_range(value):
192     """
193     Parse values like 123 +/- 4.5
194     """
195     if value.strip() == 'unknown':
196         return 0, 0
197
198     average, pm, deviation = value.split()
199     if pm != '+/-':
200         raise RuntimeError("Summary.htm file format changed")
201     return tonumber(average), tonumber(deviation)
202
203 def make_mean_range_element(parent, name, mean, deviation):
204     """
205     Make an ElementTree subelement <Name mean='mean', deviation='deviation'/>
206     """
207     element = ElementTree.SubElement(parent, name,
208                                      { 'mean': unicode(mean),
209                                        'deviation': unicode(deviation)})
210     return element
211
212 def parse_mean_range_element(element):
213     """
214     Grab mean/deviation out of element
215     """
216     return (tonumber(element.attrib['mean']), 
217             tonumber(element.attrib['deviation']))
218
219 def parse_summary_element(element):
220     """
221     Determine if we have a simple element or a mean/deviation element
222     """
223     if len(element.attrib) > 0:
224         return parse_mean_range_element(element)
225     else:
226         return element.text
227
228 class Summary(object):
229     """
230     Extract some useful information from the Summary.htm file
231     """
232     XML_VERSION = 2
233     SUMMARY = 'Summary'
234
235     class LaneResultSummary(object):
236         """
237         Parse the LaneResultSummary table out of Summary.htm
238         Mostly for the cluster number
239         """
240         LANE_RESULT_SUMMARY = 'LaneResultSummary'
241         TAGS = { 
242           'LaneYield': 'lane_yield',
243           'Cluster': 'cluster', # Raw
244           'ClusterPF': 'cluster_pass_filter',
245           'AverageFirstCycleIntensity': 'average_first_cycle_intensity',
246           'PercentIntensityAfter20Cycles': 'percent_intensity_after_20_cycles',
247           'PercentPassFilterClusters': 'percent_pass_filter_clusters',
248           'PercentPassFilterAlign': 'percent_pass_filter_align',
249           'AverageAlignmentScore': 'average_alignment_score',
250           'PercentErrorRate': 'percent_error_rate'
251         }
252                  
253         def __init__(self, html=None, xml=None):
254             self.lane = None
255             self.lane_yield = None
256             self.cluster = None
257             self.cluster_pass_filter = None
258             self.average_first_cycle_intensity = None
259             self.percent_intensity_after_20_cycles = None
260             self.percent_pass_filter_clusters = None
261             self.percent_pass_filter_align = None
262             self.average_alignment_score = None
263             self.percent_error_rate = None
264
265             if html is not None:
266                 self.set_elements_from_html(html)
267             if xml is not None:
268                 self.set_elements(xml)
269
270         def set_elements_from_html(self, data):
271             if not len(data) in (8,10):
272                 raise RuntimeError("Summary.htm file format changed")
273
274             # same in pre-0.3.0 Summary file and 0.3 summary file
275             self.lane = data[0]
276
277             if len(data) == 8:
278                 parsed_data = [ parse_mean_range(x) for x in data[1:] ]
279                 # this is the < 0.3 Pipeline version
280                 self.cluster = parsed_data[0]
281                 self.average_first_cycle_intensity = parsed_data[1]
282                 self.percent_intensity_after_20_cycles = parsed_data[2]
283                 self.percent_pass_filter_clusters = parsed_data[3]
284                 self.percent_pass_filter_align = parsed_data[4]
285                 self.average_alignment_score = parsed_data[5]
286                 self.percent_error_rate = parsed_data[6]
287             elif len(data) == 10:
288                 parsed_data = [ parse_mean_range(x) for x in data[2:] ]
289                 # this is the >= 0.3 summary file
290                 self.lane_yield = data[1]
291                 self.cluster = parsed_data[0]
292                 self.cluster_pass_filter = parsed_data[1]
293                 self.average_first_cycle_intensity = parsed_data[2]
294                 self.percent_intensity_after_20_cycles = parsed_data[3]
295                 self.percent_pass_filter_clusters = parsed_data[4]
296                 self.percent_pass_filter_align = parsed_data[5]
297                 self.average_alignment_score = parsed_data[6]
298                 self.percent_error_rate = parsed_data[7]
299
300         def get_elements(self):
301             lane_result = ElementTree.Element(
302                             Summary.LaneResultSummary.LANE_RESULT_SUMMARY, 
303                             {'lane': self.lane})
304             for tag, variable_name in Summary.LaneResultSummary.TAGS.items():
305                 value = getattr(self, variable_name)
306                 if value is None:
307                     continue
308                 # it looks like a sequence
309                 elif type(value) in (types.TupleType, types.ListType):
310                     element = make_mean_range_element(
311                       lane_result,
312                       tag,
313                       *value
314                     )
315                 else:
316                     element = ElementTree.SubElement(lane_result, tag)
317                     element.text = value
318             return lane_result
319
320         def set_elements(self, tree):
321             if tree.tag != Summary.LaneResultSummary.LANE_RESULT_SUMMARY:
322                 raise ValueError('Expected %s' % (
323                         Summary.LaneResultSummary.LANE_RESULT_SUMMARY))
324             self.lane = tree.attrib['lane']
325             tags = Summary.LaneResultSummary.TAGS
326             for element in list(tree):
327                 try:
328                     variable_name = tags[element.tag]
329                     setattr(self, variable_name, 
330                             parse_summary_element(element))
331                 except KeyError, e:
332                     logging.warn('Unrecognized tag %s' % (element.tag,))
333
334     def __init__(self, filename=None, xml=None):
335         self.lane_results = {}
336
337         if filename is not None:
338             self._extract_lane_results(filename)
339         if xml is not None:
340             self.set_elements(xml)
341
342     def __getitem__(self, key):
343         return self.lane_results[key]
344
345     def __len__(self):
346         return len(self.lane_results)
347
348     def keys(self):
349         return self.lane_results.keys()
350
351     def values(self):
352         return self.lane_results.values()
353
354     def items(self):
355         return self.lane_results.items()
356
357     def _flattened_row(self, row):
358         """
359         flatten the children of a <tr>...</tr>
360         """
361         return [flatten(x) for x in row.getchildren() ]
362     
363     def _parse_table(self, table):
364         """
365         assumes the first line is the header of a table, 
366         and that the remaining rows are data
367         """
368         rows = table.getchildren()
369         data = []
370         for r in rows:
371             data.append(self._flattened_row(r))
372         return data
373     
374     def _extract_named_tables(self, pathname):
375         """
376         extract all the 'named' tables from a Summary.htm file
377         and return as a dictionary
378         
379         Named tables are <h2>...</h2><table>...</table> pairs
380         The contents of the h2 tag is considered to the name
381         of the table.
382         """
383         tree = ElementTree.parse(pathname).getroot()
384         body = tree.find('body')
385         tables = {}
386         for i in range(len(body)):
387             if body[i].tag == 'h2' and body[i+1].tag == 'table':
388                 # we have an interesting table
389                 name = flatten(body[i])
390                 table = body[i+1]
391                 data = self._parse_table(table)
392                 tables[name] = data
393         return tables
394
395     def _extract_lane_results(self, pathname):
396         """
397         extract the Lane Results Summary table
398         """
399
400         tables = self._extract_named_tables(pathname)
401
402         # parse lane result summary
403         lane_summary = tables['Lane Results Summary']
404         # this is version 1 of the summary file
405         if len(lane_summary[-1]) == 8:
406             # strip header
407             headers = lane_summary[0]
408             # grab the lane by lane data
409             lane_summary = lane_summary[1:]
410
411         # this is version 2 of the summary file
412         if len(lane_summary[-1]) == 10:
413             # lane_summary[0] is a different less specific header row
414             headers = lane_summary[1]
415             lane_summary = lane_summary[2:10]
416             # after the last lane, there's a set of chip wide averages
417
418         for r in lane_summary:
419             lrs = Summary.LaneResultSummary(html=r)
420             self.lane_results[lrs.lane] = lrs
421
422     def get_elements(self):
423         summary = ElementTree.Element(Summary.SUMMARY, 
424                                       {'version': unicode(Summary.XML_VERSION)})
425         for lane in self.lane_results.values():
426             summary.append(lane.get_elements())
427         return summary
428
429     def set_elements(self, tree):
430         if tree.tag != Summary.SUMMARY:
431             return ValueError("Expected %s" % (Summary.SUMMARY,))
432         xml_version = int(tree.attrib.get('version', 0))
433         if xml_version > Summary.XML_VERSION:
434             logging.warn('Summary XML tree is a higher version than this class')
435         for element in list(tree):
436             lrs = Summary.LaneResultSummary()
437             lrs.set_elements(element)
438             self.lane_results[lrs.lane] = lrs
439
440     def dump(self):
441         """
442         Debugging function, report current object
443         """
444         pass
445
446
447 def build_genome_fasta_map(genome_dir):
448     # build fasta to fasta file map
449     genome = genome_dir.split(os.path.sep)[-1]
450     fasta_map = {}
451     for vld_file in glob(os.path.join(genome_dir, '*.vld')):
452         is_link = False
453         if os.path.islink(vld_file):
454             is_link = True
455         vld_file = os.path.realpath(vld_file)
456         path, vld_name = os.path.split(vld_file)
457         name, ext = os.path.splitext(vld_name)
458         if is_link:
459             fasta_map[name] = name
460         else:
461             fasta_map[name] = os.path.join(genome, name)
462     return fasta_map
463     
464 class ElandLane(object):
465     """
466     Process an eland result file
467     """
468     XML_VERSION = 1
469     LANE = 'ElandLane'
470     SAMPLE_NAME = 'SampleName'
471     LANE_ID = 'LaneID'
472     GENOME_MAP = 'GenomeMap'
473     GENOME_ITEM = 'GenomeItem'
474     MAPPED_READS = 'MappedReads'
475     MAPPED_ITEM = 'MappedItem'
476     MATCH_CODES = 'MatchCodes'
477     MATCH_ITEM = 'Code'
478     READS = 'Reads'
479
480     def __init__(self, pathname=None, genome_map=None, xml=None):
481         self.pathname = pathname
482         self._sample_name = None
483         self._lane_id = None
484         self._reads = None
485         self._mapped_reads = None
486         self._match_codes = None
487         if genome_map is None:
488             genome_map = {}
489         self.genome_map = genome_map
490         
491         if xml is not None:
492             self.set_elements(xml)
493
494     def _update(self):
495         """
496         Actually read the file and actually count the reads
497         """
498         # can't do anything if we don't have a file to process
499         if self.pathname is None:
500             return
501
502         if os.stat(self.pathname)[stat.ST_SIZE] == 0:
503             raise RuntimeError("Eland isn't done, try again later.")
504
505         reads = 0
506         mapped_reads = {}
507
508         match_codes = {'NM':0, 'QC':0, 'RM':0, 
509                        'U0':0, 'U1':0, 'U2':0,
510                        'R0':0, 'R1':0, 'R2':0,
511                       }
512         for line in autoopen(self.pathname,'r'):
513             reads += 1
514             fields = line.split()
515             # code = fields[2]
516             # match_codes[code] = match_codes.setdefault(code, 0) + 1
517             # the QC/NM etc codes are in the 3rd field and always present
518             match_codes[fields[2]] += 1
519             # ignore lines that don't have a fasta filename
520             if len(fields) < 7:
521                 continue
522             fasta = self.genome_map.get(fields[6], fields[6])
523             mapped_reads[fasta] = mapped_reads.setdefault(fasta, 0) + 1
524         self._match_codes = match_codes
525         self._mapped_reads = mapped_reads
526         self._reads = reads
527
528     def _update_name(self):
529         # extract the sample name
530         if self.pathname is None:
531             return
532
533         path, name = os.path.split(self.pathname)
534         split_name = name.split('_')
535         self._sample_name = split_name[0]
536         self._lane_id = split_name[1]
537
538     def _get_sample_name(self):
539         if self._sample_name is None:
540             self._update_name()
541         return self._sample_name
542     sample_name = property(_get_sample_name)
543
544     def _get_lane_id(self):
545         if self._lane_id is None:
546             self._update_name()
547         return self._lane_id
548     lane_id = property(_get_lane_id)
549
550     def _get_reads(self):
551         if self._reads is None:
552             self._update()
553         return self._reads
554     reads = property(_get_reads)
555
556     def _get_mapped_reads(self):
557         if self._mapped_reads is None:
558             self._update()
559         return self._mapped_reads
560     mapped_reads = property(_get_mapped_reads)
561
562     def _get_match_codes(self):
563         if self._match_codes is None:
564             self._update()
565         return self._match_codes
566     match_codes = property(_get_match_codes)
567
568     def get_elements(self):
569         lane = ElementTree.Element(ElandLane.LANE, 
570                                    {'version': 
571                                     unicode(ElandLane.XML_VERSION)})
572         sample_tag = ElementTree.SubElement(lane, ElandLane.SAMPLE_NAME)
573         sample_tag.text = self.sample_name
574         lane_tag = ElementTree.SubElement(lane, ElandLane.LANE_ID)
575         lane_tag.text = self.lane_id
576         genome_map = ElementTree.SubElement(lane, ElandLane.GENOME_MAP)
577         for k, v in self.genome_map.items():
578             item = ElementTree.SubElement(
579                 genome_map, ElandLane.GENOME_ITEM, 
580                 {'name':k, 'value':unicode(v)})
581         mapped_reads = ElementTree.SubElement(lane, ElandLane.MAPPED_READS)
582         for k, v in self.mapped_reads.items():
583             item = ElementTree.SubElement(
584                 mapped_reads, ElandLane.MAPPED_ITEM, 
585                 {'name':k, 'value':unicode(v)})
586         match_codes = ElementTree.SubElement(lane, ElandLane.MATCH_CODES)
587         for k, v in self.match_codes.items():
588             item = ElementTree.SubElement(
589                 match_codes, ElandLane.MATCH_ITEM, 
590                 {'name':k, 'value':unicode(v)})
591         reads = ElementTree.SubElement(lane, ElandLane.READS)
592         reads.text = unicode(self.reads)
593
594         return lane
595
596     def set_elements(self, tree):
597         if tree.tag != ElandLane.LANE:
598             raise ValueError('Exptecting %s' % (ElandLane.LANE,))
599
600         # reset dictionaries
601         self._mapped_reads = {}
602         self._match_codes = {}
603         
604         for element in tree:
605             tag = element.tag.lower()
606             if tag == ElandLane.SAMPLE_NAME.lower():
607                 self._sample_name = element.text
608             elif tag == ElandLane.LANE_ID.lower():
609                 self._lane_id = element.text
610             elif tag == ElandLane.GENOME_MAP.lower():
611                 for child in element:
612                     name = child.attrib['name']
613                     value = child.attrib['value']
614                     self.genome_map[name] = value
615             elif tag == ElandLane.MAPPED_READS.lower():
616                 for child in element:
617                     name = child.attrib['name']
618                     value = child.attrib['value']
619                     self._mapped_reads[name] = int(value)
620             elif tag == ElandLane.MATCH_CODES.lower():
621                 for child in element:
622                     name = child.attrib['name']
623                     value = int(child.attrib['value'])
624                     self._match_codes[name] = value
625             elif tag == ElandLane.READS.lower():
626                 self._reads = int(element.text)
627             else:
628                 logging.warn("ElandLane unrecognized tag %s" % (element.tag,))
629
630 def extract_eland_sequence(instream, outstream, start, end):
631     """
632     Extract a chunk of sequence out of an eland file
633     """
634     for line in instream:
635         record = line.split()
636         if len(record) > 1:
637             result = [record[0], record[1][start:end]]
638         else:
639             result = [record[0][start:end]]
640         outstream.write("\t".join(result))
641         outstream.write(os.linesep)
642
643 class ELAND(object):
644     """
645     Summarize information from eland files
646     """
647     XML_VERSION = 1
648
649     ELAND = 'ElandCollection'
650     LANE = 'Lane'
651     LANE_ID = 'id'
652
653     def __init__(self, xml=None):
654         # we need information from the gerald config.xml
655         self.results = {}
656         
657         if xml is not None:
658             self.set_elements(xml)
659
660     def __len__(self):
661         return len(self.results)
662
663     def keys(self):
664         return self.results.keys()
665     
666     def values(self):
667         return self.results.values()
668
669     def items(self):
670         return self.results.items()
671
672     def __getitem__(self, key):
673         return self.results[key]
674
675     def get_elements(self):
676         root = ElementTree.Element(ELAND.ELAND, 
677                                    {'version': unicode(ELAND.XML_VERSION)})
678         for lane_id, lane in self.results.items():
679             eland_lane = lane.get_elements()
680             eland_lane.attrib[ELAND.LANE_ID] = unicode(lane_id)
681             root.append(eland_lane)
682         return root
683
684     def set_elements(self, tree):
685         if tree.tag.lower() != ELAND.ELAND.lower():
686             raise ValueError('Expecting %s', ELAND.ELAND)
687         for element in list(tree):
688             lane_id = element.attrib[ELAND.LANE_ID]
689             lane = ElandLane(xml=element)
690             self.results[lane_id] = lane
691
692 def eland(basedir, gerald=None, genome_maps=None):
693     e = ELAND()
694
695     file_list = glob(os.path.join(basedir, "*_eland_result.txt"))
696     if len(file_list) == 0:
697         # lets handle compressed eland files too
698         file_list = glob(os.path.join(basedir, "*_eland_result.txt.bz2"))
699
700     for pathname in file_list:
701         # yes the lane_id is also being computed in ElandLane._update
702         # I didn't want to clutter up my constructor
703         # but I needed to persist the sample_name/lane_id for
704         # runfolder summary_report
705         path, name = os.path.split(pathname)
706         split_name = name.split('_')
707         lane_id = split_name[1]
708
709         if genome_maps is not None:
710             genome_map = genome_maps[lane_id]
711         elif gerald is not None:
712             genome_dir = gerald.lanes[lane_id].eland_genome
713             genome_map = build_genome_fasta_map(genome_dir)
714         else:
715             genome_map = {}
716
717         eland_result = ElandLane(pathname, genome_map)
718         e.results[lane_id] = eland_result
719     return e