Refactor last-dotplot
authorMartin C. Frith
Thu Nov 02 20:04:51 2017 +0900 (2017-11-02)
changeset 897f7bc7d099e2b
parent 896 2238466f36bd
child 898 f6a9c15287ea
Refactor last-dotplot
scripts/last-dotplot
     1.1 --- a/scripts/last-dotplot	Thu Nov 02 19:44:23 2017 +0900
     1.2 +++ b/scripts/last-dotplot	Thu Nov 02 20:04:51 2017 +0900
     1.3 @@ -103,7 +103,7 @@
     1.4                  yield chr1, seqlen1, chr2, seqlen2, blocks
     1.5                  mafCount = 0
     1.6  
     1.7 -def seqRangeFromText(text):
     1.8 +def seqRequestFromText(text):
     1.9      if ":" in text:
    1.10          pattern, interval = text.rsplit(":", 1)
    1.11          if "-" in interval:
    1.12 @@ -111,10 +111,10 @@
    1.13              return pattern, int(beg), int(end)  # beg may be negative
    1.14      return text, 0, sys.maxsize
    1.15  
    1.16 -def rangeFromSeqName(seqRanges, name, seqLen):
    1.17 -    if not seqRanges: return 0, seqLen
    1.18 +def rangeFromSeqName(seqRequests, name, seqLen):
    1.19 +    if not seqRequests: return 0, seqLen
    1.20      base = name.split(".")[-1]  # allow for names like hg19.chr7
    1.21 -    for pat, beg, end in seqRanges:
    1.22 +    for pat, beg, end in seqRequests:
    1.23          if fnmatchcase(name, pat) or fnmatchcase(base, pat):
    1.24              return max(beg, 0), min(end, seqLen)
    1.25      return None
    1.26 @@ -136,8 +136,8 @@
    1.27  
    1.28  def readAlignments(fileName, opts):
    1.29      '''Get alignments and sequence limits, from MAF or tabular format.'''
    1.30 -    seqRanges1 = map(seqRangeFromText, opts.seq1)
    1.31 -    seqRanges2 = map(seqRangeFromText, opts.seq2)
    1.32 +    seqRequests1 = map(seqRequestFromText, opts.seq1)
    1.33 +    seqRequests2 = map(seqRequestFromText, opts.seq2)
    1.34  
    1.35      alignments = []
    1.36      seqNames1 = []
    1.37 @@ -146,9 +146,9 @@
    1.38      seqLimits2 = {}
    1.39      lines = myOpen(fileName)
    1.40      for seqName1, seqLen1, seqName2, seqLen2, blocks in alignmentInput(lines):
    1.41 -        range1 = rangeFromSeqName(seqRanges1, seqName1, seqLen1)
    1.42 +        range1 = rangeFromSeqName(seqRequests1, seqName1, seqLen1)
    1.43          if not range1: continue
    1.44 -        range2 = rangeFromSeqName(seqRanges2, seqName2, seqLen2)
    1.45 +        range2 = rangeFromSeqName(seqRequests2, seqName2, seqLen2)
    1.46          if not range2: continue
    1.47          b = list(croppedBlocks(list(blocks), range1, range2))
    1.48          if not b: continue
    1.49 @@ -391,10 +391,10 @@
    1.50              yield 2000, opts.unbridged_color, seqName, beg, end
    1.51  
    1.52  def bedBoxes(beds, seqLimits, origins, margin, edge, isTop, bpPerPix):
    1.53 -    for layer, color, seqName, beg, end in beds:
    1.54 +    for layer, color, seqName, bedBeg, bedEnd in beds:
    1.55          cropBeg, cropEnd = seqLimits[seqName]
    1.56 -        beg = max(beg, cropBeg)
    1.57 -        end = min(end, cropEnd)
    1.58 +        beg = max(bedBeg, cropBeg)
    1.59 +        end = min(bedEnd, cropEnd)
    1.60          if beg >= end: continue
    1.61          ori = origins[seqName]
    1.62          if layer <= 1000:
    1.63 @@ -425,27 +425,28 @@
    1.64      sort_key    = text_width - range_size
    1.65      return sort_key, label_start, label_end, text, text_height
    1.66  
    1.67 -def get_nonoverlapping_labels(labels, label_space):
    1.68 +def nonoverlappingLabels(labels, minPixTweenLabels):
    1.69      '''Get a subset of non-overlapping axis labels, greedily.'''
    1.70 -    nonoverlapping_labels = []
    1.71 +    out = []
    1.72      for i in labels:
    1.73 -        if True not in [i[1] < j[2] + label_space and j[1] < i[2] + label_space
    1.74 -                        for j in nonoverlapping_labels]:
    1.75 -            nonoverlapping_labels.append(i)
    1.76 -    return nonoverlapping_labels
    1.77 +        beg = i[1] - minPixTweenLabels
    1.78 +        end = i[2] + minPixTweenLabels
    1.79 +        if all(j[2] <= beg or j[1] >= end for j in out):
    1.80 +            out.append(i)
    1.81 +    return out
    1.82  
    1.83  def axisImage(seqNames, name_sizes, seq_starts, seq_pix, textRot,
    1.84                textAln, font, image_mode, opts):
    1.85      '''Make an image of axis labels.'''
    1.86 -    min_pos = seq_starts[0]
    1.87 -    max_pos = seq_starts[-1] + seq_pix[-1]
    1.88 +    beg = seq_starts[0]
    1.89 +    end = seq_starts[-1] + seq_pix[-1]
    1.90      margin = max(i[1] for i in name_sizes)
    1.91      labels = map(make_label, seqNames, name_sizes, seq_starts, seq_pix)
    1.92 -    labels = [i for i in labels if i[1] >= min_pos and i[2] <= max_pos]
    1.93 +    labels = [i for i in labels if i[1] >= beg and i[2] <= end]
    1.94      labels.sort()
    1.95      minPixTweenLabels = 0 if textRot else opts.label_space
    1.96 -    labels = get_nonoverlapping_labels(labels, minPixTweenLabels)
    1.97 -    image_size = (margin, max_pos) if textRot else (max_pos, margin)
    1.98 +    labels = nonoverlappingLabels(labels, minPixTweenLabels)
    1.99 +    image_size = (margin, end) if textRot else (end, margin)
   1.100      im = Image.new(image_mode, image_size, opts.margin_color)
   1.101      draw = ImageDraw.Draw(im)
   1.102      for i in labels:
   1.103 @@ -477,24 +478,24 @@
   1.104      textRot1 = "vertical".startswith(opts.rot1)
   1.105      i1 = getSeqInfo(opts.sort1, seqNames1, seqLimits1,
   1.106                      font, opts.fontsize, image_mode, opts.lengths1, textRot1)
   1.107 -    seqNames1, seqSizes1, seqLabels1, labelSizes1, tMargin = i1
   1.108 +    seqNames1, rangeSizes1, seqLabels1, labelSizes1, tMargin = i1
   1.109  
   1.110      textRot2 = "horizontal".startswith(opts.rot2)
   1.111      i2 = getSeqInfo(opts.sort2, seqNames2, seqLimits2,
   1.112                      font, opts.fontsize, image_mode, opts.lengths2, textRot2)
   1.113 -    seqNames2, seqSizes2, seqLabels2, labelSizes2, lMargin = i2
   1.114 +    seqNames2, rangeSizes2, seqLabels2, labelSizes2, lMargin = i2
   1.115  
   1.116      warn("choosing bp per pixel...")
   1.117      maxPixels1 = opts.width  - lMargin
   1.118      maxPixels2 = opts.height - tMargin
   1.119 -    bpPerPix1 = get_bp_per_pix(seqSizes1, opts.border_pixels, maxPixels1)
   1.120 -    bpPerPix2 = get_bp_per_pix(seqSizes2, opts.border_pixels, maxPixels2)
   1.121 +    bpPerPix1 = get_bp_per_pix(rangeSizes1, opts.border_pixels, maxPixels1)
   1.122 +    bpPerPix2 = get_bp_per_pix(rangeSizes2, opts.border_pixels, maxPixels2)
   1.123      bpPerPix = max(bpPerPix1, bpPerPix2)
   1.124      warn("bp per pixel = " + str(bpPerPix))
   1.125  
   1.126 -    seq_pix1, seq_starts1, width  = pixelData(seqSizes1, bpPerPix,
   1.127 +    seq_pix1, seq_starts1, width  = pixelData(rangeSizes1, bpPerPix,
   1.128                                                opts.border_pixels, lMargin)
   1.129 -    seq_pix2, seq_starts2, height = pixelData(seqSizes2, bpPerPix,
   1.130 +    seq_pix2, seq_starts2, height = pixelData(rangeSizes2, bpPerPix,
   1.131                                                opts.border_pixels, tMargin)
   1.132      warn("width:  " + str(width))
   1.133      warn("height: " + str(height))