Added last-dotplot options to show sequence lengths.
authorMartin C. Frith
Tue Apr 04 12:15:35 2017 +0900 (2017-04-04)
changeset 8461f46ab956351
parent 845 16060c00b129
child 847 3ca7aa4b27a8
Added last-dotplot options to show sequence lengths.
doc/last-dotplot.txt
scripts/last-dotplot
     1.1 --- a/doc/last-dotplot.txt	Tue Apr 04 11:51:15 2017 +0900
     1.2 +++ b/doc/last-dotplot.txt	Tue Apr 04 12:15:35 2017 +0900
     1.3 @@ -79,6 +79,10 @@
     1.4        Trim unaligned sequence flanks from the 1st genome.
     1.5    --trim2
     1.6        Trim unaligned sequence flanks from the 2nd genome.
     1.7 +  --lengths1
     1.8 +      Show sequence lengths for the 1st genome.
     1.9 +  --lengths2
    1.10 +      Show sequence lengths for the 2nd genome.
    1.11    --bed1=FILE
    1.12        Read `BED-format
    1.13        <https://genome.ucsc.edu/FAQ/FAQformat.html#format1>`_
     2.1 --- a/scripts/last-dotplot	Tue Apr 04 11:51:15 2017 +0900
     2.2 +++ b/scripts/last-dotplot	Tue Apr 04 12:15:35 2017 +0900
     2.3 @@ -163,14 +163,28 @@
     2.4      draw = ImageDraw.Draw(im)
     2.5      return [draw.textsize(i, font=font) for i in my_strings]
     2.6  
     2.7 -def get_seq_info(seqLimits, font, fontsize, image_mode):
     2.8 +def sizeText(size):
     2.9 +    suffixes = "bp", "kb", "Mb", "Gb"
    2.10 +    for i, x in enumerate(suffixes):
    2.11 +        j = 10 ** (i * 3)
    2.12 +        if size < j * 10:
    2.13 +            return "%.2g" % (1.0 * size / j) + x
    2.14 +        if size < j * 1000 or i == len(suffixes) - 1:
    2.15 +            return "%.0f" % (1.0 * size / j) + x
    2.16 +
    2.17 +def seqNameAndSizeText(seqName, seqSize):
    2.18 +    return seqName + ": " + sizeText(seqSize)
    2.19 +
    2.20 +def getSeqInfo(seqLimits, font, fontsize, image_mode, isShowSize):
    2.21      '''Return miscellaneous information about the sequences.'''
    2.22      seqNames = seqLimits.keys()
    2.23      seqNames.sort(key=natural_sort_key)
    2.24      seq_sizes = [seqLimits[i][1] - seqLimits[i][0] for i in seqNames]
    2.25 -    name_sizes = get_text_sizes(seqNames, font, fontsize, image_mode)
    2.26 -    margin = max(zip(*name_sizes)[1])  # maximum text height
    2.27 -    return seqNames, seq_sizes, name_sizes, margin
    2.28 +    if isShowSize: seqLabels = map(seqNameAndSizeText, seqNames, seq_sizes)
    2.29 +    else:          seqLabels = seqNames
    2.30 +    labelSizes = get_text_sizes(seqLabels, font, fontsize, image_mode)
    2.31 +    margin = max(zip(*labelSizes)[1])  # maximum text height
    2.32 +    return seqNames, seq_sizes, seqLabels, labelSizes, margin
    2.33  
    2.34  def div_ceil(x, y):
    2.35      '''Return x / y rounded up.'''
    2.36 @@ -382,10 +396,11 @@
    2.37  
    2.38      if not alignments: raise Exception("there are no alignments")
    2.39  
    2.40 -    seq_info1 = get_seq_info(seqLimits1, font, opts.fontsize, image_mode)
    2.41 -    seq_info2 = get_seq_info(seqLimits2, font, opts.fontsize, image_mode)
    2.42 -    seqNames1, seq_sizes1, name_sizes1, margin1 = seq_info1
    2.43 -    seqNames2, seq_sizes2, name_sizes2, margin2 = seq_info2
    2.44 +    i1 = getSeqInfo(seqLimits1, font, opts.fontsize, image_mode, opts.lengths1)
    2.45 +    seqNames1, seq_sizes1, seqLabels1, labelSizes1, margin1 = i1
    2.46 +
    2.47 +    i2 = getSeqInfo(seqLimits2, font, opts.fontsize, image_mode, opts.lengths2)
    2.48 +    seqNames2, seq_sizes2, seqLabels2, labelSizes2, margin2 = i2
    2.49  
    2.50      warn("choosing bp per pixel...")
    2.51      pix_limit1 = opts.width  - margin1
    2.52 @@ -442,9 +457,9 @@
    2.53              elif store_value == 3: im.putpixel(xy, overlap_color)
    2.54  
    2.55      if opts.fontsize != 0:
    2.56 -        axis1 = get_axis_image(seqNames1, name_sizes1, seq_starts1, seq_pix1,
    2.57 +        axis1 = get_axis_image(seqLabels1, labelSizes1, seq_starts1, seq_pix1,
    2.58                                 font, image_mode, opts)
    2.59 -        axis2 = get_axis_image(seqNames2, name_sizes2, seq_starts2, seq_pix2,
    2.60 +        axis2 = get_axis_image(seqLabels2, labelSizes2, seq_starts2, seq_pix2,
    2.61                                 font, image_mode, opts)
    2.62          axis2 = axis2.transpose(Image.ROTATE_270)  # !!! bug hotspot
    2.63          im.paste(axis1, (0, 0))
    2.64 @@ -490,10 +505,15 @@
    2.65                    help="trim unaligned sequence flanks from the 1st genome")
    2.66      op.add_option("--trim2", action="store_true",
    2.67                    help="trim unaligned sequence flanks from the 2nd genome")
    2.68 +    op.add_option("--lengths1", action="store_true",
    2.69 +                  help="show sequence lengths for the 1st genome")
    2.70 +    op.add_option("--lengths2", action="store_true",
    2.71 +                  help="show sequence lengths for the 2nd genome")
    2.72      op.add_option("--bed1", metavar="FILE",
    2.73                    help="read genome1 annotations from bed file")
    2.74      op.add_option("--bed2", metavar="FILE",
    2.75                    help="read genome2 annotations from bed file")
    2.76 +
    2.77      og = optparse.OptionGroup(op, "Unsequenced gap options")
    2.78      og.add_option("--gap1", metavar="FILE",
    2.79                    help="read genome1 unsequenced gaps from agp or gap file")