Made last-dotplot faster, sometimes.
authorMartin C. Frith
Tue Apr 04 10:50:57 2017 +0900 (2017-04-04)
changeset 84423de4eb3be1d
parent 843 7872fdade23a
child 845 16060c00b129
Made last-dotplot faster, sometimes.
scripts/last-dotplot
     1.1 --- a/scripts/last-dotplot	Thu Mar 09 13:40:51 2017 +0900
     1.2 +++ b/scripts/last-dotplot	Tue Apr 04 10:50:57 2017 +0900
     1.3 @@ -9,12 +9,17 @@
     1.4  # according to the number of aligned nt-pairs within it, but the
     1.5  # result is too faint.  How can this be done better?
     1.6  
     1.7 -import fileinput, fnmatch, itertools, optparse, os, re, sys
     1.8 +import fnmatch, itertools, optparse, os, re, sys
     1.9  
    1.10  # Try to make PIL/PILLOW work:
    1.11  try: from PIL import Image, ImageDraw, ImageFont, ImageColor
    1.12  except ImportError: import Image, ImageDraw, ImageFont, ImageColor
    1.13  
    1.14 +def myOpen(fileName):  # faster than fileinput
    1.15 +    if fileName == "-":
    1.16 +        return sys.stdin
    1.17 +    return open(fileName)
    1.18 +
    1.19  def warn(message):
    1.20      prog = os.path.basename(sys.argv[0])
    1.21      sys.stderr.write(prog + ": " + message + "\n")
    1.22 @@ -107,7 +112,7 @@
    1.23      for pat, beg, end in seqRanges:
    1.24          if fnmatch.fnmatchcase(name, pat) or fnmatch.fnmatchcase(base, pat):
    1.25              return max(beg, 0), min(end, seqLen)
    1.26 -    return 0, 0
    1.27 +    return None
    1.28  
    1.29  def updateSeqLimits(isTrim, seqLimits, seqName, seqRange, blocks, index):
    1.30      if isTrim:
    1.31 @@ -130,10 +135,12 @@
    1.32      alignments = []
    1.33      seqLimits1 = {}
    1.34      seqLimits2 = {}
    1.35 -    lines = fileinput.input(fileName)
    1.36 +    lines = myOpen(fileName)
    1.37      for seqName1, seqLen1, seqName2, seqLen2, blocks in alignmentInput(lines):
    1.38          range1 = rangeFromSeqName(seqRanges1, seqName1, seqLen1)
    1.39 +        if not range1: continue
    1.40          range2 = rangeFromSeqName(seqRanges2, seqName2, seqLen2)
    1.41 +        if not range2: continue
    1.42          b = list(croppedBlocks(blocks, range1, range2))
    1.43          if not b: continue
    1.44          aln = seqName1, seqName2, b
    1.45 @@ -258,7 +265,7 @@
    1.46      '''Read locations of unsequenced gaps, from an agp or gap file.'''
    1.47      if not fileName: return
    1.48      seqLimits = expandedSeqDict(seqLimits)
    1.49 -    for line in fileinput.input(fileName):
    1.50 +    for line in myOpen(fileName):
    1.51          w = line.split()
    1.52          if not w or w[0][0] == "#": continue
    1.53          if isExtraFirstGapField(w): w = w[1:]