%prog gaps fastafile

    Print out a list of gaps in BED format (.gaps.bed).
334800
def gaps(args):
    
    from jcvi.formats.sizes import agp
    from jcvi.formats.agp import mask, build

    p = OptionParser(gaps.__doc__)
    p.add_option("--split", default=False, action="store_true",
            help="Generate .split.fasta [default: %default]")
    p.set_mingap(default=100)
    p.set_cpus()
    opts, args = p.parse_args(args)

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

    inputfasta, = args
    mingap = opts.mingap
    split = opts.split
    prefix = inputfasta.rsplit(".", 1)[0]
    bedfile = prefix + ".gaps.bed"

    if need_update(inputfasta, bedfile):
        write_gaps_bed(inputfasta, prefix, mingap, opts.cpus)

    if split:
        splitfile = prefix + ".split.fasta"
        oagpfile = prefix + ".splitobject.agp"
        cagpfile = prefix + ".splitcomponent.agp"

        if need_update((inputfasta, bedfile), splitfile):

            sizesagpfile = agp([inputfasta])

            maskedagpfile = mask([sizesagpfile, bedfile, "--splitobject"])
            shutil.move(maskedagpfile, oagpfile)
            logging.debug("AGP file written to `{0}`.".format(oagpfile))

            maskedagpfile = mask([sizesagpfile, bedfile, "--splitcomponent"])
            shutil.move(maskedagpfile, cagpfile)
            logging.debug("AGP file written to `{0}`.".format(cagpfile))

            build([oagpfile, inputfasta, splitfile])
            os.remove(sizesagpfile)

        return splitfile, oagpfile, cagpfile
632566
def estimate(args):
    
    from collections import defaultdict
    from jcvi.formats.bed import intersectBed_wao
    from jcvi.formats.posmap import MatesFile

    p = OptionParser(estimate.__doc__)
    p.add_option("--minlinks", default=3, type="int",
                 help="Minimum number of links to place [default: %default]")
    opts, args = p.parse_args(args)

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

    gapsbed, spansbed, matesfile = args
    mf = MatesFile(matesfile)
    bed = Bed(gapsbed)
    order = bed.order

    gap2mate = defaultdict(set)
    mate2gap = defaultdict(set)

    for a, b in intersectBed_wao(gapsbed, spansbed):
        gapsize = a.span
        if gapsize != 100:
            continue

        gapname = a.accn

        if b is None:
            gap2mate[gapname] = set()
            continue

        matename = b.accn
        gap2mate[gapname].add(matename)
        mate2gap[matename].add(gapname)

    omgapsbed = "gaps.linkage.bed"
    fw = open(omgapsbed, "w")
    for gapname, mates in sorted(gap2mate.items()):
        i, b = order[gapname]
        nmates = len(mates)
        if nmates < opts.minlinks:
            print("{0}\t{1}".format(b, nmates), file=fw)
            continue

        print(gapname, mates)

    fw.close()
337296
def print_menu(place, static=False):
    
    day = get_day()
    if static:
        plan = get(FILES[1])
        for meal in plan["weeks"][0]["days"][day][place]["meals"]:
            if place == "Diner":
                print(meal["category"] + " " + meal["meal"])
            else:
                print(meal["category"] + ": " + meal["meal"])
    else:
        plan = get(FILES[0])
        for meal in plan["weeks"][0]["days"][day][place]["meals"]:
            print(meal["category"] + ": " + meal["meal"])
1070367
def printout(*args, **kwargs):
    

    # TODO(Lukas): conflicts with function names
    color = kwargs.pop('color', {})
    style = kwargs.pop('style', {})
    prefx = kwargs.pop('prefix', '')
    suffx = kwargs.pop('suffix', '')
    ind = kwargs.pop('indent', 0)

    print_args = []
    for arg in args:
        arg = str(arg)
        arg = colorize(arg, **color)
        arg = stylize(arg, **style)
        arg = prefix(arg, prefx)
        arg = indent(arg, ind)
        arg += str(suffx)
        print_args.append(arg)

    print(*print_args, **kwargs)
416497
def tidy(args):
    
    p = OptionParser(tidy.__doc__)
    p.add_option("--gapsize", dest="gapsize", default=0, type="int",
            help="Set all gaps to the same size [default: %default]")
    p.add_option("--minlen", dest="minlen", default=100, type="int",
            help="Minimum component size [default: %default]")

    opts, args = p.parse_args(args)

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

    fastafile, = args
    gapsize = opts.gapsize
    minlen = opts.minlen

    tidyfastafile = fastafile.rsplit(".", 1)[0] + ".tidy.fasta"
    fw = must_open(tidyfastafile, "w")

    removed = normalized = 0
    fasta = Fasta(fastafile, lazy=True)
    for name, rec in fasta.iteritems_ordered():
        rec.seq = rec.seq.upper()
        if minlen:
            removed += remove_small_components(rec, minlen)
        trim_terminal_Ns(rec)
        if gapsize:
            normalized += normalize_gaps(rec, gapsize)

        if len(rec) == 0:
            logging.debug("Drop seq {0}".format(rec.id))
            continue
        SeqIO.write([rec], fw, "fasta")

    # Print statistics
    if removed:
        logging.debug("Total discarded bases: {0}".format(removed))
    if normalized:
        logging.debug("Gaps normalized: {0}".format(normalized))

    logging.debug("Tidy FASTA written to `{0}`.".format(tidyfastafile))
    fw.close()

    return tidyfastafile
579905
def filter(args):
    
    p = OptionParser(filter.__doc__)
    p.add_option("--less", default=False, action="store_true",
                 help="filter the sizes < certain cutoff [default: >=]")
    p.set_outfile()

    opts, args = p.parse_args(args)

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

    fastafile, cutoff = args
    try:
        cutoff = int(cutoff)
    except ValueError:
        sys.exit(not p.print_help())

    f = Fasta(fastafile, lazy=True)

    fw = must_open(opts.outfile, "w")
    for name, rec in f.iteritems_ordered():

        if opts.less and len(rec) >= cutoff:
            continue

        if (not opts.less) and len(rec) < cutoff:
            continue

        SeqIO.write([rec], fw, "fasta")
        fw.flush()

    return fw.name
1290294
def cut(args):
    
    p = OptionParser(cut.__doc__)
    p.add_option("--sep", default=".", help="Separator for splits")
    opts, args = p.parse_args(args)

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

    agpfile, bedfile = args
    sep = opts.sep

    agp = AGP(agpfile)
    bed = Bed(bedfile)
    simple_agp = agp.order
    newagpfile = agpfile.replace(".agp", ".cut.agp")
    fw = open(newagpfile, "w")

    agp_fixes = defaultdict(list)
    for component, intervals in bed.sub_beds():
        i, a = simple_agp[component]
        object = a.object
        component_span = a.component_span
        orientation = a.orientation

        assert a.component_beg, a.component_end
        cuts = set()
        for i in intervals:
            start, end = i.start, i.end
            end -= 1

            assert start <= end
            cuts.add(start)
            cuts.add(end)

        cuts.add(0)
        cuts.add(component_span)
        cuts = list(sorted(cuts))

        sum_of_spans = 0
        for i, (a, b) in enumerate(pairwise(cuts)):
            oid = object + "{0}{1}".format(sep, i + 1)
            aline = [oid, 0, 0, 0]
            cspan = b - a
            aline += ['D', component, a + 1, b, orientation]
            sum_of_spans += cspan

            aline = "\t".join(str(x) for x in aline)
            agp_fixes[component].append(aline)

        assert component_span == sum_of_spans

    # Finally write the masked agp
    for a in agp:
        if not a.is_gap and a.component_id in agp_fixes:
            print("\n".join(agp_fixes[a.component_id]), file=fw)
        else:
            print(a, file=fw)

    fw.close()
    # Reindex
    reindex([newagpfile, "--inplace"])

    return newagpfile
588384
def fillstats(args):
    
    from jcvi.utils.cbook import SummaryStats, percentage, thousands

    p = OptionParser(fillstats.__doc__)
    opts, args = p.parse_args(args)

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

    fillfile, = args
    fp = open(fillfile)
    scaffolds = 0
    gaps = []
    for row in fp:
        if row[0] == ">":
            scaffolds += 1
            continue
        fl = FillLine(row)
        gaps.append(fl)

    print("{0} scaffolds in total".format(scaffolds), file=sys.stderr)

    closed = [x for x in gaps if x.closed]
    closedbp = sum(x.before for x in closed)
    notClosed = [x for x in gaps if not x.closed]
    notClosedbp = sum(x.before for x in notClosed)

    totalgaps = len(closed) + len(notClosed)

    print("Closed gaps: {0} size: {1} bp".\
                        format(percentage(len(closed), totalgaps), thousands(closedbp)), file=sys.stderr)
    ss = SummaryStats([x.after for x in closed])
    print(ss, file=sys.stderr)

    ss = SummaryStats([x.delta for x in closed])
    print("Delta:", ss, file=sys.stderr)

    print("Remaining gaps: {0} size: {1} bp".\
                        format(percentage(len(notClosed), totalgaps), thousands(notClosedbp)), file=sys.stderr)
    ss = SummaryStats([x.after for x in notClosed])
    print(ss, file=sys.stderr)
881969
def overlapbatch(args):
    
    p = OptionParser(overlap.__doc__)
    opts, args = p.parse_args(args)
    if len(args) != 2:
        sys.exit(not p.print_help())

    ctgfasta, poolfasta = args
    f = Fasta(ctgfasta)
    for k, rec in f.iteritems_ordered():
        fastafile = k + ".fasta"
        fw = open(fastafile, "w")
        SeqIO.write([rec], fw, "fasta")
        fw.close()

        overlap([fastafile, poolfasta])
886330
def fasta(args):
    
    from jcvi.formats.sizes import Sizes

    p = OptionParser(fasta.__doc__)
    p.add_option("--extend", default=1000, type="int",
                 help="Extend seq flanking the gaps [default: %default]")
    opts, args = p.parse_args(args)

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

    mapout, sfasta = args
    Flank = opts.extend
    pf = mapout.split(".")[0]
    mapbed = pf + ".bed"
    bm = BinMap(mapout)
    bm.print_to_bed(mapbed)

    bed = Bed(mapbed, sorted=False)
    markersbed = pf + ".markers.bed"
    fw = open(markersbed, "w")
    sizes = Sizes(sfasta).mapping
    for b in bed:
        accn = b.accn
        scf, pos = accn.split(".")
        pos = int(pos)
        start = max(0, pos - Flank)
        end = min(pos + Flank, sizes[scf])
        print("\t".join(str(x) for x in \
                    (scf, start, end, accn)), file=fw)

    fw.close()

    fastaFromBed(markersbed, sfasta, name=True)
212743
def print_intervals(intervals):
        
        res = []
        for i in intervals:
            res.append(repr(i))
        print("".join(res))
908931
def gapsplit(args):
    
    p = OptionParser(gapsplit.__doc__)

    opts, args = p.parse_args(args)

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

    gffile, = args

    gff = Gff(gffile)
    for g in gff:
        if re.match("EST_match", g.type):
            """
            hacky implementation:
            since the standard urlparse.parse_qsl() replaces all "+" symbols with spaces
            we will write a regex to check either for a "-" or a " " (space)
            """
            match = re.search(r'\S+ (\d+) \d+ ([\s{1}\-])', g.attributes["Target"][0])
            if match.group(2) == "-":
                strand = match.group(2)
            else:
                strand = "+"
                g.attributes["Target"][0] = " ".join(str(x) \
                        for x in [g.attributes["Target"][0].rstrip(), strand])

            if g.strand == "?":
                g.strand = strand
        else:
            match = re.match(r'\S+ (\d+) \d+', g.attributes["Target"][0])
        target_start = int(match.group(1))

        re_cigar = re.compile(r'(\D+)(\d+)');
        cigar = g.attributes["Gap"][0].split(" ")
        g.attributes["Gap"] = None

        parts = []
        if g.strand == "+":
            for event in cigar:
                match = re_cigar.match(event)
                op, count = match.group(1), int(match.group(2))
                if op in "IHS":
                    target_start += count
                elif op in "DN":
                    g.start += count
                elif op == "P":
                    continue
                else:
                    parts.append([g.start, g.start + count - 1, \
                            target_start, target_start + count - 1])
                    g.start += count
                    target_start += count
        else:
            for event in cigar:
                match = re_cigar.match(event)
                op, count = match.group(1), int(match.group(2))
                if op in "IHS":
                    target_start += count
                elif op in "DN":
                    g.end -= count
                elif op == "P":
                    continue
                else:
                    parts.append([g.end - count + 1, g.end, \
                            target_start, target_start + count - 1])
                    g.end -= count
                    target_start += count

        g.update_attributes()
        print(g)

        parent = g.attributes["Name"][0]
        g.type = "match_part"
        g.attributes.clear()

        for part in parts:
            g.start, g.end = part[0], part[1]
            g.score, g.strand, g.phase = ".", g.strand, "."

            if re.match("EST", g.type):
                target_list = [parent, part[2], part[3], g.strand]
            else:
                target_list = [parent, part[2], part[3]]
            target = " ".join(str(x) for x in target_list)

            g.attributes["Parent"] = [parent]
            g.attributes["Target"] = [target]

            g.update_attributes()
            print(g)
1246954
def estimategaps(args):
    
    p = OptionParser(estimategaps.__doc__)
    p.add_option("--minsize", default=100, type="int",
                 help="Minimum gap size")
    p.add_option("--maxsize", default=500000, type="int",
                 help="Maximum gap size")
    p.add_option("--links", default=10, type="int",
                 help="Only use linkage grounds with matchings more than")
    p.set_verbose(help="Print details for each gap calculation")
    opts, args = p.parse_args(args)

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

    inputbed, = args
    pf = inputbed.rsplit(".", 1)[0]
    agpfile = pf + ".chr.agp"
    bedfile = pf + ".lifted.bed"

    cc = Map(bedfile, scaffold_info=True)
    agp = AGP(agpfile)
    minsize, maxsize = opts.minsize, opts.maxsize
    links = opts.links
    verbose = opts.verbose

    outagpfile = pf + ".estimategaps.agp"
    fw = must_open(outagpfile, "w")

    for ob, components in agp.iter_object():
        components = list(components)
        s = Scaffold(ob, cc)
        mlg_counts = s.mlg_counts
        gaps = [x for x in components if x.is_gap]
        gapsizes = [None] * len(gaps)   # master
        for mlg, count in mlg_counts.items():
            if count < links:
                continue
            g = GapEstimator(cc, agp, ob, mlg)
            g.compute_all_gaps(minsize=minsize, maxsize=maxsize, \
                               verbose=verbose)
            # Merge evidence from this mlg into master
            assert len(g.gapsizes) == len(gaps)
            for i, gs in enumerate(gapsizes):
                gg = g.gapsizes[i]
                if gs is None:
                    gapsizes[i] = gg
                elif gg:
                    gapsizes[i] = min(gs, gg)

        print(gapsizes)
        # Modify AGP
        i = 0
        for x in components:
            if x.is_gap:
                x.gap_length = gapsizes[i] or minsize
                x.component_type = 'U' if x.gap_length == 100 else 'N'
                i += 1
            print(x, file=fw)

    fw.close()
    reindex([outagpfile, "--inplace"])
1024172
def gaps(self) -> "IntervalList":
        
        if len(self.intervals) < 2:
            return IntervalList(None)
        gaps = []
        for i in range(len(self.intervals) - 1):
            gap = Interval(
                self.intervals[i].end,
                self.intervals[i + 1].start
            )
            gaps.append(gap)
        return IntervalList(gaps)
319485
def tofasta(args):
    
    p = OptionParser(tofasta.__doc__)
    p.add_option("--prefix", default="gbfasta",
            help="prefix of output files [default: %default]")
    filenames, accessions, idfile, opts, args = preparegb(p, args)
    prefix = opts.prefix

    GenBank(filenames=filenames, accessions=accessions, idfile=idfile).\
        write_fasta(output=prefix, individual=opts.individual)

    if opts.individual:
        logging.debug("Output written dir {0}".format(prefix))
    else:
        logging.debug("Output written to {0}.fasta".format(prefix))
426070
def blat(args):
    
    from jcvi.formats.base import is_number
    from jcvi.formats.blast import best as blast_best, bed as blast_bed
    from jcvi.apps.align import blat as blat_align

    p = OptionParser(blat.__doc__)
    opts, args = p.parse_args(args)

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

    maptxt, ref = args
    pf = maptxt.rsplit(".", 1)[0]
    register = {}
    fastafile = pf + ".fasta"
    fp = open(maptxt)
    fw = open(fastafile, "w")
    for row in fp:
        name, lg, pos, seq = row.split()
        if not is_number(pos):
            continue
        register[name] = (pf + '-' + lg, pos)
        print(">{0}\n{1}\n".format(name, seq), file=fw)
    fw.close()

    blatfile = blat_align([ref, fastafile])
    bestfile = blast_best([blatfile])
    bedfile = blast_bed([bestfile])
    b = Bed(bedfile).order

    pf = ".".join((op.basename(maptxt).split(".")[0],
                   op.basename(ref).split(".")[0]))
    csvfile = pf + ".csv"
    fp = open(maptxt)
    fw = open(csvfile, "w")
    for row in fp:
        name, lg, pos, seq = row.split()
        if name not in b:
            continue
        bbi, bb = b[name]
        scaffold, scaffold_pos = bb.seqid, bb.start
        print(",".join(str(x) for x in \
                    (scaffold, scaffold_pos, lg, pos)), file=fw)
    fw.close()
780379
def main(show_details:['-l']=False, cols:['-w', '--width']='', *files):
    '''
    
    '''
    print(files)
    print(show_details)
    print(cols)