def overlap(args):
    from jcvi.apps.console import green

    p = OptionParser(overlap.__doc__)
    p.add_option("--maxerr", default=2, type="int", help="Maximum error rate")
    p.add_option("--canvas", default=100, type="int", help="Canvas size")
    opts, args = p.parse_args(args)

    if len(args) != 2:
        sys.exit(not p.print_help())

    bestcontains, iid = args
    canvas = opts.canvas

    bestcontainscache = bestcontains + ".cache"
    if need_update(bestcontains, bestcontainscache):
        fp = open(bestcontains)
        fw = open(bestcontainscache, "w")
        exclude = set()
        for row in fp:
            if row[0] == '#':
                continue
            j = int(row.split()[0])
            exclude.add(j)
        cPickle.dump(exclude, fw)
        fw.close()

    exclude = cPickle.load(open(bestcontainscache))
    logging.debug("A total of {0} reads to exclude".format(len(exclude)))

    cmd = "overlapStore -d ../asm.ovlStore -b {0} -e {0}".format(iid)
    cmd += " -E {0}".format(opts.maxerr)
    frags = []
    for row in popen(cmd):
        r = OverlapLine(row)
        if r.bid in exclude:
            continue
        frags.append(r)

    # Also include to query fragment
    frags.append(OverlapLine("{0} {0} N 0 0 0 0".format(iid)))
    frags.sort(key=lambda x: x.ahang)

    # Determine size of the query fragment
    cmd = "gatekeeper -b {0} -e {0}".format(iid)
    cmd += " -tabular -dumpfragments ../asm.gkpStore"
    fp = popen(cmd)
    row = fp.next()
    size = int(fp.next().split()[-1])

    # Determine size of canvas
    xmin = min(x.ahang for x in frags)
    xmax = max(x.bhang for x in frags)
    xsize = -xmin + size + xmax
    ratio = xsize / canvas

    fw = sys.stdout
    for f in frags:
        fsize = -f.ahang + size + f.bhang
        a = (f.ahang - xmin) / ratio
        b = fsize / ratio
        t = '-' * b
        if f.orientation == 'N':
            t = t[:-1] + '>'
        else:
            t = '<' + t[1:]
        if f.ahang == 0 and f.bhang == 0:
            t = green(t)
        c = canvas - a - b
        fw.write(' ' * a)
        fw.write(t)
        fw.write(' ' * c)
        print >> fw, "{0} ({1})".format(str(f.bid).rjust(10), f.erate_adj)
%prog overlap best.contains iid
dataset/ETHPy150Open tanghaibao/jcvi/assembly/ca.py/overlap