/usr/bin/last-dotplot is in last-align 921-1.
This file is owned by root:root, with mode 0o755.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 | #! /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
def myOpen(fileName): # faster than fileinput
if fileName is None:
return []
if fileName == "-":
return sys.stdin
if fileName.endswith(".gz"):
return gzip.open(fileName)
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 itertools.izip(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 = map(seqRequestFromText, opts.seq1)
seqRequests2 = map(seqRequestFromText, 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.iteritems():
v.sort()
yield k, list(mergedRanges(v))
def coveredLength(mergedCoverDict):
return sum(sum(e - b for b, e in v) for v in mergedCoverDict.itervalues())
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 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:
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, [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)
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))
|