htsworkflow/pipelines/eland.py

   1 """
   2 Analyze ELAND files
   3 """
   4
   5 from glob import glob
   6 import logging
   7 import os
   8 import re
   9 import stat
  10
  11 from htsworkflow.pipelines.runfolder import ElementTree
  12 from htsworkflow.util.ethelp import indent, flatten
  13 from htsworkflow.util.opener import autoopen
  14
  15 class ElandLane(object):
  16     """
  17     Process an eland result file
  18     """
  19     XML_VERSION = 2
  20     LANE = 'ElandLane'
  21     SAMPLE_NAME = 'SampleName'
  22     LANE_ID = 'LaneID'
  23     END = 'End'
  24     GENOME_MAP = 'GenomeMap'
  25     GENOME_ITEM = 'GenomeItem'
  26     MAPPED_READS = 'MappedReads'
  27     MAPPED_ITEM = 'MappedItem'
  28     MATCH_CODES = 'MatchCodes'
  29     MATCH_ITEM = 'Code'
  30     READS = 'Reads'
  31
  32     ELAND_SINGLE = 0
  33     ELAND_MULTI = 1
  34     ELAND_EXTENDED = 2
  35     ELAND_EXPORT = 3
  36
  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
  41         self.end = end
  42         self._reads = None
  43         self._mapped_reads = None
  44         self._match_codes = None
  45         if genome_map is None:
  46             genome_map = {}
  47         self.genome_map = genome_map
  48         self.eland_type = None
  49
  50         if xml is not None:
  51             self.set_elements(xml)
  52
  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
  65           else:
  66             self.eland_type = ElandLane.ELAND_SINGLE
  67
  68     def _update(self):
  69         """
  70         Actually read the file and actually count the reads
  71         """
  72         # can't do anything if we don't have a file to process
  73         if self.pathname is None:
  74             return
  75         self._guess_eland_type(self.pathname)
  76
  77         if os.stat(self.pathname)[stat.ST_SIZE] == 0:
  78             raise RuntimeError("Eland isn't done, try again later.")
  79
  80         logging.info("summarizing results for %s" % (self.pathname))
  81
  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)
  87         else:
  88           raise NotImplementedError("Only support single/multi/extended eland files")
  89         self._match_codes, self._mapped_reads, self._reads = result
  90
  91     def _update_eland_result(self, pathname):
  92         reads = 0
  93         mapped_reads = {}
  94
  95         match_codes = {'NM':0, 'QC':0, 'RM':0,
  96                        'U0':0, 'U1':0, 'U2':0,
  97                        'R0':0, 'R1':0, 'R2':0,
  98                       }
  99         for line in autoopen(pathname,'r'):
 100             reads += 1
 101             fields = line.split()
 102             # code = fields[2]
 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
 107             if len(fields) < 7:
 108                 continue
 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
 112
 113     def _update_eland_multi(self, pathname):
 114         reads = 0
 115         mapped_reads = {}
 116
 117         match_codes = {'NM':0, 'QC':0, 'RM':0,
 118                        'U0':0, 'U1':0, 'U2':0,
 119                        'R0':0, 'R1':0, 'R2':0,
 120                       }
 121         match_counts_re = re.compile("([\d]+):([\d]+):([\d]+)")
 122         for line in autoopen(pathname,'r'):
 123             reads += 1
 124             fields = line.split()
 125             # fields[2] = QC/NM/or number of matches
 126             groups = match_counts_re.match(fields[2])
 127             if groups is None:
 128                 match_codes[fields[2]] += 1
 129             else:
 130                 # when there are too many hit, eland writes a - where
 131                 # it would have put the list of hits
 132                 if fields[3] == '-':
 133                   continue
 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
 139
 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
 145
 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
 151
 152                 chromo = None
 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]
 158
 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
 163
 164     def _update_name(self):
 165         # extract the sample name
 166         if self.pathname is None:
 167             return
 168
 169         path, name = os.path.split(self.pathname)
 170         split_name = name.split('_')
 171         self._sample_name = split_name[0]
 172
 173     def _get_sample_name(self):
 174         if self._sample_name is None:
 175             self._update_name()
 176         return self._sample_name
 177     sample_name = property(_get_sample_name)
 178
 179     def _get_reads(self):
 180         if self._reads is None:
 181             self._update()
 182         return self._reads
 183     reads = property(_get_reads)
 184
 185     def _get_mapped_reads(self):
 186         if self._mapped_reads is None:
 187             self._update()
 188         return self._mapped_reads
 189     mapped_reads = property(_get_mapped_reads)
 190
 191     def _get_match_codes(self):
 192         if self._match_codes is None:
 193             self._update()
 194         return self._match_codes
 195     match_codes = property(_get_match_codes)
 196
 197     def get_elements(self):
 198         lane = ElementTree.Element(ElandLane.LANE,
 199                                    {'version':
 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)
 225
 226         return lane
 227
 228     def set_elements(self, tree):
 229         if tree.tag != ElandLane.LANE:
 230             raise ValueError('Exptecting %s' % (ElandLane.LANE,))
 231
 232         # reset dictionaries
 233         self._mapped_reads = {}
 234         self._match_codes = {}
 235
 236         for element in tree:
 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)
 261             else:
 262                 logging.warn("ElandLane unrecognized tag %s" % (element.tag,))
 263
 264 class ELAND(object):
 265     """
 266     Summarize information from eland files
 267     """
 268     XML_VERSION = 2
 269
 270     ELAND = 'ElandCollection'
 271     LANE = 'Lane'
 272     LANE_ID = 'id'
 273     END = 'end'
 274
 275     def __init__(self, xml=None):
 276         # we need information from the gerald config.xml
 277         self.results = [{},{}]
 278
 279         if xml is not None:
 280             self.set_elements(xml)
 281
 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)
 292         return root
 293
 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
 302
 303 def check_for_eland_file(basedir, pattern, lane_id, end):
 304    if end is None:
 305       full_lane_id = lane_id
 306    else:
 307       full_lane_id = "%d_%d" % ( lane_id, end )
 308
 309    basename = pattern % (full_lane_id,)
 310    pathname = os.path.join(basedir, basename)
 311    if os.path.exists(pathname):
 312        return pathname
 313    else:
 314        return None
 315
 316 def eland(basedir, gerald=None, genome_maps=None):
 317     e = ELAND()
 318
 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"))
 323
 324     lane_ids = range(1,9)
 325     ends = [None, 1, 2]
 326
 327     # the order in patterns determines the preference for what
 328     # will be found.
 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',]
 338
 339     for end in ends:
 340         for lane_id in lane_ids:
 341             for p in patterns:
 342                 pathname = check_for_eland_file(basedir, p, lane_id, end)
 343                 if pathname is not None:
 344                   break
 345             else:
 346                 continue
 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])
 355
 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)
 361             else:
 362                 genome_map = {}
 363
 364             eland_result = ElandLane(pathname, lane_id, end, genome_map)
 365             if end is None:
 366                 effective_end =  0
 367             else:
 368                 effective_end = end - 1
 369             e.results[effective_end][lane_id] = eland_result
 370     return e
 371
 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]
 376     fasta_map = {}
 377     for vld_file in glob(os.path.join(genome_dir, '*.vld')):
 378         is_link = False
 379         if os.path.islink(vld_file):
 380             is_link = True
 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)
 384         if is_link:
 385             fasta_map[name] = name
 386         else:
 387             fasta_map[name] = os.path.join(genome, name)
 388     return fasta_map
 389
 390
 391 def extract_eland_sequence(instream, outstream, start, end):
 392     """
 393     Extract a chunk of sequence out of an eland file
 394     """
 395     for line in instream:
 396         record = line.split()
 397         if len(record) > 1:
 398             result = [record[0], record[1][start:end]]
 399         else:
 400             result = [record[0][start:end]]
 401         outstream.write("\t".join(result))
 402         outstream.write(os.linesep)