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
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_elements(self):
198 lane = ElementTree.Element(ElandLane.LANE,
200 unicode(ElandLane.XML_VERSION)})
201 sample_tag = ElementTree.SubElement(lane, ElandLane.SAMPLE_NAME)
202 sample_tag.text = self.sample_name
203 lane_tag = ElementTree.SubElement(lane, ElandLane.LANE_ID)
204 lane_tag.text = str(self.lane_id)
205 if self.end is not None:
206 end_tag = ElementTree.SubElement(lane, ElandLane.END)
207 end_tag.text = str(self.end)
208 genome_map = ElementTree.SubElement(lane, ElandLane.GENOME_MAP)
209 for k, v in self.genome_map.items():
210 item = ElementTree.SubElement(
211 genome_map, ElandLane.GENOME_ITEM,
212 {'name':k, 'value':unicode(v)})
213 mapped_reads = ElementTree.SubElement(lane, ElandLane.MAPPED_READS)
214 for k, v in self.mapped_reads.items():
215 item = ElementTree.SubElement(
216 mapped_reads, ElandLane.MAPPED_ITEM,
217 {'name':k, 'value':unicode(v)})
218 match_codes = ElementTree.SubElement(lane, ElandLane.MATCH_CODES)
219 for k, v in self.match_codes.items():
220 item = ElementTree.SubElement(
221 match_codes, ElandLane.MATCH_ITEM,
222 {'name':k, 'value':unicode(v)})
223 reads = ElementTree.SubElement(lane, ElandLane.READS)
224 reads.text = unicode(self.reads)
228 def set_elements(self, tree):
229 if tree.tag != ElandLane.LANE:
230 raise ValueError('Exptecting %s' % (ElandLane.LANE,))
233 self._mapped_reads = {}
234 self._match_codes = {}
237 tag = element.tag.lower()
238 if tag == ElandLane.SAMPLE_NAME.lower():
239 self._sample_name = element.text
240 elif tag == ElandLane.LANE_ID.lower():
241 self.lane_id = int(element.text)
242 elif tag == ElandLane.END.lower():
243 self.end = int(element.text)
244 elif tag == ElandLane.GENOME_MAP.lower():
245 for child in element:
246 name = child.attrib['name']
247 value = child.attrib['value']
248 self.genome_map[name] = value
249 elif tag == ElandLane.MAPPED_READS.lower():
250 for child in element:
251 name = child.attrib['name']
252 value = child.attrib['value']
253 self._mapped_reads[name] = int(value)
254 elif tag == ElandLane.MATCH_CODES.lower():
255 for child in element:
256 name = child.attrib['name']
257 value = int(child.attrib['value'])
258 self._match_codes[name] = value
259 elif tag == ElandLane.READS.lower():
260 self._reads = int(element.text)
262 logging.warn("ElandLane unrecognized tag %s" % (element.tag,))
266 Summarize information from eland files
270 ELAND = 'ElandCollection'
275 def __init__(self, xml=None):
276 # we need information from the gerald config.xml
277 self.results = [{},{}]
280 self.set_elements(xml)
282 def get_elements(self):
283 root = ElementTree.Element(ELAND.ELAND,
284 {'version': unicode(ELAND.XML_VERSION)})
285 for end in range(len(self.results)):
286 end_results = self.results[end]
287 for lane_id, lane in end_results.items():
288 eland_lane = lane.get_elements()
289 eland_lane.attrib[ELAND.END] = unicode (end)
290 eland_lane.attrib[ELAND.LANE_ID] = unicode(lane_id)
291 root.append(eland_lane)
294 def set_elements(self, tree):
295 if tree.tag.lower() != ELAND.ELAND.lower():
296 raise ValueError('Expecting %s', ELAND.ELAND)
297 for element in list(tree):
298 lane_id = int(element.attrib[ELAND.LANE_ID])
299 end = int(element.attrib.get(ELAND.END, 0))
300 lane = ElandLane(xml=element)
301 self.results[end][lane_id] = lane
303 def check_for_eland_file(basedir, pattern, lane_id, end):
305 full_lane_id = lane_id
307 full_lane_id = "%d_%d" % ( lane_id, end )
309 basename = pattern % (full_lane_id,)
310 pathname = os.path.join(basedir, basename)
311 if os.path.exists(pathname):
316 def eland(basedir, gerald=None, genome_maps=None):
319 #file_list = glob(os.path.join(basedir, "*_eland_result.txt"))
320 #if len(file_list) == 0:
321 # # lets handle compressed eland files too
322 # file_list = glob(os.path.join(basedir, "*_eland_result.txt.bz2"))
324 lane_ids = range(1,9)
327 # the order in patterns determines the preference for what
329 patterns = ['s_%s_eland_result.txt',
330 's_%s_eland_result.txt.bz2',
331 's_%s_eland_result.txt.gz',
332 's_%s_eland_extended.txt',
333 's_%s_eland_extended.txt.bz2',
334 's_%s_eland_extended.txt.gz',
335 's_%s_eland_multi.txt',
336 's_%s_eland_multi.txt.bz2',
337 's_%s_eland_multi.txt.gz',]
340 for lane_id in lane_ids:
342 pathname = check_for_eland_file(basedir, p, lane_id, end)
343 if pathname is not None:
347 # yes the lane_id is also being computed in ElandLane._update
348 # I didn't want to clutter up my constructor
349 # but I needed to persist the sample_name/lane_id for
350 # runfolder summary_report
351 path, name = os.path.split(pathname)
352 logging.info("Adding eland file %s" %(name,))
353 # split_name = name.split('_')
354 # lane_id = int(split_name[1])
356 if genome_maps is not None:
357 genome_map = genome_maps[lane_id]
358 elif gerald is not None:
359 genome_dir = gerald.lanes[lane_id].eland_genome
360 genome_map = build_genome_fasta_map(genome_dir)
364 eland_result = ElandLane(pathname, lane_id, end, genome_map)
368 effective_end = end - 1
369 e.results[effective_end][lane_id] = eland_result
372 def build_genome_fasta_map(genome_dir):
373 # build fasta to fasta file map
374 logging.info("Building genome map")
375 genome = genome_dir.split(os.path.sep)[-1]
377 for vld_file in glob(os.path.join(genome_dir, '*.vld')):
379 if os.path.islink(vld_file):
381 vld_file = os.path.realpath(vld_file)
382 path, vld_name = os.path.split(vld_file)
383 name, ext = os.path.splitext(vld_name)
385 fasta_map[name] = name
387 fasta_map[name] = os.path.join(genome, name)
391 def extract_eland_sequence(instream, outstream, start, end):
393 Extract a chunk of sequence out of an eland file
395 for line in instream:
396 record = line.split()
398 result = [record[0], record[1][start:end]]
400 result = [record[0][start:end]]
401 outstream.write("\t".join(result))
402 outstream.write(os.linesep)