11 from htsworkflow.pipelines.runfolder import ElementTree
12 from htsworkflow.util.ethelp import indent, flatten
13 from htsworkflow.util.opener import autoopen
15 class ElandLane(object):
17 Process an eland result file
21 SAMPLE_NAME = 'SampleName'
24 GENOME_MAP = 'GenomeMap'
25 GENOME_ITEM = 'GenomeItem'
26 MAPPED_READS = 'MappedReads'
27 MAPPED_ITEM = 'MappedItem'
28 MATCH_CODES = 'MatchCodes'
37 def __init__(self, pathname=None, lane_id=None, end=None, genome_map=None, eland_type=None, xml=None):
38 self.pathname = pathname
39 self._sample_name = None
40 self.lane_id = lane_id
43 self._mapped_reads = None
44 self._match_codes = None
45 if genome_map is None:
47 self.genome_map = genome_map
48 self.eland_type = None
51 self.set_elements(xml)
53 def _guess_eland_type(self, pathname):
54 if self.eland_type is None:
55 # attempt autodetect eland file type
56 pathn, name = os.path.split(pathname)
57 if re.search('result', name):
58 self.eland_type = ElandLane.ELAND_SINGLE
59 elif re.search('multi', name):
60 self.eland_type = ElandLane.ELAND_MULTI
61 elif re.search('extended', name):
62 self.eland_type = ElandLane.ELAND_EXTENDED
63 elif re.search('export', name):
64 self.eland_type = ElandLane.ELAND_EXPORT
66 self.eland_type = ElandLane.ELAND_SINGLE
70 Actually read the file and actually count the reads
72 # can't do anything if we don't have a file to process
73 if self.pathname is None:
75 self._guess_eland_type(self.pathname)
77 if os.stat(self.pathname)[stat.ST_SIZE] == 0:
78 raise RuntimeError("Eland isn't done, try again later.")
80 logging.info("summarizing results for %s" % (self.pathname))
82 if self.eland_type == ElandLane.ELAND_SINGLE:
83 result = self._update_eland_result(self.pathname)
84 elif self.eland_type == ElandLane.ELAND_MULTI or \
85 self.eland_type == ElandLane.ELAND_EXTENDED:
86 result = self._update_eland_multi(self.pathname)
88 raise NotImplementedError("Only support single/multi/extended eland files")
89 self._match_codes, self._mapped_reads, self._reads = result
91 def _update_eland_result(self, pathname):
95 match_codes = {'NM':0, 'QC':0, 'RM':0,
96 'U0':0, 'U1':0, 'U2':0,
97 'R0':0, 'R1':0, 'R2':0,
99 for line in autoopen(pathname,'r'):
101 fields = line.split()
103 # match_codes[code] = match_codes.setdefault(code, 0) + 1
104 # the QC/NM etc codes are in the 3rd field and always present
105 match_codes[fields[2]] += 1
106 # ignore lines that don't have a fasta filename
109 fasta = self.genome_map.get(fields[6], fields[6])
110 mapped_reads[fasta] = mapped_reads.setdefault(fasta, 0) + 1
111 return match_codes, mapped_reads, reads
113 def _update_eland_multi(self, pathname):
117 match_codes = {'NM':0, 'QC':0, 'RM':0,
118 'U0':0, 'U1':0, 'U2':0,
119 'R0':0, 'R1':0, 'R2':0,
121 match_counts_re = re.compile("([\d]+):([\d]+):([\d]+)")
122 for line in autoopen(pathname,'r'):
124 fields = line.split()
125 # fields[2] = QC/NM/or number of matches
126 groups = match_counts_re.match(fields[2])
128 match_codes[fields[2]] += 1
130 # when there are too many hit, eland writes a - where
131 # it would have put the list of hits
132 if len(fields) < 4 or fields[3] == '-':
134 zero_mismatches = int(groups.group(1))
135 if zero_mismatches == 1:
136 match_codes['U0'] += 1
137 elif zero_mismatches < 255:
138 match_codes['R0'] += zero_mismatches
140 one_mismatches = int(groups.group(2))
141 if one_mismatches == 1:
142 match_codes['U1'] += 1
143 elif one_mismatches < 255:
144 match_codes['R1'] += one_mismatches
146 two_mismatches = int(groups.group(3))
147 if two_mismatches == 1:
148 match_codes['U2'] += 1
149 elif two_mismatches < 255:
150 match_codes['R2'] += two_mismatches
153 for match in fields[3].split(','):
154 match_fragment = match.split(':')
155 if len(match_fragment) == 2:
156 chromo = match_fragment[0]
157 pos = match_fragment[1]
159 fasta = self.genome_map.get(chromo, chromo)
160 assert fasta is not None
161 mapped_reads[fasta] = mapped_reads.setdefault(fasta, 0) + 1
162 return match_codes, mapped_reads, reads
164 def _update_name(self):
165 # extract the sample name
166 if self.pathname is None:
169 path, name = os.path.split(self.pathname)
170 split_name = name.split('_')
171 self._sample_name = split_name[0]
173 def _get_sample_name(self):
174 if self._sample_name is None:
176 return self._sample_name
177 sample_name = property(_get_sample_name)
179 def _get_reads(self):
180 if self._reads is None:
183 reads = property(_get_reads)
185 def _get_mapped_reads(self):
186 if self._mapped_reads is None:
188 return self._mapped_reads
189 mapped_reads = property(_get_mapped_reads)
191 def _get_match_codes(self):
192 if self._match_codes is None:
194 return self._match_codes
195 match_codes = property(_get_match_codes)
197 def _get_no_match(self):
198 if self._mapped_reads is None:
200 return self._match_codes['NM']
201 no_match = property(_get_no_match,
202 doc="total reads that didn't match the target genome.")
204 def _get_no_match_percent(self):
205 return float(self.no_match)/self.reads * 100
206 no_match_percent = property(_get_no_match_percent,
207 doc="no match reads as percent of total")
209 def _get_qc_failed(self):
210 if self._mapped_reads is None:
212 return self._match_codes['QC']
213 qc_failed = property(_get_qc_failed,
214 doc="total reads that didn't match the target genome.")
216 def _get_qc_failed_percent(self):
217 return float(self.qc_failed)/self.reads * 100
218 qc_failed_percent = property(_get_qc_failed_percent,
219 doc="QC failed reads as percent of total")
221 def _get_unique_reads(self):
222 if self._mapped_reads is None:
225 for code in ['U0','U1','U2']:
226 sum += self._match_codes[code]
228 unique_reads = property(_get_unique_reads,
229 doc="total unique reads")
231 def _get_repeat_reads(self):
232 if self._mapped_reads is None:
235 for code in ['R0','R1','R2']:
236 sum += self._match_codes[code]
238 repeat_reads = property(_get_repeat_reads,
239 doc="total repeat reads")
241 def get_elements(self):
242 lane = ElementTree.Element(ElandLane.LANE,
244 unicode(ElandLane.XML_VERSION)})
245 sample_tag = ElementTree.SubElement(lane, ElandLane.SAMPLE_NAME)
246 sample_tag.text = self.sample_name
247 lane_tag = ElementTree.SubElement(lane, ElandLane.LANE_ID)
248 lane_tag.text = str(self.lane_id)
249 if self.end is not None:
250 end_tag = ElementTree.SubElement(lane, ElandLane.END)
251 end_tag.text = str(self.end)
252 genome_map = ElementTree.SubElement(lane, ElandLane.GENOME_MAP)
253 for k, v in self.genome_map.items():
254 item = ElementTree.SubElement(
255 genome_map, ElandLane.GENOME_ITEM,
256 {'name':k, 'value':unicode(v)})
257 mapped_reads = ElementTree.SubElement(lane, ElandLane.MAPPED_READS)
258 for k, v in self.mapped_reads.items():
259 item = ElementTree.SubElement(
260 mapped_reads, ElandLane.MAPPED_ITEM,
261 {'name':k, 'value':unicode(v)})
262 match_codes = ElementTree.SubElement(lane, ElandLane.MATCH_CODES)
263 for k, v in self.match_codes.items():
264 item = ElementTree.SubElement(
265 match_codes, ElandLane.MATCH_ITEM,
266 {'name':k, 'value':unicode(v)})
267 reads = ElementTree.SubElement(lane, ElandLane.READS)
268 reads.text = unicode(self.reads)
272 def set_elements(self, tree):
273 if tree.tag != ElandLane.LANE:
274 raise ValueError('Exptecting %s' % (ElandLane.LANE,))
277 self._mapped_reads = {}
278 self._match_codes = {}
281 tag = element.tag.lower()
282 if tag == ElandLane.SAMPLE_NAME.lower():
283 self._sample_name = element.text
284 elif tag == ElandLane.LANE_ID.lower():
285 self.lane_id = int(element.text)
286 elif tag == ElandLane.END.lower():
287 self.end = int(element.text)
288 elif tag == ElandLane.GENOME_MAP.lower():
289 for child in element:
290 name = child.attrib['name']
291 value = child.attrib['value']
292 self.genome_map[name] = value
293 elif tag == ElandLane.MAPPED_READS.lower():
294 for child in element:
295 name = child.attrib['name']
296 value = child.attrib['value']
297 self._mapped_reads[name] = int(value)
298 elif tag == ElandLane.MATCH_CODES.lower():
299 for child in element:
300 name = child.attrib['name']
301 value = int(child.attrib['value'])
302 self._match_codes[name] = value
303 elif tag == ElandLane.READS.lower():
304 self._reads = int(element.text)
306 logging.warn("ElandLane unrecognized tag %s" % (element.tag,))
310 Summarize information from eland files
314 ELAND = 'ElandCollection'
319 def __init__(self, xml=None):
320 # we need information from the gerald config.xml
321 self.results = [{},{}]
324 self.set_elements(xml)
326 def get_elements(self):
327 root = ElementTree.Element(ELAND.ELAND,
328 {'version': unicode(ELAND.XML_VERSION)})
329 for end in range(len(self.results)):
330 end_results = self.results[end]
331 for lane_id, lane in end_results.items():
332 eland_lane = lane.get_elements()
333 eland_lane.attrib[ELAND.END] = unicode (end)
334 eland_lane.attrib[ELAND.LANE_ID] = unicode(lane_id)
335 root.append(eland_lane)
338 def set_elements(self, tree):
339 if tree.tag.lower() != ELAND.ELAND.lower():
340 raise ValueError('Expecting %s', ELAND.ELAND)
341 for element in list(tree):
342 lane_id = int(element.attrib[ELAND.LANE_ID])
343 end = int(element.attrib.get(ELAND.END, 0))
344 lane = ElandLane(xml=element)
345 self.results[end][lane_id] = lane
347 def check_for_eland_file(basedir, pattern, lane_id, end):
349 full_lane_id = lane_id
351 full_lane_id = "%d_%d" % ( lane_id, end )
353 basename = pattern % (full_lane_id,)
354 pathname = os.path.join(basedir, basename)
355 if os.path.exists(pathname):
360 def eland(basedir, gerald=None, genome_maps=None):
363 # if there is a basedir/Temp change basedir to point to the temp
364 # directory, as 1.1rc1 moves most of the files we've historically
365 # cared about to that subdirectory.
366 # we should look into what the official 'result' files are.
367 basedir_temp = os.path.join(basedir, 'Temp')
368 if os.path.isdir(basedir_temp):
369 basedir = basedir_temp
371 lane_ids = range(1,9)
374 # the order in patterns determines the preference for what
376 patterns = ['s_%s_eland_result.txt',
377 's_%s_eland_result.txt.bz2',
378 's_%s_eland_result.txt.gz',
379 's_%s_eland_extended.txt',
380 's_%s_eland_extended.txt.bz2',
381 's_%s_eland_extended.txt.gz',
382 's_%s_eland_multi.txt',
383 's_%s_eland_multi.txt.bz2',
384 's_%s_eland_multi.txt.gz',]
387 for lane_id in lane_ids:
389 pathname = check_for_eland_file(basedir, p, lane_id, end)
390 if pathname is not None:
394 # yes the lane_id is also being computed in ElandLane._update
395 # I didn't want to clutter up my constructor
396 # but I needed to persist the sample_name/lane_id for
397 # runfolder summary_report
398 path, name = os.path.split(pathname)
399 logging.info("Adding eland file %s" %(name,))
400 # split_name = name.split('_')
401 # lane_id = int(split_name[1])
403 if genome_maps is not None:
404 genome_map = genome_maps[lane_id]
405 elif gerald is not None:
406 genome_dir = gerald.lanes[lane_id].eland_genome
407 genome_map = build_genome_fasta_map(genome_dir)
411 eland_result = ElandLane(pathname, lane_id, end, genome_map)
415 effective_end = end - 1
416 e.results[effective_end][lane_id] = eland_result
419 def build_genome_fasta_map(genome_dir):
420 # build fasta to fasta file map
421 logging.info("Building genome map")
422 genome = genome_dir.split(os.path.sep)[-1]
424 for vld_file in glob(os.path.join(genome_dir, '*.vld')):
426 if os.path.islink(vld_file):
428 vld_file = os.path.realpath(vld_file)
429 path, vld_name = os.path.split(vld_file)
430 name, ext = os.path.splitext(vld_name)
432 fasta_map[name] = name
434 fasta_map[name] = os.path.join(genome, name)
438 def extract_eland_sequence(instream, outstream, start, end):
440 Extract a chunk of sequence out of an eland file
442 for line in instream:
443 record = line.split()
445 result = [record[0], record[1][start:end]]
447 result = [record[0][start:end]]
448 outstream.write("\t".join(result))
449 outstream.write(os.linesep)