scripts/last-dotplot
changeset 845 16060c00b129
parent 844 23de4eb3be1d
child 846 1f46ab956351
     1.1 --- a/scripts/last-dotplot	Tue Apr 04 10:50:57 2017 +0900
     1.2 +++ b/scripts/last-dotplot	Tue Apr 04 11:51:15 2017 +0900
     1.3 @@ -258,13 +258,35 @@
     1.4          newDict[base] = x
     1.5      return newDict
     1.6  
     1.7 +def readBed(fileName, seqLimits):
     1.8 +    if not fileName: return
     1.9 +    for line in myOpen(fileName):
    1.10 +        w = line.split()
    1.11 +        if not w or w[0][0] == "#": continue
    1.12 +        seqName = w[0]
    1.13 +        if seqName not in seqLimits: continue
    1.14 +        cropBeg, cropEnd = seqLimits[seqName]
    1.15 +        beg = int(w[1])
    1.16 +        end = int(w[2])
    1.17 +        b = max(beg, cropBeg)
    1.18 +        e = min(end, cropEnd)
    1.19 +        if b >= e: continue
    1.20 +        if len(w) > 8:
    1.21 +            color = tuple(map(int, w[8].split(",")))
    1.22 +        elif len(w) > 5:
    1.23 +            if   w[5] == "+": color = 255, 244, 244
    1.24 +            elif w[5] == "-": color = 244, 244, 255
    1.25 +            else:             color = 255, 228, 255
    1.26 +        else:
    1.27 +            color = 255, 228, 255
    1.28 +        yield seqName, b, e, color
    1.29 +
    1.30  def isExtraFirstGapField(fields):
    1.31      return fields[4].isdigit()
    1.32  
    1.33  def readGaps(fileName, seqLimits):
    1.34      '''Read locations of unsequenced gaps, from an agp or gap file.'''
    1.35      if not fileName: return
    1.36 -    seqLimits = expandedSeqDict(seqLimits)
    1.37      for line in myOpen(fileName):
    1.38          w = line.split()
    1.39          if not w or w[0][0] == "#": continue
    1.40 @@ -281,6 +303,16 @@
    1.41          bridgedText = w[7]
    1.42          yield seqName, b, e, bridgedText
    1.43  
    1.44 +def drawAnnotations(im, beds, origins, margin, limit, isTop, bp_per_pix):
    1.45 +    # XXX no consideration of different-color overlaps
    1.46 +    for seqName, beg, end, color in beds:
    1.47 +        ori = origins[seqName]
    1.48 +        b = (ori + beg) // bp_per_pix
    1.49 +        e = div_ceil(ori + end, bp_per_pix)
    1.50 +        if isTop: box = b, margin, e, limit
    1.51 +        else:     box = margin, b, limit, e
    1.52 +        im.paste(color, box)
    1.53 +
    1.54  def drawUnsequencedGaps(im, gaps, origins, margin, limit, isTop, bridgedText,
    1.55                          bp_per_pix, color):
    1.56      '''Draw rectangles representing unsequenced gaps.'''
    1.57 @@ -379,8 +411,16 @@
    1.58      image_size = width, height
    1.59      im = Image.new(image_mode, image_size, opts.background_color)
    1.60  
    1.61 +    seqLimits1 = expandedSeqDict(seqLimits1)
    1.62 +    seqLimits2 = expandedSeqDict(seqLimits2)
    1.63      origins1 = expandedSeqDict(origins1)
    1.64      origins2 = expandedSeqDict(origins2)
    1.65 +
    1.66 +    beds1 = list(readBed(opts.bed1, seqLimits1))
    1.67 +    beds2 = list(readBed(opts.bed2, seqLimits2))
    1.68 +    drawAnnotations(im, beds1, origins1, margin2, height, True, bp_per_pix)
    1.69 +    drawAnnotations(im, beds2, origins2, margin1, width, False, bp_per_pix)
    1.70 +
    1.71      gaps1 = list(readGaps(opts.gap1, seqLimits1))
    1.72      gaps2 = list(readGaps(opts.gap2, seqLimits2))
    1.73      # draw bridged gaps first, then unbridged gaps on top:
    1.74 @@ -450,6 +490,10 @@
    1.75                    help="trim unaligned sequence flanks from the 1st genome")
    1.76      op.add_option("--trim2", action="store_true",
    1.77                    help="trim unaligned sequence flanks from the 2nd genome")
    1.78 +    op.add_option("--bed1", metavar="FILE",
    1.79 +                  help="read genome1 annotations from bed file")
    1.80 +    op.add_option("--bed2", metavar="FILE",
    1.81 +                  help="read genome2 annotations from bed file")
    1.82      og = optparse.OptionGroup(op, "Unsequenced gap options")
    1.83      og.add_option("--gap1", metavar="FILE",
    1.84                    help="read genome1 unsequenced gaps from agp or gap file")