gaworkflow/util/makebed.py

   1 """
   2 Utility functions to make bedfiles.
   3 """
   4 import os
   5
   6 def make_bed_from_eland_stream(instream, outstream, name, description, chromosome_prefix='chr'):
   7   """
   8   read an eland result file from instream and write a bedfile to outstream
   9   """
  10   # indexes into fields in eland_result.txt file
  11   SEQ = 1
  12   CHR = 6
  13   START = 7
  14   SENSE = 8
  15   # map eland_result.txt sense
  16   sense_map = { 'F': '+', 'R': '-'}
  17   sense_color = { 'F': '0,0,255', 'R': '255,255,0' }
  18   # provide default track names
  19   if name is None: name = "track"
  20   if description is None: description = "eland result file"
  21   bed_header = 'track name="%s" description="%s" visibility=4 itemRgb="ON"'
  22   bed_header += os.linesep
  23   outstream.write(bed_header % (name, description))
  24
  25   for line in instream:
  26     fields = line.split()
  27     # we need more than the CHR field, and it needs to match a chromosome
  28     if len(fields) <= CHR or \
  29           (chromosome_prefix is not None and \
  30              fields[CHR][:3] != chromosome_prefix):
  31       continue
  32     start = fields[START]
  33     stop = int(start) + len(fields[SEQ])
  34     chromosome, extension = fields[CHR].split('.')
  35     assert extension == "fa"
  36     outstream.write('%s %s %d read 0 %s - - %s%s' % (
  37       chromosome,
  38       start,
  39       stop,
  40       sense_map[fields[SENSE]],
  41       sense_color[fields[SENSE]],
  42       os.linesep
  43     ))
  44
  45 def make_description(database, flowcell_id, lane):
  46     """
  47     compute a bedfile name and description from the fctracker database
  48     """
  49     from gaworkflow.util.fctracker import fctracker
  50
  51     fc = fctracker(database)
  52     cells = fc._get_flowcells("where flowcell_id='%s'" % (flowcell_id))
  53     if len(cells) != 1:
  54       raise RuntimeError("couldn't find flowcell id %s" % (flowcell_id))
  55     lane = int(lane)
  56     if lane < 1 or lane > 8:
  57       raise RuntimeError("flowcells only have lanes 1-8")
  58
  59     name = "%s-%s" % (flowcell_id, lane)
  60
  61     cell_id, cell = cells.items()[0]
  62     assert cell_id == flowcell_id
  63
  64     cell_library_id = cell['lane_%d_library_id' %(lane,)]
  65     cell_library = cell['lane_%d_library' %(lane,)]
  66     description = "%s-%s" % (cell_library['library_name'], cell_library_id)
  67     return name, description