%prog bed btabfile

    Convert btab to bed format.
15114
def bed(args):
    
    from jcvi.formats.blast import BlastLine
    p = OptionParser(bed.__doc__)

    opts, args = p.parse_args(args)

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

    btabfile, = args
    btab = Btab(btabfile)
    for b in btab:
        Bline = BlastLine(b.blastline)
        print(Bline.bedline)
878989
def bins(args):
    
    from jcvi.formats.sizes import Sizes

    p = OptionParser(bins.__doc__)
    p.add_option("--binsize", default=100000, type="int",
                 help="Size of the bins [default: %default]")
    p.add_option("--subtract",
                 help="Subtract bases from window [default: %default]")
    p.add_option("--mode", default="span", choices=("span", "count", "score"),
                 help="Accumulate feature based on [default: %default]")
    p.add_option("--nomerge", default=False, action="store_true",
                 help="Do not merge features")
    opts, args = p.parse_args(args)

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

    bedfile, fastafile = args
    subtract = opts.subtract
    mode = opts.mode
    assert op.exists(bedfile), "File `{0}` not found".format(bedfile)

    binsize = opts.binsize
    binfile = bedfile + ".{0}".format(binsize)
    binfile += ".{0}.bins".format(mode)

    if not need_update(bedfile, binfile):
        return binfile

    sz = Sizes(fastafile)
    sizesfile = sz.filename
    sizes = sz.mapping
    fw = open(binfile, "w")
    scores = "median" if mode == "score" else None
    if not opts.nomerge:
        bedfile = mergeBed(bedfile, nms=True, scores=scores)
    if subtract:
        subtractmerge = mergeBed(subtract)
        subtract_complement = complementBed(subtractmerge, sizesfile)
        bedfile = intersectBed(bedfile, subtract_complement)

    bedfile = sort([bedfile, "-i"])

    bed = Bed(bedfile)
    sbdict = dict(bed.sub_beds())
    for chr, chr_len in sorted(sizes.items()):
        chr_len = sizes[chr]
        subbeds = sbdict.get(chr, [])
        nbins = chr_len / binsize
        last_bin = chr_len % binsize
        if last_bin:
            nbins += 1

        a = np.zeros(nbins)  # values
        b = np.zeros(nbins, dtype="int")  # bases
        c = np.zeros(nbins, dtype="int")  # count
        b[:-1] = binsize
        b[-1] = last_bin

        for bb in subbeds:

            start, end = bb.start, bb.end
            startbin = start / binsize
            endbin = end / binsize

            assert startbin <= endbin
            c[startbin:endbin + 1] += 1

            if mode == "score":
                a[startbin:endbin + 1] += float(bb.score)

            elif mode == "span":
                if startbin == endbin:
                    a[startbin] += end - start + 1

                if startbin < endbin:
                    firstsize = (startbin + 1) * binsize - start + 1
                    lastsize = end - endbin * binsize
                    a[startbin] += firstsize
                    if startbin + 1 < endbin:
                        a[startbin + 1:endbin] += binsize
                    a[endbin] += lastsize

        if mode == "count":
            a = c

        for xa, xb in zip(a, b):
            print("\t".join(str(x) for x in (chr, xa, xb)), file=fw)

    fw.close()

    if subtract:
        subtractbinfile = bins([subtract, fastafile, "--binsize={0}".format(binsize)])
        binfile = subtractbins(binfile, subtractbinfile)

    return binfile
48798
def deletion(args):
    
    p = OptionParser(deletion.__doc__)
    p.add_option("--mindepth", default=3, type="int",
                 help="Minimum depth to call a deletion")
    p.add_option("--minspan", default=30, type="int",
                 help="Minimum span to call a deletion")
    p.add_option("--split", default=False, action="store_true",
                 help="Break at cigar N into separate parts")
    p.set_tmpdir()
    opts, args = p.parse_args(args)

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

    bedfile, gapsbedfile = args
    if bedfile.endswith(".bam"):
        bamfile = bedfile
        bedfile = bamfile.replace(".sorted.", ".").replace(".bam", ".bed")
        if need_update(bamfile, bedfile):
            cmd = "bamToBed -i {0}".format(bamfile)
            if opts.split:
                cmd += " -split"
            cmd += " | cut -f1-4"
            sh(cmd, outfile=bedfile)

    sort_tmpdir = "--tmpdir={0}".format(opts.tmpdir)
    if bedfile.endswith(".sorted.bed"):
        pf = bedfile.rsplit(".", 2)[0]
        sortedbedfile = bedfile
    else:
        pf = bedfile.rsplit(".", 1)[0]
        sortedbedfile = pf + ".sorted.bed"
        if need_update(bedfile, sortedbedfile):
            sort([bedfile, "-u", "--accn", sort_tmpdir])

    # Find reads that contain multiple matches
    ibedfile = pf + ".d.bed"
    if need_update(sortedbedfile, ibedfile):
        bed = Bed(sortedbedfile, sorted=False)
        fw = open(ibedfile, "w")
        logging.debug("Write deletions to `{0}`.".format(ibedfile))
        for accn, bb in groupby(bed, key=lambda x: x.accn):
            bb = list(bb)
            branges = [(x.seqid, x.start, x.end) for x in bb]
            iranges = range_interleave(branges)
            for seqid, start, end in iranges:
                if end - start + 1 < opts.minspan:
                    continue
                print("\t".join(str(x) for x in \
                            (seqid, start - 1, end, accn + '-d')), file=fw)
        fw.close()

    # Uniqify the insertions and count occurrences
    countbedfile = pf + ".uniq.bed"
    if need_update(ibedfile, countbedfile):
        bed = Bed(ibedfile)
        fw = open(countbedfile, "w")
        logging.debug("Write counts to `{0}`.".format(countbedfile))
        registry = Counter((x.seqid, x.start, x.end) for x in bed)
        ies_id = 1
        for (seqid, start, end), count in registry.items():
            ies_name = "{0:05d}-r{1}".format(ies_id, count)
            if count < opts.mindepth:
                continue
            print("\t".join(str(x) for x in \
                            (seqid, start - 1, end, ies_name)), file=fw)
            ies_id += 1
        fw.close()
        sort([countbedfile, "-i", sort_tmpdir])

    # Remove deletions that contain some read depth
    depthbedfile = pf + ".depth.bed"
    if need_update((sortedbedfile, countbedfile), depthbedfile):
        depth([sortedbedfile, countbedfile, "--outfile={0}".format(depthbedfile)])

    validbedfile = pf + ".valid.bed"
    if need_update(depthbedfile, validbedfile):
        fw = open(validbedfile, "w")
        logging.debug("Filter valid deletions to `{0}`.".format(validbedfile))
        bed = Bed(depthbedfile)
        all_scores = [float(b.score) for b in bed]
        lb, ub = outlier_cutoff(all_scores)
        logging.debug("Bounds for depths: LB={0:.2f} (ignored)  UB={1:.2f}".format(lb, ub))
        for b in bed:
            if float(b.score) > ub:
                continue
            print(b, file=fw)
        fw.close()

    # Remove deletions that contain sequencing gaps on its flanks
    selectedbedfile = pf + ".selected.bed"
    if need_update(validbedfile, selectedbedfile):
        flanksbedfile = pf + ".flanks.bed"
        fw = open(flanksbedfile, "w")
        bed = Bed(validbedfile)
        flank = 100
        logging.debug("Write deletion flanks to `{0}`.".format(flanksbedfile))
        for b in bed:
            start, end = b.start, b.end
            b.start, b.end = start, min(start + flank - 1, end)
            print(b, file=fw)
            b.start, b.end = max(start, end - flank + 1), end
            print(b, file=fw)
        fw.close()

        intersectidsfile = pf + ".intersect.ids"
        cmd = "intersectBed -a {0} -b {1}".format(flanksbedfile, gapsbedfile)
        cmd += " | cut -f4 | sort -u"
        sh(cmd, outfile=intersectidsfile)
        some([validbedfile, intersectidsfile, "-v",
                "--outfile={0}".format(selectedbedfile)])

    # Find best-scoring non-overlapping set
    iesbedfile = pf + ".ies.bed"
    if need_update(selectedbedfile, iesbedfile):
        bed = Bed(selectedbedfile)
        fw = open(iesbedfile, "w")
        logging.debug("Write IES to `{0}`.".format(iesbedfile))
        branges = [Range(x.seqid, x.start, x.end, int(x.accn.rsplit("r")[-1]), i) \
                        for i, x in enumerate(bed)]
        iranges, iscore = range_chain(branges)
        logging.debug("Best chain score: {0} ({1} IES)".\
                        format(iscore, len(iranges)))
        ies_id = 1
        for seqid, start, end, score, id in iranges:
            ies_name = "IES-{0:05d}-r{1}".format(ies_id, score)
            span = end - start + 1
            print("\t".join(str(x) for x in \
                            (seqid, start - 1, end, ies_name, span)), file=fw)
            ies_id += 1
        fw.close()
677722
def gb(args):
    
    from Bio.Alphabet import generic_dna
    try:
        from BCBio import GFF
    except ImportError:
        print("You need to install dep first: $ easy_install bcbio-gff", file=sys.stderr)

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

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

    gff_file, fasta_file = args
    pf = op.splitext(gff_file)[0]
    out_file = pf + ".gb"
    fasta_input = SeqIO.to_dict(SeqIO.parse(fasta_file, "fasta", generic_dna))
    gff_iter = GFF.parse(gff_file, fasta_input)
    SeqIO.write(gff_iter, out_file, "genbank")
1169809
def mutect_to_bed(df):
    
    s = (df.contig.astype(str) + '\t' + 
         (df.position - 1).astype(int).astype(str) + '\t' + 
         df.position.astype(int).astype(str) + '\t' + 
         df.tumor_name)
    s = '\n'.join(s.values) + '\n'
    bt = pbt.BedTool(s, from_string=True)
    return bt
571329
def pairs(args):
    
    import jcvi.formats.bed

    p = OptionParser(pairs.__doc__)
    p.set_pairs()
    opts, targs = p.parse_args(args)

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

    samfile, = targs
    bedfile = samfile.rsplit(".", 1)[0] + ".bed"
    if need_update(samfile, bedfile):
        cmd = "bamToBed -i {0}".format(samfile)
        sh(cmd, outfile=bedfile)

    args[args.index(samfile)] = bedfile

    return jcvi.formats.bed.pairs(args)
572038
def mcscanx(args):
    
    p = OptionParser(mcscanx.__doc__)
    opts, args = p.parse_args(args)

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

    blastfile = args[0]
    bedfiles = args[1:]
    prefix = "_".join(op.basename(x)[:2] for x in bedfiles)
    symlink(blastfile, prefix + ".blast")
    allbedfile = prefix + ".gff"
    fw = open(allbedfile, "w")
    for i, bedfile in enumerate(bedfiles):
        prefix = chr(ord('A') + i)
        make_gff(bedfile, prefix, fw)
    fw.close()
245273
def binary_dumps(obj, alt_format=False):
    
    return b''.join(_binary_dump_gen(obj, alt_format=alt_format))
324916
def bedpe(args):
    
    from jcvi.assembly.coverage import bed_to_bedpe

    p = OptionParser(bedpe.__doc__)
    p.add_option("--span", default=False, action="store_true",
                 help="Write span bed file [default: %default]")
    p.add_option("--strand", default=False, action="store_true",
            help="Write the strand columns [default: %default]")
    p.add_option("--mates", help="Check the library stats from .mates file")
    opts, args = p.parse_args(args)

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

    bedfile, = args
    pf = bedfile.rsplit(".", 1)[0]
    bedpefile = pf + ".bedpe"
    bedspanfile = pf + ".spans.bed" if opts.span else None
    bed_to_bedpe(bedfile, bedpefile, \
                 pairsbedfile=bedspanfile, matesfile=opts.mates, \
                 strand=opts.strand)
    return bedpefile, bedspanfile
978779
def dbfpack(fn, usecols=None):
    """"""

    loadfunc = pulldbf.dbf_asdict
    cp = ChannelPack(loadfunc)
    cp.load(fn, usecols)
    names = pulldbf.channel_names(fn, usecols)
    cp.set_channel_names(names)
    # cp.set_basefilemtime()
    return cp
35955
def _bam_to_outbase(bam_file, work_dir, data):
    
    batch = dd.get_batch(data) or dd.get_sample_name(data)
    out_base = os.path.splitext(os.path.basename(bam_file))[0].split(".")[0]
    base = os.path.join(work_dir, out_base)
    return "%s-%s" % (base, batch), base
1164013
def _bed_to_bed6(orig_file, out_dir):
    
    bed6_file = os.path.join(out_dir, "%s-bed6%s" % os.path.splitext(os.path.basename(orig_file)))
    if not utils.file_exists(bed6_file):
        with open(bed6_file, "w") as out_handle:
            for i, region in enumerate(list(x) for x in pybedtools.BedTool(orig_file)):
                region = [x for x in list(region) if x]
                fillers = [str(i), "1.0", "+"]
                full = region + fillers[:6 - len(region)]
                out_handle.write("\t".join(full) + "\n")
    return bed6_file
1369968
def clean_file(in_file, data, prefix="", bedprep_dir=None, simple=None):
    
    # Remove non-ascii characters. Used in coverage analysis, to support JSON code in one column
    #   and be happy with sambamba:
    simple = "iconv -c -f utf-8 -t ascii | sed 's/ //g' |" if simple else ""
    if in_file:
        if not bedprep_dir:
            bedprep_dir = utils.safe_makedir(os.path.join(data["dirs"]["work"], "bedprep"))
        # Avoid running multiple times with same prefix
        if prefix and os.path.basename(in_file).startswith(prefix):
            return in_file
        out_file = os.path.join(bedprep_dir, "%s%s" % (prefix, os.path.basename(in_file)))
        out_file = out_file.replace(".interval_list", ".bed")
        if out_file.endswith(".gz"):
            out_file = out_file[:-3]
        if not utils.file_uptodate(out_file, in_file):
            check_bed_contigs(in_file, data)
            check_bed_coords(in_file, data)
            with file_transaction(data, out_file) as tx_out_file:
                bcbio_py = sys.executable
                cat_cmd = "zcat" if in_file.endswith(".gz") else "cat"
                sort_cmd = get_sort_cmd(os.path.dirname(tx_out_file))
                cmd = ("{cat_cmd} {in_file} | grep -v ^track | grep -v ^browser | grep -v ^@ | "
                       "grep -v ^# | {simple} "
                       "{bcbio_py} -c 'from bcbio.variation import bedutils; bedutils.remove_bad()' | "
                       "{sort_cmd} -k1,1 -k2,2n > {tx_out_file}")
                do.run(cmd.format(**locals()), "Prepare cleaned BED file", data)
        vcfutils.bgzip_and_index(out_file, data.get("config", {}), remove_orig=False)
        return out_file
511061
def fasta(args):
    
    from jcvi.formats.sizes import Sizes
    from jcvi.formats.agp import OO, build

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

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

    bedfile, scffasta, pmolfasta = args
    pf = bedfile.rsplit(".", 1)[0]
    bed = Bed(bedfile)
    selected = select_bed(bed)
    oo = OO()
    seen = set()
    sizes = Sizes(scffasta).mapping
    agpfile = pf + ".agp"
    agp = open(agpfile, "w")
    for b in selected:
        scf = range_parse(b.accn).seqid
        chr = b.seqid
        cs = (chr, scf)
        if cs not in seen:
            oo.add(chr, scf, sizes[scf], b.strand)
            seen.add(cs)
        else:
            logging.debug("Seen {0}, ignored.".format(cs))

    oo.write_AGP(agp, gaptype="contig")
    agp.close()
    build([agpfile, scffasta, pmolfasta])
1332559
def mergebydepth(args):
    
    p = OptionParser(mergebydepth.__doc__)
    p.add_option("--mindepth", default=3, type="int",
                 help="Minimum depth required")
    opts, args = p.parse_args(args)

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

    bedfile, fastafile = args
    mindepth = opts.mindepth
    bedgraph = make_bedgraph(bedfile)

    bedgraphfiltered = bedgraph + ".d{0}".format(mindepth)
    if need_update(bedgraph, bedgraphfiltered):
        filter([bedgraph, "--minaccn={0}".format(mindepth),
                "--outfile={0}".format(bedgraphfiltered)])

    merged = bedgraphfiltered + ".merge.fasta"
    if need_update(bedgraphfiltered, merged):
        mergeBed(bedgraphfiltered, sorted=True)
723000
def bed(args):
    
    p = OptionParser(bed.__doc__)
    opts, args = p.parse_args(args)

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

    idsfile, = args
    ids = set(x.strip() for x in open(idsfile))
    data = get_bed_from_phytozome(list(ids))

    pf = idsfile.rsplit(".", 1)[0]
    bedfile = pf + ".bed"
    fw = open(bedfile, "w")
    for i, row in enumerate(data):
        row = row.strip()
        if row == "":
            continue

        print(row, file=fw)

    logging.debug("A total of {0} records written to `{1}`.".format(i + 1, bedfile))
991339
def _prep_callable_bed(in_file, work_dir, stats, data):
    
    out_file = os.path.join(work_dir, "%s-merge.bed.gz" % utils.splitext_plus(os.path.basename(in_file))[0])
    gsort = config_utils.get_program("gsort", data)
    if not utils.file_uptodate(out_file, in_file):
        with file_transaction(data, out_file) as tx_out_file:
            fai_file = ref.fasta_idx(dd.get_ref_file(data))
            cmd = ("{gsort} {in_file} {fai_file} | bedtools merge -i - -d {stats[merge_size]} | "
                   "bgzip -c > {tx_out_file}")
            do.run(cmd.format(**locals()), "Prepare SV callable BED regions")
    return vcfutils.bgzip_and_index(out_file, data["config"])