gaworkflow/util/makebed.py

   1 """
   2 Utility functions to make bedfiles.
   3 """
   4 import os
   5 import re
   6
   7 # map eland_result.txt sense
   8 sense_map = { 'F': '+', 'R': '-'}
   9 sense_color = { 'F': '0,0,255', 'R': '255,255,0' }
  10
  11 def write_bed_header(outstream, name, description):
  12   """
  13   Produce the headerline for a bedfile
  14   """
  15   # provide default track names
  16   if name is None: name = "track"
  17   if description is None: description = "eland result file"
  18   bed_header = 'track name="%s" description="%s" visibility=4 itemRgb="ON"'
  19   bed_header += os.linesep
  20   outstream.write(bed_header % (name, description))
  21
  22 def make_bed_from_eland_stream(instream, outstream, name, description, chromosome_prefix='chr'):
  23   """
  24   read an eland result file from instream and write a bedfile to outstream
  25   """
  26   # indexes into fields in eland_result.txt file
  27   SEQ = 1
  28   CHR = 6
  29   START = 7
  30   SENSE = 8
  31
  32   write_bed_header(outstream, name, description)
  33
  34   for line in instream:
  35     fields = line.split()
  36     # we need more than the CHR field, and it needs to match a chromosome
  37     if len(fields) <= CHR or \
  38           (chromosome_prefix is not None and \
  39              fields[CHR][:3] != chromosome_prefix):
  40       continue
  41     start = fields[START]
  42     stop = int(start) + len(fields[SEQ])
  43     chromosome, extension = fields[CHR].split('.')
  44     assert extension == "fa"
  45     outstream.write('%s %s %d read 0 %s - - %s%s' % (
  46       chromosome,
  47       start,
  48       stop,
  49       sense_map[fields[SENSE]],
  50       sense_color[fields[SENSE]],
  51       os.linesep
  52     ))
  53
  54
  55 def make_bed_from_multi_eland_stream(
  56   instream,
  57   outstream,
  58   name,
  59   description,
  60   chr_prefix='chr',
  61   max_reads=255
  62   ):
  63   """
  64   read a multi eland stream and write a bedfile
  65   """
  66   write_bed_header(outstream, name, description)
  67   parse_multi_eland(instream, outstream, chr_prefix, max_reads)
  68
  69 def parse_multi_eland(instream, outstream, chr_prefix, max_reads=255):
  70
  71   loc_pattern = '(?P<fullloc>(?P<start>[0-9]+)(?P<dir>[FR])(?P<count>[0-9]+))'
  72   other_pattern = '(?P<chr>[^:,]+)'
  73   split_re = re.compile('(%s|%s)' % (loc_pattern, other_pattern))
  74
  75   for line in instream:
  76     rec = line.split()
  77     if len(rec) > 3:
  78       # colony_id = rec[0]
  79       seq = rec[1]
  80       # number of matches for 0, 1, and 2 mismatches
  81       # m0, m1, m2 = [int(x) for x in rec[2].split(':')]
  82       compressed_reads = rec[3]
  83       cur_chr = ""
  84       reads = {0: [], 1: [], 2:[]}
  85
  86       for token in split_re.finditer(compressed_reads):
  87         if token.group('chr') is not None:
  88           cur_chr =  token.group('chr')[:-3] # strip off .fa
  89         elif token.group('fullloc') is not None:
  90           matches = int(token.group('count'))
  91           # only emit a bed line if
  92           #  our current chromosome starts with chromosome pattern
  93           if chr_prefix is None or cur_chr.startswith(chr_prefix):
  94             start = int(token.group('start'))
  95             stop = start + len(seq)
  96             orientation = token.group('dir')
  97             strand = sense_map[orientation]
  98             color = sense_color[orientation]
  99             # build up list of reads for this record
 100             reads[matches].append((cur_chr, start, stop, strand, color))
 101
 102       # report up to our max_read threshold reporting the fewer-mismatch
 103       # matches first
 104       reported_reads = 0
 105       keys = [0,1,2]
 106       for mismatch, read_list in ((k, reads[k]) for k in keys):
 107         reported_reads += len(read_list)
 108         if reported_reads <= max_reads:
 109           for cur_chr, start, stop, strand, color in read_list:
 110             reported_reads += 1
 111             outstream.write('%s %d %d read 0 %s - - %s%s' % (
 112                 cur_chr,
 113                 start,
 114                 stop,
 115                 sense_map[orientation],
 116                 sense_color[orientation],
 117                 os.linesep
 118             ))
 119
 120 def make_description(database, flowcell_id, lane):
 121     """
 122     compute a bedfile name and description from the fctracker database
 123     """
 124     from gaworkflow.util.fctracker import fctracker
 125
 126     fc = fctracker(database)
 127     cells = fc._get_flowcells("where flowcell_id='%s'" % (flowcell_id))
 128     if len(cells) != 1:
 129       raise RuntimeError("couldn't find flowcell id %s" % (flowcell_id))
 130     lane = int(lane)
 131     if lane < 1 or lane > 8:
 132       raise RuntimeError("flowcells only have lanes 1-8")
 133
 134     name = "%s-%s" % (flowcell_id, lane)
 135
 136     cell_id, cell = cells.items()[0]
 137     assert cell_id == flowcell_id
 138
 139     cell_library_id = cell['lane_%d_library_id' %(lane,)]
 140     cell_library = cell['lane_%d_library' %(lane,)]
 141     description = "%s-%s" % (cell_library['library_name'], cell_library_id)
 142     return name, description