scripts/last-dotplot
author Martin C. Frith
Fri Jul 31 15:24:44 2020 +0900 (21 months ago)
changeset 1077 5209deca661c
parent 1076 3b8d229f8359
child 1078 579d74b5066f
permissions -rwxr-xr-x
last-dotplot: show annotation names!
     1 #! /usr/bin/env python
     2 # Author: Martin C. Frith 2008
     3 # SPDX-License-Identifier: GPL-3.0-or-later
     4 
     5 # Read pair-wise alignments in MAF or LAST tabular format: write an
     6 # "Oxford grid", a.k.a. dotplot.
     7 
     8 # TODO: Currently, pixels with zero aligned nt-pairs are white, and
     9 # pixels with one or more aligned nt-pairs are black.  This can look
    10 # too crowded for large genome alignments.  I tried shading each pixel
    11 # according to the number of aligned nt-pairs within it, but the
    12 # result is too faint.  How can this be done better?
    13 
    14 import collections
    15 import functools
    16 import gzip
    17 from fnmatch import fnmatchcase
    18 import logging
    19 from operator import itemgetter
    20 import subprocess
    21 import itertools, optparse, os, re, sys
    22 
    23 # Try to make PIL/PILLOW work:
    24 try:
    25     from PIL import Image, ImageDraw, ImageFont, ImageColor
    26 except ImportError:
    27     import Image, ImageDraw, ImageFont, ImageColor
    28 
    29 try:
    30     from future_builtins import zip
    31 except ImportError:
    32     pass
    33 
    34 def myOpen(fileName):  # faster than fileinput
    35     if fileName is None:
    36         return []
    37     if fileName == "-":
    38         return sys.stdin
    39     if fileName.endswith(".gz"):
    40         return gzip.open(fileName, "rt")  # xxx dubious for Python2
    41     return open(fileName)
    42 
    43 def groupByFirstItem(things):
    44     for k, v in itertools.groupby(things, itemgetter(0)):
    45         yield k, [i[1:] for i in v]
    46 
    47 def croppedBlocks(blocks, ranges1, ranges2):
    48     headBeg1, headBeg2, headSize = blocks[0]
    49     for r1 in ranges1:
    50         for r2 in ranges2:
    51             cropBeg1, cropEnd1 = r1
    52             if headBeg1 < 0:
    53                 cropBeg1, cropEnd1 = -cropEnd1, -cropBeg1
    54             cropBeg2, cropEnd2 = r2
    55             if headBeg2 < 0:
    56                 cropBeg2, cropEnd2 = -cropEnd2, -cropBeg2
    57             for beg1, beg2, size in blocks:
    58                 b1 = max(cropBeg1, beg1)
    59                 e1 = min(cropEnd1, beg1 + size)
    60                 if b1 >= e1: continue
    61                 offset = beg2 - beg1
    62                 b2 = max(cropBeg2, b1 + offset)
    63                 e2 = min(cropEnd2, e1 + offset)
    64                 if b2 >= e2: continue
    65                 yield b2 - offset, b2, e2 - b2
    66 
    67 def tabBlocks(beg1, beg2, blocks):
    68     '''Get the gapless blocks of an alignment, from LAST tabular format.'''
    69     for i in blocks.split(","):
    70         if ":" in i:
    71             x, y = i.split(":")
    72             beg1 += int(x)
    73             beg2 += int(y)
    74         else:
    75             size = int(i)
    76             yield beg1, beg2, size
    77             beg1 += size
    78             beg2 += size
    79 
    80 def mafBlocks(beg1, beg2, seq1, seq2):
    81     '''Get the gapless blocks of an alignment, from MAF format.'''
    82     size = 0
    83     for x, y in zip(seq1, seq2):
    84         if x == "-":
    85             if size:
    86                 yield beg1, beg2, size
    87                 beg1 += size
    88                 beg2 += size
    89                 size = 0
    90             beg2 += 1
    91         elif y == "-":
    92             if size:
    93                 yield beg1, beg2, size
    94                 beg1 += size
    95                 beg2 += size
    96                 size = 0
    97             beg1 += 1
    98         else:
    99             size += 1
   100     if size: yield beg1, beg2, size
   101 
   102 def alignmentFromSegment(qrySeqName, qrySeqLen, segment):
   103     refSeqLen = sys.maxsize  # XXX
   104     refSeqName, refSeqBeg, qrySeqBeg, size = segment
   105     block = refSeqBeg, qrySeqBeg, size
   106     return refSeqName, refSeqLen, qrySeqName, qrySeqLen, [block]
   107 
   108 def alignmentInput(lines):
   109     '''Get alignments and sequence lengths, from MAF or tabular format.'''
   110     mafCount = 0
   111     qrySeqName = ""
   112     segments = []
   113     for line in lines:
   114         w = line.split()
   115         if line[0] == "#":
   116             pass
   117         elif len(w) == 1:
   118             for i in segments:
   119                 yield alignmentFromSegment(qrySeqName, qrySeqLen, i)
   120             qrySeqName = w[0]
   121             qrySeqLen = 0
   122             segments = []
   123         elif len(w) == 2 and qrySeqName and w[1].isdigit():
   124             qrySeqLen += int(w[1])
   125         elif len(w) == 4 and qrySeqName and w[1].isdigit() and w[3].isdigit():
   126             refSeqName, refSeqBeg, refSeqEnd = w[0], int(w[1]), int(w[3])
   127             size = abs(refSeqEnd - refSeqBeg)
   128             if refSeqBeg > refSeqEnd:
   129                 refSeqBeg = -refSeqBeg
   130             segments.append((refSeqName, refSeqBeg, qrySeqLen, size))
   131             qrySeqLen += size
   132         elif line[0].isdigit():  # tabular format
   133             chr1, beg1, seqlen1 = w[1], int(w[2]), int(w[5])
   134             if w[4] == "-": beg1 -= seqlen1
   135             chr2, beg2, seqlen2 = w[6], int(w[7]), int(w[10])
   136             if w[9] == "-": beg2 -= seqlen2
   137             blocks = tabBlocks(beg1, beg2, w[11])
   138             yield chr1, seqlen1, chr2, seqlen2, blocks
   139         elif line[0] == "s":  # MAF format
   140             if mafCount == 0:
   141                 chr1, beg1, seqlen1, seq1 = w[1], int(w[2]), int(w[5]), w[6]
   142                 if w[4] == "-": beg1 -= seqlen1
   143                 mafCount = 1
   144             else:
   145                 chr2, beg2, seqlen2, seq2 = w[1], int(w[2]), int(w[5]), w[6]
   146                 if w[4] == "-": beg2 -= seqlen2
   147                 blocks = mafBlocks(beg1, beg2, seq1, seq2)
   148                 yield chr1, seqlen1, chr2, seqlen2, blocks
   149                 mafCount = 0
   150     for i in segments:
   151         yield alignmentFromSegment(qrySeqName, qrySeqLen, i)
   152 
   153 def seqRequestFromText(text):
   154     if ":" in text:
   155         pattern, interval = text.rsplit(":", 1)
   156         if "-" in interval:
   157             beg, end = interval.rsplit("-", 1)
   158             return pattern, int(beg), int(end)  # beg may be negative
   159     return text, 0, sys.maxsize
   160 
   161 def rangesFromSeqName(seqRequests, name, seqLen):
   162     if seqRequests:
   163         base = name.split(".")[-1]  # allow for names like hg19.chr7
   164         for pat, beg, end in seqRequests:
   165             if fnmatchcase(name, pat) or fnmatchcase(base, pat):
   166                 yield max(beg, 0), min(end, seqLen)
   167     else:
   168         yield 0, seqLen
   169 
   170 def updateSeqs(coverDict, seqRanges, seqName, ranges, coveredRange):
   171     beg, end = coveredRange
   172     if beg < 0:
   173         coveredRange = -end, -beg
   174     if seqName in coverDict:
   175         coverDict[seqName].append(coveredRange)
   176     else:
   177         coverDict[seqName] = [coveredRange]
   178         for beg, end in ranges:
   179             r = seqName, beg, end
   180             seqRanges.append(r)
   181 
   182 def readAlignments(fileName, opts):
   183     '''Get alignments and sequence limits, from MAF or tabular format.'''
   184     seqRequests1 = [seqRequestFromText(i) for i in opts.seq1]
   185     seqRequests2 = [seqRequestFromText(i) for i in opts.seq2]
   186 
   187     alignments = []
   188     seqRanges1 = []
   189     seqRanges2 = []
   190     coverDict1 = {}
   191     coverDict2 = {}
   192     lines = myOpen(fileName)
   193     for seqName1, seqLen1, seqName2, seqLen2, blocks in alignmentInput(lines):
   194         ranges1 = sorted(rangesFromSeqName(seqRequests1, seqName1, seqLen1))
   195         if not ranges1: continue
   196         ranges2 = sorted(rangesFromSeqName(seqRequests2, seqName2, seqLen2))
   197         if not ranges2: continue
   198         b = list(croppedBlocks(list(blocks), ranges1, ranges2))
   199         if not b: continue
   200         aln = seqName1, seqName2, b
   201         alignments.append(aln)
   202         coveredRange1 = b[0][0], b[-1][0] + b[-1][2]
   203         updateSeqs(coverDict1, seqRanges1, seqName1, ranges1, coveredRange1)
   204         coveredRange2 = b[0][1], b[-1][1] + b[-1][2]
   205         updateSeqs(coverDict2, seqRanges2, seqName2, ranges2, coveredRange2)
   206     return alignments, seqRanges1, coverDict1, seqRanges2, coverDict2
   207 
   208 def nameAndRangesFromDict(cropDict, seqName):
   209     if seqName in cropDict:
   210         return seqName, cropDict[seqName]
   211     n = seqName.split(".")[-1]
   212     if n in cropDict:
   213         return n, cropDict[n]
   214     return seqName, []
   215 
   216 def rangesForSecondaryAlignments(primaryRanges, seqLen):
   217     if primaryRanges:
   218         return primaryRanges
   219     return [(0, seqLen)]
   220 
   221 def readSecondaryAlignments(opts, cropRanges1, cropRanges2):
   222     cropDict1 = dict(groupByFirstItem(cropRanges1))
   223     cropDict2 = dict(groupByFirstItem(cropRanges2))
   224 
   225     alignments = []
   226     seqRanges1 = []
   227     seqRanges2 = []
   228     coverDict1 = {}
   229     coverDict2 = {}
   230     lines = myOpen(opts.alignments)
   231     for seqName1, seqLen1, seqName2, seqLen2, blocks in alignmentInput(lines):
   232         seqName1, ranges1 = nameAndRangesFromDict(cropDict1, seqName1)
   233         seqName2, ranges2 = nameAndRangesFromDict(cropDict2, seqName2)
   234         if not ranges1 and not ranges2:
   235             continue
   236         r1 = rangesForSecondaryAlignments(ranges1, seqLen1)
   237         r2 = rangesForSecondaryAlignments(ranges2, seqLen2)
   238         b = list(croppedBlocks(list(blocks), r1, r2))
   239         if not b: continue
   240         aln = seqName1, seqName2, b
   241         alignments.append(aln)
   242         if not ranges1:
   243             coveredRange1 = b[0][0], b[-1][0] + b[-1][2]
   244             updateSeqs(coverDict1, seqRanges1, seqName1, r1, coveredRange1)
   245         if not ranges2:
   246             coveredRange2 = b[0][1], b[-1][1] + b[-1][2]
   247             updateSeqs(coverDict2, seqRanges2, seqName2, r2, coveredRange2)
   248     return alignments, seqRanges1, coverDict1, seqRanges2, coverDict2
   249 
   250 def twoValuesFromOption(text, separator):
   251     if separator in text:
   252         return text.split(separator)
   253     return text, text
   254 
   255 def mergedRanges(ranges):
   256     oldBeg, maxEnd = ranges[0]
   257     for beg, end in ranges:
   258         if beg > maxEnd:
   259             yield oldBeg, maxEnd
   260             oldBeg = beg
   261             maxEnd = end
   262         elif end > maxEnd:
   263             maxEnd = end
   264     yield oldBeg, maxEnd
   265 
   266 def mergedRangesPerSeq(coverDict):
   267     for k, v in coverDict.items():
   268         v.sort()
   269         yield k, list(mergedRanges(v))
   270 
   271 def coveredLength(mergedCoverDict):
   272     return sum(sum(e - b for b, e in v) for v in mergedCoverDict.values())
   273 
   274 def trimmed(seqRanges, coverDict, minAlignedBases, maxGapFrac, endPad, midPad):
   275     maxEndGapFrac, maxMidGapFrac = twoValuesFromOption(maxGapFrac, ",")
   276     maxEndGap = max(float(maxEndGapFrac) * minAlignedBases, endPad * 1.0)
   277     maxMidGap = max(float(maxMidGapFrac) * minAlignedBases, midPad * 2.0)
   278 
   279     for seqName, rangeBeg, rangeEnd in seqRanges:
   280         seqBlocks = coverDict[seqName]
   281         blocks = [i for i in seqBlocks if i[0] < rangeEnd and i[1] > rangeBeg]
   282         if blocks[0][0] - rangeBeg > maxEndGap:
   283             rangeBeg = blocks[0][0] - endPad
   284         for j, y in enumerate(blocks):
   285             if j:
   286                 x = blocks[j - 1]
   287                 if y[0] - x[1] > maxMidGap:
   288                     yield seqName, rangeBeg, x[1] + midPad
   289                     rangeBeg = y[0] - midPad
   290         if rangeEnd - blocks[-1][1] > maxEndGap:
   291             rangeEnd = blocks[-1][1] + endPad
   292         yield seqName, rangeBeg, rangeEnd
   293 
   294 def rangesWithStrandInfo(seqRanges, strandOpt, alignments, seqIndex):
   295     if strandOpt == "1":
   296         forwardMinusReverse = collections.defaultdict(int)
   297         for i in alignments:
   298             blocks = i[2]
   299             beg1, beg2, size = blocks[0]
   300             numOfAlignedLetterPairs = sum(i[2] for i in blocks)
   301             if (beg1 < 0) != (beg2 < 0):  # opposite-strand alignment
   302                 numOfAlignedLetterPairs *= -1
   303             forwardMinusReverse[i[seqIndex]] += numOfAlignedLetterPairs
   304     strandNum = 0
   305     for seqName, beg, end in seqRanges:
   306         if strandOpt == "1":
   307             strandNum = 1 if forwardMinusReverse[seqName] >= 0 else 2
   308         yield seqName, beg, end, strandNum
   309 
   310 def natural_sort_key(my_string):
   311     '''Return a sort key for "natural" ordering, e.g. chr9 < chr10.'''
   312     parts = re.split(r'(\d+)', my_string)
   313     parts[1::2] = map(int, parts[1::2])
   314     return parts
   315 
   316 def nameKey(oneSeqRanges):
   317     return natural_sort_key(oneSeqRanges[0][0])
   318 
   319 def sizeKey(oneSeqRanges):
   320     return sum(b - e for n, b, e, s in oneSeqRanges), nameKey(oneSeqRanges)
   321 
   322 def alignmentKey(seqNamesToLists, oneSeqRanges):
   323     seqName = oneSeqRanges[0][0]
   324     alignmentsOfThisSequence = seqNamesToLists[seqName]
   325     numOfAlignedLetterPairs = sum(i[3] for i in alignmentsOfThisSequence)
   326     toMiddle = numOfAlignedLetterPairs // 2
   327     for i in alignmentsOfThisSequence:
   328         toMiddle -= i[3]
   329         if toMiddle < 0:
   330             return i[1:3]  # sequence-rank and "position" of this alignment
   331 
   332 def rankAndFlipPerSeq(seqRanges):
   333     rangesGroupedBySeqName = itertools.groupby(seqRanges, itemgetter(0))
   334     for rank, group in enumerate(rangesGroupedBySeqName):
   335         seqName, ranges = group
   336         strandNum = next(ranges)[3]
   337         flip = 1 if strandNum < 2 else -1
   338         yield seqName, (rank, flip)
   339 
   340 def alignmentSortData(alignments, seqIndex, otherNamesToRanksAndFlips):
   341     otherIndex = 1 - seqIndex
   342     for i in alignments:
   343         blocks = i[2]
   344         otherRank, otherFlip = otherNamesToRanksAndFlips[i[otherIndex]]
   345         otherPos = otherFlip * abs(blocks[0][otherIndex] +
   346                                    blocks[-1][otherIndex] + blocks[-1][2])
   347         numOfAlignedLetterPairs = sum(i[2] for i in blocks)
   348         yield i[seqIndex], otherRank, otherPos, numOfAlignedLetterPairs
   349 
   350 def mySortedRanges(seqRanges, sortOpt, seqIndex, alignments, otherRanges):
   351     rangesGroupedBySeqName = itertools.groupby(seqRanges, itemgetter(0))
   352     g = [list(ranges) for seqName, ranges in rangesGroupedBySeqName]
   353     for i in g:
   354         if i[0][3] > 1:
   355             i.reverse()
   356     if sortOpt == "1":
   357         g.sort(key=nameKey)
   358     if sortOpt == "2":
   359         g.sort(key=sizeKey)
   360     if sortOpt == "3":
   361         otherNamesToRanksAndFlips = dict(rankAndFlipPerSeq(otherRanges))
   362         alns = sorted(alignmentSortData(alignments, seqIndex,
   363                                         otherNamesToRanksAndFlips))
   364         alnsGroupedBySeqName = itertools.groupby(alns, itemgetter(0))
   365         seqNamesToLists = dict((k, list(v)) for k, v in alnsGroupedBySeqName)
   366         g.sort(key=functools.partial(alignmentKey, seqNamesToLists))
   367     return [j for i in g for j in i]
   368 
   369 def allSortedRanges(opts, alignments, alignmentsB,
   370                     seqRanges1, seqRangesB1, seqRanges2, seqRangesB2):
   371     o1, oB1 = twoValuesFromOption(opts.strands1, ":")
   372     o2, oB2 = twoValuesFromOption(opts.strands2, ":")
   373     if o1 == "1" and o2 == "1":
   374         raise Exception("the strand options have circular dependency")
   375     seqRanges1 = list(rangesWithStrandInfo(seqRanges1, o1, alignments, 0))
   376     seqRanges2 = list(rangesWithStrandInfo(seqRanges2, o2, alignments, 1))
   377     seqRangesB1 = list(rangesWithStrandInfo(seqRangesB1, oB1, alignmentsB, 0))
   378     seqRangesB2 = list(rangesWithStrandInfo(seqRangesB2, oB2, alignmentsB, 1))
   379 
   380     o1, oB1 = twoValuesFromOption(opts.sort1, ":")
   381     o2, oB2 = twoValuesFromOption(opts.sort2, ":")
   382     if o1 == "3" and o2 == "3":
   383         raise Exception("the sort options have circular dependency")
   384     if o1 != "3":
   385         s1 = mySortedRanges(seqRanges1, o1, None, None, None)
   386     if o2 != "3":
   387         s2 = mySortedRanges(seqRanges2, o2, None, None, None)
   388     if o1 == "3":
   389         s1 = mySortedRanges(seqRanges1, o1, 0, alignments, s2)
   390     if o2 == "3":
   391         s2 = mySortedRanges(seqRanges2, o2, 1, alignments, s1)
   392     t1 = mySortedRanges(seqRangesB1, oB1, 0, alignmentsB, s2)
   393     t2 = mySortedRanges(seqRangesB2, oB2, 1, alignmentsB, s1)
   394     return s1 + t1, s2 + t2
   395 
   396 def sizesPerText(texts, font, textDraw):
   397     sizes = 0, 0
   398     for t in texts:
   399         if textDraw is not None:
   400             sizes = textDraw.textsize(t, font=font)
   401         yield t, sizes
   402 
   403 def prettyNum(n):
   404     t = str(n)
   405     groups = []
   406     while t:
   407         groups.append(t[-3:])
   408         t = t[:-3]
   409     return ",".join(reversed(groups))
   410 
   411 def sizeText(size):
   412     suffixes = "bp", "kb", "Mb", "Gb"
   413     for i, x in enumerate(suffixes):
   414         j = 10 ** (i * 3)
   415         if size < j * 10:
   416             return "%.2g" % (1.0 * size / j) + x
   417         if size < j * 1000 or i == len(suffixes) - 1:
   418             return "%.0f" % (1.0 * size / j) + x
   419 
   420 def labelText(seqRange, labelOpt):
   421     seqName, beg, end, strandNum = seqRange
   422     if labelOpt == 1:
   423         return seqName + ": " + sizeText(end - beg)
   424     if labelOpt == 2:
   425         return seqName + ":" + prettyNum(beg) + ": " + sizeText(end - beg)
   426     if labelOpt == 3:
   427         return seqName + ":" + prettyNum(beg) + "-" + prettyNum(end)
   428     return seqName
   429 
   430 def rangeLabels(seqRanges, labelOpt, font, textDraw, textRot):
   431     x = y = 0
   432     for r in seqRanges:
   433         text = labelText(r, labelOpt)
   434         if textDraw is not None:
   435             x, y = textDraw.textsize(text, font=font)
   436             if textRot:
   437                 x, y = y, x
   438         yield text, x, y, r[3]
   439 
   440 def dataFromRanges(sortedRanges, font, textDraw, labelOpt, textRot):
   441     for seqName, rangeBeg, rangeEnd, strandNum in sortedRanges:
   442         out = [seqName, str(rangeBeg), str(rangeEnd)]
   443         if strandNum > 0:
   444             out.append(".+-"[strandNum])
   445         logging.info("\t".join(out))
   446     logging.info("")
   447     rangeSizes = [e - b for n, b, e, s in sortedRanges]
   448     labs = list(rangeLabels(sortedRanges, labelOpt, font, textDraw, textRot))
   449     margin = max(i[2] for i in labs)
   450     # xxx the margin may be too big, because some labels may get omitted
   451     return rangeSizes, labs, margin
   452 
   453 def div_ceil(x, y):
   454     '''Return x / y rounded up.'''
   455     q, r = divmod(x, y)
   456     return q + (r != 0)
   457 
   458 def get_bp_per_pix(rangeSizes, pixTweenRanges, maxPixels):
   459     '''Get the minimum bp-per-pixel that fits in the size limit.'''
   460     logging.info("choosing bp per pixel...")
   461     numOfRanges = len(rangeSizes)
   462     maxPixelsInRanges = maxPixels - pixTweenRanges * (numOfRanges - 1)
   463     if maxPixelsInRanges < numOfRanges:
   464         raise Exception("can't fit the image: too many sequences?")
   465     negLimit = -maxPixelsInRanges
   466     negBpPerPix = sum(rangeSizes) // negLimit
   467     while True:
   468         if sum(i // negBpPerPix for i in rangeSizes) >= negLimit:
   469             return -negBpPerPix
   470         negBpPerPix -= 1
   471 
   472 def getRangePixBegs(rangePixLens, pixTweenRanges, margin):
   473     '''Get the start pixel for each range.'''
   474     rangePixBegs = []
   475     pix_tot = margin - pixTweenRanges
   476     for i in rangePixLens:
   477         pix_tot += pixTweenRanges
   478         rangePixBegs.append(pix_tot)
   479         pix_tot += i
   480     return rangePixBegs
   481 
   482 def pixelData(rangeSizes, bp_per_pix, pixTweenRanges, margin):
   483     '''Return pixel information about the ranges.'''
   484     rangePixLens = [div_ceil(i, bp_per_pix) for i in rangeSizes]
   485     rangePixBegs = getRangePixBegs(rangePixLens, pixTweenRanges, margin)
   486     tot_pix = rangePixBegs[-1] + rangePixLens[-1]
   487     return rangePixBegs, rangePixLens, tot_pix
   488 
   489 def drawLineForward(hits, width, bp_per_pix, beg1, beg2, size):
   490     while True:
   491         q1, r1 = divmod(beg1, bp_per_pix)
   492         q2, r2 = divmod(beg2, bp_per_pix)
   493         hits[q2 * width + q1] |= 1
   494         next_pix = min(bp_per_pix - r1, bp_per_pix - r2)
   495         if next_pix >= size: break
   496         beg1 += next_pix
   497         beg2 += next_pix
   498         size -= next_pix
   499 
   500 def drawLineReverse(hits, width, bp_per_pix, beg1, beg2, size):
   501     while True:
   502         q1, r1 = divmod(beg1, bp_per_pix)
   503         q2, r2 = divmod(beg2, bp_per_pix)
   504         hits[q2 * width + q1] |= 2
   505         next_pix = min(bp_per_pix - r1, r2 + 1)
   506         if next_pix >= size: break
   507         beg1 += next_pix
   508         beg2 -= next_pix
   509         size -= next_pix
   510 
   511 def strandAndOrigin(ranges, beg, size):
   512     isReverseStrand = (beg < 0)
   513     if isReverseStrand:
   514         beg = -(beg + size)
   515     for rangeBeg, rangeEnd, isReverseRange, origin in ranges:
   516         if rangeEnd > beg:  # assumes the ranges are sorted
   517             return (isReverseStrand != isReverseRange), origin
   518 
   519 def alignmentPixels(width, height, alignments, bp_per_pix,
   520                     rangeDict1, rangeDict2):
   521     hits = [0] * (width * height)  # the image data
   522     for seq1, seq2, blocks in alignments:
   523         beg1, beg2, size = blocks[0]
   524         isReverse1, ori1 = strandAndOrigin(rangeDict1[seq1], beg1, size)
   525         isReverse2, ori2 = strandAndOrigin(rangeDict2[seq2], beg2, size)
   526         for beg1, beg2, size in blocks:
   527             if isReverse1:
   528                 beg1 = -(beg1 + size)
   529                 beg2 = -(beg2 + size)
   530             if isReverse1 == isReverse2:
   531                 drawLineForward(hits, width, bp_per_pix,
   532                                 ori1 + beg1, ori2 + beg2, size)
   533             else:
   534                 drawLineReverse(hits, width, bp_per_pix,
   535                                 ori1 + beg1, ori2 - beg2 - 1, size)
   536     return hits
   537 
   538 def orientedBlocks(alignments, seqIndex):
   539     otherIndex = 1 - seqIndex
   540     for a in alignments:
   541         seq1, seq2, blocks = a
   542         for b in blocks:
   543             beg1, beg2, size = b
   544             if b[seqIndex] < 0:
   545                 b = -(beg1 + size), -(beg2 + size), size
   546             yield a[seqIndex], b[seqIndex], a[otherIndex], b[otherIndex], size
   547 
   548 def drawJoins(im, alignments, bpPerPix, seqIndex, rangeDict1, rangeDict2):
   549     blocks = orientedBlocks(alignments, seqIndex)
   550     oldSeq1 = ""
   551     for seq1, beg1, seq2, beg2, size in sorted(blocks):
   552         isReverse1, ori1 = strandAndOrigin(rangeDict1[seq1], beg1, size)
   553         isReverse2, ori2 = strandAndOrigin(rangeDict2[seq2], beg2, size)
   554         end1 = beg1 + size - 1
   555         end2 = beg2 + size - 1
   556         if isReverse1:
   557             beg1 = -(beg1 + 1)
   558             end1 = -(end1 + 1)
   559         if isReverse2:
   560             beg2 = -(beg2 + 1)
   561             end2 = -(end2 + 1)
   562         newPix1 = (ori1 + beg1) // bpPerPix
   563         newPix2 = (ori2 + beg2) // bpPerPix
   564         if seq1 == oldSeq1:
   565             lowerPix2 = min(oldPix2, newPix2)
   566             upperPix2 = max(oldPix2, newPix2)
   567             midPix1 = (oldPix1 + newPix1) // 2
   568             if isReverse1:
   569                 midPix1 = (oldPix1 + newPix1 + 1) // 2
   570                 oldPix1, newPix1 = newPix1, oldPix1
   571             if upperPix2 - lowerPix2 > 1 and oldPix1 <= newPix1 <= oldPix1 + 1:
   572                 if seqIndex == 0:
   573                     box = midPix1, lowerPix2, midPix1 + 1, upperPix2 + 1
   574                 else:
   575                     box = lowerPix2, midPix1, upperPix2 + 1, midPix1 + 1
   576                 im.paste("lightgray", box)
   577         oldPix1 = (ori1 + end1) // bpPerPix
   578         oldPix2 = (ori2 + end2) // bpPerPix
   579         oldSeq1 = seq1
   580 
   581 def expandedSeqDict(seqDict):
   582     '''Allow lookup by short sequence names, e.g. chr7 as well as hg19.chr7.'''
   583     newDict = seqDict.copy()
   584     for name, x in seqDict.items():
   585         if "." in name:
   586             base = name.split(".")[-1]
   587             if base in newDict:  # an ambiguous case was found:
   588                 return seqDict   # so give up completely
   589             newDict[base] = x
   590     return newDict
   591 
   592 def readBed(fileName, rangeDict):
   593     for line in myOpen(fileName):
   594         w = line.split()
   595         if not w: continue
   596         seqName = w[0]
   597         if seqName not in rangeDict: continue
   598         seqRanges = rangeDict[seqName]
   599         beg = int(w[1])
   600         end = int(w[2])
   601         if all(beg >= i[2] or end <= i[1] for i in seqRanges):
   602             continue
   603         itemName = w[3] if len(w) > 3 and w[3] != "." else ""
   604         layer = 900
   605         color = "#fbf"
   606         if len(w) > 4:
   607             if w[4] != ".":
   608                 layer = float(w[4])
   609             if len(w) > 5:
   610                 if len(w) > 8 and w[8].count(",") == 2:
   611                     color = "rgb(" + w[8] + ")"
   612                 else:
   613                     strand = w[5]
   614                     isRev = (seqRanges[0][3] > 1)
   615                     if strand == "+" and not isRev or strand == "-" and isRev:
   616                         color = "#ffe8e8"
   617                     if strand == "-" and not isRev or strand == "+" and isRev:
   618                         color = "#e8e8ff"
   619         yield layer, color, seqName, beg, end, itemName
   620 
   621 def commaSeparatedInts(text):
   622     return map(int, text.rstrip(",").split(","))
   623 
   624 def readGenePred(opts, fileName, rangeDict):
   625     for line in myOpen(fileName):
   626         fields = line.split()  # xxx tab ???
   627         if not fields: continue
   628         geneName = fields[12 if len(fields) > 12 else 0]  # XXX ???
   629         if fields[2] not in "+-":
   630             fields = fields[1:]
   631         seqName = fields[1]
   632         if seqName not in rangeDict: continue
   633         seqRanges = rangeDict[seqName]
   634         #strand = fields[2]
   635         cdsBeg = int(fields[5])
   636         cdsEnd = int(fields[6])
   637         exonBegs = commaSeparatedInts(fields[8])
   638         exonEnds = commaSeparatedInts(fields[9])
   639         for beg, end in zip(exonBegs, exonEnds):
   640             if all(beg >= i[2] or end <= i[1] for i in seqRanges):
   641                 continue
   642             yield 300, opts.exon_color, seqName, beg, end, geneName
   643             b = max(beg, cdsBeg)
   644             e = min(end, cdsEnd)
   645             if b < e: yield 400, opts.cds_color, seqName, b, e, ""
   646 
   647 def readRmsk(fileName, rangeDict):
   648     for line in myOpen(fileName):
   649         fields = line.split()
   650         if len(fields) == 17:  # rmsk.txt
   651             seqName = fields[5]
   652             if seqName not in rangeDict: continue  # do this ASAP for speed
   653             beg = int(fields[6])
   654             end = int(fields[7])
   655             strand = fields[9]
   656             repeatName = fields[10]
   657             repeatClass = fields[11]
   658         elif len(fields) == 15:  # .out
   659             seqName = fields[4]
   660             if seqName not in rangeDict: continue
   661             beg = int(fields[5]) - 1
   662             end = int(fields[6])
   663             strand = fields[8]
   664             repeatName = fields[9]
   665             repeatClass = fields[10].split("/")[0]
   666         else:
   667             continue
   668         seqRanges = rangeDict[seqName]
   669         if all(beg >= i[2] or end <= i[1] for i in seqRanges):
   670             continue
   671         if repeatClass in ("Low_complexity", "Simple_repeat", "Satellite"):
   672             yield 200, "#fbf", seqName, beg, end, repeatName
   673         elif (strand == "+") != (seqRanges[0][3] > 1):
   674             yield 100, "#ffe8e8", seqName, beg, end, repeatName
   675         else:
   676             yield 100, "#e8e8ff", seqName, beg, end, repeatName
   677 
   678 def isExtraFirstGapField(fields):
   679     return fields[4].isdigit()
   680 
   681 def readGaps(opts, fileName, rangeDict):
   682     '''Read locations of unsequenced gaps, from an agp or gap file.'''
   683     for line in myOpen(fileName):
   684         w = line.split()
   685         if not w or w[0][0] == "#": continue
   686         if isExtraFirstGapField(w): w = w[1:]
   687         if w[4] not in "NU": continue
   688         seqName = w[0]
   689         if seqName not in rangeDict: continue
   690         end = int(w[2])
   691         beg = end - int(w[5])  # zero-based coordinate
   692         if w[7] == "yes":
   693             yield 3000, opts.bridged_color, seqName, beg, end, ""
   694         else:
   695             yield 2000, opts.unbridged_color, seqName, beg, end, ""
   696 
   697 def bedBoxes(beds, rangeDict, edge, isTop, bpPerPix, textSizes):
   698     cover = [(edge, edge)]
   699     for layer, color, seqName, bedBeg, bedEnd, name in reversed(beds):
   700         for rangeBeg, rangeEnd, isReverseRange, origin in rangeDict[seqName]:
   701             beg = max(bedBeg, rangeBeg)
   702             end = min(bedEnd, rangeEnd)
   703             if beg >= end: continue
   704             if isReverseRange:
   705                 beg, end = -end, -beg
   706             if layer <= 1000:
   707                 # include partly-covered pixels
   708                 pixBeg = (origin + beg) // bpPerPix
   709                 pixEnd = div_ceil(origin + end, bpPerPix)
   710             else:
   711                 # exclude partly-covered pixels
   712                 pixBeg = div_ceil(origin + beg, bpPerPix)
   713                 pixEnd = (origin + end) // bpPerPix
   714                 if pixEnd <= pixBeg: continue
   715                 if bedEnd >= rangeEnd:  # include partly-covered end pixels
   716                     if isReverseRange:
   717                         pixBeg = (origin + beg) // bpPerPix
   718                     else:
   719                         pixEnd = div_ceil(origin + end, bpPerPix)
   720             textWidth, textHeight = textSizes[name]
   721             nameBeg = (pixBeg + pixEnd - textHeight) // 2
   722             nameEnd = nameBeg + textHeight
   723             if name and all(e <= nameBeg or b >= nameEnd for b, e in cover):
   724                 cover.append((nameBeg, nameEnd))
   725             else:
   726                 name = ""
   727             yield layer, color, isTop, pixBeg, pixEnd, name, nameBeg, textWidth
   728 
   729 def drawAnnotations(im, boxes, tMargin, bMarginBeg, lMargin, rMarginBeg):
   730     # xxx use partial transparency for different-color overlaps?
   731     for layer, color, isTop, beg, end, name, nameBeg, nameLen in boxes:
   732         if isTop:
   733             box = beg, tMargin, end, bMarginBeg
   734         else:
   735             box = lMargin, beg, rMarginBeg, end
   736         im.paste(color, box)
   737 
   738 def placedLabels(labels, rangePixBegs, rangePixLens, beg, end):
   739     '''Return axis labels with endpoint & sort-order information.'''
   740     maxWidth = end - beg
   741     for i, j, k in zip(labels, rangePixBegs, rangePixLens):
   742         text, textWidth, textHeight, strandNum = i
   743         if textWidth > maxWidth:
   744             continue
   745         labelBeg = j + (k - textWidth) // 2
   746         labelEnd = labelBeg + textWidth
   747         sortKey = textWidth - k
   748         if labelBeg < beg:
   749             sortKey += maxWidth * (beg - labelBeg)
   750             labelBeg = beg
   751             labelEnd = beg + textWidth
   752         if labelEnd > end:
   753             sortKey += maxWidth * (labelEnd - end)
   754             labelEnd = end
   755             labelBeg = end - textWidth
   756         yield sortKey, labelBeg, labelEnd, text, textHeight, strandNum
   757 
   758 def nonoverlappingLabels(labels, minPixTweenLabels):
   759     '''Get a subset of non-overlapping axis labels, greedily.'''
   760     out = []
   761     for i in labels:
   762         beg = i[1] - minPixTweenLabels
   763         end = i[2] + minPixTweenLabels
   764         if all(j[2] <= beg or j[1] >= end for j in out):
   765             out.append(i)
   766     return out
   767 
   768 def axisImage(labels, rangePixBegs, rangePixLens, textRot,
   769               textAln, font, image_mode, opts):
   770     '''Make an image of axis labels.'''
   771     beg = rangePixBegs[0]
   772     end = rangePixBegs[-1] + rangePixLens[-1]
   773     margin = max(i[2] for i in labels)
   774     labels = sorted(placedLabels(labels, rangePixBegs, rangePixLens, beg, end))
   775     minPixTweenLabels = 0 if textRot else opts.label_space
   776     labels = nonoverlappingLabels(labels, minPixTweenLabels)
   777     image_size = (margin, end) if textRot else (end, margin)
   778     im = Image.new(image_mode, image_size, opts.margin_color)
   779     draw = ImageDraw.Draw(im)
   780     for sortKey, labelBeg, labelEnd, text, textHeight, strandNum in labels:
   781         base = margin - textHeight if textAln else 0
   782         position = (base, labelBeg) if textRot else (labelBeg, base)
   783         fill = ("black", opts.forwardcolor, opts.reversecolor)[strandNum]
   784         draw.text(position, text, font=font, fill=fill)
   785     return im
   786 
   787 def annoTextImage(opts, image_mode, font, margin, length, boxes, isLeftAlign):
   788     image_size = margin, length
   789     im = Image.new(image_mode, image_size, opts.margin_color)
   790     draw = ImageDraw.Draw(im)
   791     for layer, color, isTop, beg, end, name, nameBeg, nameLen in boxes:
   792         xPosition = 0 if isLeftAlign else margin - nameLen
   793         position = xPosition, nameBeg
   794         draw.text(position, name, font=font, fill="black")
   795     return im
   796 
   797 def rangesPerSeq(sortedRanges):
   798     for seqName, group in itertools.groupby(sortedRanges, itemgetter(0)):
   799         yield seqName, sorted(group)
   800 
   801 def rangesWithOrigins(sortedRanges, rangePixBegs, rangePixLens, bpPerPix):
   802     for i, j, k in zip(sortedRanges, rangePixBegs, rangePixLens):
   803         seqName, rangeBeg, rangeEnd, strandNum = i
   804         isReverseRange = (strandNum > 1)
   805         if isReverseRange:
   806             origin = bpPerPix * (j + k) + rangeBeg
   807         else:
   808             origin = bpPerPix * j - rangeBeg
   809         yield seqName, (rangeBeg, rangeEnd, isReverseRange, origin)
   810 
   811 def rangesAndOriginsPerSeq(sortedRanges, rangePixBegs, rangePixLens, bpPerPix):
   812     a = rangesWithOrigins(sortedRanges, rangePixBegs, rangePixLens, bpPerPix)
   813     for seqName, group in itertools.groupby(a, itemgetter(0)):
   814         yield seqName, sorted(i[1] for i in group)
   815 
   816 def getFont(opts):
   817     if opts.fontfile:
   818         return ImageFont.truetype(opts.fontfile, opts.fontsize)
   819     fileNames = []
   820     try:
   821         x = ["fc-match", "-f%{file}", "arial"]
   822         p = subprocess.Popen(x, stdout=subprocess.PIPE, stderr=subprocess.PIPE,
   823                              universal_newlines=True)
   824         out, err = p.communicate()
   825         fileNames.append(out)
   826     except OSError as e:
   827         logging.info("fc-match error: " + str(e))
   828     fileNames.append("/Library/Fonts/Arial.ttf")  # for Mac
   829     for i in fileNames:
   830         try:
   831             font = ImageFont.truetype(i, opts.fontsize)
   832             logging.info("font: " + i)
   833             return font
   834         except IOError as e:
   835             logging.info("font load error: " + str(e))
   836     return ImageFont.load_default()
   837 
   838 def sequenceSizesAndNames(seqRanges):
   839     for seqName, ranges in itertools.groupby(seqRanges, itemgetter(0)):
   840         size = sum(e - b for n, b, e in ranges)
   841         yield size, seqName
   842 
   843 def biggestSequences(seqRanges, maxNumOfSequences):
   844     s = sorted(sequenceSizesAndNames(seqRanges), reverse=True)
   845     if len(s) > maxNumOfSequences:
   846         logging.warning("too many sequences - discarding the smallest ones")
   847         s = s[:maxNumOfSequences]
   848     return set(i[1] for i in s)
   849 
   850 def remainingSequenceRanges(seqRanges, alignments, seqIndex):
   851     remainingSequences = set(i[seqIndex] for i in alignments)
   852     return [i for i in seqRanges if i[0] in remainingSequences]
   853 
   854 def readAnnotations(opts, font, textDraw, sortedRanges,
   855                     bedFile, repFile, geneFile, gapFile):
   856     rangeDict = expandedSeqDict(dict(rangesPerSeq(sortedRanges)))
   857     annots = sorted(itertools.chain(readBed(bedFile, rangeDict),
   858                                     readRmsk(repFile, rangeDict),
   859                                     readGenePred(opts, geneFile, rangeDict),
   860                                     readGaps(opts, gapFile, rangeDict)))
   861     names = set(i[5] for i in annots)
   862     textSizes = dict(sizesPerText(names, font, textDraw))
   863     margin = max(i[0] for i in textSizes.values()) if textSizes else 0
   864     return annots, textSizes, margin
   865 
   866 def lastDotplot(opts, args):
   867     logLevel = logging.INFO if opts.verbose else logging.WARNING
   868     logging.basicConfig(format="%(filename)s: %(message)s", level=logLevel)
   869 
   870     font = getFont(opts)
   871     image_mode = 'RGB'
   872     forward_color = ImageColor.getcolor(opts.forwardcolor, image_mode)
   873     reverse_color = ImageColor.getcolor(opts.reversecolor, image_mode)
   874     zipped_colors = zip(forward_color, reverse_color)
   875     overlap_color = tuple([(i + j) // 2 for i, j in zipped_colors])
   876 
   877     maxGap1, maxGapB1 = twoValuesFromOption(opts.max_gap1, ":")
   878     maxGap2, maxGapB2 = twoValuesFromOption(opts.max_gap2, ":")
   879 
   880     logging.info("reading alignments...")
   881     alnData = readAlignments(args[0], opts)
   882     alignments, seqRanges1, coverDict1, seqRanges2, coverDict2 = alnData
   883     if not alignments: raise Exception("there are no alignments")
   884     logging.info("cutting...")
   885     coverDict1 = dict(mergedRangesPerSeq(coverDict1))
   886     coverDict2 = dict(mergedRangesPerSeq(coverDict2))
   887     minAlignedBases = min(coveredLength(coverDict1), coveredLength(coverDict2))
   888     pad = int(opts.pad * minAlignedBases)
   889     cutRanges1 = list(trimmed(seqRanges1, coverDict1, minAlignedBases,
   890                               maxGap1, pad, pad))
   891     cutRanges2 = list(trimmed(seqRanges2, coverDict2, minAlignedBases,
   892                               maxGap2, pad, pad))
   893 
   894     biggestSeqs1 = biggestSequences(cutRanges1, opts.maxseqs)
   895     cutRanges1 = [i for i in cutRanges1 if i[0] in biggestSeqs1]
   896     alignments = [i for i in alignments if i[0] in biggestSeqs1]
   897     cutRanges2 = remainingSequenceRanges(cutRanges2, alignments, 1)
   898 
   899     biggestSeqs2 = biggestSequences(cutRanges2, opts.maxseqs)
   900     cutRanges2 = [i for i in cutRanges2 if i[0] in biggestSeqs2]
   901     alignments = [i for i in alignments if i[1] in biggestSeqs2]
   902     cutRanges1 = remainingSequenceRanges(cutRanges1, alignments, 0)
   903 
   904     logging.info("reading secondary alignments...")
   905     alnDataB = readSecondaryAlignments(opts, cutRanges1, cutRanges2)
   906     alignmentsB, seqRangesB1, coverDictB1, seqRangesB2, coverDictB2 = alnDataB
   907     logging.info("cutting...")
   908     coverDictB1 = dict(mergedRangesPerSeq(coverDictB1))
   909     coverDictB2 = dict(mergedRangesPerSeq(coverDictB2))
   910     cutRangesB1 = trimmed(seqRangesB1, coverDictB1, minAlignedBases,
   911                           maxGapB1, 0, 0)
   912     cutRangesB2 = trimmed(seqRangesB2, coverDictB2, minAlignedBases,
   913                           maxGapB2, 0, 0)
   914 
   915     logging.info("sorting...")
   916     sortOut = allSortedRanges(opts, alignments, alignmentsB,
   917                               cutRanges1, cutRangesB1, cutRanges2, cutRangesB2)
   918     sortedRanges1, sortedRanges2 = sortOut
   919 
   920     textDraw = None
   921     if opts.fontsize:
   922         textDraw = ImageDraw.Draw(Image.new(image_mode, (1, 1)))
   923 
   924     textRot1 = "vertical".startswith(opts.rot1)
   925     i1 = dataFromRanges(sortedRanges1, font, textDraw, opts.labels1, textRot1)
   926     rangeSizes1, labelData1, tMargin = i1
   927 
   928     textRot2 = "horizontal".startswith(opts.rot2)
   929     i2 = dataFromRanges(sortedRanges2, font, textDraw, opts.labels2, textRot2)
   930     rangeSizes2, labelData2, lMargin = i2
   931 
   932     logging.info("reading annotations...")
   933 
   934     a1 = readAnnotations(opts, font, textDraw, sortedRanges1,
   935                          opts.bed1, opts.rmsk1, opts.genePred1, opts.gap1)
   936     annots1, annoTextSizes1, bMargin = a1
   937 
   938     a2 = readAnnotations(opts, font, textDraw, sortedRanges2,
   939                          opts.bed2, opts.rmsk2, opts.genePred2, opts.gap2)
   940     annots2, annoTextSizes2, rMargin = a2
   941 
   942     maxPixels1 = opts.width  - lMargin - rMargin
   943     maxPixels2 = opts.height - tMargin - bMargin
   944     bpPerPix1 = get_bp_per_pix(rangeSizes1, opts.border_pixels, maxPixels1)
   945     bpPerPix2 = get_bp_per_pix(rangeSizes2, opts.border_pixels, maxPixels2)
   946     bpPerPix = max(bpPerPix1, bpPerPix2)
   947     logging.info("bp per pixel = " + str(bpPerPix))
   948 
   949     p1 = pixelData(rangeSizes1, bpPerPix, opts.border_pixels, lMargin)
   950     rangePixBegs1, rangePixLens1, rMarginBeg = p1
   951     width = rMarginBeg + rMargin
   952     rangeDict1 = dict(rangesAndOriginsPerSeq(sortedRanges1, rangePixBegs1,
   953                                              rangePixLens1, bpPerPix))
   954 
   955     p2 = pixelData(rangeSizes2, bpPerPix, opts.border_pixels, tMargin)
   956     rangePixBegs2, rangePixLens2, bMarginBeg = p2
   957     height = bMarginBeg + bMargin
   958     rangeDict2 = dict(rangesAndOriginsPerSeq(sortedRanges2, rangePixBegs2,
   959                                              rangePixLens2, bpPerPix))
   960 
   961     logging.info("width:  " + str(width))
   962     logging.info("height: " + str(height))
   963 
   964     logging.info("processing alignments...")
   965     allAlignments = alignments + alignmentsB
   966     hits = alignmentPixels(width, height, allAlignments, bpPerPix,
   967                            rangeDict1, rangeDict2)
   968 
   969     rangeDict1 = expandedSeqDict(rangeDict1)
   970     rangeDict2 = expandedSeqDict(rangeDict2)
   971 
   972     boxes1 = list(bedBoxes(annots1, rangeDict1, bMarginBeg, True, bpPerPix,
   973                            annoTextSizes1))
   974     boxes2 = list(bedBoxes(annots2, rangeDict2, rMarginBeg, False, bpPerPix,
   975                            annoTextSizes2))
   976     boxes = sorted(itertools.chain(boxes1, boxes2))
   977 
   978     logging.info("drawing...")
   979 
   980     image_size = width, height
   981     im = Image.new(image_mode, image_size, opts.background_color)
   982 
   983     drawAnnotations(im, boxes, tMargin, bMarginBeg, lMargin, rMarginBeg)
   984 
   985     joinA, joinB = twoValuesFromOption(opts.join, ":")
   986     if joinA in "13":
   987         drawJoins(im, alignments, bpPerPix, 0, rangeDict1, rangeDict2)
   988     if joinB in "13":
   989         drawJoins(im, alignmentsB, bpPerPix, 0, rangeDict1, rangeDict2)
   990     if joinA in "23":
   991         drawJoins(im, alignments, bpPerPix, 1, rangeDict2, rangeDict1)
   992     if joinB in "23":
   993         drawJoins(im, alignmentsB, bpPerPix, 1, rangeDict2, rangeDict1)
   994 
   995     for i in range(height):
   996         for j in range(width):
   997             store_value = hits[i * width + j]
   998             xy = j, i
   999             if   store_value == 1: im.putpixel(xy, forward_color)
  1000             elif store_value == 2: im.putpixel(xy, reverse_color)
  1001             elif store_value == 3: im.putpixel(xy, overlap_color)
  1002 
  1003     if opts.fontsize != 0:
  1004         axis1 = axisImage(labelData1, rangePixBegs1, rangePixLens1,
  1005                           textRot1, False, font, image_mode, opts)
  1006         if textRot1:
  1007             axis1 = axis1.transpose(Image.ROTATE_90)
  1008         axis2 = axisImage(labelData2, rangePixBegs2, rangePixLens2,
  1009                           textRot2, textRot2, font, image_mode, opts)
  1010         if not textRot2:
  1011             axis2 = axis2.transpose(Image.ROTATE_270)
  1012         im.paste(axis1, (0, 0))
  1013         im.paste(axis2, (0, 0))
  1014 
  1015         annoImage1 = annoTextImage(opts, image_mode, font, bMargin, width,
  1016                                    boxes1, False)
  1017         annoImage1 = annoImage1.transpose(Image.ROTATE_90)
  1018         annoImage2 = annoTextImage(opts, image_mode, font, rMargin, height,
  1019                                    boxes2, True)
  1020         im.paste(annoImage1, (0, bMarginBeg))
  1021         im.paste(annoImage2, (rMarginBeg, 0))
  1022 
  1023     for i in rangePixBegs1[1:]:
  1024         box = i - opts.border_pixels, tMargin, i, bMarginBeg
  1025         im.paste(opts.border_color, box)
  1026 
  1027     for i in rangePixBegs2[1:]:
  1028         box = lMargin, i - opts.border_pixels, rMarginBeg, i
  1029         im.paste(opts.border_color, box)
  1030 
  1031     im.save(args[1])
  1032 
  1033 if __name__ == "__main__":
  1034     usage = """%prog --help
  1035    or: %prog [options] maf-or-tab-alignments dotplot.png
  1036    or: %prog [options] maf-or-tab-alignments dotplot.gif
  1037    or: ..."""
  1038     description = "Draw a dotplot of pair-wise sequence alignments in MAF or tabular format."
  1039     op = optparse.OptionParser(usage=usage, description=description)
  1040     op.add_option("-v", "--verbose", action="count",
  1041                   help="show progress messages & data about the plot")
  1042     # Replace "width" & "height" with a single "length" option?
  1043     op.add_option("-x", "--width", metavar="INT", type="int", default=1000,
  1044                   help="maximum width in pixels (default: %default)")
  1045     op.add_option("-y", "--height", metavar="INT", type="int", default=1000,
  1046                   help="maximum height in pixels (default: %default)")
  1047     op.add_option("-m", "--maxseqs", type="int", default=100, metavar="M",
  1048                   help="maximum number of horizontal or vertical sequences "
  1049                   "(default=%default)")
  1050     op.add_option("-1", "--seq1", metavar="PATTERN", action="append",
  1051                   default=[],
  1052                   help="which sequences to show from the 1st genome")
  1053     op.add_option("-2", "--seq2", metavar="PATTERN", action="append",
  1054                   default=[],
  1055                   help="which sequences to show from the 2nd genome")
  1056     op.add_option("-c", "--forwardcolor", metavar="COLOR", default="red",
  1057                   help="color for forward alignments (default: %default)")
  1058     op.add_option("-r", "--reversecolor", metavar="COLOR", default="blue",
  1059                   help="color for reverse alignments (default: %default)")
  1060     op.add_option("--alignments", metavar="FILE", help="secondary alignments")
  1061     op.add_option("--sort1", default="1", metavar="N",
  1062                   help="genome1 sequence order: 0=input order, 1=name order, "
  1063                   "2=length order, 3=alignment order (default=%default)")
  1064     op.add_option("--sort2", default="1", metavar="N",
  1065                   help="genome2 sequence order: 0=input order, 1=name order, "
  1066                   "2=length order, 3=alignment order (default=%default)")
  1067     op.add_option("--strands1", default="0", metavar="N", help=
  1068                   "genome1 sequence orientation: 0=forward orientation, "
  1069                   "1=alignment orientation (default=%default)")
  1070     op.add_option("--strands2", default="0", metavar="N", help=
  1071                   "genome2 sequence orientation: 0=forward orientation, "
  1072                   "1=alignment orientation (default=%default)")
  1073     op.add_option("--max-gap1", metavar="FRAC", default="0.5,2", help=
  1074                   "maximum unaligned (end,mid) gap in genome1: "
  1075                   "fraction of aligned length (default=%default)")
  1076     op.add_option("--max-gap2", metavar="FRAC", default="0.5,2", help=
  1077                   "maximum unaligned (end,mid) gap in genome2: "
  1078                   "fraction of aligned length (default=%default)")
  1079     op.add_option("--pad", metavar="FRAC", type="float", default=0.04, help=
  1080                   "pad length when cutting unaligned gaps: "
  1081                   "fraction of aligned length (default=%default)")
  1082     op.add_option("-j", "--join", default="0", metavar="N", help=
  1083                   "join: 0=nothing, 1=alignments adjacent in genome1, "
  1084                   "2=alignments adjacent in genome2 (default=%default)")
  1085     op.add_option("--border-pixels", metavar="INT", type="int", default=1,
  1086                   help="number of pixels between sequences (default=%default)")
  1087     op.add_option("--border-color", metavar="COLOR", default="black",
  1088                   help="color for pixels between sequences (default=%default)")
  1089     # --break-color and/or --break-pixels for intra-sequence breaks?
  1090     op.add_option("--margin-color", metavar="COLOR", default="#dcdcdc",
  1091                   help="margin color")
  1092 
  1093     og = optparse.OptionGroup(op, "Text options")
  1094     og.add_option("-f", "--fontfile", metavar="FILE",
  1095                   help="TrueType or OpenType font file")
  1096     og.add_option("-s", "--fontsize", metavar="SIZE", type="int", default=14,
  1097                   help="TrueType or OpenType font size (default: %default)")
  1098     og.add_option("--labels1", type="int", default=0, metavar="N", help=
  1099                   "genome1 labels: 0=name, 1=name:length, "
  1100                   "2=name:start:length, 3=name:start-end (default=%default)")
  1101     og.add_option("--labels2", type="int", default=0, metavar="N", help=
  1102                   "genome2 labels: 0=name, 1=name:length, "
  1103                   "2=name:start:length, 3=name:start-end (default=%default)")
  1104     og.add_option("--rot1", metavar="ROT", default="h",
  1105                   help="text rotation for the 1st genome (default=%default)")
  1106     og.add_option("--rot2", metavar="ROT", default="v",
  1107                   help="text rotation for the 2nd genome (default=%default)")
  1108     op.add_option_group(og)
  1109 
  1110     og = optparse.OptionGroup(op, "Annotation options")
  1111     og.add_option("--bed1", metavar="FILE",
  1112                   help="read genome1 annotations from BED file")
  1113     og.add_option("--bed2", metavar="FILE",
  1114                   help="read genome2 annotations from BED file")
  1115     og.add_option("--rmsk1", metavar="FILE", help="read genome1 repeats from "
  1116                   "RepeatMasker .out or rmsk.txt file")
  1117     og.add_option("--rmsk2", metavar="FILE", help="read genome2 repeats from "
  1118                   "RepeatMasker .out or rmsk.txt file")
  1119     op.add_option_group(og)
  1120 
  1121     og = optparse.OptionGroup(op, "Gene options")
  1122     og.add_option("--genePred1", metavar="FILE",
  1123                   help="read genome1 genes from genePred file")
  1124     og.add_option("--genePred2", metavar="FILE",
  1125                   help="read genome2 genes from genePred file")
  1126     og.add_option("--exon-color", metavar="COLOR", default="PaleGreen",
  1127                   help="color for exons (default=%default)")
  1128     og.add_option("--cds-color", metavar="COLOR", default="LimeGreen",
  1129                   help="color for protein-coding regions (default=%default)")
  1130     op.add_option_group(og)
  1131 
  1132     og = optparse.OptionGroup(op, "Unsequenced gap options")
  1133     og.add_option("--gap1", metavar="FILE",
  1134                   help="read genome1 unsequenced gaps from agp or gap file")
  1135     og.add_option("--gap2", metavar="FILE",
  1136                   help="read genome2 unsequenced gaps from agp or gap file")
  1137     og.add_option("--bridged-color", metavar="COLOR", default="yellow",
  1138                   help="color for bridged gaps (default: %default)")
  1139     og.add_option("--unbridged-color", metavar="COLOR", default="orange",
  1140                   help="color for unbridged gaps (default: %default)")
  1141     op.add_option_group(og)
  1142     (opts, args) = op.parse_args()
  1143     if len(args) != 2: op.error("2 arguments needed")
  1144 
  1145     opts.background_color = "white"
  1146     opts.label_space = 5     # minimum number of pixels between axis labels
  1147 
  1148     try: lastDotplot(opts, args)
  1149     except KeyboardInterrupt: pass  # avoid silly error message
  1150     except Exception as e:
  1151         prog = os.path.basename(sys.argv[0])
  1152         sys.exit(prog + ": error: " + str(e))