%prog oropetium mcscan.out all.bed layout switch.ids

    Build a composite figure that calls graphis.synteny.
359532
def oropetium(args):
    
    p = OptionParser(oropetium.__doc__)
    p.add_option("--extra", help="Extra features in BED format")
    opts, args, iopts = p.set_image_options(args, figsize="9x6")

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

    datafile, bedfile, slayout, switch = args
    fig = plt.figure(1, (iopts.w, iopts.h))
    root = fig.add_axes([0, 0, 1, 1])

    Synteny(fig, root, datafile, bedfile, slayout,
            switch=switch, extra_features=opts.extra)

    # legend showing the orientation of the genes
    draw_gene_legend(root, .4, .57, .74, text=True, repeat=True)

    # On the left panel, make a species tree
    fc = 'lightslategrey'

    coords = {}
    xs, xp = .16, .03
    coords["oropetium"] = (xs, .7)
    coords["setaria"] = (xs, .6)
    coords["sorghum"] = (xs, .5)
    coords["rice"] = (xs, .4)
    coords["brachypodium"] = (xs, .3)
    xs -= xp
    coords["Panicoideae"] = join_nodes(root, coords, "setaria", "sorghum", xs)
    xs -= xp
    coords["BEP"] = join_nodes(root, coords, "rice", "brachypodium", xs)
    coords["PACMAD"] = join_nodes(root, coords, "oropetium", "Panicoideae", xs)
    xs -= xp
    coords["Poaceae"] = join_nodes(root, coords, "BEP", "PACMAD", xs)

    # Names of the internal nodes
    for tag in ("BEP", "Poaceae"):
        nx, ny = coords[tag]
        nx, ny = nx - .005, ny - .02
        root.text(nx, ny, tag, rotation=90, ha="right", va="top", color=fc)
    for tag in ("PACMAD",):
        nx, ny = coords[tag]
        nx, ny = nx - .005, ny + .02
        root.text(nx, ny, tag, rotation=90, ha="right", va="bottom", color=fc)

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

    pf = "oropetium"
    image_name = pf + "." + iopts.format
    savefig(image_name, dpi=iopts.dpi, iopts=iopts)
1345794
def omgprepare(args):
    
    from jcvi.formats.blast import cscore
    from jcvi.formats.base import DictFile

    p = OptionParser(omgprepare.__doc__)
    p.add_option("--norbh", action="store_true",
                 help="Disable RBH hits [default: %default]")
    p.add_option("--pctid", default=0, type="int",
                 help="Percent id cutoff for RBH hits [default: %default]")
    p.add_option("--cscore", default=90, type="int",
                 help="C-score cutoff for RBH hits [default: %default]")
    p.set_stripnames()
    p.set_beds()

    opts, args = p.parse_args(args)

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

    ploidy, anchorfile, blastfile = args
    norbh = opts.norbh
    pctid = opts.pctid
    cs = opts.cscore
    qbed, sbed, qorder, sorder, is_self = check_beds(anchorfile, p, opts)

    fp = open(ploidy)
    genomeidx = dict((x.split()[0], i) for i, x in enumerate(fp))
    fp.close()

    ploidy = DictFile(ploidy)

    geneinfo(qbed, qorder, genomeidx, ploidy)
    geneinfo(sbed, sorder, genomeidx, ploidy)

    pf = blastfile.rsplit(".", 1)[0]
    cscorefile = pf + ".cscore"
    cscore([blastfile, "-o", cscorefile, "--cutoff=0", "--pct"])
    ac = AnchorFile(anchorfile)
    pairs = set((a, b) for a, b, i in ac.iter_pairs())
    logging.debug("Imported {0} pairs from `{1}`.".format(len(pairs), anchorfile))

    weightsfile = pf + ".weights"
    fp = open(cscorefile)
    fw = open(weightsfile, "w")
    npairs = 0
    for row in fp:
        a, b, c, pct = row.split()
        c, pct = float(c), float(pct)
        c = int(c * 100)
        if (a, b) not in pairs:
            if norbh:
                continue
            if c < cs:
                continue
            if pct < pctid:
                continue
            c /= 10  # This severely penalizes RBH against synteny

        print("\t".join((a, b, str(c))), file=fw)
        npairs += 1
    fw.close()

    logging.debug("Write {0} pairs to `{1}`.".format(npairs, weightsfile))
1077647
public void startMcGregorIteration(int largestMappingSize, Map<Integer, Integer> presentMapping)
            throws IOException {

        this.globalMCSSize = (largestMappingSize / 2);
        List<String> cTab1Copy = McGregorChecks.generateCTabCopy(source);
        List<String> cTab2Copy = McGregorChecks.generateCTabCopy(target);

        //find mapped atoms of both molecules and store these in mappedAtoms
        List<Integer> mappedAtoms = new ArrayList<Integer>();
        //        System.out.println("\nMapped Atoms");
        for (Map.Entry<Integer, Integer> map : presentMapping.entrySet()) {
            //            System.out.println("i:" + map.getKey() + " j:" + map.getValue());
            mappedAtoms.add(map.getKey());
            mappedAtoms.add(map.getValue());
        }
        int mappingSize = presentMapping.size();

        List<Integer> iBondNeighborsA = new ArrayList<Integer>();
        List<String> cBondNeighborsA = new ArrayList<String>();

        List<Integer> iBondSetA = new ArrayList<Integer>();
        List<String> cBondSetA = new ArrayList<String>();

        List<Integer> iBondNeighborsB = new ArrayList<Integer>();
        List<Integer> iBondSetB = new ArrayList<Integer>();
        List<String> cBondNeighborsB = new ArrayList<String>();
        List<String> cBondSetB = new ArrayList<String>();

        //find unmapped atoms of molecule A

        List<Integer> unmappedAtomsMolA = McGregorChecks.markUnMappedAtoms(true, source, presentMapping);
        int counter = 0;
        int gSetBondNumA = 0;
        int gSetBondNumB = 0;
        int gNeighborBondnumA = 0; //number of remaining molecule A bonds after the clique search, which are neighbors of the MCS_1
        int gNeighborBondNumB = 0; //number of remaining molecule B bonds after the clique search, which are neighbors of the MCS_1

        QueryProcessor queryProcess = new QueryProcessor(cTab1Copy, cTab2Copy, SIGNS, gNeighborBondnumA,
                gSetBondNumA, iBondNeighborsA, cBondNeighborsA, mappingSize, iBondSetA, cBondSetA);

        if (!(source instanceof IQueryAtomContainer)) {
            queryProcess.process(source, target, unmappedAtomsMolA, mappedAtoms, counter);
        } else {
            queryProcess.process((IQueryAtomContainer) source, target, unmappedAtomsMolA, mappedAtoms, counter);
        }

        cTab1Copy = queryProcess.getCTab1();
        cTab2Copy = queryProcess.getCTab2();
        gSetBondNumA = queryProcess.getBondNumA();
        gNeighborBondnumA = queryProcess.getNeighborBondNumA();
        iBondNeighborsA = queryProcess.getIBondNeighboursA();
        cBondNeighborsA = queryProcess.getCBondNeighborsA();

        //find unmapped atoms of molecule B
        List<Integer> unmappedAtomsMolB = McGregorChecks.markUnMappedAtoms(false, target, presentMapping);

        //        System.out.println("unmappedAtomsMolB: " + unmappedAtomsMolB.size());

        //Extract bonds which are related with unmapped atoms of molecule B.
        //In case that unmapped atoms are connected with already mapped atoms, the mapped atoms are labelled with
        //new special signs -> the result are two vectors: cBondNeighborsA and int_bonds_molB, which contain those
        //bonds of molecule B, which are relevant for the McGregorBondTypeInSensitive algorithm.
        //The special signs must be transfered to the corresponding atoms of molecule A

        TargetProcessor targetProcess = new TargetProcessor(cTab1Copy, cTab2Copy, SIGNS, gNeighborBondNumB,
                gSetBondNumB, iBondNeighborsB, cBondNeighborsB, gNeighborBondnumA, iBondNeighborsA,
                cBondNeighborsA);

        targetProcess.process(target, unmappedAtomsMolB, mappingSize, iBondSetB, cBondSetB, mappedAtoms,
                counter);

        cTab1Copy = targetProcess.getCTab1();
        cTab2Copy = targetProcess.getCTab2();
        gSetBondNumB = targetProcess.getBondNumB();
        gNeighborBondNumB = targetProcess.getNeighborBondNumB();
        iBondNeighborsB = targetProcess.getIBondNeighboursB();
        cBondNeighborsB = targetProcess.getCBondNeighborsB();

        boolean dummy = false;

        McgregorHelper mcGregorHelper = new McgregorHelper(dummy, presentMapping.size(), mappedAtoms,
                gNeighborBondnumA, gNeighborBondNumB, iBondNeighborsA, iBondNeighborsB, cBondNeighborsA,
                cBondNeighborsB, gSetBondNumA, gSetBondNumB, iBondSetA, iBondSetB, cBondSetA, cBondSetB);
        iterator(mcGregorHelper);
    }
1046376
def benchmark(args):
    
    p = OptionParser(benchmark.__doc__)
    opts, args = p.parse_args(args)

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

    pf, bedfile = args
    truth = pf + ".truth"
    synfind = pf + ".synfind"
    mcscanx = pf + ".mcscanx"
    iadhore = pf + ".iadhore"
    orthofinder = pf + ".orthofinder"
    pivots = set([x.accn for x in Bed(bedfile)])

    fp = open(truth)
    truth = set()
    for row in fp:
        a, b = row.strip().split("\t")[:2]
        pivots.add(a)
        truth.add(tuple(sorted((a, b))))
    logging.debug("Truth: {0} pairs".format(len(truth)))

    fp = open(synfind)
    benchmarkfile = pf + ".benchmark"
    fw = must_open(benchmarkfile, "w")
    synfind = set()
    for row in fp:
        atoms = row.strip().split("\t")
        query, hit, tag = atoms[:3]
        if tag != "S":
            continue
        synfind.add(tuple(sorted((query, hit))))
    calc_sensitivity_specificity(synfind, truth, "SynFind", fw)

    fp = open(mcscanx)
    mcscanx = set()
    for row in fp:
        if row[0] == '#':
            continue
        atoms = row.strip().split(":")[1].split()
        query, hit = atoms[:2]
        mcscanx.add(tuple(sorted((query, hit))))
    calc_sensitivity_specificity(mcscanx, truth, "MCScanX", fw)

    fp = open(iadhore)
    iadhore = set()
    next(fp)
    for row in fp:
        atoms = row.strip().split("\t")
        query, hit = atoms[3:5]
        iadhore.add(tuple(sorted((query, hit))))
    calc_sensitivity_specificity(iadhore, truth, "iADHoRe", fw)

    fp = open(orthofinder)
    orthofinder = set()
    next(fp)
    for row in fp:
        row = row.replace('"', "")
        atoms = row.replace(",", " ").split()
        genes = [x.strip() for x in atoms if not x.startswith("OG")]
        genes = [gene_name(x) for x in genes]
        pps = [x for x in genes if x in pivots]
        for p in pps:
            for g in genes:
                if p == g:
                    continue
                orthofinder.add(tuple(sorted((p, g))))
    #write_pairs(orthofinder, "orthofinder.pairs")
    calc_sensitivity_specificity(orthofinder, truth, "OrthoFinder", fw)
    fw.close()
1119724
def fig3(args):
    
    from jcvi.formats.bed import Bed

    p = OptionParser(fig3.__doc__)
    p.add_option("--gauge_step", default=10000000, type="int",
                help="Step size for the base scale")
    opts, args, iopts = p.set_image_options(args, figsize="12x9")

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

    chrs, sizes, bedfile, datadir = args
    gauge_step = opts.gauge_step
    diverge = iopts.diverge
    rr, gg = diverge
    chrs = [[x] for x in chrs.split(",")]
    sizes = Sizes(sizes).mapping

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

    chr_sizes, chr_sum_sizes, ratio = calc_ratio(chrs, sizes)

    # Synteny panel
    seqidsfile = make_seqids(chrs)
    klayout = make_layout(chrs, chr_sum_sizes, ratio, template_f3a, shift=.05)
    height = .07
    r = height / 4
    K = Karyotype(fig, root, seqidsfile, klayout, gap=gap,
                  height=height, lw=2, generank=False, sizes=sizes,
                  heightpad=r, roundrect=True, plot_label=False)

    # Chromosome labels
    for kl in K.layout:
        if kl.empty:
            continue
        lx, ly = kl.xstart, kl.y
        if lx < .11:
            lx += .1
            ly += .06
        label = kl.label
        root.text(lx - .015, ly, label, fontsize=15,
                  ha="right", va="center")

    # Inset with datafiles
    datafiles = ("chrA02.bzh.forxmgr", "parent.A02.per10kb.forxmgr",
                 "parent.C2.per10kb.forxmgr", "chrC02.bzh.forxmgr")
    datafiles = [op.join(datadir, x) for x in datafiles]
    tracks = K.tracks
    hlfile = op.join(datadir, "bzh.regions.forhaibao")
    xy_axes = []
    for t, datafile in zip(tracks, datafiles):
        ax = make_affix_axis(fig, t, -r, height=2 * r)
        xy_axes.append(ax)
        chr = t.seqids[0]
        xy = XYtrack(ax, datafile, color="lightslategray")
        start, end = 0, t.total
        xy.interpolate(end)
        xy.cap(ymax=40)
        xy.import_hlfile(hlfile, chr, diverge=diverge)
        xy.draw()
        ax.set_xlim(start, end)
        gauge_ax = make_affix_axis(fig, t, -r)
        adjust_spines(gauge_ax, ["bottom"])
        setup_gauge_ax(gauge_ax, start, end, gauge_step)

    # Converted gene tracks
    ax_Ar = make_affix_axis(fig, tracks[1], r, height=r/2)
    ax_Co = make_affix_axis(fig, tracks[2], r, height=r/2)

    order = Bed(bedfile).order
    for asterisk in (False, True):
        conversion_track(order, "data/Genes.Converted.seuil.0.6.AtoC.txt",
                         0, "A02", ax_Ar, rr, asterisk=asterisk)
        conversion_track(order, "data/Genes.Converted.seuil.0.6.AtoC.txt",
                         1, "C2", ax_Co, gg, asterisk=asterisk)
        conversion_track(order, "data/Genes.Converted.seuil.0.6.CtoA.txt",
                         0, "A02", ax_Ar, gg, ypos=1, asterisk=asterisk)
        conversion_track(order, "data/Genes.Converted.seuil.0.6.CtoA.txt",
                         1, "C2", ax_Co, rr, ypos=1, asterisk=asterisk)

    Ar, Co = xy_axes[1:3]
    annotations = ((Ar, "Bra028920 Bra028897", "center", "1DAn2+"),
                   (Ar, "Bra020081 Bra020171", "right", "2DAn2+"),
                   (Ar, "Bra020218 Bra020286", "left", "3DAn2+"),
                   (Ar, "Bra008143 Bra008167", "left", "4DAn2-"),
                   (Ar, "Bra029317 Bra029251", "right", "5DAn2+ (GSL)"),
                   (Co, "Bo2g001000 Bo2g001300", "left", "1DCn2-"),
                   (Co, "Bo2g018560 Bo2g023700", "right", "2DCn2-"),
                   (Co, "Bo2g024450 Bo2g025390", "left", "3DCn2-"),
                   (Co, "Bo2g081060 Bo2g082340", "left", "4DCn2+"),
                   (Co, "Bo2g161510 Bo2g164260", "right", "5DCn2-"))

    for ax, genes, ha, label in annotations:
        g1, g2 = genes.split()
        x1, x2 = order[g1][1].start, order[g2][1].start
        if ha == "center":
            x = (x1 + x2) / 2 * .8
        elif ha == "left":
            x = x2
        else:
            x = x1
        label = r"\textit{{{0}}}".format(label)
        color = rr if "+" in label else gg
        ax.text(x, 30, label, color=color, fontsize=9, ha=ha, va="center")

    ax_Ar.set_xlim(0, tracks[1].total)
    ax_Ar.set_ylim(-1, 1)
    ax_Co.set_xlim(0, tracks[2].total)
    ax_Co.set_ylim(-1, 1)

    # Plot coverage in resequencing lines
    gstep = 5000000
    order = "swede,kale,h165,yudal,aviso,abu,bristol".split(",")
    labels_dict = {"h165": "Resynthesized (H165)", "abu": "Aburamasari"}
    hlsuffix = "regions.forhaibao"
    chr1, chr2 = "chrA02", "chrC02"
    t1, t2 = tracks[0], tracks[-1]
    s1, s2 = sizes[chr1], sizes[chr2]

    canvas1 = (t1.xstart, .75, t1.xend - t1.xstart, .2)
    c = Coverage(fig, root, canvas1, chr1, (0, s1), datadir,
                 order=order, gauge=None, plot_chr_label=False,
                 gauge_step=gstep, palette="gray",
                 cap=40, hlsuffix=hlsuffix, labels_dict=labels_dict,
                 diverge=diverge)
    yys = c.yys
    x1, x2 = .37, .72
    tip = .02
    annotations = ((x1, yys[2] + .3 * tip, tip, tip / 2, "FLC"),
                   (x1, yys[3] + .6 * tip, tip, tip / 2, "FLC"),
                   (x1, yys[5] + .6 * tip, tip, tip / 2, "FLC"),
                   (x2, yys[0] + .9 * tip, -1.2 * tip, 0, "GSL"),
                   (x2, yys[4] + .9 * tip, -1.2 * tip, 0, "GSL"),
                   (x2, yys[6] + .9 * tip, -1.2 * tip, 0, "GSL"))

    arrowprops=dict(facecolor='black', shrink=.05, frac=.5,
                    width=1, headwidth=4)
    for x, y, dx, dy, label in annotations:
        label = r"\textit{{{0}}}".format(label)
        root.annotate(label, xy=(x, y), xytext=(x + dx, y + dy),
                      arrowprops=arrowprops, color=rr, fontsize=9,
                      ha="center", va="center")

    canvas2 = (t2.xstart, .05, t2.xend - t2.xstart, .2)
    Coverage(fig, root, canvas2, chr2, (0, s2), datadir,
                 order=order, gauge=None, plot_chr_label=False,
                 gauge_step=gstep, palette="gray",
                 cap=40, hlsuffix=hlsuffix, labels_dict=labels_dict,
                 diverge=diverge)

    pad = .03
    labels = ((.1, .67, "A"), (t1.xstart - 3 * pad, .95 + pad, "B"),
              (t2.xstart - 3 * pad, .25 + pad, "C"))
    panel_labels(root, labels)
    normalize_axes(root)

    image_name = "napus-fig3." + iopts.format
    savefig(image_name, dpi=iopts.dpi, iopts=iopts)
561361
def :
    
    if procSteps is not None:
        procSteps.addStep('Blot')

    blot_name = util.getSectionName(configObj, _blot_step_num_)

    # This can be called directly from MultiDrizle, so only execute if
    # switch has been turned on (no guarantee MD will check before calling).
    if configObj[blot_name]['blot']:
        paramDict = buildBlotParamDict(configObj)

        log.info('USER INPUT PARAMETERS for Blot Step:')
        util.printParams(paramDict, log=log)

        run_blot(imageObjectList, output_wcs.single_wcs, paramDict,
                 wcsmap=wcsmap)
    else:
        log.info('Blot step not performed.')

    if procSteps is not None:
        procSteps.endStep('Blot')
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)))
779929
def __get_node_text(self, diagram, node, fmt):
        '''
        
        '''
        fmt_proc = fmt

        # IF blocks
        while (1):
            if_block = re.search('<%if ([^%]*): ([^%]*)%>', fmt_proc)
            if (if_block == None):
                break

            # evaluate condition
            if_cond = if_block[1]
            if_val  = if_block[2]
            if (self.__eval_if_block(if_cond, node) == 0):
                if_val = ''
            fmt_proc = fmt_proc[:if_block.span()[0]] + if_val + fmt_proc[if_block.span()[1]:]

        # {node.ip} = best IP
        ip = node.get_ipaddr()
        fmt_proc = fmt_proc.replace('{node.ip}', ip)

        # stackwise
        stack_block = re.search('<%stack ([^%]*)%>', fmt_proc)
        if (stack_block != None):
            if (node.stack.count == 0):
                # no stackwise, remove this
                fmt_proc = fmt_proc[:stack_block.span()[0]] + fmt_proc[stack_block.span()[1]:]
            else:
                val = ''
                if (self.config.diagram.expand_stackwise == 0):
                    if (self.config.diagram.get_stack_members):
                        for smem in node.stack.members:
                            nval = stack_block[1]
                            nval = nval.replace('{stack.num}',    str(smem.num))
                            nval = nval.replace('{stack.plat}',   smem.plat)
                            nval = nval.replace('{stack.serial}', smem.serial)
                            nval = nval.replace('{stack.role}',   smem.role)
                            val += nval
                fmt_proc = fmt_proc[:stack_block.span()[0]] + val + fmt_proc[stack_block.span()[1]:]
        
        # loopbacks
        loopback_block = re.search('<%loopback ([^%]*)%>', fmt_proc)
        if (loopback_block != None):
            val = ''
            for lo in node.loopbacks:
                for lo_ip in lo.ips:
                    nval = loopback_block[1]
                    nval = nval.replace('{lo.name}', lo.name)
                    nval = nval.replace('{lo.ip}', lo_ip)
                    val += nval
            fmt_proc = fmt_proc[:loopback_block.span()[0]] + val + fmt_proc[loopback_block.span()[1]:]

        # SVIs
        svi_block = re.search('<%svi ([^%]*)%>', fmt_proc)
        if (svi_block != None):
            val = ''
            for svi in node.svis:
                for svi_ip in svi.ip:
                    nval = svi_block[1]
                    nval = nval.replace('{svi.vlan}', svi.vlan)
                    nval = nval.replace('{svi.ip}', svi_ip)
                    val += nval
            fmt_proc = fmt_proc[:svi_block.span()[0]] + val + fmt_proc[svi_block.span()[1]:]

        # replace {stack.} with magic
        fmt_proc = re.sub('{stack\.(([a-zA-Z])*)}', '$stack2354$\g<1>$stack2354$', fmt_proc)
        fmt_proc = re.sub('{vss\.(([a-zA-Z])*)}', '$vss2354$\g<1>$vss2354$', fmt_proc)

        # {node.} variables
        fmt_proc = fmt_proc.format(node=node)
        
        # replace magics
        fmt_proc = re.sub('\$stack2354\$(([a-zA-Z])*)\$stack2354\$', '{stack.\g<1>}', fmt_proc)
        fmt_proc = re.sub('\$vss2354\$(([a-zA-Z])*)\$vss2354\$', '{vss.\g<1>}', fmt_proc)

        return fmt_proc
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
929977
def cross(self,p_i,p_j, max_depth = 2):
        """"""
        # only choose crossover points for out_types available in both programs
        # pdb.set_trace()
        # determine possible outttypes
        types_p_i = [t for t in [p.out_type for p in p_i]]
        types_p_j = [t for t in [p.out_type for p in p_j]]
        types = set(types_p_i).intersection(types_p_j)

        # grab subtree of p_i
        p_i_sub = [i for i,n in enumerate(p_i) if n.out_type in types]
        x_i_end = self.random_state.choice(p_i_sub)
        x_i_begin = x_i_end
        arity_sum = p_i[x_i_end].arity[p_i[x_i_end].in_type]
        # print("x_i_end:",x_i_end)
        # i = 0
        while (arity_sum > 0): #and i < 1000:
            if x_i_begin == 0:
                print("arity_sum:",arity_sum,"x_i_begin:",x_i_begin,"x_i_end:",x_i_end)
            x_i_begin -= 1
            arity_sum += p_i[x_i_begin].arity[p_i[x_i_begin].in_type]-1
        #     i += 1
        # if i == 1000:
        #     print("in variation")
        #     pdb.set_trace()

        # grab subtree of p_j with matching out_type to p_i[x_i_end]
        p_j_sub = [i for i,n in enumerate(p_j) if n.out_type == p_i[x_i_end].out_type]
        x_j_end = self.random_state.choice(p_j_sub)
        x_j_begin = x_j_end
        arity_sum = p_j[x_j_end].arity[p_j[x_j_end].in_type]
        # i = 0
        while (arity_sum > 0): #and i < 1000:
            if x_j_begin == 0:
                print("arity_sum:",arity_sum,"x_j_begin:",x_j_begin,"x_j_end:",x_j_end)
                print("p_j:",p_j)
            x_j_begin -= 1
            arity_sum += p_j[x_j_begin].arity[p_j[x_j_begin].in_type]-1
        #     i += 1
        # if i == 1000:
        #     print("in variation")
        #     pdb.set_trace()
        #swap subtrees
        tmpi = p_i[:]
        tmpj = p_j[:]
        tmpi[x_i_begin:x_i_end+1:],tmpj[x_j_begin:x_j_end+1:] = \
                        tmpj[x_j_begin:x_j_end+1:],tmpi[x_i_begin:x_i_end+1:]

        if not self.is_valid_program(p_i) or not self.is_valid_program(p_j):
            # pdb.set_trace()
            print("parent 1:",p_i,"x_i_begin:",x_i_begin,"x_i_end:",x_i_end)
            print("parent 2:",p_j,"x_j_begin:",x_j_begin,"x_j_end:",x_j_end)
            print("child 1:",tmpi)
            print("child 2:",tmpj)
            raise ValueError('Crossover produced an invalid program.')

        # size check, then assignment
        if len(tmpi) <= 2**max_depth-1:
            p_i[:] = tmpi
        if len(tmpj) <= 2**max_depth-1:
            p_j[:] = tmpj
873152
def main():
    
    allowed_format = ("emf", "eps", "pdf", "png", "ps", \
                      "raw", "rgba", "svg", "svgz")
    p = OptionParser(main.__doc__)
    p.add_option("--skip", default=0, type="int",
            help="skip the first several lines [default: %default]")
    p.add_option("--col", default=0, type="int", help="Get the n-th column")
    p.set_histogram()
    p.add_option("--tags", dest="tags", default=None,
            help="tags for data if multiple input files, comma sep")
    p.add_option("--ascii", default=False, action="store_true",
            help="print ASCII text stem-leaf plot [default: %default]")
    p.add_option("--base", default="0", choices=("0", "2", "10"),
            help="use logarithm axis with base, 0 to disable [default: %default]")
    p.add_option("--facet", default=False, action="store_true",
            help="place multiple histograms side-by-side [default: %default]")
    p.add_option("--fill", default="white",
            help="color of the bin [default: %default]")
    p.add_option("--format", default="pdf", choices=allowed_format,
            help="Generate image of format [default: %default]")
    p.add_option("--quick", default=False, action="store_true",
            help="Use quick plot, assuming bins are already counted")
    p.add_option("--noprintstats", default=False, action="store_true",
            help="Write basic stats when using --quick")
    opts, args = p.parse_args()

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

    skip = opts.skip
    vmin, vmax = opts.vmin, opts.vmax
    bins = opts.bins
    xlabel, title = opts.xlabel, opts.title
    title = title or args[0]
    base = int(opts.base)
    fileno = len(args)

    if opts.quick:
        assert fileno == 1, "Single input file expected using --quick"
        filename = args[0]
        figname = filename.rsplit(".", 1)[0] + ".pdf"
        data = DictFile(filename, keycast=int, cast=int)
        quickplot(data, vmin, vmax, xlabel, title, figname=figname,
                  print_stats=(not opts.noprintstats))
        return

    if fileno == 1:
        histogram(args[0], vmin, vmax, xlabel, title, outfmt=opts.format,
                bins=bins, skip=skip, ascii=opts.ascii,
                base=base, fill=opts.fill, col=opts.col)
    else:
        histogram_multiple(args, vmin, vmax, xlabel, title, outfmt=opts.format,
                tags=opts.tags, bins=bins, skip=skip, ascii=opts.ascii,
                facet=opts.facet, fill=opts.fill)
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)
442648
def _importSNPs_AgnosticSNP(setName, species, genomeSource, snpsFile) :
	""
	printf('importing SNP set %s for species %s...' % (setName, species))

	snpData = CSVFile()
	snpData.parse(snpsFile, separator = "\t")

	AgnosticSNP.dropIndex(('start', 'chromosomeNumber', 'setName'))
	conf.db.beginTransaction()
	
	pBar = ProgressBar(len(snpData))
	pLabel = ''
	currChrNumber = None
	for snpEntry in snpData :
		tmpChr = snpEntry['chromosomeNumber']
		if tmpChr != currChrNumber :
			currChrNumber = tmpChr
			pLabel = 'Chr %s...' % currChrNumber

		snp = AgnosticSNP()
		snp.species = species
		snp.setName = setName
		for f in snp.getFields() :
			try :
				setattr(snp, f, snpEntry[f])
			except KeyError :
				if f != 'species' and f != 'setName' :
					printf("Warning filetype as no key %s", f)
		snp.quality = float(snp.quality)
		snp.start = int(snp.start)
		snp.end = int(snp.end)
		snp.save()
		pBar.update(label = pLabel)

	pBar.close()
	
	snpMaster = SNPMaster()
	snpMaster.set(setName = setName, SNPType = 'AgnosticSNP', species = species)
	snpMaster.save()

	printf('saving...')
	conf.db.endTransaction()
	printf('creating indexes...')
	AgnosticSNP.ensureGlobalIndex(('start', 'chromosomeNumber', 'setName'))
	printf('importation of SNP set %s for species %s done.' %(setName, species))
	
	return True
55629
def Main_Frame(self):
        

        #--------------------------------------------------------------------
        # initialize first specimen in list as current specimen
        #--------------------------------------------------------------------

        try:
            self.s = self.specimens[0]
        except:
            self.s = ""
            print("No specimens during UI build")

        #--------------------------------------------------------------------
        # create main panel in the right size
        #--------------------------------------------------------------------

        dw, dh = wx.DisplaySize()
        w, h = self.GetSize()
        r1 = dw / 1250.
        r2 = dw / 750.

        GUI_RESOLUTION = min(r1, r2, 1.3)
        if 'gui_resolution' in list(self.preferences.keys()):
            if float(self.preferences['gui_resolution']) != 1:
                self.GUI_RESOLUTION = float(
                    self.preferences['gui_resolution']) / 100
            else:
                self.GUI_RESOLUTION = min(r1, r2, 1.3)
        else:
            self.GUI_RESOLUTION = min(r1, r2, 1.3)

        #--------------------------------------------------------------------
        # adjust font size
        #--------------------------------------------------------------------

        self.font_type = "Arial"
        if sys.platform.startswith("linux"):
            self.font_type = "Liberation Serif"

        if self.GUI_RESOLUTION >= 1.1 and self.GUI_RESOLUTION <= 1.3:
            font2 = wx.Font(13, wx.SWISS, wx.NORMAL,
                            wx.NORMAL, False, self.font_type)
        elif self.GUI_RESOLUTION <= 0.9 and self.GUI_RESOLUTION < 1.0:
            font2 = wx.Font(11, wx.SWISS, wx.NORMAL,
                            wx.NORMAL, False, self.font_type)
        elif self.GUI_RESOLUTION <= 0.9:
            font2 = wx.Font(10, wx.SWISS, wx.NORMAL,
                            wx.NORMAL, False, self.font_type)
        else:
            font2 = wx.Font(12, wx.SWISS, wx.NORMAL,
                            wx.NORMAL, False, self.font_type)
        print("    self.GUI_RESOLUTION", self.GUI_RESOLUTION)

        h_space = 4
        v_space = 4

        # set font size and style
        #font1 = wx.Font(10, wx.SWISS, wx.NORMAL, wx.NORMAL, False, u'Comic Sans MS')
        FONT_RATIO = self.GUI_RESOLUTION + (self.GUI_RESOLUTION - 1) * 5
        font1 = wx.Font(9 + FONT_RATIO, wx.SWISS, wx.NORMAL,
                        wx.NORMAL, False, self.font_type)
        # GUI headers

        font3 = wx.Font(11 + FONT_RATIO, wx.SWISS, wx.NORMAL,
                        wx.NORMAL, False, self.font_type)
        font = wx.SystemSettings.GetFont(wx.SYS_SYSTEM_FONT)
        font.SetPointSize(10 + FONT_RATIO)

        #--------------------------------------------------------------------
        # Create Figures and FigCanvas objects.
        #--------------------------------------------------------------------

        self.fig1 = Figure((5. * self.GUI_RESOLUTION, 5. *
                            self.GUI_RESOLUTION), dpi=self.dpi)
        self.canvas1 = FigCanvas(self.plot_panel, wx.ID_ANY, self.fig1)
        self.fig1.text(0.01, 0.98, "Arai plot", {
                       'family': self.font_type, 'fontsize': 10, 'style': 'normal', 'va': 'center', 'ha': 'left'})
        self.toolbar1 = NavigationToolbar(self.canvas1)
        self.toolbar1.Hide()
        self.fig1_setting = "Zoom"
        self.toolbar1.zoom()
        self.canvas1.Bind(wx.EVT_RIGHT_DOWN, self.on_right_click_fig)
        self.canvas1.Bind(wx.EVT_MIDDLE_DOWN, self.on_home_fig)

        self.fig2 = Figure((2.5 * self.GUI_RESOLUTION, 2.5 *
                            self.GUI_RESOLUTION), dpi=self.dpi)
        self.canvas2 = FigCanvas(self.plot_panel, wx.ID_ANY, self.fig2)
        self.fig2.text(0.02, 0.96, "Zijderveld", {
                       'family': self.font_type, 'fontsize': 10, 'style': 'normal', 'va': 'center', 'ha': 'left'})
        self.toolbar2 = NavigationToolbar(self.canvas2)
        self.toolbar2.Hide()
        self.fig2_setting = "Zoom"
        self.toolbar2.zoom()
        self.canvas2.Bind(wx.EVT_RIGHT_DOWN, self.on_right_click_fig)
        self.canvas2.Bind(wx.EVT_MIDDLE_DOWN, self.on_home_fig)

        self.fig3 = Figure((2.5 * self.GUI_RESOLUTION, 2.5 *
                            self.GUI_RESOLUTION), dpi=self.dpi)
        self.canvas3 = FigCanvas(self.plot_panel, wx.ID_ANY, self.fig3)
        #self.fig3.text(0.02,0.96,"Equal area",{'family':self.font_type, 'fontsize':10*self.GUI_RESOLUTION, 'style':'normal','va':'center', 'ha':'left' })
        self.toolbar3 = NavigationToolbar(self.canvas3)
        self.toolbar3.Hide()
        self.fig3_setting = "Zoom"
        self.toolbar3.zoom()
        self.canvas3.Bind(wx.EVT_RIGHT_DOWN, self.on_right_click_fig)
        self.canvas3.Bind(wx.EVT_MIDDLE_DOWN, self.on_home_fig)

        self.fig4 = Figure((2.5 * self.GUI_RESOLUTION, 2.5 *
                            self.GUI_RESOLUTION), dpi=self.dpi)
        self.canvas4 = FigCanvas(self.plot_panel, wx.ID_ANY, self.fig4)
        if self.acceptance_criteria['average_by_sample_or_site']['value'] == 'site':
            TEXT = "Site data"
        else:
            TEXT = "Sample data"
        self.fig4.text(0.02, 0.96, TEXT, {
                       'family': self.font_type, 'fontsize': 10, 'style': 'normal', 'va': 'center', 'ha': 'left'})
        self.toolbar4 = NavigationToolbar(self.canvas4)
        self.toolbar4.Hide()
        self.fig4_setting = "Zoom"
        self.toolbar4.zoom()
        self.canvas4.Bind(wx.EVT_RIGHT_DOWN, self.on_right_click_fig)
        self.canvas4.Bind(wx.EVT_MIDDLE_DOWN, self.on_home_fig)

        self.fig5 = Figure((2.5 * self.GUI_RESOLUTION, 2.5 *
                            self.GUI_RESOLUTION), dpi=self.dpi)
        self.canvas5 = FigCanvas(self.plot_panel, wx.ID_ANY, self.fig5)
        #self.fig5.text(0.02,0.96,"M/M0",{'family':self.font_type, 'fontsize':10, 'style':'normal','va':'center', 'ha':'left' })
        self.toolbar5 = NavigationToolbar(self.canvas5)
        self.toolbar5.Hide()
        self.fig5_setting = "Zoom"
        self.toolbar5.zoom()
        self.canvas5.Bind(wx.EVT_RIGHT_DOWN, self.on_right_click_fig)
        self.canvas5.Bind(wx.EVT_MIDDLE_DOWN, self.on_home_fig)

        # make axes of the figures
        self.araiplot = self.fig1.add_axes([0.1, 0.1, 0.8, 0.8])
        self.zijplot = self.fig2.add_subplot(111)
        self.eqplot = self.fig3.add_subplot(111)
        self.sampleplot = self.fig4.add_axes(
            [0.2, 0.3, 0.7, 0.6], frameon=True, facecolor='None')
        self.mplot = self.fig5.add_axes(
            [0.2, 0.15, 0.7, 0.7], frameon=True, facecolor='None')

        #--------------------------------------------------------------------
        # text box displaying measurement data
        #--------------------------------------------------------------------

        self.logger = wx.ListCtrl(self.side_panel, id=wx.ID_ANY, size=(
            100 * self.GUI_RESOLUTION, 100 * self.GUI_RESOLUTION), style=wx.LC_REPORT)
        self.logger.SetFont(font1)
        self.logger.InsertColumn(0, 'i', width=45 * self.GUI_RESOLUTION)
        self.logger.InsertColumn(1, 'Step', width=45 * self.GUI_RESOLUTION)
        self.logger.InsertColumn(2, 'Tr', width=65 * self.GUI_RESOLUTION)
        self.logger.InsertColumn(3, 'Dec', width=65 * self.GUI_RESOLUTION)
        self.logger.InsertColumn(4, 'Inc', width=65 * self.GUI_RESOLUTION)
        self.logger.InsertColumn(5, 'M', width=75 * self.GUI_RESOLUTION)
        self.Bind(wx.EVT_LIST_ITEM_ACTIVATED,
                  self.on_click_listctrl, self.logger)
        self.Bind(wx.EVT_LIST_ITEM_RIGHT_CLICK,
                  self.on_right_click_listctrl, self.logger)

        #--------------------------------------------------------------------
        # select a specimen box
        #--------------------------------------------------------------------

        # Combo-box with a list of specimen
        self.specimens_box = wx.ComboBox(self.top_panel, wx.ID_ANY, self.s, (250 * self.GUI_RESOLUTION, 25),
                                         wx.DefaultSize, self.specimens, wx.CB_DROPDOWN | wx.TE_PROCESS_ENTER, name="specimen")
        self.specimens_box.SetFont(font2)
        self.Bind(wx.EVT_COMBOBOX, self.onSelect_specimen, self.specimens_box)
        self.Bind(wx.EVT_TEXT_ENTER, self.onSelect_specimen,
                  self.specimens_box)

        # buttons to move forward and backwards from specimens
        nextbutton = wx.Button(self.top_panel, id=wx.ID_ANY, label='next', size=(
            75 * self.GUI_RESOLUTION, 25))  # ,style=wx.BU_EXACTFIT)#, size=(175, 28))
        self.Bind(wx.EVT_BUTTON, self.on_next_button, nextbutton)
        nextbutton.SetFont(font2)

        prevbutton = wx.Button(self.top_panel, id=wx.ID_ANY, label='previous', size=(
            75 * self.GUI_RESOLUTION, 25))  # ,style=wx.BU_EXACTFIT)#, size=(175, 28))
        prevbutton.SetFont(font2)
        self.Bind(wx.EVT_BUTTON, self.on_prev_button, prevbutton)

        #--------------------------------------------------------------------
        # select temperature bounds
        #--------------------------------------------------------------------

        try:
            if self.Data[self.s]['T_or_MW'] == "T":
                self.temperatures = np.array(
                    self.Data[self.s]['t_Arai']) - 273.
                self.T_list = ["%.0f" % T for T in self.temperatures]
            elif self.Data[self.s]['T_or_MW'] == "MW":
                self.temperatures = np.array(self.Data[self.s]['t_Arai'])
                self.T_list = ["%.0f" % T for T in self.temperatures]
        except (ValueError, TypeError, KeyError) as e:
            self.T_list = []

        self.tmin_box = wx.ComboBox(self.top_panel, wx.ID_ANY, size=(
            100 * self.GUI_RESOLUTION, 25), choices=self.T_list, style=wx.CB_DROPDOWN | wx.TE_READONLY)
        self.Bind(wx.EVT_COMBOBOX, self.get_new_T_PI_parameters, self.tmin_box)

        self.tmax_box = wx.ComboBox(self.top_panel, -1, size=(100 * self.GUI_RESOLUTION, 25),
                                    choices=self.T_list, style=wx.CB_DROPDOWN | wx.TE_READONLY)
        self.Bind(wx.EVT_COMBOBOX, self.get_new_T_PI_parameters, self.tmax_box)

        #--------------------------------------------------------------------
        # save/delete buttons
        #--------------------------------------------------------------------

        # save/delete interpretation buttons
        self.save_interpretation_button = wx.Button(self.top_panel, id=-1, label='save', size=(
            75 * self.GUI_RESOLUTION, 25))  # ,style=wx.BU_EXACTFIT)#, size=(175, 28))
        self.save_interpretation_button.SetFont(font2)
        self.delete_interpretation_button = wx.Button(self.top_panel, id=-1, label='delete', size=(
            75 * self.GUI_RESOLUTION, 25))  # ,style=wx.BU_EXACTFIT)#, size=(175, 28))
        self.delete_interpretation_button.SetFont(font2)
        self.Bind(wx.EVT_BUTTON, self.on_save_interpretation_button,
                  self.save_interpretation_button)
        self.Bind(wx.EVT_BUTTON, self.on_delete_interpretation_button,
                  self.delete_interpretation_button)

        self.auto_save = wx.CheckBox(self.top_panel, wx.ID_ANY, 'auto-save')
        self.auto_save_info = wx.Button(self.top_panel, wx.ID_ANY, "?")
        self.Bind(wx.EVT_BUTTON, self.on_info_click, self.auto_save_info)
        #self.auto_save_text = wx.StaticText(self.top_panel, wx.ID_ANY, label="(saves with 'next')")


        #--------------------------------------------------------------------
        # specimen interpretation and statistics window (Blab; Banc, Dec, Inc, correction factors etc.)
        #--------------------------------------------------------------------

        self.Blab_window = wx.TextCtrl(
            self.top_panel, style=wx.TE_CENTER | wx.TE_READONLY, size=(50 * self.GUI_RESOLUTION, 25))
        self.Banc_window = wx.TextCtrl(
            self.top_panel, style=wx.TE_CENTER | wx.TE_READONLY, size=(50 * self.GUI_RESOLUTION, 25))
        self.Aniso_factor_window = wx.TextCtrl(
            self.top_panel, style=wx.TE_CENTER | wx.TE_READONLY, size=(50 * self.GUI_RESOLUTION, 25))
        self.NLT_factor_window = wx.TextCtrl(
            self.top_panel, style=wx.TE_CENTER | wx.TE_READONLY, size=(50 * self.GUI_RESOLUTION, 25))
        self.CR_factor_window = wx.TextCtrl(
            self.top_panel, style=wx.TE_CENTER | wx.TE_READONLY, size=(50 * self.GUI_RESOLUTION, 25))
        self.declination_window = wx.TextCtrl(
            self.top_panel, style=wx.TE_CENTER | wx.TE_READONLY, size=(50 * self.GUI_RESOLUTION, 25))
        self.inclination_window = wx.TextCtrl(
            self.top_panel, style=wx.TE_CENTER | wx.TE_READONLY, size=(50 * self.GUI_RESOLUTION, 25))

        for stat in ['Blab', 'Banc', 'Aniso_factor', 'NLT_factor', 'CR_factor', 'declination', 'inclination']:
            exec("self.%s_window.SetBackgroundColour(wx.WHITE)" % stat)

        self.Blab_label = wx.StaticText(
            self.top_panel, label="\nB_lab", style=wx.ALIGN_CENTRE)
        self.Blab_label.SetFont(font2)
        self.Banc_label = wx.StaticText(
            self.top_panel, label="\nB_anc", style=wx.ALIGN_CENTRE)
        self.Banc_label.SetFont(font2)
        self.aniso_corr_label = wx.StaticText(
            self.top_panel, label="Aniso\ncorr", style=wx.ALIGN_CENTRE)
        self.aniso_corr_label.SetFont(font2)
        self.nlt_corr_label = wx.StaticText(
            self.top_panel, label="NLT\ncorr", style=wx.ALIGN_CENTRE)
        self.nlt_corr_label.SetFont(font2)
        self.cr_corr_label = wx.StaticText(
            self.top_panel, label="CR\ncorr", style=wx.ALIGN_CENTRE)
        self.cr_corr_label.SetFont(font2)
        self.dec_label = wx.StaticText(
            self.top_panel, label="\nDec", style=wx.ALIGN_CENTRE)
        self.dec_label.SetFont(font2)
        self.inc_label = wx.StaticText(
            self.top_panel, label="\nInc", style=wx.ALIGN_CENTRE)
        self.inc_label.SetFont(font2)

        # handle Specimen Results Sizer
        sizer_specimen_results = wx.StaticBoxSizer(wx.StaticBox(
            self.top_panel, wx.ID_ANY, "specimen results"), wx.HORIZONTAL)
        specimen_stat_window = wx.GridSizer(2, 7, h_space, v_space)
        specimen_stat_window.AddMany([(self.Blab_label, 1, wx.ALIGN_BOTTOM),
                                      ((self.Banc_label), 1, wx.ALIGN_BOTTOM),
                                      ((self.aniso_corr_label), 1, wx.ALIGN_BOTTOM),
                                      ((self.nlt_corr_label), 1, wx.ALIGN_BOTTOM),
                                      ((self.cr_corr_label), 1, wx.ALIGN_BOTTOM),
                                      ((self.dec_label), 1,
                                       wx.TE_CENTER | wx.ALIGN_BOTTOM),
                                      ((self.inc_label), 1, wx.ALIGN_BOTTOM),
                                      (self.Blab_window, 1, wx.EXPAND),
                                      (self.Banc_window, 1, wx.EXPAND),
                                      (self.Aniso_factor_window, 1, wx.EXPAND),
                                      (self.NLT_factor_window, 1, wx.EXPAND),
                                      (self.CR_factor_window, 1, wx.EXPAND),
                                      (self.declination_window, 1, wx.EXPAND),
                                      (self.inclination_window, 1, wx.EXPAND)])
        sizer_specimen_results.Add(
            specimen_stat_window, 1, wx.EXPAND | wx.ALIGN_LEFT, 0)

        #--------------------------------------------------------------------
        # Sample interpretation window
        #--------------------------------------------------------------------

        for key in ["sample_int_n", "sample_int_uT", "sample_int_sigma", "sample_int_sigma_perc"]:
            command = "self.%s_window=wx.TextCtrl(self.top_panel,style=wx.TE_CENTER|wx.TE_READONLY,size=(50*self.GUI_RESOLUTION,25))" % key
            exec(command)
            exec("self.%s_window.SetBackgroundColour(wx.WHITE)" % key)

        sample_mean_label = wx.StaticText(
            self.top_panel, label="\nmean", style=wx.TE_CENTER)
        sample_mean_label.SetFont(font2)
        sample_N_label = wx.StaticText(
            self.top_panel, label="\nN ", style=wx.TE_CENTER)
        sample_N_label.SetFont(font2)
        sample_std_label = wx.StaticText(
            self.top_panel, label="\nstd uT", style=wx.TE_CENTER)
        sample_std_label.SetFont(font2)
        sample_std_per_label = wx.StaticText(
            self.top_panel, label="\nstd %", style=wx.TE_CENTER)
        sample_std_per_label.SetFont(font2)

        # handle samples/sites results sizers
        sizer_sample_results = wx.StaticBoxSizer(wx.StaticBox(
            self.top_panel, wx.ID_ANY, "sample/site results"), wx.HORIZONTAL)
        sample_stat_window = wx.GridSizer(2, 4, h_space, v_space)
        sample_stat_window.AddMany([(sample_mean_label, 1, wx.ALIGN_BOTTOM),
                                    (sample_N_label, 1, wx.ALIGN_BOTTOM),
                                    (sample_std_label, 1, wx.ALIGN_BOTTOM),
                                    (sample_std_per_label, 1, wx.ALIGN_BOTTOM),
                                    (self.sample_int_uT_window, 1, wx.EXPAND),
                                    (self.sample_int_n_window, 1, wx.EXPAND),
                                    (self.sample_int_sigma_window, 1, wx.EXPAND),
                                    (self.sample_int_sigma_perc_window, 1, wx.EXPAND)])
        sizer_sample_results.Add(
            sample_stat_window, 1, wx.EXPAND | wx.ALIGN_LEFT, 0)

        #--------------------------------------------------------------------

        label_0 = wx.StaticText(
            self.bottom_panel, label="                    ", style=wx.ALIGN_CENTER, size=(180, 25))
        label_1 = wx.StaticText(
            self.bottom_panel, label="Acceptance criteria:", style=wx.ALIGN_CENTER, size=(180, 25))
        label_2 = wx.StaticText(
            self.bottom_panel, label="Specimen statistics:", style=wx.ALIGN_CENTER, size=(180, 25))

        for statistic in self.preferences['show_statistics_on_gui']:
            self.stat_windows[statistic] = wx.TextCtrl(
                self.bottom_panel, style=wx.TE_CENTER | wx.TE_READONLY, size=(50 * self.GUI_RESOLUTION, 25))
            self.stat_windows[statistic].SetBackgroundColour(wx.WHITE)
            self.stat_windows[statistic].SetFont(font2)
            self.threshold_windows[statistic] = wx.TextCtrl(
                self.bottom_panel, style=wx.TE_CENTER | wx.TE_READONLY, size=(50 * self.GUI_RESOLUTION, 25))
            self.threshold_windows[statistic].SetFont(font2)
            self.threshold_windows[statistic].SetBackgroundColour(wx.WHITE)
            label = statistic.replace("specimen_", "").replace("int_", "")
            self.stat_labels[statistic] = wx.StaticText(
                self.bottom_panel, label=label, style=wx.ALIGN_CENTRE_HORIZONTAL | wx.ALIGN_BOTTOM)
            self.stat_labels[statistic].SetFont(font2)

        #-------------------------------------------------------------------
        # Design the panels
        #-------------------------------------------------------------------

        # Plots Panel--------------------------------------------------------
        sizer_grid_plots = wx.GridSizer(2, 2, 0, 0)
        sizer_grid_plots.AddMany([(self.canvas2, 1, wx.EXPAND),
                                  (self.canvas4, 1, wx.EXPAND),
                                  (self.canvas3, 1, wx.EXPAND),
                                  (self.canvas5, 1, wx.EXPAND)])

        sizer_plots_outer = wx.BoxSizer(wx.HORIZONTAL)
        sizer_plots_outer.Add(self.canvas1, 1, wx.EXPAND)
        sizer_plots_outer.Add(sizer_grid_plots, 1, wx.EXPAND)

        # Top Bar Sizer-------------------------------------------------------
        #-------------Specimens Sizer----------------------------------------
        sizer_prev_next_btns = wx.BoxSizer(wx.HORIZONTAL)
        sizer_prev_next_btns.Add(prevbutton, 1, wx.EXPAND | wx.RIGHT, h_space)
        sizer_prev_next_btns.Add(nextbutton, 1, wx.EXPAND | wx.LEFT, h_space)

        sizer_select_specimen = wx.StaticBoxSizer(wx.StaticBox(
            self.top_panel, wx.ID_ANY, "specimen"), wx.VERTICAL)
        sizer_select_specimen.Add(
            self.specimens_box, 1, wx.EXPAND | wx.BOTTOM, v_space)
        sizer_select_specimen.Add(
            sizer_prev_next_btns, 1, wx.EXPAND | wx.TOP, v_space)

        #-------------Bounds Sizer----------------------------------------
        sizer_grid_bounds_btns = wx.GridSizer(2, 3, 2 * h_space, 2 * v_space)
        sizer_grid_bounds_btns.AddMany([(self.tmin_box, 1, wx.EXPAND),
                                        (self.save_interpretation_button,
                                         1, wx.EXPAND), (self.auto_save, 1, wx.EXPAND),
                                        (self.tmax_box, 1, wx.EXPAND),
                                        (self.delete_interpretation_button, 1, wx.EXPAND),
                                        (self.auto_save_info, 1, wx.EXPAND)])


        if self.s in list(self.Data.keys()) and self.Data[self.s]['T_or_MW'] == "T":
            sizer_select_temp = wx.StaticBoxSizer(wx.StaticBox(
                self.top_panel, wx.ID_ANY, "temperatures"), wx.HORIZONTAL)
        else:
            sizer_select_temp = wx.StaticBoxSizer(wx.StaticBox(
                self.top_panel, wx.ID_ANY, "MW power"), wx.HORIZONTAL)
        sizer_select_temp.Add(sizer_grid_bounds_btns, 1, wx.EXPAND)

        #-------------Top Bar Outer Sizer------------------------------------
        sizer_top_bar = wx.BoxSizer(wx.HORIZONTAL)
        sizer_top_bar.AddMany([(sizer_select_specimen, 1, wx.EXPAND | wx.ALIGN_LEFT | wx.RIGHT, 2 * h_space),
                               (sizer_select_temp, 1, wx.EXPAND |
                                wx.ALIGN_LEFT | wx.RIGHT, 2 * h_space),
                               (sizer_specimen_results, 2, wx.EXPAND |
                                wx.ALIGN_LEFT | wx.RIGHT, 2 * h_space),
                               (sizer_sample_results, 1, wx.EXPAND | wx.ALIGN_LEFT | wx.RIGHT, 0)])

        # Bottom Bar Sizer----------------------------------------------------
        #----------------Criteria Labels Sizer-------------------------------
        sizer_criteria_labels = wx.BoxSizer(wx.HORIZONTAL)
        sizer_criteria_labels.Add(label_0, 3, wx.EXPAND | wx.LEFT, 2 * h_space)
        sizer_criteria_boxes = wx.BoxSizer(wx.HORIZONTAL)
        sizer_criteria_boxes.Add(label_1, 3, wx.EXPAND | wx.LEFT, 2 * h_space)
        sizer_stats_boxes = wx.BoxSizer(wx.HORIZONTAL)
        sizer_stats_boxes.Add(label_2, 3, wx.EXPAND | wx.LEFT, 2 * h_space)
        for statistic in self.preferences['show_statistics_on_gui']:
            sizer_criteria_labels.Add(
                self.stat_labels[statistic], 1, wx.ALIGN_BOTTOM, 0)

        #----------------Acceptance Criteria Boxes---------------------------
            sizer_criteria_boxes.Add(
                self.threshold_windows[statistic], 1, wx.EXPAND | wx.LEFT, h_space)

        #----------------Specimen Statistics Boxes---------------------------
            sizer_stats_boxes.Add(
                self.stat_windows[statistic], 1, wx.EXPAND | wx.LEFT, h_space)

        #----------------Bottom Outer Sizer----------------------------------
        sizer_bottom_bar = wx.BoxSizer(wx.VERTICAL)
        sizer_bottom_bar.AddMany([(sizer_criteria_labels, 1, wx.EXPAND | wx.ALIGN_BOTTOM | wx.BOTTOM, v_space),
                                  (sizer_criteria_boxes, 1, wx.EXPAND |
                                   wx.BOTTOM | wx.ALIGN_TOP, v_space),
                                  (sizer_stats_boxes, 1, wx.EXPAND | wx.ALIGN_TOP)])

        # Logger Sizer--------------------------------------------------------
        sizer_logger = wx.BoxSizer(wx.HORIZONTAL)
        sizer_logger.Add(self.logger, 1, wx.EXPAND)

        # Set Panel Sizers----------------------------------------------------
        self.plot_panel.SetSizer(sizer_plots_outer)
        self.side_panel.SetSizerAndFit(sizer_logger)
        self.top_panel.SetSizerAndFit(sizer_top_bar)
        self.bottom_panel.SetSizerAndFit(sizer_bottom_bar)

        # Outer Sizer for Frame-----------------------------------------------
        sizer_logger_plots = wx.BoxSizer(wx.HORIZONTAL)
        sizer_logger_plots.Add(self.side_panel, 1, wx.EXPAND | wx.ALIGN_LEFT)
        sizer_logger_plots.Add(self.plot_panel, 3, wx.EXPAND | wx.ALIGN_LEFT)

        sizer_outer = wx.BoxSizer(wx.VERTICAL)
        sizer_outer.AddMany([(self.top_panel, 1, wx.EXPAND | wx.ALIGN_TOP | wx.BOTTOM, v_space / 2),
                             (sizer_logger_plots, 4, wx.EXPAND |
                              wx.ALIGN_TOP | wx.BOTTOM, v_space / 2),
                             (self.bottom_panel, 1, wx.EXPAND | wx.ALIGN_TOP)])

        self.SetSizer(sizer_outer)
        sizer_outer.Fit(self)
        self.Layout()
1011581
def ploidy(args):
    
    p = OptionParser(ploidy.__doc__)
    opts, args, iopts = p.set_image_options(args, figsize="8x7")

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

    seqidsfile, klayout = args

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

    Karyotype(fig, root, seqidsfile, klayout)

    fc = "darkslategrey"
    radius = .012
    ot = -.05  # use this to adjust vertical position of the left panel
    TextCircle(root, .1, .9 + ot, r'$\gamma$', radius=radius, fc=fc)
    root.text(.1, .88 + ot, r"$\times3$", ha="center", va="top", color=fc)
    TextCircle(root, .08, .79 + ot, r'$\alpha$', radius=radius, fc=fc)
    TextCircle(root, .12, .79 + ot, r'$\beta$', radius=radius, fc=fc)
    root.text(.1, .77 + ot, r"$\times3\times2\times2$", ha="center", va="top", color=fc)
    root.text(.1, .67 + ot, r"Brassica triplication", ha="center",
                va="top", color=fc, size=11)
    root.text(.1, .65 + ot, r"$\times3\times2\times2\times3$", ha="center", va="top", color=fc)
    root.text(.1, .42 + ot, r"Allo-tetraploidy", ha="center",
                va="top", color=fc, size=11)
    root.text(.1, .4 + ot, r"$\times3\times2\times2\times3\times2$", ha="center", va="top", color=fc)

    bb = dict(boxstyle="round,pad=.5", fc="w", ec="0.5", alpha=0.5)
    root.text(.5, .2 + ot, r"\noindent\textit{Brassica napus}\\"
                "(A$\mathsf{_n}$C$\mathsf{_n}$ genome)", ha="center",
                size=16, color="k", bbox=bb)

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

    pf = "napus"
    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)
1152758
def flanking(args):
    
    from numpy import array, argsort

    p = OptionParser(flanking.__doc__)
    p.add_option("--chrom", default=None, type="string",
            help="chrom name of the position in query. Make sure it matches bedfile.")
    p.add_option("--coord", default=None, type="int",
            help="coordinate of the position in query.")
    p.add_option("-n", default=10, type="int",
            help="number of flanking features to get [default: %default]")
    p.add_option("--side", default="both", choices=("upstream", "downstream", "both"),
            help="which side to get flanking features [default: %default]")
    p.add_option("--max_d", default=None, type="int",
            help="features <= max_d away from position [default: %default]")
    p.set_outfile()

    opts, args = p.parse_args(args)

    if any([len(args) != 1, opts.chrom is None, opts.coord is None]):
        sys.exit(not p.print_help())

    bedfile, = args
    position = (opts.chrom, opts.coord)
    n, side, maxd = opts.n, opts.side, opts.max_d

    chrombed = Bed(bedfile).sub_bed(position[0])

    if side == "upstream":
        data = [(abs(f.start-position[1]), f) for f in chrombed \
            if f.start <= position[1]]
    elif side == "downstream":
        data = [(abs(f.start-position[1]), f) for f in chrombed \
            if f.start >= position[1]]
    else:
        data = [(abs(f.start-position[1]), f) for f in chrombed]

    if maxd:
        data = [f for f in data if f[0]<=maxd]

    n += 1 # not counting self
    n = min(n, len(data))
    distances, subbed = zip(*data)
    distances = array(distances)
    idx = argsort(distances)[:n]
    flankingbed = [f for (i, f) in enumerate(subbed) if i in idx]

    fw = must_open(opts.outfile, "w")
    for atom in flankingbed:
        print(str(atom), file=fw)

    return (position, flankingbed)