scripts/last-dotplot
changeset 851 5b2acb7fdb3e
parent 850 db0ddc92f2ec
child 852 f3b6c666afad
     1.1 --- a/scripts/last-dotplot	Tue May 02 13:18:23 2017 +0900
     1.2 +++ b/scripts/last-dotplot	Tue May 02 13:35:07 2017 +0900
     1.3 @@ -114,7 +114,9 @@
     1.4              return max(beg, 0), min(end, seqLen)
     1.5      return None
     1.6  
     1.7 -def updateSeqLimits(isTrim, seqLimits, seqName, seqRange, blocks, index):
     1.8 +def updateSeqs(isTrim, seqNames, seqLimits, seqName, seqRange, blocks, index):
     1.9 +    if seqName not in seqLimits:
    1.10 +        seqNames.append(seqName)
    1.11      if isTrim:
    1.12          beg = blocks[0][index]
    1.13          end = blocks[-1][index] + blocks[-1][2]
    1.14 @@ -133,6 +135,8 @@
    1.15      seqRanges2 = map(seqRangeFromText, opts.seq2)
    1.16  
    1.17      alignments = []
    1.18 +    seqNames1 = []
    1.19 +    seqNames2 = []
    1.20      seqLimits1 = {}
    1.21      seqLimits2 = {}
    1.22      lines = myOpen(fileName)
    1.23 @@ -145,9 +149,9 @@
    1.24          if not b: continue
    1.25          aln = seqName1, seqName2, b
    1.26          alignments.append(aln)
    1.27 -        updateSeqLimits(opts.trim1, seqLimits1, seqName1, range1, b, 0)
    1.28 -        updateSeqLimits(opts.trim2, seqLimits2, seqName2, range2, b, 1)
    1.29 -    return alignments, seqLimits1, seqLimits2
    1.30 +        updateSeqs(opts.trim1, seqNames1, seqLimits1, seqName1, range1, b, 0)
    1.31 +        updateSeqs(opts.trim2, seqNames2, seqLimits2, seqName2, range2, b, 1)
    1.32 +    return alignments, seqNames1, seqNames2, seqLimits1, seqLimits2
    1.33  
    1.34  def natural_sort_key(my_string):
    1.35      '''Return a sort key for "natural" ordering, e.g. chr9 < chr10.'''
    1.36 @@ -175,11 +179,16 @@
    1.37  def seqNameAndSizeText(seqName, seqSize):
    1.38      return seqName + ": " + sizeText(seqSize)
    1.39  
    1.40 -def getSeqInfo(seqLimits, font, fontsize, image_mode, isShowSize):
    1.41 +def getSeqInfo(sortOpt, seqNames, seqLimits,
    1.42 +               font, fontsize, image_mode, isShowSize):
    1.43      '''Return miscellaneous information about the sequences.'''
    1.44 -    seqNames = seqLimits.keys()
    1.45 -    seqNames.sort(key=natural_sort_key)
    1.46 +    if sortOpt == 1:
    1.47 +        seqNames.sort(key=natural_sort_key)
    1.48      seqSizes = [seqLimits[i][1] - seqLimits[i][0] for i in seqNames]
    1.49 +    if sortOpt == 2:
    1.50 +        seqRecords = sorted(zip(seqSizes, seqNames), reverse=True)
    1.51 +        seqSizes = [i[0] for i in seqRecords]
    1.52 +        seqNames = [i[1] for i in seqRecords]
    1.53      if isShowSize:
    1.54          seqLabels = map(seqNameAndSizeText, seqNames, seqSizes)
    1.55      else:
    1.56 @@ -392,15 +401,17 @@
    1.57      overlap_color = tuple([(i + j) // 2 for i, j in zipped_colors])
    1.58  
    1.59      warn("reading alignments...")
    1.60 -    alignments, seqLimits1, seqLimits2 = readAlignments(args[0], opts)
    1.61 +    alignmentInfo = readAlignments(args[0], opts)
    1.62 +    alignments, seqNames1, seqNames2, seqLimits1, seqLimits2 = alignmentInfo
    1.63      warn("done")
    1.64 -
    1.65      if not alignments: raise Exception("there are no alignments")
    1.66  
    1.67 -    i1 = getSeqInfo(seqLimits1, font, opts.fontsize, image_mode, opts.lengths1)
    1.68 +    i1 = getSeqInfo(opts.sort1, seqNames1, seqLimits1,
    1.69 +                    font, opts.fontsize, image_mode, opts.lengths1)
    1.70      seqNames1, seqSizes1, seqLabels1, labelSizes1, margin1 = i1
    1.71  
    1.72 -    i2 = getSeqInfo(seqLimits2, font, opts.fontsize, image_mode, opts.lengths2)
    1.73 +    i2 = getSeqInfo(opts.sort2, seqNames2, seqLimits2,
    1.74 +                    font, opts.fontsize, image_mode, opts.lengths2)
    1.75      seqNames2, seqSizes2, seqLabels2, labelSizes2, margin2 = i2
    1.76  
    1.77      warn("choosing bp per pixel...")
    1.78 @@ -500,6 +511,12 @@
    1.79                    help="color for forward alignments (default: %default)")
    1.80      op.add_option("-r", "--reversecolor", metavar="COLOR", default="blue",
    1.81                    help="color for reverse alignments (default: %default)")
    1.82 +    op.add_option("--sort1", type="int", default=1, metavar="N",
    1.83 +                  help="genome1 sequence order: 0=input order, 1=name order, "
    1.84 +                  "2=length order (default=%default)")
    1.85 +    op.add_option("--sort2", type="int", default=1, metavar="N",
    1.86 +                  help="genome2 sequence order: 0=input order, 1=name order, "
    1.87 +                  "2=length order (default=%default)")
    1.88      op.add_option("--trim1", action="store_true",
    1.89                    help="trim unaligned sequence flanks from the 1st genome")
    1.90      op.add_option("--trim2", action="store_true",