Added unsequenced gap options to last-dotplot.
authorMartin C. Frith
Tue Oct 20 16:26:13 2015 +0900 (2015-10-20)
changeset 650157301c5c6fd
parent 649 94f7cda3507b
child 651 fbb4c3b60998
Added unsequenced gap options to last-dotplot.
scripts/last-dotplot
     1.1 --- a/scripts/last-dotplot	Mon Oct 19 13:54:50 2015 +0900
     1.2 +++ b/scripts/last-dotplot	Tue Oct 20 16:26:13 2015 +0900
     1.3 @@ -192,6 +192,46 @@
     1.4                                  beg1, beg2, size)
     1.5      return hits
     1.6  
     1.7 +def expandedSeqDict(seqDict):
     1.8 +    '''Allow lookup by short sequence names, e.g. chr7 as well as hg19.chr7.'''
     1.9 +    newDict = {}
    1.10 +    for name, x in seqDict.items():
    1.11 +        base = name.split(".")[-1]
    1.12 +        newDict[name] = x
    1.13 +        newDict[base] = x
    1.14 +    return newDict
    1.15 +
    1.16 +def isExtraFirstGapField(fields):
    1.17 +    return fields[4].isdigit()
    1.18 +
    1.19 +def readGaps(fileName):
    1.20 +    '''Read locations of unsequenced gaps, from an agp or gap file.'''
    1.21 +    if not fileName: return
    1.22 +    for line in fileinput.input(fileName):
    1.23 +        w = line.split()
    1.24 +        if not w or w[0][0] == "#": continue
    1.25 +        if isExtraFirstGapField(w): w = w[1:]
    1.26 +        if w[4] not in "NU": continue
    1.27 +        seqName = w[0]
    1.28 +        end = int(w[2])
    1.29 +        beg = end - int(w[5])  # zero-based coordinate
    1.30 +        bridgedText = w[7]
    1.31 +        yield seqName, beg, end, bridgedText
    1.32 +
    1.33 +def drawUnsequencedGaps(im, gaps, start_dic, margin, limit, isTop, bridgedText,
    1.34 +                        bp_per_pix, color):
    1.35 +    '''Draw rectangles representing unsequenced gaps.'''
    1.36 +    for seqName, beg, end, b in gaps:
    1.37 +        if b != bridgedText: continue
    1.38 +        if seqName not in start_dic: continue
    1.39 +        origin = start_dic[seqName]
    1.40 +        b = div_ceil(beg, bp_per_pix)  # use fully-covered pixels only
    1.41 +        e = end // bp_per_pix
    1.42 +        if e <= b: continue
    1.43 +        if isTop: box = origin + b, margin, origin + e, limit
    1.44 +        else:     box = margin, origin + b, limit, origin + e
    1.45 +        im.paste(color, box)
    1.46 +
    1.47  def make_label(text, text_size, range_start, range_size):
    1.48      '''Return an axis label with endpoint & sort-order information.'''
    1.49      text_width  = text_size[0]
    1.50 @@ -273,6 +313,20 @@
    1.51      image_size = width, height
    1.52      im = Image.new(image_mode, image_size, opts.background_color)
    1.53  
    1.54 +    start_dic1 = expandedSeqDict(seq_start_dic1)
    1.55 +    start_dic2 = expandedSeqDict(seq_start_dic2)
    1.56 +    gaps1 = list(readGaps(opts.gap1))
    1.57 +    gaps2 = list(readGaps(opts.gap2))
    1.58 +    # draw bridged gaps first, then unbridged gaps on top:
    1.59 +    drawUnsequencedGaps(im, gaps1, start_dic1, margin2, height, True, "yes",
    1.60 +                        bp_per_pix, opts.bridged_color)
    1.61 +    drawUnsequencedGaps(im, gaps2, start_dic2, margin1, width, False, "yes",
    1.62 +                        bp_per_pix, opts.bridged_color)
    1.63 +    drawUnsequencedGaps(im, gaps1, start_dic1, margin2, height, True, "no",
    1.64 +                        bp_per_pix, opts.unbridged_color)
    1.65 +    drawUnsequencedGaps(im, gaps2, start_dic2, margin1, width, False, "no",
    1.66 +                        bp_per_pix, opts.unbridged_color)
    1.67 +
    1.68      for i in range(height):
    1.69          for j in range(width):
    1.70              store_value = hits[i * width + j]
    1.71 @@ -320,6 +374,16 @@
    1.72                    help="color for forward alignments (default: %default)")
    1.73      op.add_option("-r", "--reversecolor", metavar="COLOR", default="blue",
    1.74                    help="color for reverse alignments (default: %default)")
    1.75 +    og = optparse.OptionGroup(op, "Unsequenced gap options")
    1.76 +    og.add_option("--gap1", metavar="FILE",
    1.77 +                  help="read genome1 unsequenced gaps from agp or gap file")
    1.78 +    og.add_option("--gap2", metavar="FILE",
    1.79 +                  help="read genome2 unsequenced gaps from agp or gap file")
    1.80 +    og.add_option("--bridged-color", metavar="COLOR", default="yellow",
    1.81 +                  help="color for bridged gaps (default: %default)")
    1.82 +    og.add_option("--unbridged-color", metavar="COLOR", default="pink",
    1.83 +                  help="color for unbridged gaps (default: %default)")
    1.84 +    op.add_option_group(og)
    1.85      (opts, args) = op.parse_args()
    1.86      if len(args) != 2: op.error("2 arguments needed")
    1.87