def sample(args):
    """
    %prog sample bedfile sizesfile

    Sample bed file and remove high-coverage regions.

    When option --targetsize is used, this program uses a differnent mode. It
    first calculates the current total bases from all ranges and then compare to
    targetsize, if more, then sample down as close to targetsize as possible.

    Selection via --raindrop has the effect of making coverage even. Selected
    reads have the property that their end points are not within a certain
    window from one another. One sweep goes from left to right, the other in
    the reverse direction.
    """
    import random
    from jcvi.assembly.coverage import Coverage

    p = OptionParser(sample.__doc__)
    p.add_option("--raindrop", default=0, type="int",
                 help="Raindrop selection, ignores all other options")
    p.add_option("--max", default=10, type="int",
                 help="Max depth allowed")
    p.add_option("--targetsize", type="int",
                 help="Sample bed file to get target base number")
    p.set_outfile()
    opts, args = p.parse_args(args)

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

    bedfile, sizesfile = args
    pf = bedfile.rsplit(".", 1)[0]
    raindrop = opts.raindrop

    # Raindrop method
    if raindrop:
        bed = Bed(bedfile)
        forward = []
        for b in bed:
            if not forward or abs(b.start - forward[-1].start) >= raindrop:
                forward.append(b)

        reverse = []
        bed.sort(key=lambda x: -x.end)
        for b in bed:
            if not reverse or abs(b.end - reverse[-1].end) >= raindrop:
                reverse.append(b)

        for tag, L in zip(("forward", "reverse"), (forward, reverse)):
            logging.debug("Selected {0} features in {1} direction, span: {2}".\
                        format(len(L), tag, sum(x.span for x in L)))

        selected = Bed()
        selected.extend(set(forward + reverse))
        selected.print_to_file(opts.outfile, sorted=True)
        return

    targetsize = opts.targetsize
    if targetsize:
        bed = Bed(bedfile)
        samplebed = pf + ".sample.bed"
        fw = open(samplebed, "w")
        nfeats = len(bed)
        nbases = bed.sum(unique=False)
        targetfeats = int(round(nfeats * targetsize / nbases))
        sub_bed = random.sample(bed, targetfeats)
        for b in sub_bed:
            print >> fw, b

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

    c = Coverage(bedfile, sizesfile)
    coveragefile = c.filename
    samplecoveragefile = pf + ".sample.coverage"
    fw = open(samplecoveragefile, "w")
    fp = open(coveragefile)
    for row in fp:
        seqid, start, end, cov = row.split()
        cov = int(cov)
        if cov <= opts.max:
            fw.write(row)
    fw.close()

    samplebedfile = pf + ".sample.bed"
    cmd = "intersectBed -a {0} -b {1} -wa -u".format(bedfile, samplecoveragefile)
    sh(cmd, outfile=samplebedfile)
    logging.debug("Sampled bedfile written to `{0}`.".format(samplebedfile))
dataset/ETHPy150Open tanghaibao/jcvi/formats/bed.py/sample/original