last-dotplot: start implementing "layers"
authorMartin C. Frith
Thu May 18 14:01:30 2017 +0900 (2017-05-18)
changeset 8575e4e3f2cf68b
parent 856 d086361eabec
child 858 93d47bb8527f
last-dotplot: start implementing "layers"
scripts/last-dotplot
     1.1 --- a/scripts/last-dotplot	Thu May 18 13:21:30 2017 +0900
     1.2 +++ b/scripts/last-dotplot	Thu May 18 14:01:30 2017 +0900
     1.3 @@ -296,6 +296,7 @@
     1.4          b = max(beg, cropBeg)
     1.5          e = min(end, cropEnd)
     1.6          if b >= e: continue
     1.7 +        layer = 900
     1.8          color = "#ffe4ff"
     1.9          if len(w) > 5:
    1.10              if len(w) > 8 and w[8].count(",") == 2:
    1.11 @@ -304,12 +305,12 @@
    1.12                  color = "#fff4f4"
    1.13              elif w[5] == "-":
    1.14                  color = "#f4f4ff"
    1.15 -        yield seqName, b, e, color
    1.16 +        yield layer, color, seqName, b, e
    1.17  
    1.18  def isExtraFirstGapField(fields):
    1.19      return fields[4].isdigit()
    1.20  
    1.21 -def readGaps(fileName, seqLimits):
    1.22 +def readGaps(opts, fileName, seqLimits):
    1.23      '''Read locations of unsequenced gaps, from an agp or gap file.'''
    1.24      if not fileName: return
    1.25      for line in myOpen(fileName):
    1.26 @@ -325,30 +326,32 @@
    1.27          b = max(beg, cropBeg)
    1.28          e = min(end, cropEnd)
    1.29          if b >= e: continue
    1.30 -        bridgedText = w[7]
    1.31 -        yield seqName, b, e, bridgedText
    1.32 +        if w[7] == "yes":
    1.33 +            yield 3000, opts.bridged_color, seqName, b, e
    1.34 +        else:
    1.35 +            yield 2000, opts.unbridged_color, seqName, b, e
    1.36  
    1.37 -def drawAnnotations(im, beds, origins, margin, limit, isTop, bp_per_pix):
    1.38 -    # XXX no consideration of different-color overlaps
    1.39 -    for seqName, beg, end, color in beds:
    1.40 +def bedBoxes(beds, origins, margin, edge, isTop, bpPerPix):
    1.41 +    for layer, color, seqName, beg, end in beds:
    1.42          ori = origins[seqName]
    1.43 -        b = (ori + beg) // bp_per_pix
    1.44 -        e = div_ceil(ori + end, bp_per_pix)
    1.45 -        if isTop: box = b, margin, e, limit
    1.46 -        else:     box = margin, b, limit, e
    1.47 -        im.paste(color, box)
    1.48 +        if layer <= 1000:
    1.49 +            # include partly-covered pixels
    1.50 +            b = (ori + beg) // bpPerPix
    1.51 +            e = div_ceil(ori + end, bpPerPix)
    1.52 +        else:
    1.53 +            # exclude partly-covered pixels
    1.54 +            b = div_ceil(ori + beg, bpPerPix)
    1.55 +            e = (ori + end) // bpPerPix
    1.56 +            if e <= b: continue
    1.57 +        if isTop:
    1.58 +            box = b, margin, e, edge
    1.59 +        else:
    1.60 +            box = margin, b, edge, e
    1.61 +        yield layer, color, box
    1.62  
    1.63 -def drawUnsequencedGaps(im, gaps, origins, margin, limit, isTop, bridgedText,
    1.64 -                        bp_per_pix, color):
    1.65 -    '''Draw rectangles representing unsequenced gaps.'''
    1.66 -    for seqName, beg, end, b in gaps:
    1.67 -        if b != bridgedText: continue
    1.68 -        ori = origins[seqName]
    1.69 -        b = div_ceil(ori + beg, bp_per_pix)  # use fully-covered pixels only
    1.70 -        e = (ori + end) // bp_per_pix
    1.71 -        if e <= b: continue
    1.72 -        if isTop: box = b, margin, e, limit
    1.73 -        else:     box = margin, b, limit, e
    1.74 +def drawAnnotations(im, boxes):
    1.75 +    # xxx use partial transparency for different-color overlaps?
    1.76 +    for layer, color, box in boxes:
    1.77          im.paste(color, box)
    1.78  
    1.79  def make_label(text, text_size, range_start, range_size):
    1.80 @@ -445,22 +448,16 @@
    1.81      origins1 = expandedSeqDict(origins1)
    1.82      origins2 = expandedSeqDict(origins2)
    1.83  
    1.84 -    beds1 = list(readBed(opts.bed1, seqLimits1))
    1.85 -    beds2 = list(readBed(opts.bed2, seqLimits2))
    1.86 -    drawAnnotations(im, beds1, origins1, margin2, height, True, bpPerPix)
    1.87 -    drawAnnotations(im, beds2, origins2, margin1, width, False, bpPerPix)
    1.88 +    beds1 = itertools.chain(readBed(opts.bed1, seqLimits1),
    1.89 +                            readGaps(opts, opts.gap1, seqLimits1))
    1.90 +    b1 = bedBoxes(beds1, origins1, margin2, height, True, bpPerPix)
    1.91  
    1.92 -    gaps1 = list(readGaps(opts.gap1, seqLimits1))
    1.93 -    gaps2 = list(readGaps(opts.gap2, seqLimits2))
    1.94 -    # draw bridged gaps first, then unbridged gaps on top:
    1.95 -    drawUnsequencedGaps(im, gaps1, origins1, margin2, height, True, "yes",
    1.96 -                        bpPerPix, opts.bridged_color)
    1.97 -    drawUnsequencedGaps(im, gaps2, origins2, margin1, width, False, "yes",
    1.98 -                        bpPerPix, opts.bridged_color)
    1.99 -    drawUnsequencedGaps(im, gaps1, origins1, margin2, height, True, "no",
   1.100 -                        bpPerPix, opts.unbridged_color)
   1.101 -    drawUnsequencedGaps(im, gaps2, origins2, margin1, width, False, "no",
   1.102 -                        bpPerPix, opts.unbridged_color)
   1.103 +    beds2 = itertools.chain(readBed(opts.bed2, seqLimits2),
   1.104 +                            readGaps(opts, opts.gap2, seqLimits2))
   1.105 +    b2 = bedBoxes(beds2, origins2, margin1, width, False, bpPerPix)
   1.106 +
   1.107 +    boxes = sorted(itertools.chain(b1, b2))
   1.108 +    drawAnnotations(im, boxes)
   1.109  
   1.110      for i in range(height):
   1.111          for j in range(width):