scripts/last-dotplot
changeset 908 6840398323e1
parent 907 87991b57ed04
child 909 92dfad507286
     1.1 --- a/scripts/last-dotplot	Wed Nov 08 14:22:00 2017 +0900
     1.2 +++ b/scripts/last-dotplot	Wed Nov 08 15:45:08 2017 +0900
     1.3 @@ -33,20 +33,25 @@
     1.4          prog = os.path.basename(sys.argv[0])
     1.5          sys.stderr.write(prog + ": " + message + "\n")
     1.6  
     1.7 -def croppedBlocks(blocks, range1, range2):
     1.8 -    cropBeg1, cropEnd1 = range1
     1.9 -    cropBeg2, cropEnd2 = range2
    1.10 -    if blocks[0][0] < 0: cropBeg1, cropEnd1 = -cropEnd1, -cropBeg1
    1.11 -    if blocks[0][1] < 0: cropBeg2, cropEnd2 = -cropEnd2, -cropBeg2
    1.12 -    for beg1, beg2, size in blocks:
    1.13 -        b1 = max(cropBeg1, beg1)
    1.14 -        e1 = min(cropEnd1, beg1 + size)
    1.15 -        if b1 >= e1: continue
    1.16 -        offset = beg2 - beg1
    1.17 -        b2 = max(cropBeg2, b1 + offset)
    1.18 -        e2 = min(cropEnd2, e1 + offset)
    1.19 -        if b2 >= e2: continue
    1.20 -        yield b2 - offset, b2, e2 - b2
    1.21 +def croppedBlocks(blocks, ranges1, ranges2):
    1.22 +    headBeg1, headBeg2, headSize = blocks[0]
    1.23 +    for r1 in ranges1:
    1.24 +        for r2 in ranges2:
    1.25 +            cropBeg1, cropEnd1 = r1
    1.26 +            if headBeg1 < 0:
    1.27 +                cropBeg1, cropEnd1 = -cropEnd1, -cropBeg1
    1.28 +            cropBeg2, cropEnd2 = r2
    1.29 +            if headBeg2 < 0:
    1.30 +                cropBeg2, cropEnd2 = -cropEnd2, -cropBeg2
    1.31 +            for beg1, beg2, size in blocks:
    1.32 +                b1 = max(cropBeg1, beg1)
    1.33 +                e1 = min(cropEnd1, beg1 + size)
    1.34 +                if b1 >= e1: continue
    1.35 +                offset = beg2 - beg1
    1.36 +                b2 = max(cropBeg2, b1 + offset)
    1.37 +                e2 = min(cropEnd2, e1 + offset)
    1.38 +                if b2 >= e2: continue
    1.39 +                yield b2 - offset, b2, e2 - b2
    1.40  
    1.41  def tabBlocks(beg1, beg2, blocks):
    1.42      '''Get the gapless blocks of an alignment, from LAST tabular format.'''
    1.43 @@ -116,12 +121,14 @@
    1.44      return text, 0, sys.maxsize
    1.45  
    1.46  def rangeFromSeqName(seqRequests, name, seqLen):
    1.47 -    if not seqRequests: return 0, seqLen
    1.48 -    base = name.split(".")[-1]  # allow for names like hg19.chr7
    1.49 -    for pat, beg, end in seqRequests:
    1.50 -        if fnmatchcase(name, pat) or fnmatchcase(base, pat):
    1.51 -            return max(beg, 0), min(end, seqLen)
    1.52 -    return None
    1.53 +    if seqRequests:
    1.54 +        base = name.split(".")[-1]  # allow for names like hg19.chr7
    1.55 +        for pat, beg, end in seqRequests:
    1.56 +            if fnmatchcase(name, pat) or fnmatchcase(base, pat):
    1.57 +                return max(beg, 0), min(end, seqLen)
    1.58 +        return None
    1.59 +    else:
    1.60 +        return 0, seqLen
    1.61  
    1.62  def updateSeqs(isTrim, seqNames, seqLimits, seqName, seqRange, blocks, index):
    1.63      if seqName not in seqLimits:
    1.64 @@ -154,13 +161,15 @@
    1.65          if not range1: continue
    1.66          range2 = rangeFromSeqName(seqRequests2, seqName2, seqLen2)
    1.67          if not range2: continue
    1.68 -        b = list(croppedBlocks(list(blocks), range1, range2))
    1.69 +        b = list(croppedBlocks(list(blocks), [range1], [range2]))
    1.70          if not b: continue
    1.71          aln = seqName1, seqName2, b
    1.72          alignments.append(aln)
    1.73          updateSeqs(opts.trim1, seqNames1, seqLimits1, seqName1, range1, b, 0)
    1.74          updateSeqs(opts.trim2, seqNames2, seqLimits2, seqName2, range2, b, 1)
    1.75 -    return alignments, seqNames1, seqNames2, seqLimits1, seqLimits2
    1.76 +    seqRanges1 = [(i, seqLimits1[i][0], seqLimits1[i][1]) for i in seqNames1]
    1.77 +    seqRanges2 = [(i, seqLimits2[i][0], seqLimits2[i][1]) for i in seqNames2]
    1.78 +    return alignments, seqRanges1, seqRanges2
    1.79  
    1.80  def natural_sort_key(my_string):
    1.81      '''Return a sort key for "natural" ordering, e.g. chr9 < chr10.'''
    1.82 @@ -223,9 +232,8 @@
    1.83                  x, y = y, x
    1.84          yield text, x, y
    1.85  
    1.86 -def getSeqInfo(sortOpt, seqNames, seqLimits,
    1.87 -               font, fontsize, image_mode, labelOpt, textRot):
    1.88 -    seqRanges = [(i, seqLimits[i][0], seqLimits[i][1]) for i in seqNames]
    1.89 +def dataFromRanges(seqRanges, sortOpt, font,
    1.90 +                   fontsize, image_mode, labelOpt, textRot):
    1.91      s = getSortedRanges(seqRanges, sortOpt)
    1.92      for i in s:
    1.93          warn("\t".join(map(str, i)))
    1.94 @@ -536,18 +544,18 @@
    1.95  
    1.96      warn("reading alignments...")
    1.97      alnData = readAlignments(args[0], opts)
    1.98 -    alignments, seqNames1, seqNames2, seqLimits1, seqLimits2 = alnData
    1.99 +    alignments, seqRanges1, seqRanges2 = alnData
   1.100      warn("done")
   1.101      if not alignments: raise Exception("there are no alignments")
   1.102  
   1.103      textRot1 = "vertical".startswith(opts.rot1)
   1.104 -    i1 = getSeqInfo(opts.sort1, seqNames1, seqLimits1,
   1.105 -                    font, opts.fontsize, image_mode, opts.labels1, textRot1)
   1.106 +    i1 = dataFromRanges(seqRanges1, opts.sort1, font,
   1.107 +                        opts.fontsize, image_mode, opts.labels1, textRot1)
   1.108      sortedRanges1, rangeSizes1, labelData1, 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.labels2, textRot2)
   1.113 +    i2 = dataFromRanges(seqRanges2, opts.sort2, font,
   1.114 +                        opts.fontsize, image_mode, opts.labels2, textRot2)
   1.115      sortedRanges2, rangeSizes2, labelData2, lMargin = i2
   1.116  
   1.117      maxPixels1 = opts.width  - lMargin