#! /usr/bin/env python

 # Read pair-wise alignments in MAF or LAST tabular format: write an

 # "Oxford grid", a.k.a. dotplot.

 # TODO: Currently, pixels with zero aligned nt-pairs are white, and

 # pixels with one or more aligned nt-pairs are black.  This can look

 # too crowded for large genome alignments.  I tried shading each pixel

 # according to the number of aligned nt-pairs within it, but the

 # result is too faint.  How can this be done better?

 import collections

 import functools

 import gzip

 from fnmatch import fnmatchcase

 from operator import itemgetter

 import subprocess

 import itertools, optparse, os, re, sys

 # Try to make PIL/PILLOW work:

 try:

     from PIL import Image, ImageDraw, ImageFont, ImageColor

 except ImportError:

     import Image, ImageDraw, ImageFont, ImageColor

 try:

     from future_builtins import zip

 except ImportError:

     pass

 def myOpen(fileName):  # faster than fileinput

     if fileName is None:

         return []

     if fileName == "-":

         return sys.stdin

     if fileName.endswith(".gz"):

         return gzip.open(fileName, "rt")  # xxx dubious for Python2

     return open(fileName)

 def warn(message):

     if opts.verbose:

         prog = os.path.basename(sys.argv[0])

         sys.stderr.write(prog + ": " + message + "\n")

 def groupByFirstItem(things):

     for k, v in itertools.groupby(things, itemgetter(0)):

         yield k, [i[1:] for i in v]

 def croppedBlocks(blocks, ranges1, ranges2):

     headBeg1, headBeg2, headSize = blocks[0]

     for r1 in ranges1:

         for r2 in ranges2:

             cropBeg1, cropEnd1 = r1

             if headBeg1 < 0:

                 cropBeg1, cropEnd1 = -cropEnd1, -cropBeg1

             cropBeg2, cropEnd2 = r2

             if headBeg2 < 0:

                 cropBeg2, cropEnd2 = -cropEnd2, -cropBeg2

             for beg1, beg2, size in blocks:

                 b1 = max(cropBeg1, beg1)

                 e1 = min(cropEnd1, beg1 + size)

                 if b1 >= e1: continue

                 offset = beg2 - beg1

                 b2 = max(cropBeg2, b1 + offset)

                 e2 = min(cropEnd2, e1 + offset)

                 if b2 >= e2: continue

                 yield b2 - offset, b2, e2 - b2

 def tabBlocks(beg1, beg2, blocks):

     '''Get the gapless blocks of an alignment, from LAST tabular format.'''

     for i in blocks.split(","):

         if ":" in i:

             x, y = i.split(":")

             beg1 += int(x)

             beg2 += int(y)

         else:

             size = int(i)

             yield beg1, beg2, size

             beg1 += size

             beg2 += size

 def mafBlocks(beg1, beg2, seq1, seq2):

     '''Get the gapless blocks of an alignment, from MAF format.'''

     size = 0

     for x, y in zip(seq1, seq2):

         if x == "-":

             if size:

                 yield beg1, beg2, size

                 beg1 += size

                 beg2 += size

                 size = 0

             beg2 += 1

         elif y == "-":

             if size:

                 yield beg1, beg2, size

                 beg1 += size

                 beg2 += size

                 size = 0

             beg1 += 1

         else:

             size += 1

     if size: yield beg1, beg2, size

 def alignmentInput(lines):

     '''Get alignments and sequence lengths, from MAF or tabular format.'''

     mafCount = 0

     for line in lines:

         w = line.split()

         if line[0].isdigit():  # tabular format

             chr1, beg1, seqlen1 = w[1], int(w[2]), int(w[5])

             if w[4] == "-": beg1 -= seqlen1

             chr2, beg2, seqlen2 = w[6], int(w[7]), int(w[10])

             if w[9] == "-": beg2 -= seqlen2

             blocks = tabBlocks(beg1, beg2, w[11])

             yield chr1, seqlen1, chr2, seqlen2, blocks

         elif line[0] == "s":  # MAF format

             if mafCount == 0:

                 chr1, beg1, seqlen1, seq1 = w[1], int(w[2]), int(w[5]), w[6]

                 if w[4] == "-": beg1 -= seqlen1

                 mafCount = 1

             else:

                 chr2, beg2, seqlen2, seq2 = w[1], int(w[2]), int(w[5]), w[6]

                 if w[4] == "-": beg2 -= seqlen2

                 blocks = mafBlocks(beg1, beg2, seq1, seq2)

                 yield chr1, seqlen1, chr2, seqlen2, blocks

                 mafCount = 0

 def seqRequestFromText(text):

     if ":" in text:

         pattern, interval = text.rsplit(":", 1)

         if "-" in interval:

             beg, end = interval.rsplit("-", 1)

             return pattern, int(beg), int(end)  # beg may be negative

     return text, 0, sys.maxsize

 def rangesFromSeqName(seqRequests, name, seqLen):

     if seqRequests:

         base = name.split(".")[-1]  # allow for names like hg19.chr7

         for pat, beg, end in seqRequests:

             if fnmatchcase(name, pat) or fnmatchcase(base, pat):

                 yield max(beg, 0), min(end, seqLen)

     else:

         yield 0, seqLen

 def updateSeqs(coverDict, seqRanges, seqName, ranges, coveredRange):

     beg, end = coveredRange

     if beg < 0:

         coveredRange = -end, -beg

     if seqName in coverDict:

         coverDict[seqName].append(coveredRange)

     else:

         coverDict[seqName] = [coveredRange]

         for beg, end in ranges:

             r = seqName, beg, end

             seqRanges.append(r)

 def readAlignments(fileName, opts):

     '''Get alignments and sequence limits, from MAF or tabular format.'''

     seqRequests1 = [seqRequestFromText(i) for i in opts.seq1]

     seqRequests2 = [seqRequestFromText(i) for i in opts.seq2]

     alignments = []

     seqRanges1 = []

     seqRanges2 = []

     coverDict1 = {}

     coverDict2 = {}

     lines = myOpen(fileName)

     for seqName1, seqLen1, seqName2, seqLen2, blocks in alignmentInput(lines):

         ranges1 = sorted(rangesFromSeqName(seqRequests1, seqName1, seqLen1))

         if not ranges1: continue

         ranges2 = sorted(rangesFromSeqName(seqRequests2, seqName2, seqLen2))

         if not ranges2: continue

         b = list(croppedBlocks(list(blocks), ranges1, ranges2))

         if not b: continue

         aln = seqName1, seqName2, b

         alignments.append(aln)

         coveredRange1 = b[0][0], b[-1][0] + b[-1][2]

         updateSeqs(coverDict1, seqRanges1, seqName1, ranges1, coveredRange1)

         coveredRange2 = b[0][1], b[-1][1] + b[-1][2]

         updateSeqs(coverDict2, seqRanges2, seqName2, ranges2, coveredRange2)

     return alignments, seqRanges1, coverDict1, seqRanges2, coverDict2

 def nameAndRangesFromDict(cropDict, seqName):

     if seqName in cropDict:

         return seqName, cropDict[seqName]

     n = seqName.split(".")[-1]

     if n in cropDict:

         return n, cropDict[n]

     return seqName, []

 def rangesForSecondaryAlignments(primaryRanges, seqLen):

     if primaryRanges:

         return primaryRanges

     return [(0, seqLen)]

 def readSecondaryAlignments(opts, cropRanges1, cropRanges2):

     cropDict1 = dict(groupByFirstItem(cropRanges1))

     cropDict2 = dict(groupByFirstItem(cropRanges2))

     alignments = []

     seqRanges1 = []

     seqRanges2 = []

     coverDict1 = {}

     coverDict2 = {}

     lines = myOpen(opts.alignments)

     for seqName1, seqLen1, seqName2, seqLen2, blocks in alignmentInput(lines):

         seqName1, ranges1 = nameAndRangesFromDict(cropDict1, seqName1)

         seqName2, ranges2 = nameAndRangesFromDict(cropDict2, seqName2)

         if not ranges1 and not ranges2:

             continue

         r1 = rangesForSecondaryAlignments(ranges1, seqLen1)

         r2 = rangesForSecondaryAlignments(ranges2, seqLen2)

         b = list(croppedBlocks(list(blocks), r1, r2))

         if not b: continue

         aln = seqName1, seqName2, b

         alignments.append(aln)

         if not ranges1:

             coveredRange1 = b[0][0], b[-1][0] + b[-1][2]

             updateSeqs(coverDict1, seqRanges1, seqName1, r1, coveredRange1)

         if not ranges2:

             coveredRange2 = b[0][1], b[-1][1] + b[-1][2]

             updateSeqs(coverDict2, seqRanges2, seqName2, r2, coveredRange2)

     return alignments, seqRanges1, coverDict1, seqRanges2, coverDict2

 def twoValuesFromOption(text, separator):

     if separator in text:

         return text.split(separator)

     return text, text

 def mergedRanges(ranges):

     oldBeg, maxEnd = ranges[0]

     for beg, end in ranges:

         if beg > maxEnd:

             yield oldBeg, maxEnd

             oldBeg = beg

             maxEnd = end

         elif end > maxEnd:

             maxEnd = end

     yield oldBeg, maxEnd

 def mergedRangesPerSeq(coverDict):

     for k, v in coverDict.items():

         v.sort()

         yield k, list(mergedRanges(v))

 def coveredLength(mergedCoverDict):

     return sum(sum(e - b for b, e in v) for v in mergedCoverDict.values())

 def trimmed(seqRanges, coverDict, minAlignedBases, maxGapFrac, endPad, midPad):

     maxEndGapFrac, maxMidGapFrac = twoValuesFromOption(maxGapFrac, ",")

     maxEndGap = max(float(maxEndGapFrac) * minAlignedBases, endPad * 1.0)

     maxMidGap = max(float(maxMidGapFrac) * minAlignedBases, midPad * 2.0)

     for seqName, rangeBeg, rangeEnd in seqRanges:

         seqBlocks = coverDict[seqName]

         blocks = [i for i in seqBlocks if i[0] < rangeEnd and i[1] > rangeBeg]

         if blocks[0][0] - rangeBeg > maxEndGap:

             rangeBeg = blocks[0][0] - endPad

         for j, y in enumerate(blocks):

             if j:

                 x = blocks[j - 1]

                 if y[0] - x[1] > maxMidGap:

                     yield seqName, rangeBeg, x[1] + midPad

                     rangeBeg = y[0] - midPad

         if rangeEnd - blocks[-1][1] > maxEndGap:

             rangeEnd = blocks[-1][1] + endPad

         yield seqName, rangeBeg, rangeEnd

 def rangesWithStrandInfo(seqRanges, strandOpt, alignments, seqIndex):

     if strandOpt == "1":

         forwardMinusReverse = collections.defaultdict(int)

         for i in alignments:

             blocks = i[2]

             beg1, beg2, size = blocks[0]

             numOfAlignedLetterPairs = sum(i[2] for i in blocks)

             if (beg1 < 0) != (beg2 < 0):  # opposite-strand alignment

                 numOfAlignedLetterPairs *= -1

             forwardMinusReverse[i[seqIndex]] += numOfAlignedLetterPairs

     strandNum = 0

     for seqName, beg, end in seqRanges:

         if strandOpt == "1":

             strandNum = 1 if forwardMinusReverse[seqName] >= 0 else 2

         yield seqName, beg, end, strandNum

 def natural_sort_key(my_string):

     '''Return a sort key for "natural" ordering, e.g. chr9 < chr10.'''

     parts = re.split(r'(\d+)', my_string)

     parts[1::2] = map(int, parts[1::2])

     return parts

 def nameKey(oneSeqRanges):

     return natural_sort_key(oneSeqRanges[0][0])

 def sizeKey(oneSeqRanges):

     return sum(b - e for n, b, e, s in oneSeqRanges), nameKey(oneSeqRanges)

 def alignmentKey(seqNamesToLists, oneSeqRanges):

     seqName = oneSeqRanges[0][0]

     alignmentsOfThisSequence = seqNamesToLists[seqName]

     numOfAlignedLetterPairs = sum(i[3] for i in alignmentsOfThisSequence)

     toMiddle = numOfAlignedLetterPairs // 2

     for i in alignmentsOfThisSequence:

         toMiddle -= i[3]

         if toMiddle < 0:

             return i[1:3]  # sequence-rank and "position" of this alignment

 def rankAndFlipPerSeq(seqRanges):

     rangesGroupedBySeqName = itertools.groupby(seqRanges, itemgetter(0))

     for rank, group in enumerate(rangesGroupedBySeqName):

         seqName, ranges = group

         strandNum = next(ranges)[3]

         flip = 1 if strandNum < 2 else -1

         yield seqName, (rank, flip)

 def alignmentSortData(alignments, seqIndex, otherNamesToRanksAndFlips):

     otherIndex = 1 - seqIndex

     for i in alignments:

         blocks = i[2]

         otherRank, otherFlip = otherNamesToRanksAndFlips[i[otherIndex]]

         otherPos = otherFlip * abs(blocks[0][otherIndex] +

                                    blocks[-1][otherIndex] + blocks[-1][2])

         numOfAlignedLetterPairs = sum(i[2] for i in blocks)

         yield i[seqIndex], otherRank, otherPos, numOfAlignedLetterPairs

 def mySortedRanges(seqRanges, sortOpt, seqIndex, alignments, otherRanges):

     rangesGroupedBySeqName = itertools.groupby(seqRanges, itemgetter(0))

     g = [list(ranges) for seqName, ranges in rangesGroupedBySeqName]

     for i in g:

         if i[0][3] > 1:

             i.reverse()

     if sortOpt == "1":

         g.sort(key=nameKey)

     if sortOpt == "2":

         g.sort(key=sizeKey)

     if sortOpt == "3":

         otherNamesToRanksAndFlips = dict(rankAndFlipPerSeq(otherRanges))

         alns = sorted(alignmentSortData(alignments, seqIndex,

                                         otherNamesToRanksAndFlips))

         alnsGroupedBySeqName = itertools.groupby(alns, itemgetter(0))

         seqNamesToLists = dict((k, list(v)) for k, v in alnsGroupedBySeqName)

         g.sort(key=functools.partial(alignmentKey, seqNamesToLists))

     return [j for i in g for j in i]

 def allSortedRanges(opts, alignments, alignmentsB,

                     seqRanges1, seqRangesB1, seqRanges2, seqRangesB2):

     o1, oB1 = twoValuesFromOption(opts.strands1, ":")

     o2, oB2 = twoValuesFromOption(opts.strands2, ":")

     if o1 == "1" and o2 == "1":

         raise Exception("the strand options have circular dependency")

     seqRanges1 = list(rangesWithStrandInfo(seqRanges1, o1, alignments, 0))

     seqRanges2 = list(rangesWithStrandInfo(seqRanges2, o2, alignments, 1))

     seqRangesB1 = list(rangesWithStrandInfo(seqRangesB1, oB1, alignmentsB, 0))

     seqRangesB2 = list(rangesWithStrandInfo(seqRangesB2, oB2, alignmentsB, 1))

     o1, oB1 = twoValuesFromOption(opts.sort1, ":")

     o2, oB2 = twoValuesFromOption(opts.sort2, ":")

     if o1 == "3" and o2 == "3":

         raise Exception("the sort options have circular dependency")

     if o1 != "3":

         s1 = mySortedRanges(seqRanges1, o1, None, None, None)

     if o2 != "3":

         s2 = mySortedRanges(seqRanges2, o2, None, None, None)

     if o1 == "3":

         s1 = mySortedRanges(seqRanges1, o1, 0, alignments, s2)

     if o2 == "3":

         s2 = mySortedRanges(seqRanges2, o2, 1, alignments, s1)

     t1 = mySortedRanges(seqRangesB1, oB1, 0, alignmentsB, s2)

     t2 = mySortedRanges(seqRangesB2, oB2, 1, alignmentsB, s1)

     return s1 + t1, s2 + t2

 def prettyNum(n):

     t = str(n)

     groups = []

     while t:

         groups.append(t[-3:])

         t = t[:-3]

     return ",".join(reversed(groups))

 def sizeText(size):

     suffixes = "bp", "kb", "Mb", "Gb"

     for i, x in enumerate(suffixes):

         j = 10 ** (i * 3)

         if size < j * 10:

             return "%.2g" % (1.0 * size / j) + x

         if size < j * 1000 or i == len(suffixes) - 1:

             return "%.0f" % (1.0 * size / j) + x

 def labelText(seqRange, labelOpt):

     seqName, beg, end, strandNum = seqRange

     if labelOpt == 1:

         return seqName + ": " + sizeText(end - beg)

     if labelOpt == 2:

         return seqName + ":" + prettyNum(beg) + ": " + sizeText(end - beg)

     if labelOpt == 3:

         return seqName + ":" + prettyNum(beg) + "-" + prettyNum(end)

     return seqName

 def rangeLabels(seqRanges, labelOpt, font, fontsize, image_mode, textRot):

     if fontsize:

         image_size = 1, 1

         im = Image.new(image_mode, image_size)

         draw = ImageDraw.Draw(im)

     x = y = 0

     for r in seqRanges:

         text = labelText(r, labelOpt)

         if fontsize:

             x, y = draw.textsize(text, font=font)

             if textRot:

                 x, y = y, x

         yield text, x, y, r[3]

 def dataFromRanges(sortedRanges, font, fontSize, imageMode, labelOpt, textRot):

     for seqName, rangeBeg, rangeEnd, strandNum in sortedRanges:

         out = [seqName, str(rangeBeg), str(rangeEnd)]

         if strandNum > 0:

             out.append(".+-"[strandNum])

         warn("\t".join(out))

     warn("")

     rangeSizes = [e - b for n, b, e, s in sortedRanges]

     labs = list(rangeLabels(sortedRanges, labelOpt, font, fontSize,

                             imageMode, textRot))

     margin = max(i[2] for i in labs)

     # xxx the margin may be too big, because some labels may get omitted

     return rangeSizes, labs, margin

 def div_ceil(x, y):

     '''Return x / y rounded up.'''

     q, r = divmod(x, y)

     return q + (r != 0)

 def get_bp_per_pix(rangeSizes, pixTweenRanges, maxPixels):

     '''Get the minimum bp-per-pixel that fits in the size limit.'''

     warn("choosing bp per pixel...")

     numOfRanges = len(rangeSizes)

     maxPixelsInRanges = maxPixels - pixTweenRanges * (numOfRanges - 1)

     if maxPixelsInRanges < numOfRanges:

         raise Exception("can't fit the image: too many sequences?")

     negLimit = -maxPixelsInRanges

     negBpPerPix = sum(rangeSizes) // negLimit

     while True:

         if sum(i // negBpPerPix for i in rangeSizes) >= negLimit:

             return -negBpPerPix

         negBpPerPix -= 1

 def getRangePixBegs(rangePixLens, pixTweenRanges, margin):

     '''Get the start pixel for each range.'''

     rangePixBegs = []

     pix_tot = margin - pixTweenRanges

     for i in rangePixLens:

         pix_tot += pixTweenRanges

         rangePixBegs.append(pix_tot)

         pix_tot += i

     return rangePixBegs

 def pixelData(rangeSizes, bp_per_pix, pixTweenRanges, margin):

     '''Return pixel information about the ranges.'''

     rangePixLens = [div_ceil(i, bp_per_pix) for i in rangeSizes]

     rangePixBegs = getRangePixBegs(rangePixLens, pixTweenRanges, margin)

     tot_pix = rangePixBegs[-1] + rangePixLens[-1]

     return rangePixBegs, rangePixLens, tot_pix

 def drawLineForward(hits, width, bp_per_pix, beg1, beg2, size):

     while True:

         q1, r1 = divmod(beg1, bp_per_pix)

         q2, r2 = divmod(beg2, bp_per_pix)

         hits[q2 * width + q1] |= 1

         next_pix = min(bp_per_pix - r1, bp_per_pix - r2)

         if next_pix >= size: break

         beg1 += next_pix

         beg2 += next_pix

         size -= next_pix

 def drawLineReverse(hits, width, bp_per_pix, beg1, beg2, size):

     while True:

         q1, r1 = divmod(beg1, bp_per_pix)

         q2, r2 = divmod(beg2, bp_per_pix)

         hits[q2 * width + q1] |= 2

         next_pix = min(bp_per_pix - r1, r2 + 1)

         if next_pix >= size: break

         beg1 += next_pix

         beg2 -= next_pix

         size -= next_pix

 def strandAndOrigin(ranges, beg, size):

     isReverseStrand = (beg < 0)

     if isReverseStrand:

         beg = -(beg + size)

     for rangeBeg, rangeEnd, isReverseRange, origin in ranges:

         if rangeEnd > beg:  # assumes the ranges are sorted

             return (isReverseStrand != isReverseRange), origin

 def alignmentPixels(width, height, alignments, bp_per_pix,

                     rangeDict1, rangeDict2):

     hits = [0] * (width * height)  # the image data

     for seq1, seq2, blocks in alignments:

         beg1, beg2, size = blocks[0]

         isReverse1, ori1 = strandAndOrigin(rangeDict1[seq1], beg1, size)

         isReverse2, ori2 = strandAndOrigin(rangeDict2[seq2], beg2, size)

         for beg1, beg2, size in blocks:

             if isReverse1:

                 beg1 = -(beg1 + size)

                 beg2 = -(beg2 + size)

             if isReverse1 == isReverse2:

                 drawLineForward(hits, width, bp_per_pix,

                                 ori1 + beg1, ori2 + beg2, size)

             else:

                 drawLineReverse(hits, width, bp_per_pix,

                                 ori1 + beg1, ori2 - beg2 - 1, size)

     return hits

 def expandedSeqDict(seqDict):

     '''Allow lookup by short sequence names, e.g. chr7 as well as hg19.chr7.'''

     newDict = seqDict.copy()

     for name, x in seqDict.items():

         if "." in name:

             base = name.split(".")[-1]

             if base in newDict:  # an ambiguous case was found:

                 return seqDict   # so give up completely

             newDict[base] = x

     return newDict

 def readBed(fileName, rangeDict):

     for line in myOpen(fileName):

         w = line.split()

         if not w: continue

         seqName = w[0]

         if seqName not in rangeDict: continue

         beg = int(w[1])

         end = int(w[2])

         layer = 900

         color = "#fbf"

         if len(w) > 4:

             if w[4] != ".":

                 layer = float(w[4])

             if len(w) > 5:

                 if len(w) > 8 and w[8].count(",") == 2:

                     color = "rgb(" + w[8] + ")"

                 else:

                     strand = w[5]

                     isRev = rangeDict[seqName][0][2]

                     if strand == "+" and not isRev or strand == "-" and isRev:

                         color = "#ffe8e8"

                     if strand == "-" and not isRev or strand == "+" and isRev:

                         color = "#e8e8ff"

         yield layer, color, seqName, beg, end

 def commaSeparatedInts(text):

     return map(int, text.rstrip(",").split(","))

 def readGenePred(opts, fileName, rangeDict):

     for line in myOpen(fileName):

         fields = line.split()

         if not fields: continue

         if fields[2] not in "+-": fields = fields[1:]

         seqName = fields[1]

         if seqName not in rangeDict: continue

         #strand = fields[2]

         cdsBeg = int(fields[5])

         cdsEnd = int(fields[6])

         exonBegs = commaSeparatedInts(fields[8])

         exonEnds = commaSeparatedInts(fields[9])

         for beg, end in zip(exonBegs, exonEnds):

             yield 300, opts.exon_color, seqName, beg, end

             b = max(beg, cdsBeg)

             e = min(end, cdsEnd)

             if b < e: yield 400, opts.cds_color, seqName, b, e

 def readRmsk(fileName, rangeDict):

     for line in myOpen(fileName):

         fields = line.split()

         if len(fields) == 17:  # rmsk.txt

             seqName = fields[5]

             if seqName not in rangeDict: continue  # do this ASAP for speed

             beg = int(fields[6])

             end = int(fields[7])

             strand = fields[9]

             repeatClass = fields[11]

         elif len(fields) == 15:  # .out

             seqName = fields[4]

             if seqName not in rangeDict: continue

             beg = int(fields[5]) - 1

             end = int(fields[6])

             strand = fields[8]

             repeatClass = fields[10]

         else:

             continue

         if repeatClass in ("Low_complexity", "Simple_repeat"):

             yield 200, "#fbf", seqName, beg, end

         elif (strand == "+") != rangeDict[seqName][0][2]:

             yield 100, "#ffe8e8", seqName, beg, end

         else:

             yield 100, "#e8e8ff", seqName, beg, end

 def isExtraFirstGapField(fields):

     return fields[4].isdigit()

 def readGaps(opts, fileName, rangeDict):

     '''Read locations of unsequenced gaps, from an agp or gap file.'''

     for line in myOpen(fileName):

         w = line.split()

         if not w or w[0][0] == "#": continue

         if isExtraFirstGapField(w): w = w[1:]

         if w[4] not in "NU": continue

         seqName = w[0]

         if seqName not in rangeDict: continue

         end = int(w[2])

         beg = end - int(w[5])  # zero-based coordinate

         if w[7] == "yes":

             yield 3000, opts.bridged_color, seqName, beg, end

         else:

             yield 2000, opts.unbridged_color, seqName, beg, end

 def bedBoxes(beds, rangeDict, margin, edge, isTop, bpPerPix):

     for layer, color, seqName, bedBeg, bedEnd in beds:

         for rangeBeg, rangeEnd, isReverseRange, origin in rangeDict[seqName]:

             beg = max(bedBeg, rangeBeg)

             end = min(bedEnd, rangeEnd)

             if beg >= end: continue

             if isReverseRange:

                 beg, end = -end, -beg

             if layer <= 1000:

                 # include partly-covered pixels

                 b = (origin + beg) // bpPerPix

                 e = div_ceil(origin + end, bpPerPix)

             else:

                 # exclude partly-covered pixels

                 b = div_ceil(origin + beg, bpPerPix)

                 e = (origin + end) // bpPerPix

                 if e <= b: continue

                 if bedEnd >= rangeEnd:  # include partly-covered end pixels

                     if isReverseRange:

                         b = (origin + beg) // bpPerPix

                     else:

                         e = div_ceil(origin + end, bpPerPix)

             if isTop:

                 box = b, margin, e, edge

             else:

                 box = margin, b, edge, e

             yield layer, color, box

 def drawAnnotations(im, boxes):

     # xxx use partial transparency for different-color overlaps?

     for layer, color, box in boxes:

         im.paste(color, box)

 def placedLabels(labels, rangePixBegs, rangePixLens, beg, end):

     '''Return axis labels with endpoint & sort-order information.'''

     maxWidth = end - beg

     for i, j, k in zip(labels, rangePixBegs, rangePixLens):

         text, textWidth, textHeight, strandNum = i

         if textWidth > maxWidth:

             continue

         labelBeg = j + (k - textWidth) // 2

         labelEnd = labelBeg + textWidth

         sortKey = textWidth - k

         if labelBeg < beg:

             sortKey += maxWidth * (beg - labelBeg)

             labelBeg = beg

             labelEnd = beg + textWidth

         if labelEnd > end:

             sortKey += maxWidth * (labelEnd - end)

             labelEnd = end

             labelBeg = end - textWidth

         yield sortKey, labelBeg, labelEnd, text, textHeight, strandNum

 def nonoverlappingLabels(labels, minPixTweenLabels):

     '''Get a subset of non-overlapping axis labels, greedily.'''

     out = []

     for i in labels:

         beg = i[1] - minPixTweenLabels

         end = i[2] + minPixTweenLabels

         if all(j[2] <= beg or j[1] >= end for j in out):

             out.append(i)

     return out

 def axisImage(labels, rangePixBegs, rangePixLens, textRot,

               textAln, font, image_mode, opts):

     '''Make an image of axis labels.'''

     beg = rangePixBegs[0]

     end = rangePixBegs[-1] + rangePixLens[-1]

     margin = max(i[2] for i in labels)

     labels = sorted(placedLabels(labels, rangePixBegs, rangePixLens, beg, end))

     minPixTweenLabels = 0 if textRot else opts.label_space

     labels = nonoverlappingLabels(labels, minPixTweenLabels)

     image_size = (margin, end) if textRot else (end, margin)

     im = Image.new(image_mode, image_size, opts.margin_color)

     draw = ImageDraw.Draw(im)

     for sortKey, labelBeg, labelEnd, text, textHeight, strandNum in labels:

         base = margin - textHeight if textAln else 0

         position = (base, labelBeg) if textRot else (labelBeg, base)

         fill = ("black", opts.forwardcolor, opts.reversecolor)[strandNum]

         draw.text(position, text, font=font, fill=fill)

     return im

 def rangesWithOrigins(sortedRanges, rangePixBegs, rangePixLens, bpPerPix):

     for i, j, k in zip(sortedRanges, rangePixBegs, rangePixLens):

         seqName, rangeBeg, rangeEnd, strandNum = i

         isReverseRange = (strandNum > 1)

         if isReverseRange:

             origin = bpPerPix * (j + k) + rangeBeg

         else:

             origin = bpPerPix * j - rangeBeg

         yield seqName, (rangeBeg, rangeEnd, isReverseRange, origin)

 def rangesPerSeq(sortedRanges, rangePixBegs, rangePixLens, bpPerPix):

     a = rangesWithOrigins(sortedRanges, rangePixBegs, rangePixLens, bpPerPix)

     for k, v in itertools.groupby(a, itemgetter(0)):

         yield k, sorted(i[1] for i in v)

 def getFont(opts):

     if opts.fontfile:

         return ImageFont.truetype(opts.fontfile, opts.fontsize)

     fileNames = []

     try:

         x = ["fc-match", "-f%{file}", "arial"]

         p = subprocess.Popen(x, stdout=subprocess.PIPE, stderr=subprocess.PIPE,

                              universal_newlines=True)

         out, err = p.communicate()

         fileNames.append(out)

     except OSError as e:

         warn("fc-match error: " + str(e))

     fileNames.append("/Library/Fonts/Arial.ttf")  # for Mac

     for i in fileNames:

         try:

             font = ImageFont.truetype(i, opts.fontsize)

             warn("font: " + i)

             return font

         except IOError as e:

             warn("font load error: " + str(e))

     return ImageFont.load_default()

 def lastDotplot(opts, args):

     font = getFont(opts)

     image_mode = 'RGB'

     forward_color = ImageColor.getcolor(opts.forwardcolor, image_mode)

     reverse_color = ImageColor.getcolor(opts.reversecolor, image_mode)

     zipped_colors = zip(forward_color, reverse_color)

     overlap_color = tuple([(i + j) // 2 for i, j in zipped_colors])

     maxGap1, maxGapB1 = twoValuesFromOption(opts.max_gap1, ":")

     maxGap2, maxGapB2 = twoValuesFromOption(opts.max_gap2, ":")

     warn("reading alignments...")

     alnData = readAlignments(args[0], opts)

     alignments, seqRanges1, coverDict1, seqRanges2, coverDict2 = alnData

     if not alignments: raise Exception("there are no alignments")

     warn("cutting...")

     coverDict1 = dict(mergedRangesPerSeq(coverDict1))

     coverDict2 = dict(mergedRangesPerSeq(coverDict2))

     minAlignedBases = min(coveredLength(coverDict1), coveredLength(coverDict2))

     pad = int(opts.pad * minAlignedBases)

     cutRanges1 = list(trimmed(seqRanges1, coverDict1, minAlignedBases,

                               maxGap1, pad, pad))

     cutRanges2 = list(trimmed(seqRanges2, coverDict2, minAlignedBases,

                               maxGap2, pad, pad))

     warn("reading secondary alignments...")

     alnDataB = readSecondaryAlignments(opts, cutRanges1, cutRanges2)

     alignmentsB, seqRangesB1, coverDictB1, seqRangesB2, coverDictB2 = alnDataB

     warn("cutting...")

     coverDictB1 = dict(mergedRangesPerSeq(coverDictB1))

     coverDictB2 = dict(mergedRangesPerSeq(coverDictB2))

     cutRangesB1 = trimmed(seqRangesB1, coverDictB1, minAlignedBases,

                           maxGapB1, 0, 0)

     cutRangesB2 = trimmed(seqRangesB2, coverDictB2, minAlignedBases,

                           maxGapB2, 0, 0)

     warn("sorting...")

     sortOut = allSortedRanges(opts, alignments, alignmentsB,

                               cutRanges1, cutRangesB1, cutRanges2, cutRangesB2)

     sortedRanges1, sortedRanges2 = sortOut

     textRot1 = "vertical".startswith(opts.rot1)

     i1 = dataFromRanges(sortedRanges1, font,

                         opts.fontsize, image_mode, opts.labels1, textRot1)

     rangeSizes1, labelData1, tMargin = i1

     textRot2 = "horizontal".startswith(opts.rot2)

     i2 = dataFromRanges(sortedRanges2, font,

                         opts.fontsize, image_mode, opts.labels2, textRot2)

     rangeSizes2, labelData2, lMargin = i2

     maxPixels1 = opts.width  - lMargin

     maxPixels2 = opts.height - tMargin

     bpPerPix1 = get_bp_per_pix(rangeSizes1, opts.border_pixels, maxPixels1)

     bpPerPix2 = get_bp_per_pix(rangeSizes2, opts.border_pixels, maxPixels2)

     bpPerPix = max(bpPerPix1, bpPerPix2)

     warn("bp per pixel = " + str(bpPerPix))

     p1 = pixelData(rangeSizes1, bpPerPix, opts.border_pixels, lMargin)

     rangePixBegs1, rangePixLens1, width = p1

     rangeDict1 = dict(rangesPerSeq(sortedRanges1, rangePixBegs1,

                                    rangePixLens1, bpPerPix))

     p2 = pixelData(rangeSizes2, bpPerPix, opts.border_pixels, tMargin)

     rangePixBegs2, rangePixLens2, height = p2

     rangeDict2 = dict(rangesPerSeq(sortedRanges2, rangePixBegs2,

                                    rangePixLens2, bpPerPix))

     warn("width:  " + str(width))

     warn("height: " + str(height))

     warn("processing alignments...")

     hits = alignmentPixels(width, height, alignments + alignmentsB, bpPerPix,

                            rangeDict1, rangeDict2)

     warn("reading annotations...")

     rangeDict1 = expandedSeqDict(rangeDict1)

     beds1 = itertools.chain(readBed(opts.bed1, rangeDict1),

                             readRmsk(opts.rmsk1, rangeDict1),

                             readGenePred(opts, opts.genePred1, rangeDict1),

                             readGaps(opts, opts.gap1, rangeDict1))

     b1 = bedBoxes(beds1, rangeDict1, tMargin, height, True, bpPerPix)

     rangeDict2 = expandedSeqDict(rangeDict2)

     beds2 = itertools.chain(readBed(opts.bed2, rangeDict2),

                             readRmsk(opts.rmsk2, rangeDict2),

                             readGenePred(opts, opts.genePred2, rangeDict2),

                             readGaps(opts, opts.gap2, rangeDict2))

     b2 = bedBoxes(beds2, rangeDict2, lMargin, width, False, bpPerPix)

     boxes = sorted(itertools.chain(b1, b2))

     warn("drawing...")

     image_size = width, height

     im = Image.new(image_mode, image_size, opts.background_color)

     drawAnnotations(im, boxes)

     for i in range(height):

         for j in range(width):

             store_value = hits[i * width + j]

             xy = j, i

             if   store_value == 1: im.putpixel(xy, forward_color)

             elif store_value == 2: im.putpixel(xy, reverse_color)

             elif store_value == 3: im.putpixel(xy, overlap_color)

     if opts.fontsize != 0:

         axis1 = axisImage(labelData1, rangePixBegs1, rangePixLens1,

                           textRot1, False, font, image_mode, opts)

         if textRot1:

             axis1 = axis1.transpose(Image.ROTATE_90)

         axis2 = axisImage(labelData2, rangePixBegs2, rangePixLens2,

                           textRot2, textRot2, font, image_mode, opts)

         if not textRot2:

             axis2 = axis2.transpose(Image.ROTATE_270)

         im.paste(axis1, (0, 0))

         im.paste(axis2, (0, 0))

     for i in rangePixBegs1[1:]:

         box = i - opts.border_pixels, tMargin, i, height

         im.paste(opts.border_color, box)

     for i in rangePixBegs2[1:]:

         box = lMargin, i - opts.border_pixels, width, i

         im.paste(opts.border_color, box)

     im.save(args[1])

 if __name__ == "__main__":

     usage = """%prog --help

    or: %prog [options] maf-or-tab-alignments dotplot.png

    or: %prog [options] maf-or-tab-alignments dotplot.gif

    or: ..."""

     description = "Draw a dotplot of pair-wise sequence alignments in MAF or tabular format."

     op = optparse.OptionParser(usage=usage, description=description)

     op.add_option("-v", "--verbose", action="count",

                   help="show progress messages & data about the plot")

     op.add_option("-1", "--seq1", metavar="PATTERN", action="append",

                   default=[],

                   help="which sequences to show from the 1st genome")

     op.add_option("-2", "--seq2", metavar="PATTERN", action="append",

                   default=[],

                   help="which sequences to show from the 2nd genome")

     # Replace "width" & "height" with a single "length" option?

     op.add_option("-x", "--width", type="int", default=1000,

                   help="maximum width in pixels (default: %default)")

     op.add_option("-y", "--height", type="int", default=1000,

                   help="maximum height in pixels (default: %default)")

     op.add_option("-c", "--forwardcolor", metavar="COLOR", default="red",

                   help="color for forward alignments (default: %default)")

     op.add_option("-r", "--reversecolor", metavar="COLOR", default="blue",

                   help="color for reverse alignments (default: %default)")

     op.add_option("--alignments", metavar="FILE", help="secondary alignments")

     op.add_option("--sort1", default="1", metavar="N",

                   help="genome1 sequence order: 0=input order, 1=name order, "

                   "2=length order, 3=alignment order (default=%default)")

     op.add_option("--sort2", default="1", metavar="N",

                   help="genome2 sequence order: 0=input order, 1=name order, "

                   "2=length order, 3=alignment order (default=%default)")

     op.add_option("--strands1", default="0", metavar="N", help=

                   "genome1 sequence orientation: 0=forward orientation, "

                   "1=alignment orientation (default=%default)")

     op.add_option("--strands2", default="0", metavar="N", help=

                   "genome2 sequence orientation: 0=forward orientation, "

                   "1=alignment orientation (default=%default)")

     op.add_option("--max-gap1", metavar="FRAC", default="0.5,2", help=

                   "maximum unaligned (end,mid) gap in genome1: "

                   "fraction of aligned length (default=%default)")

     op.add_option("--max-gap2", metavar="FRAC", default="0.5,2", help=

                   "maximum unaligned (end,mid) gap in genome2: "

                   "fraction of aligned length (default=%default)")

     op.add_option("--pad", metavar="FRAC", type="float", default=0.04, help=

                   "pad length when cutting unaligned gaps: "

                   "fraction of aligned length (default=%default)")

     op.add_option("--border-pixels", metavar="INT", type="int", default=1,

                   help="number of pixels between sequences (default=%default)")

     op.add_option("--border-color", metavar="COLOR", default="black",

                   help="color for pixels between sequences (default=%default)")

     # --break-color and/or --break-pixels for intra-sequence breaks?

     op.add_option("--margin-color", metavar="COLOR", default="#dcdcdc",

                   help="margin color")

     og = optparse.OptionGroup(op, "Text options")

     og.add_option("-f", "--fontfile", metavar="FILE",

                   help="TrueType or OpenType font file")

     og.add_option("-s", "--fontsize", metavar="SIZE", type="int", default=14,

                   help="TrueType or OpenType font size (default: %default)")

     og.add_option("--labels1", type="int", default=0, metavar="N", help=

                   "genome1 labels: 0=name, 1=name:length, "

                   "2=name:start:length, 3=name:start-end (default=%default)")

     og.add_option("--labels2", type="int", default=0, metavar="N", help=

                   "genome2 labels: 0=name, 1=name:length, "

                   "2=name:start:length, 3=name:start-end (default=%default)")

     og.add_option("--rot1", metavar="ROT", default="h",

                   help="text rotation for the 1st genome (default=%default)")

     og.add_option("--rot2", metavar="ROT", default="v",

                   help="text rotation for the 2nd genome (default=%default)")

     op.add_option_group(og)

     og = optparse.OptionGroup(op, "Annotation options")

     og.add_option("--bed1", metavar="FILE",

                   help="read genome1 annotations from BED file")

     og.add_option("--bed2", metavar="FILE",

                   help="read genome2 annotations from BED file")

     og.add_option("--rmsk1", metavar="FILE", help="read genome1 repeats from "

                   "RepeatMasker .out or rmsk.txt file")

     og.add_option("--rmsk2", metavar="FILE", help="read genome2 repeats from "

                   "RepeatMasker .out or rmsk.txt file")

     op.add_option_group(og)

     og = optparse.OptionGroup(op, "Gene options")

     og.add_option("--genePred1", metavar="FILE",

                   help="read genome1 genes from genePred file")

     og.add_option("--genePred2", metavar="FILE",

                   help="read genome2 genes from genePred file")

     og.add_option("--exon-color", metavar="COLOR", default="PaleGreen",

                   help="color for exons (default=%default)")

     og.add_option("--cds-color", metavar="COLOR", default="LimeGreen",

                   help="color for protein-coding regions (default=%default)")

     op.add_option_group(og)

     og = optparse.OptionGroup(op, "Unsequenced gap options")

     og.add_option("--gap1", metavar="FILE",

                   help="read genome1 unsequenced gaps from agp or gap file")

     og.add_option("--gap2", metavar="FILE",

                   help="read genome2 unsequenced gaps from agp or gap file")

     og.add_option("--bridged-color", metavar="COLOR", default="yellow",

                   help="color for bridged gaps (default: %default)")

     og.add_option("--unbridged-color", metavar="COLOR", default="orange",

                   help="color for unbridged gaps (default: %default)")

     op.add_option_group(og)

     (opts, args) = op.parse_args()

     if len(args) != 2: op.error("2 arguments needed")

     opts.background_color = "white"

     opts.label_space = 5     # minimum number of pixels between axis labels

     try: lastDotplot(opts, args)

     except KeyboardInterrupt: pass  # avoid silly error message

     except Exception as e:

         prog = os.path.basename(sys.argv[0])

         sys.exit(prog + ": error: " + str(e))