%prog astat coverage.log

    Create coverage-rho scatter plot.
775503
def astat(args):
    
    p = OptionParser(astat.__doc__)
    p.add_option("--cutoff", default=1000, type="int",
                 help="Length cutoff [default: %default]")
    p.add_option("--genome", default="",
                 help="Genome name [default: %default]")
    p.add_option("--arrDist", default=False, action="store_true",
                 help="Use arrDist instead [default: %default]")
    opts, args = p.parse_args(args)

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

    covfile, = args
    cutoff = opts.cutoff
    genome = opts.genome
    plot_arrDist = opts.arrDist

    suffix = ".{0}".format(cutoff)
    small_covfile = covfile + suffix
    update_covfile = need_update(covfile, small_covfile)
    if update_covfile:
        fw = open(small_covfile, "w")
    else:
        logging.debug("Found `{0}`, will use this one".format(small_covfile))
        covfile = small_covfile

    fp = open(covfile)
    header = next(fp)
    if update_covfile:
        fw.write(header)

    data = []
    msg = "{0} tigs scanned ..."
    for row in fp:
        tigID, rho, covStat, arrDist = row.split()
        tigID = int(tigID)
        if tigID % 1000000 == 0:
            sys.stderr.write(msg.format(tigID) + "\r")

        rho, covStat, arrDist = [float(x) for x in (rho, covStat, arrDist)]
        if rho < cutoff:
            continue

        if update_covfile:
            fw.write(row)
        data.append((tigID, rho, covStat, arrDist))

    print(msg.format(tigID), file=sys.stderr)

    from jcvi.graphics.base import plt, savefig

    logging.debug("Plotting {0} data points.".format(len(data)))
    tigID, rho, covStat, arrDist = zip(*data)

    y = arrDist if plot_arrDist else covStat
    ytag = "arrDist" if plot_arrDist else "covStat"

    fig = plt.figure(1, (7, 7))
    ax = fig.add_axes([.12, .1, .8, .8])
    ax.plot(rho, y, ".", color="lightslategrey")

    xtag = "rho"
    info = (genome, xtag, ytag)
    title = "{0} {1} vs. {2}".format(*info)
    ax.set_title(title)
    ax.set_xlabel(xtag)
    ax.set_ylabel(ytag)

    if plot_arrDist:
        ax.set_yscale('log')

    imagename = "{0}.png".format(".".join(info))
    savefig(imagename, dpi=150)
719576
def postprocess(options):
    

    resdir = options.resdir
    out_file = options.outfile
    tol = options.tol

    print('.. load permutation results')
    file_name = os.path.join(resdir,'perm*','*.res')
    files = glob.glob(file_name)
    LLR0 = []
    for _file in files:
        print(_file)
        LLR0.append(NP.loadtxt(_file,usecols=[6]))
    LLR0 = NP.concatenate(LLR0)

    print('.. fit test statistics')
    t0 = time.time()
    c2m = C2M.Chi2mixture(tol=4e-3)
    c2m.estimate_chi2mixture(LLR0)
    pv0 = c2m.sf(LLR0)
    t1 = time.time()
    print(('finished in %s seconds'%(t1-t0)))

    print('.. export permutation results')
    perm_file = out_file+'.perm'
    RV = NP.array([LLR0,pv0]).T
    NP.savetxt(perm_file,RV,delimiter='\t',fmt='%.6f %.6e')

    print('.. load test results')
    file_name = os.path.join(resdir,'test','*.res')
    files = glob.glob(file_name)
    RV_test = []
    for _file in files:
        print(_file)
        RV_test.append(NP.loadtxt(_file))
    RV_test = NP.concatenate(RV_test)

    print('.. calc pvalues')
    pv = c2m.sf(RV_test[:,-1])[:,NP.newaxis]

    print('.. export test results')
    perm_file = out_file+'.test'
    RV_test = NP.hstack([RV_test,pv])
    NP.savetxt(perm_file,RV_test,delimiter='\t',fmt='%d %d %d %d %d %d %.6e %.6e')

    if options.manhattan:
        manhattan_file = out_file+'.manhattan.jpg'
        plot_manhattan(pv,manhattan_file)
1050005
def _bam_coverage(name, bam_input, data):
    """"""
    cmd = ("{bam_coverage} -b {bam_input} -o {bw_output} "
          "--binSize 20 --effectiveGenomeSize {size} "
          "--smoothLength 60 --extendReads 150 --centerReads -p {cores}")
    size = bam.fasta.total_sequence_length(dd.get_ref_file(data))
    cores = dd.get_num_cores(data)
    try:
        bam_coverage = config_utils.get_program("bamCoverage", data)
    except config_utils.CmdNotFound:
        logger.info("No bamCoverage found, skipping bamCoverage.")
        return None
    resources = config_utils.get_resources("bamCoverage", data["config"])
    if resources:
        options = resources.get("options")
        if options:
            cmd += " %s" % " ".join([str(x) for x in options])
    bw_output = os.path.join(os.path.dirname(bam_input), "%s.bw" % name)
    if utils.file_exists(bw_output):
        return bw_output
    with file_transaction(bw_output) as out_tx:
        do.run(cmd.format(**locals()), "Run bamCoverage in %s" % name)
    return bw_output
443013
def prepare(args):
    
    p = OptionParser(prepare.__doc__)
    opts, args = p.parse_args(args)

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

    counts, families = args
    countfiles = glob(op.join(counts, "*.count"))
    countsdb = defaultdict(list)
    for c in countfiles:
        rs = RiceSample(c)
        countsdb[(rs.tissue, rs.ind)].append(rs)

    # Merge duplicates - data sequenced in different batches
    key = lambda x: (x.label, x.rep)
    for (tissue, ind), rs in sorted(countsdb.items()):
        rs.sort(key=key)
        nrs = len(rs)
        for i in xrange(nrs):
            ri = rs[i]
            if not ri.working:
                continue
            for j in xrange(i + 1, nrs):
                rj = rs[j]
                if key(ri) != key(rj):
                    continue
                ri.merge(rj)
                rj.working = False
        countsdb[(tissue, ind)] = [x for x in rs if x.working]

    # Group into families
    mkdir("families")
    for (tissue, ind), r in sorted(countsdb.items()):
        r = list(r)
        if r[0].label != "F1":
            continue
        P1, P2 = r[0].P1, r[0].P2
        P1, P2 = countsdb[(tissue, P1)], countsdb[(tissue, P2)]
        rs = P1 + P2 + r
        groups = [1] * len(P1) + [2] * len(P2) + [3] * len(r)
        assert len(rs) == len(groups)

        outfile = "-".join((tissue, ind))
        merge_counts(rs, op.join(families, outfile))
        groupsfile = outfile + ".groups"
        fw = open(op.join(families, groupsfile), "w")
        print(",".join(str(x) for x in groups), file=fw)
        fw.close()
558665
def create_report(self):
        """"""
        data = self.create_data()
        try:
            json_string = json.dumps(data)
        except UnicodeDecodeError as e:
            log.error('ERROR: While preparing JSON:', exc_info=e)
            self.debug_bad_encoding(data)
            raise

        log_string = re.sub(r'"repo_token": "(.+?)"',
                            '"repo_token": "[secure]"', json_string)
        log.debug(log_string)
        log.debug('==\nReporting %s files\n==\n', len(data['source_files']))
        for source_file in data['source_files']:
            log.debug('%s - %s/%s', source_file['name'],
                      sum(filter(None, source_file['coverage'])),
                      len(source_file['coverage']))
        return json_string
538220
def plot(self):
        ''''''
        import matplotlib.pyplot as plt
        import matplotlib.patches as patches

        # Constants to use throughout drawing
        n = len(self.results)
        nbases = len(self.aligned_reference)
        barheight = 0.4

        # Vary height of figure based on number of results
        figheight = 3 + 3 * (n - 1)
        fig = plt.figure(figsize=(9, figheight))
        ax1 = fig.add_subplot(111)

        # Plot bars to represent coverage area
        # Reference sequence
        ax1.add_patch(patches.Rectangle((0, 0), nbases, barheight,
                                        facecolor='black'))
        # Results
        for i, report in enumerate(self.nonmatches()):
            j = i + 1
            start, stop = report['coverage']
            patch = patches.Rectangle((start, j), stop - start, barheight,
                                      facecolor='darkgray')
            ax1.add_patch(patch)

            # Draw a vertical line for each type of result
            plt.vlines(report['mismatches'], j, j + barheight,
                       colors='b')
            plt.vlines(report['insertions'], j, j + barheight,
                       colors='r')

            # Terminal trailing deletions shouldn't be added
            deletions = []
            crange = range(*report['coverage'])
            deletions = [idx for idx in report['deletions'] if idx in crange]
            plt.vlines(deletions, j, j + barheight,
                       colors='g')

        ax1.set_xlim((0, nbases))
        ax1.set_ylim((-0.3, n + 1))
        ax1.set_yticks([i + barheight / 2 for i in range(n + 1)])
        ax1.set_yticklabels(['Reference'] + self.names)

        # Add legend
        mismatch_patch = patches.Patch(color='blue', label='Mismatch')
        insertion_patch = patches.Patch(color='red', label='Insertion')
        deletion_patch = patches.Patch(color='green', label='Deletion')
        plt.legend(handles=[mismatch_patch, insertion_patch, deletion_patch],
                   loc=1, ncol=3, mode='expand', borderaxespad=0.)

        plt.show()
547969
def statstable(args):
    
    p = OptionParser(statstable.__doc__)
    opts, args = p.parse_args(args)

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

    gff_files = args
    for metric in metrics:
        logging.debug("Parsing files in `{0}`..".format(metric))

        table = {}
        for x in gff_files:
            pf = op.basename(x).split(".")[0]
            numberfile = op.join(metric, pf + ".txt")
            ar = [int(x.strip()) for x in open(numberfile)]
            sum = SummaryStats(ar).todict().items()
            keys, vals = zip(*sum)
            keys = [(pf, x) for x in keys]
            table.update(dict(zip(keys, vals)))

        print(tabulate(table), file=sys.stderr)
190286
function(app, options) {
  app.post(WRITE_COVERAGE,
    bodyParser,
    (req, res) => {
      let map = istanbul.libCoverage.createCoverageMap();
      coverageHandler(map, options, req, res);
    },
    logError);
}
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))
1285056
def main(argv=None):  # pragma: no coverage
    
    try:
        colorama.init()
        if argv is None:
            argv = sys.argv[1:]
        _main(argv)
    except RuntimeError as e:
        print(colorama.Fore.RED + 'ERROR: ' +
              str(e) + colorama.Style.RESET_ALL)
        sys.exit(1)
    else:
        sys.exit(0)
262743
def CheckProg(context, prog_name):
    
    res = SCons.Conftest.CheckProg(context, prog_name)
    context.did_show_result = 1
    return res
325981
def plot(args):
    
    from jcvi.graphics.base import plt, savefig
    from jcvi.graphics.chromosome import ChromosomeMap

    p = OptionParser(plot.__doc__)
    p.add_option("--firstn", type="int", help="Only plot the first N genes")
    p.add_option("--ymax", type="int", help="Y-axis max value")
    p.add_option("--log", action="store_true",
                help="Write plotting data [default: %default]")
    opts, args, iopts = p.set_image_options(args, figsize="6x4")

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

    taggedbed, chr = args
    bed = Bed(taggedbed)
    beds = list(bed.sub_bed(chr))
    old, new = [], []
    i = 0
    for b in beds:
        accn = b.extra[0]
        if "te" in accn:
            continue

        accn, tag = accn.split("|")
        if tag == "OVERLAP":
            continue

        c, r = atg_name(accn)
        if tag == "NEW":
            new.append((i, r))
        else:
            old.append((i, r))
        i += 1

    ngenes = i
    assert ngenes == len(new) + len(old)

    logging.debug("Imported {0} ranks on {1}.".format(ngenes, chr))
    fig = plt.figure(1, (iopts.w, iopts.h))
    root = fig.add_axes([0, 0, 1, 1])

    xstart, xend = .2, .8
    ystart, yend = .2, .8
    pad = .02

    ngenes = opts.firstn or ngenes
    ymax = opts.ymax or 500000

    title = "Assignment of Medtr identifiers"
    if opts.ymax:
        subtitle = "{0}, first {1} genes".format(chr, ngenes)
    else:
        subtitle = "{0}, {1} genes ({2} new)".format(chr, ngenes, len(new))

    chr_map = ChromosomeMap(fig, root, xstart, xend, ystart, yend, pad, 0,
                        ymax, 5, title, subtitle)

    ax = chr_map.axes

    if opts.log:
        from jcvi.utils.table import write_csv
        header = ["x", "y"]
        write_csv(header, new, filename=chr + ".new")
        write_csv(header, old, filename=chr + ".old")

    x, y = zip(*new)
    ax.plot(x, y, "b,")
    x, y = zip(*old)
    ax.plot(x, y, "r,")

    # Legends
    ymid = (ystart + yend) / 2
    y = ymid + pad
    root.plot([.2], [y], "r.", lw=2)
    root.text(.2 + pad, y, "Existing Medtr ids", va="center", size=10)
    y = ymid - pad
    root.plot([.2], [y], "b.", lw=2)
    root.text(.2 + pad, y, "Newly instantiated ids", va="center", size=10)

    ax.set_xlim(0, ngenes)
    ax.set_ylim(0, ymax)
    ax.set_axis_off()

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

    image_name = chr + ".identifiers." + iopts.format
    savefig(image_name, dpi=iopts.dpi, iopts=iopts)
1131440
def plot_prep_methods(df, prep, prepi, out_file_base, outtype, title=None,
                      size=None):
    
    samples = df[(df["bamprep"] == prep)]["sample"].unique()
    assert len(samples) >= 1, samples
    out_file = "%s-%s.%s" % (out_file_base, samples[0], outtype)
    df = df[df["category"].isin(cat_labels)]
    _seaborn(df, prep, prepi, out_file, title, size)
    return out_file
303123
public static void main(String[] args) {
        // Parse arguments.
        String coverageDirName = null;
        if (args.length == 0) {
            System.out.println("Incorrect arguments.  Directory with coverage has to be specified.");
            System.exit(1);
        }
        coverageDirName = args[0];
        String mode = null;
        if (args.length > 1) {
            mode = args[1];
        }
        boolean forceCacheUse = false;
        if (args.length > 2) {
            forceCacheUse = args[2].equals(FORCE_CACHE_USE);
        }

        Set<String> allClasses = new HashSet<String>();
        Set<String> affectedClasses = new HashSet<String>();
        if (args.length > 3) {
            String options = args[3];
            Config.loadConfig(options, true);
        } else {
            Config.loadConfig();
        }

        List<String> nonAffectedClasses = findNonAffectedClasses(coverageDirName, forceCacheUse, allClasses, affectedClasses);
        
        // Print non affected classes.
        printNonAffectedClasses(allClasses, affectedClasses, nonAffectedClasses, mode);
    }
640937
def coverage(args):
    
    p = OptionParser(coverage.__doc__)
    p.add_option("-c", dest="cutoff", default=0.5, type="float",
            help="only report query with coverage greater than [default: %default]")

    opts, args = p.parse_args(args)

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

    coordsfile, = args
    fp = open(coordsfile)

    coords = []
    for row in fp:
        try:
            c = CoordsLine(row)
        except AssertionError:
            continue
        coords.append(c)

    coords.sort(key=lambda x: x.query)

    coverages = []
    for query, lines in groupby(coords, key=lambda x: x.query):
        cumulative_cutoff = sum(x.querycov for x in lines)
        coverages.append((query, cumulative_cutoff))

    coverages.sort(key=lambda x: (-x[1], x[0]))
    for query, cumulative_cutoff in coverages:
        if cumulative_cutoff < opts.cutoff:
            break
        print("{0}\t{1:.2f}".format(query, cumulative_cutoff))
246430
def count(args):
    
    p = OptionParser(count.__doc__)
    opts, args = p.parse_args(args)

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

    coveragefile, fastafile = args

    countsfile = coveragefile.split(".")[0] + ".bin"
    if op.exists(countsfile):
        logging.error("`{0}` file exists. Remove before proceed."\
                .format(countsfile))
        return

    fastasize, sizes, offsets = get_offsets(fastafile)
    logging.debug("Initialize array of uint8 with size {0}".format(fastasize))
    ar = np.zeros(fastasize, dtype=np.uint8)

    update_array(ar, coveragefile, sizes, offsets)

    ar.tofile(countsfile)
    logging.debug("Array written to `{0}`".format(countsfile))
635109
def stats(args):
    
    from jcvi.utils.cbook import percentage

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

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

    blocksfile, = args
    fp = open(blocksfile)
    counts = defaultdict(int)
    total = orthologous = 0
    for row in fp:
        atoms = row.rstrip().split("\t")
        hits = [x for x in atoms[1:] if x != '.']
        counts[len(hits)] += 1
        total += 1
        if atoms[1] != '.':
            orthologous += 1

    print("Total lines: {0}".format(total), file=sys.stderr)
    for i, n in sorted(counts.items()):
        print("Count {0}: {1}".format(i, percentage(n, total)), file=sys.stderr)

    print(file=sys.stderr)

    matches = sum(n for i, n in counts.items() if i != 0)
    print("Total lines with matches: {0}".\
                format(percentage(matches, total)), file=sys.stderr)
    for i, n in sorted(counts.items()):
        if i == 0:
            continue

        print("Count {0}: {1}".format(i, percentage(n, matches)), file=sys.stderr)

    print(file=sys.stderr)
    print("Orthologous matches: {0}".\
                format(percentage(orthologous, matches)), file=sys.stderr)