%prog depth reads.bed features.bed

    Calculate depth depth per feature using coverageBed.
342567
def depth(args):
    
    p = OptionParser(depth.__doc__)
    p.set_outfile()
    opts, args = p.parse_args(args)

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

    readsbed, featsbed = args
    fp = open(featsbed)
    nargs = len(fp.readline().split("\t"))
    keepcols = ",".join(str(x) for x in range(1, nargs + 1))
    cmd = "coverageBed -a {0} -b {1} -d".format(readsbed, featsbed)
    cmd += " | groupBy -g {0} -c {1} -o mean".format(keepcols, nargs + 2)
    sh(cmd, outfile=opts.outfile)
1209385
def frombed(args):
    
    from jcvi.formats.fasta import Fasta
    from jcvi.formats.bed import Bed
    from jcvi.utils.cbook import fill

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

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

    bedfile, contigfasta, readfasta = args
    prefix = bedfile.rsplit(".", 1)[0]
    contigfile = prefix + ".contig"
    idsfile = prefix + ".ids"

    contigfasta = Fasta(contigfasta)
    readfasta = Fasta(readfasta)

    bed = Bed(bedfile)
    checksum = "00000000 checksum."
    fw_ids = open(idsfile, "w")
    fw = open(contigfile, "w")

    for ctg, reads in bed.sub_beds():
        ctgseq = contigfasta[ctg]
        ctgline = "##{0} {1} {2} bases, {3}".format(\
                ctg, len(reads), len(ctgseq), checksum)

        print(ctg, file=fw_ids)
        print(ctgline, file=fw)
        print(fill(ctgseq.seq), file=fw)

        for b in reads:
            read = b.accn
            strand = b.strand
            readseq = readfasta[read]
            rc = " [RC]" if strand == "-" else ""
            readlen = len(readseq)
            rstart, rend = 1, readlen
            if strand == "-":
                rstart, rend = rend, rstart

            readrange = "{{{0} {1}}}".format(rstart, rend)
            conrange = "<{0} {1}>".format(b.start, b.end)
            readline = "#{0}(0){1} {2} bases, {3} {4} {5}".format(\
                    read, rc, readlen, checksum, readrange, conrange)
            print(readline, file=fw)
            print(fill(readseq.seq), file=fw)

    logging.debug("Mapped contigs written to `{0}`.".format(contigfile))
    logging.debug("Contig IDs written to `{0}`.".format(idsfile))
1086400
def coverage(args):
    
    from jcvi.graphics.base import savefig

    p = OptionParser(coverage.__doc__)
    opts, args, iopts = p.set_image_options(args, format="png")

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

    covfile, = args
    df = pd.read_csv(covfile, sep='\t', names=["Ref", "Position", "Depth"])

    xlabel, ylabel = "Position", "Depth"
    df.plot(xlabel, ylabel, color='g')

    image_name = covfile + "." + iopts.format
    savefig(image_name)
1361386
function readRawFeature(pbf, layer, i) {
  pbf.pos = layer.features[i];
  const end = pbf.readVarint() + pbf.pos;

  const feature = {
    layer: layer,
    type: 0,
    properties: {}
  };
  pbf.readFields(featurePBFReader, feature, end);
  return feature;
}
135922
def bed(args):
    
    p = OptionParser(bed.__doc__)

    opts, args = p.parse_args(args)

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

    flist = args
    prefix = flist[0].split(".")[0]

    j = 0
    for f in flist:
        reader = Maf(f).reader
        for rec in reader:
            a, b = rec.components

            for a, tag in zip((a, b), "ab"):
                name = "{0}_{1:07d}{2}".format(prefix, j, tag)
                print("\t".join(str(x) for x in (a.src, a.forward_strand_start, \
                        a.forward_strand_end, name)))

            j += 1
1048187
def v_depth(d, depth):
        
        if depth == 0:
            yield d
        else:
            for node in DictTree.v(d):
                for node1 in DictTree.v_depth(node, depth-1):
                    yield node1
231663
def bed(args):
    '''
    
    '''
    from jcvi.utils.cbook import gene_name

    p = OptionParser(bed.__doc__)
    p.add_option("--type", dest="type", default="gene",
            help="Feature type to extract, use comma for multiple [default: %default]")
    p.add_option("--key", default="ID", help="Key in the attributes to extract")
    p.add_option("--source",
            help="Source to extract from, use comma for multiple [default: %default]")
    p.add_option("--span", default=False, action="store_true",
            help="Use feature span in the score column")
    p.add_option("--score_attrib", dest="score_attrib", default=False,
            help="Attribute whose value is to be used as score in `bedline` [default: %default]")
    p.add_option("--append_source", default=False, action="store_true",
            help="Append GFF source name to extracted key value")
    p.add_option("--append_ftype", default=False, action="store_true",
            help="Append GFF feature type to extracted key value")
    p.add_option("--nosort", default=False, action="store_true",
            help="Do not sort the output bed file [default: %default]")
    p.set_stripnames(default=False)
    p.set_outfile()

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

    gffile, = args
    key = opts.key or None
    type = opts.type or set()
    source = opts.source or set()
    strip_names = opts.strip_names
    span = opts.span

    if opts.type:
        type = set(x.strip() for x in opts.type.split(","))
    if opts.source:
        source = set(x.strip() for x in opts.source.split(","))

    gff = Gff(gffile, key=key, append_source=opts.append_source, \
        append_ftype=opts.append_ftype, score_attrib=opts.score_attrib)
    b = Bed()

    for g in gff:
        if type and g.type not in type:
            continue
        if source and g.source not in source:
            continue

        bl = g.bedline
        if strip_names:
            bl.accn = gene_name(bl.accn)
        if span:
            bl.score = bl.span
        b.append(bl)

    sorted = not opts.nosort
    b.print_to_file(opts.outfile, sorted=sorted)
    logging.debug("Extracted {0} features (type={1} id={2})".\
                    format(len(b), ",".join(type), key))
1283315
def _calc_regional_coverage(in_bam, chrom, start, end, samplename, work_dir, data):
    
    region_bt = pybedtools.BedTool("%s\t%s\t%s\n" % (chrom, start, end), from_string=True).saveas()
    region_file = region_bt.fn
    coords = "%s:%s-%s" % (chrom, start, end)
    tx_tmp_file = os.path.join(work_dir, "coverage-%s-%s.txt" % (samplename, coords.replace(":", "_")))
    samtools = config_utils.get_program("samtools", data)
    bedtools = config_utils.get_program("bedtools", data)
    cmd = ("{samtools} view -b {in_bam} {coords} | "
           "{bedtools} coverage -a {region_file} -b - -d > {tx_tmp_file}")
    do.run(cmd.format(**locals()), "Plotting coverage for %s %s" % (samplename, coords))
    names = ["chom", "start", "end", "offset", "coverage"]
    df = pd.io.parsers.read_table(tx_tmp_file, sep="\t", header=None,
                                  names=names).dropna()
    os.remove(tx_tmp_file)
    df["sample"] = samplename
    df["chrom"] = chrom
    df["position"] = df["start"] + df["offset"] - 1
    return df[["chrom", "position", "coverage", "sample"]]
1119558
function getFeature(args) {
    var item = args.item,
      params = args.params,
      propFunc = args.propFunc;

    var feature = { "type": "Feature" };

    feature.geometry = buildGeom(item, params);
    feature.properties = propFunc.call(item);

    return feature;
  }
242093
def len_on_depth(d, depth):
        
        counter = 0
        for node in DictTree.v_depth(d, depth-1):
            counter += DictTree.length(node)
        return counter
204639
def frombed(args):
    
    p = OptionParser(frombed.__doc__)
    p.add_option("--type", default="match",
                 help="GFF feature type [default: %default]")
    p.add_option("--source", default="default",
                help="GFF source qualifier [default: %default]")
    opts, args = p.parse_args(args)

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

    bedfile, = args
    bed = Bed(bedfile)

    for b in bed:
        print(b.gffline(type=opts.type, source=opts.source))
1149615
def multilineplot(args):
    
    p = OptionParser(multilineplot.__doc__)
    p.add_option("--lines",
                 help="Features to plot in lineplot [default: %default]")
    p.add_option("--colors",
                 help="List of colors matching number of input bed files")
    p.add_option("--mode", default="span", choices=("span", "count", "score"),
                 help="Accumulate feature based on [default: %default]")
    p.add_option("--binned", default=False, action="store_true",
                 help="Specify whether the input is already binned; " +
                 "if True, input files are considered to be binfiles")
    p.add_option("--ymax", type="int", help="Set Y-axis max")
    add_window_options(p)
    opts, args, iopts = p.set_image_options(args, figsize="8x5")

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

    fastafile, chr = args
    window, shift, subtract, merge = check_window_options(opts)
    linebeds = []
    colors = opts.colors
    if opts.lines:
        lines = opts.lines.split(",")
        assert len(colors) == len(lines), "Number of chosen colors must match" + \
                " number of input bed files"
        linebeds = get_beds(lines, binned=opts.binned)

    linebins = get_binfiles(linebeds, fastafile, shift, mode=opts.mode,
                            binned=opts.binned, merge=merge)

    clen = Sizes(fastafile).mapping[chr]
    nbins = get_nbins(clen, shift)

    plt.rcParams["xtick.major.size"] = 0
    plt.rcParams["ytick.major.size"] = 0
    plt.rcParams["figure.figsize"] = iopts.w, iopts.h

    fig, axarr = plt.subplots(nrows=len(lines))
    if len(linebeds) == 1:
        axarr = (axarr, )
    fig.suptitle(latex(chr), color="darkslategray")

    for i, ax in enumerate(axarr):
        lineplot(ax, [linebins[i]], nbins, chr, window, shift, \
                color="{0}{1}".format(colors[i], 'r'))

    if opts.ymax:
        ax.set_ylim(0, opts.ymax)

    plt.subplots_adjust(hspace=0.5)

    image_name = chr + "." + iopts.format
    savefig(image_name, dpi=iopts.dpi, iopts=iopts)
741976
def _get_basin_depth_term(self, C, sites, period):
        
        return np.zeros(len(sites.vs30), dtype=float)
834742
def bed(args):
    
    p = OptionParser(bed.__doc__)
    opts, args = p.parse_args(args)

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

    frgscffile, = args
    bedfile = frgscffile.rsplit(".", 1)[0] + ".bed"
    fw = open(bedfile, "w")

    fp = open(frgscffile)
    for row in fp:
        f = FrgScfLine(row)
        print(f.bedline, file=fw)

    logging.debug("File written to `{0}`.".format(bedfile))

    return bedfile
1121809
def calculate_sv_coverage(data):
    
    calcfns = {"cnvkit": _calculate_sv_coverage_cnvkit, "gatk-cnv": _calculate_sv_coverage_gatk}
    from bcbio.structural import cnvkit
    data = utils.to_single_data(data)
    if not cnvkit.use_general_sv_bins(data):
        out_target_file, out_anti_file = (None, None)
    else:
        work_dir = utils.safe_makedir(os.path.join(dd.get_work_dir(data), "structural",
                                                   dd.get_sample_name(data), "bins"))
        out_target_file, out_anti_file = calcfns[cnvkit.bin_approach(data)](data, work_dir)
        if not os.path.exists(out_target_file):
            out_target_file, out_anti_file = (None, None)
    if "seq2c" in dd.get_svcaller(data):
        from bcbio.structural import seq2c
        seq2c_target = seq2c.precall(data)
    else:
        seq2c_target = None

    if not tz.get_in(["depth", "bins"], data):
        data = tz.update_in(data, ["depth", "bins"], lambda x: {})
    data["depth"]["bins"] = {"target": out_target_file, "antitarget": out_anti_file, "seq2c": seq2c_target}
    return [[data]]
280236
def closest(args):
    
    p = OptionParser(closest.__doc__)
    p.add_option("--om", default=False, action="store_true",
                 help="The bedfile is OM blocks [default: %default]")
    opts, args = p.parse_args(args)

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

    candidates, gapsbed, fastafile = args
    sizes = Sizes(fastafile).mapping
    bed = Bed(candidates)
    ranges = []
    for b in bed:
        r = range_parse(b.accn) if opts.om else b
        ranges.append([r.seqid, r.start, r.end])

    gapsbed = Bed(gapsbed)
    granges = [(x.seqid, x.start, x.end) for x in gapsbed]

    ranges = range_merge(ranges)
    for r in ranges:
        a = range_closest(granges, r)
        b = range_closest(granges, r, left=False)
        seqid = r[0]

        if a is not None and a[0] != seqid:
            a = None
        if b is not None and b[0] != seqid:
            b = None

        mmin = 1 if a is None else a[1]
        mmax = sizes[seqid] if b is None else b[2]

        print("\t".join(str(x) for x in (seqid, mmin - 1, mmax)))
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()