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 fields[CHR][:3] != chromosome_prefix:
  29       continue
  30     start = fields[START]
  31     stop = int(start) + len(fields[SEQ])
  32     chromosome, extension = fields[CHR].split('.')
  33     assert extension == "fa"
  34     outstream.write('%s %s %d read 0 %s - - %s%s' % (
  35       chromosome,
  36       start,
  37       stop,
  38       sense_map[fields[SENSE]],
  39       sense_color[fields[SENSE]],
  40       os.linesep
  41     ))
  42
  43 def make_description(database, flowcell_id, lane):
  44     """
  45     compute a bedfile name and description from the fctracker database
  46     """
  47     from gaworkflow.util.fctracker import fctracker
  48
  49     fc = fctracker(database)
  50     cells = fc._get_flowcells("where flowcell_id='%s'" % (flowcell_id))
  51     if len(cells) != 1:
  52       raise RuntimeError("couldn't find flowcell id %s" % (flowcell_id))
  53     lane = int(lane)
  54     if lane < 1 or lane > 8:
  55       raise RuntimeError("flowcells only have lanes 1-8")
  56
  57     name = "%s-%s" % (flowcell_id, lane)
  58
  59     cell_id, cell = cells.items()[0]
  60     assert cell_id == flowcell_id
  61
  62     cell_library_id = cell['lane_%d_library_id' %(lane,)]
  63     cell_library = cell['lane_%d_library' %(lane,)]
  64     description = "%s-%s" % (cell_library['library_name'], cell_library_id)
  65     return name, description