%prog bed agpfile

    print out the tiling paths in bed/gff3 format
161490
def bed(args):
    
    from jcvi.formats.obo import validate_term

    p = OptionParser(bed.__doc__)
    p.add_option("--gaps", default=False, action="store_true",
            help="Only print bed lines for gaps [default: %default]")
    p.add_option("--nogaps", default=False, action="store_true",
            help="Do not print bed lines for gaps [default: %default]")
    p.add_option("--bed12", default=False, action="store_true",
            help="Produce bed12 formatted output [default: %default]")
    p.add_option("--component", default=False, action="store_true",
            help="Generate bed file for components [default: %default]")
    p.set_outfile()
    g1 = OptionGroup(p, "GFF specific parameters",
            "Note: If not specified, output will be in `bed` format")
    g1.add_option("--gff", default=False, action="store_true",
            help="Produce gff3 formatted output. By default, ignores " +\
                 "AGP gap lines. [default: %default]")
    g1.add_option("--source", default="MGSC",
            help="Specify a gff3 source [default: `%default`]")
    g1.add_option("--feature", default="golden_path_fragment",
            help="Specify a gff3 feature type [default: `%default`]")
    p.add_option_group(g1)
    p.set_SO_opts()

    opts, args = p.parse_args(args)

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

    if opts.component:
        opts.nogaps = True

    # If output format is gff3 and 'verifySO' option is invoked, validate the SO term
    if opts.gff and opts.verifySO:
        validate_term(opts.feature, method=opts.verifySO)

    agpfile, = args
    agp = AGP(agpfile)
    fw = must_open(opts.outfile, "w")
    if opts.gff:
        print("##gff-version 3", file=fw)

    for a in agp:
        if opts.nogaps and a.is_gap:
            continue
        if opts.gaps and not a.is_gap:
            continue
        if opts.bed12:
            print(a.bed12line, file=fw)
        elif opts.gff:
            print(a.gffline(gff_source=opts.source, gff_feat_type=opts.feature), file=fw)
        elif opts.component:
            name = "{0}:{1}-{2}".\
                        format(a.component_id, a.component_beg, a.component_end)
            print("\t".join(str(x) for x in (a.component_id, a.component_beg - 1,
                                a.component_end, name,
                                a.component_type, a.orientation)), file=fw)
        else:
            print(a.bedline, file=fw)
    fw.close()

    return fw.name
948389
def gtf_to_bed(gtf, alt_out_dir=None):
    
    out_file = os.path.splitext(gtf)[0] + '.bed'
    if file_exists(out_file):
        return out_file
    if not os.access(os.path.dirname(out_file), os.W_OK | os.X_OK):
        if not alt_out_dir:
            raise IOError('Cannot write transcript BED output file %s' % out_file)
        else:
            out_file = os.path.join(alt_out_dir, os.path.basename(out_file))
    with open(out_file, "w") as out_handle:
        db = get_gtf_db(gtf)
        for feature in db.features_of_type('transcript', order_by=("seqid", "start", "end")):
            chrom = feature.chrom
            start = feature.start
            end = feature.end
            attributes = feature.attributes.keys()
            strand = feature.strand
            name = (feature['gene_name'][0] if 'gene_name' in attributes else
                    feature['gene_id'][0])
            line = "\t".join([str(x) for x in [chrom, start, end, name, ".",
                                               strand]])
            out_handle.write(line + "\n")
    return out_file
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"])
13945
def write_equalwidth_bedfile(bedfile, width, outfile):
    """"""

    BUFSIZE = 10000
    f = open(bedfile)
    out = open(outfile, "w")
    lines = f.readlines(BUFSIZE)
    line_count = 0
    while lines:
        for line in lines:
            line_count += 1
            if not line.startswith("#") and not line.startswith("track") and not line.startswith("browser"):
                vals = line.strip().split("\t")
                try:
                    start, end = int(vals[1]), int(vals[2])
                except ValueError:
                    print("Error on line %s while reading %s. Is the file in BED or WIG format?" % (line_count, bedfile))
                    sys.exit(1)

                start = (start + end) // 2 - (width // 2)
                # This shifts the center, but ensures the width is identical... maybe not ideal
                if start < 0:
                    start = 0
                end = start + width
                # Keep all the other information in the bedfile if it's there
                if len(vals) > 3:
                    out.write("%s\t%s\t%s\t%s\n" % (vals[0], start, end, "\t".join(vals[3:])))
                else:
                    out.write("%s\t%s\t%s\n" % (vals[0], start, end))
        lines = f.readlines(BUFSIZE)
    
    out.close()
    f.close()
1143992
def estimategaps(args):
    
    p = OptionParser(estimategaps.__doc__)
    opts, args, iopts = p.set_image_options(args, figsize="6x6", dpi=300)

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

    pf, seqid, mlg = args
    bedfile = pf + ".lifted.bed"
    agpfile = pf + ".agp"

    function = lambda x: x.cm
    cc = Map(bedfile, scaffold_info=True, function=function)
    agp = AGP(agpfile)

    g = GapEstimator(cc, agp, seqid, mlg, function=function)
    pp, chrsize, mlgsize = g.pp, g.chrsize, g.mlgsize
    spl, spld = g.spl, g.spld
    g.compute_all_gaps(verbose=False)

    fig = plt.figure(1, (iopts.w, iopts.h))
    root = fig.add_axes([0, 0, 1, 1])

    # Panel A
    xstart, ystart = .15, .65
    w, h = .7, .3
    t = np.linspace(0, chrsize, 1000)
    ax = fig.add_axes([xstart, ystart, w, h])
    mx, my = zip(*g.scatter_data)
    rho = spearmanr(mx, my)

    dsg = "g"
    ax.vlines(pp, 0, mlgsize, colors="beige")
    ax.plot(mx, my, ".", color=set2[3])
    ax.plot(t, spl(t), "-", color=dsg)
    ax.text(.05, .95, mlg, va="top", transform=ax.transAxes)
    normalize_lms_axis(ax, xlim=chrsize, ylim=mlgsize,
                       ylabel="Genetic distance (cM)")
    if rho < 0:
        ax.invert_yaxis()

    # Panel B
    ystart -= .28
    h = .25
    ax = fig.add_axes([xstart, ystart, w, h])
    ax.vlines(pp, 0, mlgsize, colors="beige")
    ax.plot(t, spld(t), "-", lw=2, color=dsg)
    ax.plot(pp, spld(pp), "o", mfc="w", mec=dsg, ms=5)
    normalize_lms_axis(ax, xlim=chrsize, ylim=25 * 1e-6,
                       xfactor=1e-6, xlabel="Physical position (Mb)",
                       yfactor=1000000, ylabel="Recomb. rate\n(cM / Mb)")
    ax.xaxis.grid(False)

    # Panel C (specific to JMMale-1)
    a, b = "scaffold_1076", "scaffold_861"
    sizes = dict((x.component_id, (x.object_beg, x.object_end,
                                   x.component_span, x.orientation)) \
                                   for x in g.agp if not x.is_gap)
    a_beg, a_end, asize, ao = sizes[a]
    b_beg, b_end, bsize, bo = sizes[b]
    gapsize = g.get_gapsize(a)
    total_size = asize + gapsize + bsize
    ratio = .6 / total_size
    y = .16
    pad = .03
    pb_ratio = w / chrsize

    # Zoom
    lsg = "lightslategray"
    root.plot((.15 + pb_ratio * a_beg, .2),
              (ystart, ystart - .14), ":", color=lsg)
    root.plot((.15 + pb_ratio * b_end, .3),
              (ystart, ystart - .08), ":", color=lsg)
    ends = []
    for tag, size, marker, beg in zip((a, b), (asize, bsize), (49213, 81277),
                              (.2, .2 + (asize + gapsize) * ratio)):
        end = beg + size * ratio
        marker = beg + marker * ratio
        ends.append((beg, end, marker))
        root.plot((marker,), (y,), "o", color=lsg)
        root.text((beg + end) / 2, y + pad, latex(tag),
                  ha="center", va="center")
        HorizontalChromosome(root, beg, end, y, height=.025, fc='gainsboro')

    begs, ends, markers = zip(*ends)
    fontprop = dict(color=lsg, ha="center", va="center")
    ypos = y + pad * 2
    root.plot(markers, (ypos, ypos), "-", lw=2, color=lsg)
    root.text(sum(markers) / 2, ypos + pad,
              "Distance: 1.29cM $\Leftrightarrow$ 211,824bp (6.1 cM/Mb)", **fontprop)

    ypos = y - pad
    xx = markers[0], ends[0]
    root.plot(xx, (ypos, ypos), "-", lw=2, color=lsg)
    root.text(sum(xx) / 2, ypos - pad, "34,115bp", **fontprop)
    xx = markers[1], begs[1]
    root.plot(xx, (ypos, ypos), "-", lw=2, color=lsg)
    root.text(sum(xx) / 2, ypos - pad, "81,276bp", **fontprop)

    root.plot((ends[0], begs[1]), (y, y), ":", lw=2, color=lsg)
    root.text(sum(markers) / 2, ypos - 3 * pad, r"$\textit{Estimated gap size: 96,433bp}$",
                                  color="r", ha="center", va="center")

    labels = ((.05, .95, 'A'), (.05, .6, 'B'), (.05, .27, 'C'))
    panel_labels(root, labels)
    normalize_axes(root)

    pf = "estimategaps"
    image_name = pf + "." + iopts.format
    savefig(image_name, dpi=iopts.dpi, iopts=iopts)
1267497
def main():
    
    fmt,plot='svg',0
    if '-h' in sys.argv: # check if help is needed
        print(main.__doc__)
        sys.exit() # graceful quit
    elif '-f' in sys.argv: # ask for filename
        ind=sys.argv.index('-f')
        file=sys.argv[ind+1]
        f=open(file,'r')
        input=f.readlines()
    Data=[]
    for line in input:
        line.replace('\n','')
        if '\t' in line:   # read in the data from standard input
            rec=line.split('\t') # split each line on space to get records
        else:
            rec=line.split() # split each line on space to get records
        Data.append(float(rec[0]))
    
#
    if len(Data) >=10: 
        QQ={'unf1':1}
        pmagplotlib.plot_init(QQ['unf1'],5,5)
        pmagplotlib.plot_qq_unf(QQ['unf1'],Data,'QQ-Uniform') # make plot
    else:
        print('you need N> 10')
        sys.exit()
    pmagplotlib.draw_figs(QQ) 
    files={}
    for key in list(QQ.keys()):
        files[key]=key+'.'+fmt 
    if pmagplotlib.isServer:
        black     = '#000000'
        purple    = '#800080'
        titles={}
        titles['eq']='Equal Area Plot'
        EQ = pmagplotlib.add_borders(EQ,titles,black,purple)
        pmagplotlib.save_plots(QQ,files)
    elif plot==1:
        files['qq']=file+'.'+fmt 
        pmagplotlib.save_plots(QQ,files)
    else:
        ans=input(" S[a]ve to save plot, [q]uit without saving:  ")
        if ans=="a": pmagplotlib.save_plots(QQ,files)
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))
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
953076
def stats(args):
    
    from jcvi.utils.table import tabulate

    p = OptionParser(stats.__doc__)
    p.add_option("--warn", default=False, action="store_true",
                 help="Warnings on small component spans [default: %default]")
    opts, args = p.parse_args(args)

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

    agpfile, = args

    agp = AGP(agpfile)
    gap_lengths = []
    component_lengths = []
    for a in agp:
        span = a.object_span
        if a.is_gap:
            label = a.gap_type
            gap_lengths.append((span, label))
        else:
            label = "{0}:{1}-{2}".format(a.component_id, a.component_beg, \
                   a.component_end)
            component_lengths.append((span, label))
            if opts.warn and span < 50:
                logging.error("component span too small ({0}):\n{1}".\
                    format(span, a))

    table = dict()
    for label, lengths in zip(("Gaps", "Components"),
            (gap_lengths, component_lengths)):

        if not lengths:
            table[(label, "Min")] = table[(label, "Max")] \
                                  = table[(label, "Sum")] = "n.a."
            continue

        table[(label, "Min")] = "{0} ({1})".format(*min(lengths))
        table[(label, "Max")] = "{0} ({1})".format(*max(lengths))
        table[(label, "Sum")] = sum(x[0] for x in lengths)

    print(tabulate(table), file=sys.stderr)
11993
def merge(args):
    
    p = OptionParser(merge.__doc__)
    p.set_outfile()
    opts, args = p.parse_args(args)

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

    bedfiles = args
    fw = must_open(opts.outfile, "w")
    for bedfile in bedfiles:
        bed = Bed(bedfile)
        pf = op.basename(bedfile).split(".")[0]
        for b in bed:
            b.seqid = "_".join((pf, b.seqid))
            print(b, file=fw)
1099247
def agp(args):
    
    p = OptionParser(agp.__doc__)
    opts, args = p.parse_args(args)

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

    evidencefile, contigs = args
    ef = EvidenceFile(evidencefile, contigs)

    agpfile = evidencefile.replace(".evidence", ".agp")
    ef.write_agp(agpfile)
306271
def _printf(fmt, *args, **print3opts):
    '''
    '''
    if print3opts:  # like Python 3.0
        f = print3opts.get('file', None) or sys.stdout
        if args:
            f.write(fmt % args)
        else:
            f.write(fmt)
        f.write(print3opts.get('end', linesep))
    elif args:
        print(fmt % args)
    else:
        print(fmt)
3718
def gff(args):
    
    align_choices = ("left", "center", "right")
    p = OptionParser(gff.__doc__)
    p.add_option("--align", default="left", choices=align_choices,
                 help="Horizontal alignment [default: %default]")
    p.add_option("--noUTR", default=False, action="store_true",
                 help="Do not plot UTRs [default: %default]")
    opts, args = p.parse_args(args)

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

    fig = plt.figure(1, (8, 5))
    root = fig.add_axes([0, 0, 1, 1])

    gffiles = args
    ngenes = len(gffiles)

    canvas = .6
    setups, ratio = get_setups(gffiles, canvas=canvas, noUTR=opts.noUTR)
    align = opts.align
    xs = .2 if align == "left" else .8
    yinterval = canvas / ngenes
    ys = .8
    tip = .01
    for genename, mrnabed, cdsbeds in setups:
        ExonGlyph(root, xs, ys, mrnabed, cdsbeds, ratio=ratio, align=align)
        if align == "left":
            root.text(xs - tip, ys, genename, ha="right", va="center")
        elif align == "right":
            root.text(xs + tip, ys, genename, ha="left", va="center")
        ys -= yinterval

    root.set_xlim(0, 1)
    root.set_ylim(0, 1)
    root.set_axis_off()

    figname = "exons.pdf"
    savefig(figname, dpi=300)
169530
public static void main(String[] argv) throws IOException {
		File dir = new File(".");
		dir = new File("C:\\eclipse\\grobot-workspace\\anon\\res\\tiles\\indoor1");
		
		ArrayList list = new ArrayList();
		File[] files = dir.listFiles();
		for (int i=0;i<files.length;i++) {
			if (files[i].getName().endsWith(".png")) {
				if (!files[i].getName().startsWith("output")) {
					list.add(files[i]);
				}
			}
		}
		
		Pack packer = new Pack();
		packer.pack(list, 512, 512, 1, new File(dir, "output.png"));
		System.out.println("Output Generated.");
	}
553147
def write_agp(self, obj, sizes, fw=sys.stdout, gapsize=100,
                  gaptype="contig", evidence="map"):
        '''
        '''
        contigorder = [(x.contig_name, x.strand) for x in self]
        order_to_agp(obj, contigorder, sizes, fw,
                     gapsize=gapsize, gaptype=gaptype, evidence=evidence)
1162698
def allpaths(args):
    
    p = OptionParser(allpaths.__doc__)
    p.add_option("--ploidy", default="1", choices=("1", "2"),
                 help="Ploidy [default: %default]")
    opts, args = p.parse_args(args)

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

    folders = args
    for pf in folders:
        if not op.isdir(pf):
            continue
        assemble_dir(pf, target=["final.contigs.fasta", "final.assembly.fasta"],
                     ploidy=opts.ploidy)
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